/*
 * dcn.c		Jim Piper		04-09-86
 *
 * find and display dicentric centromere candidates, perform box classification
 *
 * dcn [-d <stoxxxx>] [-p] [-m <PC1> <PC2>] [-q <QUANTUM>] [file]
 */

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <wstruct.h>
#include <chromanal.h>
#include <dicen.h>
#include <otypes.h>

#define SCALE	24


static struct pframe f;
static int wh[2] = {0,1};

/* which features (numbering from zero - beware, dicen.h numbers from 1) */
int nfeat = 9;
int selfeat[] = {0,2,4,5,7,8,12,13,14};


main(argc,argv)
char **argv;
{
	struct object *spobj, *mpol, *cvh, *rect;
	struct chromosome *obj, *readobj();
	struct chromplist *plist, *makechromplist();
	struct object *fipconvhull(), *spinsqueeze(), *minwrect();
	struct object *cpol, *centpoly(), *chromaxis2();
	struct object *prof, *p[2], *cprof, *polyprof();
	struct intervaldomain *idom;
	struct polygondomain *poly, *uspaxis8();
	struct frect *rdom;
	struct chromcands *chcp,chc[MAXOBJS];
	struct cand *cndp,cnd[MAXCANDS];
	int i, j, number, numberused, Cbtype[MAXOBJS];
	int PC1, PC2, QUANTUM;
	char *cell,s[256],*timestring;
	FILE *infile, *pf;
	long daytime;
	int DISPLAY = 0;
	int PRINT = 0;

	daytime = time();
	timestring = ctime(&daytime);
	timestring[24] = ' ';

	/*
	 * defaults for PC1, PC2, QUANTUM
	 */
	PC1 = 90;
	PC2 = 120;
	QUANTUM = 12;

	while (argc > 1 && argv[1][0] == '-') {
		switch (argv[1][1]) {
		case 'd':
			DISPLAY++;
			setout(argv[2],'p');
			argc--;
			argv++;
			break;
		case 'p':
			PRINT++;
			break;
		case 'm':
			PC1 = atoi(argv[2]);
			PC2 = atoi(argv[3]);
			argc -= 2;
			argv += 2;
			break;
		case 'q':
			QUANTUM = atoi(argv[2]);
			argc--;
			argv++;
			break;
		default:
			fprintf(stderr,"unrecogised switch -%c\n",argv[1][1]);
			break;
		}
		argc--;
		argv++;
	}
	if (argc > 1) {
		cell = argv[1];
		infile = fopen(cell,"r");
		if (infile == NULL) {
			fprintf(stderr,"Can't open %s\n",cell);
			exit(1);
		}
	} else {
		cell = "<stdin";
		infile = stdin;
	}
	pf = stdout;
	if (DISPLAY != 0) {
		dinit();
		lwidth(3);
	}
	fprintf(pf,"CELL FILE %s PC1 %d PC2 %d QUANTUM %d\n",cell,PC1,PC2,QUANTUM);
	if (PRINT >= 1) {
		fprintf(pf,"FEATURES USED BY BOX CLASSIFIER: ");
		for (i=0; i<nfeat; i++)
			fprintf(pf,"%2d ",selfeat[i]+1);
		fprintf(pf,"\n\n");
	}
	numberused = number = 0;
	chcp = chc;
	cndp = cnd;
	while ((obj = readobj(infile)) != NULL) {
		number++;
		if (obj->plist == NULL) {
			obj->plist = makechromplist();
		}
		obj->plist->area = area(obj);
		obj->plist->mass = mass(obj);
		if (obj->plist->area < 100) {
			if (PRINT >= 1)
				fprintf(pf,"tiny object %d ignored\n",number);
		}
		else if (obj->plist->otype == 2 && (obj->plist->Cbtype & DICENTRIC))
			fprintf(pf,"MISSED3 DICENTRIC (false composite) obj %d\n",number);
		else if (obj->plist->otype <= 1) {
			Cbtype[numberused] = obj->plist->Cbtype;
			numberused++;
			f.type = 60;
			f.scale = SCALE/8;
			f.ix = f.iy = 8;
			f.dx = 2048;
			f.dy = 2048;
			cvh = fipconvhull(obj);
			rect = minwrect(cvh);
			rdom = (struct frect *)rect->idom;
			/*
			 * rotate to vertical, find axis,
			 * find density and 1st moment profile
			 */
			spobj = spinsqueeze(obj,rdom->rangle,XSCALE,YSCALE);
			mpol = chromaxis2(spobj);
			multiprofiles(spobj,mpol,p,wh,2);
			shorten(mpol,p[0],p[1]);
			prof = p[0];

			/*
			 * find centromere candidates from density profile,
			 * measure features, display
			 */
			plist = obj->plist;
			chcp->number = number;
			chcp->ncands = 0;
			chcp->otype = plist->otype;
			chcp->Cotype = plist->Cotype;
			chcp->area = plist->area;
			chcp->density = plist->mass / plist->area;
			chcp->cands = cndp;
			dcmerecands(prof,QUANTUM,chcp);
			if (DISPLAY >= 1) {
				clear();
				logo();
				sprintf(s,"cell-file %s object %d",cell,number);
				dispstring(s,100,200);
				dispstring(timestring,100,100);
				idom = spobj->idom;
				f.ox = (idom->kol1 + idom->lastkl)/2;
				f.oy = (idom->line1 + idom->lastln)/2;
				picframe(spobj,&f);
				picframe(mpol,&f);
				f.dx += 1024;
				f.ix = 1;
				histo256normalise(prof);
				picframe(prof,&f);
				f.dx -= 1024;
				sprintf(s,"%d candidates from profile; peak %d",chcp->ncands,chcp->profmax);
				dispstring(s,100,3900);
			}
			for (i=0; i<chcp->ncands; i++) {
				cpol = centpoly(mpol,cndp->profpos);
				poly = (struct polygondomain *)cpol->idom;
				cpol->idom = (struct intervaldomain *)uspaxis8(poly);
				free(poly);
				cprof = polyprof(spobj,cpol);
				dcfeatures(chcp,cndp,cprof,prof);
				if (DISPLAY >= 1) {
					sprintf(s,"cand %d pos %d val %d",cndp->candnum,cndp->profpos,cndp->u.sy.profval);
					dispstring(s,100,3800-100*i);
					f.ix = f.iy = 1;
					f.ox *= 8;
					f.oy *= 8;
					picframe(cpol,&f);
					f.ix = 8;
					f.ox /= 8;
					f.oy /= 8;
					f.dx -= 1024;
					f.dy = 1024 + i*512;
					picframe(cprof,&f);
					f.iy = 8;
					f.dx += 1024;
					f.dy = 2048;
					update();
				}
				cndp++;
				freeobj(cpol);
				freeobj(cprof);
			}
			if (PRINT >= 1)
				dcprint(chcp,pf,nfeat,selfeat);

			chcp++;

			freeobj(spobj);
			freeobj(mpol);
			freeobj(p[0]);
			freeobj(p[1]);
			freeobj(cvh);
			freeobj(rect);
			if (DISPLAY >= 1)
				sleep(1);
		}
		freeobj(obj);
	}
	/*
	 * sumary about candidates
	 */
	dccandsummary(chc,numberused,pf);
	/*
	 * normalise features
	 */
	dcnormalise(chc,numberused);
	/*
	 * do the classification
	 */
	dcboxclass(chc,numberused,nfeat,selfeat,PC1,PC2,PRINT,pf);
	dcresprint(chc,numberused,Cbtype,pf);
	/*
	 * display the found true and false positives
	 * and the false negatives
	 */
	if (DISPLAY >= 1)
		finish();
}

logo()
{
	moveto(2,2);
	lwidth(3);
	lineby(4092,0);
	lineby(0,4092);
	lineby(-4092,0);
	lineby(0,-4092);
	lwidth(1);
	user(9,0);
}