/*
 * dcsk.c		02-10-89		Jim Piper
 *
 * Extract centromere data for Simon Kirby.
 * Program based on "dcdemo.c"
 *
 * dcdemo [-p] [-q <QUANTUM>] [-P] [-B] [-M] [-o ofile] [file]
 *
 * FLAGS:
 * -----
 * -p		turn on printing
 * -q quantum	quantum (factor) for minimum detection in profile
 * -P		extract profile generated candidates
 * -B		extract boundary generated candidates
 * -M		refine position of boundary candidate crossing profile
 *		if there is an nearby profile local minimum.
 * -o ofile	send feature measurements to ofile
 */

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

#define RADIANS 57.296
#define FEATTEMP "/tmp/dcdemo.temp"

int QUANTUM, PROF_CANDS, BOUND_CANDS;
int DISPLAY,PRINT,BPOLMOD;
int POSPRINT = 0;
int CANDSEL = 1;
int OTHNEARDIST = 10;
int SELNEARDIST = 8;
int ENDDIST = 6;
int BOXMETHOD = 3;
int REVIEWP = 0;
int Cbtype[MAXOBJS];
static struct pframe f;
int PC1 = 60, PC2 = 120;
/* which features (numbering from one, as in dicen.h) */
int nfeat = 4;
int selfeat[MAXFEATURES] = { 9,13,18,19};
int selllim[MAXFEATURES] = { 1, 1, 1, 1};
static FILE *fin, *outfile;
char *outfilename;

main(argc,argv)
char **argv;
{
	int i,objnum,number,nonum;
	int scale,cellw,celll;
	int numberused, cellkol1,cellline1,celllastkl,celllastln;
	long daytime;
	struct chromosome *objlist[200], *obj, *readobj();
	struct chromcands chc[MAXOBJS];
	struct cand cnd[MAXCANDS];
	char s[512], *cell, *timestring;
	FILE *infile;
	int FIP=1;

	/*
	 * defaults
	 */
	DISPLAY = PRINT = BPOLMOD = PROF_CANDS = BOUND_CANDS = 0;
	QUANTUM = 12;
	scale = -1;
	nonum = 1;
	outfile = NULL;
	outfilename = "<no output>";
	setout("stosun -p270,0 -d1 -r870,870",'p');
	setbuf(stderr,NULL);

	while (argc > 1 && argv[1][0] == '-') {
		switch (argv[1][1]) {
		case 'p':
			PRINT++;
			break;
		case 'q':
			QUANTUM = atoi(argv[2]);
			argc--;
			argv++;
			break;
		case 'o':
			outfilename = argv[2];
			outfile = fopen(outfilename,"w");
			if (outfile == NULL) {
				fprintf(stderr,"Cannot open %s for output\n",outfilename);
				exit(1);
			}
			argc--;
			argv++;
			break;
		case 'B':
			BOUND_CANDS++;
			break;
		case 'P':
			PROF_CANDS++;
			break;
		case 'M':
			BPOLMOD++;
			break;

		case 'R':
			REVIEWP++;
			break;
		case 'm':
			PC1 = atoi(argv[2]);
			PC2 = atoi(argv[3]);
			argc -= 2;
			argv += 2;
			break;
		case 'C':
			CANDSEL = atoi(argv[2]);
			argc--;
			argv++;
			if (CANDSEL != 1 && CANDSEL != 2) {
				fprintf(stderr,"Illegal -C value\n");
				exit(1);
			}
			break;
		case 'D':
			OTHNEARDIST = atoi(argv[2]);
			argc--;
			argv++;
			break;
		case 'E':
			ENDDIST = atoi(argv[2]);
			argc--;
			argv++;
			break;
		case 'd':
			SELNEARDIST = atoi(argv[2]);
			argc--;
			argv++;
			break;
		case 'F':
			nfeat = atoi(argv[2]);
			argc--;
			argv++;
			for (i=0; i<nfeat; i++) {
				selfeat[i] = atoi(argv[2]);
				argc--;
				argv++;
			}
			break;
		case 'X':
			BOXMETHOD = atoi(argv[2]);
			argc--;
			argv++;
			if (BOXMETHOD < 0 || BOXMETHOD > 7) {
				fprintf(stderr,"Unknown box method\n");
				exit(1);
			}
			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;
	}
	number = 0;
	while((obj = readobj(infile)) != NULL) {
		if (obj->plist != NULL && (obj->plist->number <= 0 || obj->plist->number > MAXOBJS))
			swabplist(obj->plist);
		if (obj->type == 1 && (obj->plist==NULL || obj->plist->otype <= 1)) {
			objlist[number++] = obj;
			if (number == 1) {
				cellline1 = obj->idom->line1;
				cellkol1 = obj->idom->kol1;
				celllastln = obj->idom->lastln;
				celllastkl = obj->idom->lastkl;
			} else {
				cellline1 = min(cellline1,obj->idom->line1);
				cellkol1 = min(cellkol1,obj->idom->kol1);
				celllastln = max(celllastln,obj->idom->lastln);
				celllastkl = max(celllastkl,obj->idom->lastkl);
			}
		}
	}
	normaliseset(objlist,number);

	f.type = 60;
	f.dx = f.dy = 2048;
	f.ix = f.iy = 1;
	f.ox = (cellkol1 + celllastkl)/2;
	f.oy = (cellline1 + celllastln)/2;

	if (scale < 0) {
		celll = celllastln - cellline1 + 1;
		if (FIP)
			celll = 3*celll/4;
		cellw = celllastkl - cellkol1 + 1;
		f.scale = 3840/max(celll,cellw);
		if (f.scale > 8)
		f.scale = 8;
	}
	if (FIP) {
		f.ix = f.scale;
		f.iy = (3*f.scale+2)/4;
		f.scale = 1;
	}

	if (number > 0) {
		fin = fopen("/dev/tty","r");
		dinit();
		lwidth(2);
		colourlut();
		colour(248);
		for (objnum=1; objnum<=number; objnum++)
			picframe(objlist[objnum-1],&f);
		tsize(50,80);
		
		/*
		 * measure features
		 */
		colour(2);
		tsize(70,100);
		dispstring("Profile candidate metacentric centromeres",32,4000);
		colour(4);
		dispstring("Boundary candidate metacentric centromeres",32,3900);
		numberused = dccfmain(objlist,number,chc,cnd);
		update();
		fprintf(stderr,"Type return to exit\n");
		getc(fin);
		setrgb(0,255,255,180);	/* restore LUT */
		update();
		finish();
	}
}


max(i,j)
{
	return ( i>j? i: j);
}


min(i,j)
{
	return ( i<j? i: j);
}


#define SCALE	32
#define FIP 1			/* FIP digitisation assumed */

/* which profiles to be obtained from from multiprofiles */
static int wh[2] = {0,1};


dccfmain(objl,totalnumber,chc,cnd)
struct chromcands *chc;
struct cand *cnd;
struct chromosome **objl;
{
	struct object *spobj, *mpol, *oli[20];
	struct chromosome *obj;
	struct chromplist *plist, *makechromplist();
	struct object *dcvertical();
	struct object *cpol, *centpoly(), *boundpoly(), *chromaxis2();
	struct object *prof, *p[2], *cprof, *polyprof();
	struct polygondomain *poly, *uspaxis8();
	struct chromcands *chcp;
	struct cand *cndp,*cndp1;
	int c, i, j, dist, number, numberused, nlab, ao, amin, imin;
	int totalcands, sernum, otype;
	char *cell,s[256];
	FILE *pf;
	struct polygondomain *boundary;
	struct histogramdomain *hist, *hist_diff, *bound_curv();

	totalcands = 0;
	pf = stdout;

	if (PROF_CANDS == 0 && BOUND_CANDS == 0)
		PROF_CANDS = BOUND_CANDS = 1;
	if (PRINT) {
		fprintf(pf,"CELL FILE %s DATA TO %s QUANTUM %d\n",cell,outfilename,QUANTUM);
		fprintf(pf,"Including  ");
		if (PROF_CANDS)
			fprintf(pf,"profile candidates  ");
		if (BOUND_CANDS)
			fprintf(pf,"boundary candidates  ");
		fprintf(pf,"\n");
	}
	numberused = number = sernum = 0;
	chcp = chc;
	cndp = cnd;
	while (sernum < totalnumber) {
		obj = objl[sernum++];
		/*
		 * give it a propertylist and serial number if necessary
		 */
		if (obj->plist == NULL) {
			obj->plist = makechromplist();
			++number;
		}
		else {
			number = obj->plist->number;
		}
		/*
		 * compute essential features, reject tiny objects
		 */
		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);
		}

		/*----------------------------------------------------------------------*/
		/*
		 * only analyse objects believed to be chromosome (i.e.
		 * NOT composite, noise, overlap, etc..) by any previous pass.
		 */
		else if (obj->plist->otype <= 1) {
			if (PRINT >= 1)
				fprintf(pf,"object %d area %d mass %d\n",number,obj->plist->area,obj->plist->mass);
			chcp->Cbtype = obj->plist->Cbtype;
			numberused++;

			/*
			 * rotate to vertical, find axis,
			 * find density and 1st moment profile
			 */
			spobj = dcvertical(obj,FIP);
			/*
			 * if disconnected, find biggest piece
			 */
			label(spobj,&nlab,oli,20,5);
			amin = imin = 0;
			if (nlab > 1) {
				fprintf(stderr,"disconnected object\n");
				for (i=0; i<nlab; i++) {
					ao = area(oli[i]);
					if (ao > amin) {
						amin = ao;
						imin = i;
					}
				}
			}
			else if (nlab <= 0) {
				fprintf(stderr,"object atomised !!\n");
			}
			else {
				oli[imin]->vdom = spobj->vdom;
				oli[imin]->plist = spobj->plist;
				spobj = oli[imin];
			}

			mpol = chromaxis2(spobj);
			multiprofiles(spobj,mpol,p,wh,2);
			shorten(mpol,p[0],p[1]);
			prof = p[0];

			/*
			 * find linear piece-wise boundary
			 * and associated histograms
			 */
			hist = bound_curv(spobj, &boundary, &hist_diff);
			
			/*
			 * identify composite/overlap objects
			 */
		/*	bound_comp(obj, spobj, boundary, hist, hist_diff);	*/
			
			/*
			 * fill up top of chromcands structure -
			 * WATCH OUT FOR SWABs !!!
			 */
			plist = obj->plist;
			chcp->number = number;
			chcp->ncands = 0;
			chcp->otype = plist->otype;
			chcp->Cotype = plist->Cotype;
			if (plist->pnumber != 0)
				chcp->porigin = plist->history[0];
			chcp->area = plist->area;
			chcp->axislen = ((struct histogramdomain *)prof->idom)->npoints;
			chcp->density = plist->mass / plist->area;
			chcp->cands = cndp;
			prof_basics(prof,QUANTUM,chcp);

			/*
			 * find centromere candidates from density profile
			 */
			if (PROF_CANDS)
				prof_cands(prof,chcp);
			/*
			 * find centromere candidates from boundary profile
			 */
			if (BOUND_CANDS) {
				boundary_cands(spobj, boundary, hist, hist_diff, chcp);
			}
			/*
			 * Check on total number of candidates in this cell
			 */
			totalcands += chcp->ncands;
			if (totalcands > MAXCANDS - 20) {
				fprintf(stderr,"Too many candidates for comfort\n");
				break;	/* from outside while */
			}
			/*--------------------------------------------------------------*/
			/*
			 * For each candidate:
			 * If a profile candidate, find "crossing-line".
			 * If a boundary candidate, find "modified crossing-line".
			 * In each case, the relevant routine returns an
			 * integer-scaled-by-8 polygon.
			 */
			for (i=0; i<chcp->ncands; i++, cndp++) {
  			    switch (cndp->candtype) {
			    case 1:
			        cpol = centpoly(mpol,cndp->profpos);
			        break;
			    case 2:
			        cpol = boundpoly(boundary, cndp, spobj);
				if (cpol == NULL) {
					cndp->candtype = -1;
					fprintf(stderr,"obj %d cand %d no cpol\n",number,i+1);
					continue;
				}
			        if (BPOLMOD)
			        	dcperturb(cpol,mpol,prof);
			        break;
			    default:
			        continue;
			    }
			/*
			 * unit-space the crossing-line
			 * polygon (still integer-scaled-by-8)
			 */
			    poly = (struct polygondomain *)cpol->idom;
			    cpol->idom = (struct intervaldomain *)uspaxis8(poly);
			    free(poly);
			/*
			 * Extract density values along the crossing polygon -
			 * the "crossing profile".
			 */
			    cprof = polyprof(spobj,cpol);
			/*
			 * Fill up density features for boundary candidates,
			 * and boundary features for density profile candidates.
			 */
  			    switch (cndp->candtype) {
			    case 1:
				xtra_b_feats(cndp,cprof,cpol,boundary,hist_diff);
				break;
			    case 2:
				xtra_p_feats(chcp,cndp,cpol,mpol,prof);
				break;
			    default:
				fprintf(stderr,"Help - extra feature switch\n");
				continue;
			    }
			/*
			 * measure features derived from the crossing profile
			 */
			    dcfeatures(chcp,cndp,cprof,prof);

			/*
			 * display
			 */
			    cndp->cpol = cpol;
			    chcp->obj = (struct object *) obj;
			    dispcand(obj,cpol,cndp,&f,1,1);
			/*
			 * Centromere training data
			 */
ASK:
			fprintf(stderr,"true centromere ?? (y or n) : ");
			if ((c=getc(fin)) == 'y') {
				dispcand(obj,cpol,cndp,&f,1,0);
				dispcand(obj,cpol,cndp,&f,2,1);
				cndp->Cctype = 1;
			}
			else if (c == 'n') {
				dispcand(obj,cpol,cndp,&f,1,0);
				cndp->Cctype = 0;
			}
			else {
				fprintf(stderr," ..... try again ..... : ");
				goto ASK;
			}
			while (getc(fin) != '\n')
				;
			    freeobj(cprof);
			}

			/*--------------------------------------------------------------*/
			/*
			 * for each candidate, find nearest candidate of
			 * alternative type and fill in "nearest" and
			 * "neardist" slots
			 */
			cndp = chcp->cands;
			for (i=0; i<chcp->ncands; i++, cndp++) {
			    cndp->nearest = -1;
			    cndp->neardist = 1000000;
			    cndp->u.sy.neardistfeat = 20;
			    if (cndp->candtype == -1)
			    	continue;
			    cndp1 = chcp->cands;
			    for (j=0; j<chcp->ncands; j++, cndp1++) {
			    	if (j == i || cndp1->candtype == -1 || cndp1->candtype == cndp->candtype)
			    		continue;
			    	dist = cndp->profpos - cndp1->profpos;
			    	if (dist < 0)
			    		dist = -dist;
			    	if (dist < cndp->neardist) {
			    		cndp->neardist = dist;
			    		cndp->u.sy.neardistfeat = dist;
			    		cndp->nearest = cndp1->candnum;
			    	}
			    }
			}
			/*--------------------------------------------------------------*/
			freeobj(spobj);
			freeobj(mpol);
			freeobj(p[0]);
			freeobj(p[1]);
			free(boundary);
			free(hist);
			free(hist_diff);

			chcp++;
		}
		/*----------------------------------------------------------------------*/
	}

	/*
	 * print summary about candidates
	 */
	dccandsummary(chc,numberused,pf);
	if (PRINT) {
		chcp = chc;
		for (i=0; i<numberused; i++)
			dcprint(chcp++,pf);
	}

	/*
	 * write to disc and exit
	 */
	if (outfile != NULL)
		dcwrite(chc,numberused,outfile);
	fclose(outfile);
	return(numberused);
}


dispcand(obj,cpol,cnd,f,c,i)
struct chromosome *obj;
struct object *cpol;
struct cand *cnd;
struct pframe *f;
{
	int korig,lorig;
	struct ivertex *vtx;
	struct ivector v;
	struct polygondomain *pdom = (struct polygondomain *)cpol->idom;
	double angle = - obj->plist->rangle/RADIANS;
	korig = (obj->idom->kol1 + obj->idom->lastkl)/2;
	lorig = (obj->idom->line1 + obj->idom->lastln)/2;
	vtx = pdom->vtx;
	v.type = 30;
	v.style = 1;	
	relrotate(vtx->vtX,vtx->vtY,korig,lorig,angle,1.0,1.333,&v.k1,&v.l1);
	vtx += (pdom->nvertices - 1);
	relrotate(vtx->vtX,vtx->vtY,korig,lorig,angle,1.0,1.333,&v.k2,&v.l2);
	intens(i);
	colour(c);
	picframe(&v,f);
}

/*
 * k, l are integer scaled by 8 !!
 */
relrotate(k,l,korig,lorig,theta,xscale,yscale,knew,lnew)
int *knew,*lnew;
double theta,xscale,yscale;
{
	double sin(),cos(),s1,s2,c1,c2;
	l -= lorig*8;
	k -= korig*8;
	s1 = sin(theta);
	c1 = cos(theta);
	/* modify theta */
	theta = atan(yscale*s1/(xscale*c1));
	s1 = sin(theta);
	c1 = cos(theta);
	s2 = xscale * c1;
	c2 = - xscale * s1;
	s1 *= yscale;
	c1 *= yscale;
	*lnew = lorig + iround((l*c1 + k*s1)/8.0);
	*knew = korig + iround((l*c2 + k*s2)/8.0);
}

iround(d)
double d;
{
	if (d >= 0.0)
		return((int)(d+0.4999));
	else
		return((int)(d-0.4999));
}


dcselcanddisp(chc,n,f)
struct chromcands *chc;
struct pframe *f;
{
	struct cand *cndp;
	int i,j,nc;
	for (i=0; i<n; i++,chc++) {
		cndp = chc->cands;
		nc = 0;
		for (j=0; j<chc->ncands; j++,cndp++) {
			if (cndp->candtype >= 1 && cndp->ctype == 1) {
				nc++;
				dispcand(chc->obj,cndp->cpol,cndp,f,6,0);
				dispcand(chc->obj,cndp->cpol,cndp,f,1,1);
			}
		}
		if (nc > 1) {
			picframe(convhull(chc->obj),f);
		}
	}
}

swab4(i)
int i;
{
	int b1,b2,b3,b4;
	b1 = (i>>24)&255;
	b2 = (i>>16)&255;
	b3 = (i>>8)&255;
	b4 = i&255;
	return((b4<<24)|(b3<<16)|(b2<<8)|b1);
}

/*
 * normalise grey-tables of set of objects jointly to lie in range 0-255
 */
normaliseset(objl,unum)
struct object **objl;
{
	int v,min,max,range, *g;
	int i,n;
	struct object *obj;
	struct iwspace iwsp;
	struct gwspace gwsp;
	min = 1;
	max = 0;
	for (n=0; n<unum; n++) {
		obj = objl[n];
		initgreyscan(obj,&iwsp,&gwsp);
		while (nextgreyinterval(&iwsp) == 0) {
			g = gwsp.grintptr;
			for (i=0; i<iwsp.colrmn; i++) {
				v = *g++;
				if (min > max)
					min = max = v;
				else if (v > max)
					max = v;
				else if (v < min)
					min = v;
			}
		}
	}
	range = max - min + (int) 1;
	for (n=0; n<unum; n++) {
		obj = objl[n];
		initgreyscan(obj,&iwsp,&gwsp);
		while (nextgreyinterval(&iwsp) == 0) {
			g = gwsp.grintptr;
			for (i=0; i<iwsp.colrmn; i++) {
				*g = ((int) 256) * (*g - min) / range;
				g++;
			}
		}
	}
}