/*
 * dccf.c		Jim Piper		21-01-87
 *
 * find and display dicentric centromere candidates
 *
 * dccf [-d stoxxx] [-p] [-q <QUANTUM>] [-P] [-B] [-o ofile] [-M] [-S] [file]
 *
 * FLAGS:
 * -----
 * -d stoxxx	display all srts of things on stoxxx
 * -p		turn on printing
 * -q quantum	quantum (factor) for minimum detection in profile
 * -P		extract profile generated candidates
 * -B		extract boundary generated candidates
 * -o ofile	send feature measurements to ofile
 * -M		refine position of boundary candidate crossing profile
 *		if there is an nearby profile local minimum.
 * -S		straighten the chromosome
 *
 * MODIFICATIONS
 * -------------
 * 19-06-90	jimp	"Type 3" candidate generated for every chromosome,
 *			not just those with single type 1 candidate.
 * 07-06-90	jimp	Added "type 3" candidates in switch -Y: if just
 *			one type 1 candidate, then generate a type 3
 *			candidate somewhere between the type 1 candidate
 *			and the further extremity of the chromosome.
 * 11-04-88	jimp	Too many candidates ??
 * 24-03-88	jimp	check for connectedness, remove small disconnected bits
 * 23-03-88	simon	Added optional straightening of chromosome
 * 23-03-88	jimp	Added porigin for parental origin of derived objects
 * 01-03-88	jimp	Perturb position of boundary candidate crossing
 *			polygon if close to profile minimum (flag -M).
 * 29-02-88	jimp	Change call parameters to "display1"
 * 25-02-88	jimp	Added chcp->axislen, cndp->u.sy.neardistfeat.
 * 20-01-88	jimp	Measure distance between nearest candidates of
 *			types 1 and 2.
 * 25-02-87	JP	"prof_cands" split into "prof_basics" and "prof_cands".
 *			This allows optional runs with only boundary candidates
 *			(argument -B) or profile candidates (argument -P).
 * 10-02-87	JP	Some comments, renamed "dcmerecands" to be "prof_cands",
 *			colour look-up table for displays.  Boundary display
 *			moved from "display2" to "display1" (so only drawn once)
 */

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

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

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


main(argc,argv)
char **argv;
{
	struct object *spobj, *mpol, *oli[20];
	struct chromosome *obj, *readobj();
	struct chromplist *plist, *makechromplist();
	struct object *dcvertical();
	struct object *cpol, *centpoly(), *boundpoly(), *chromaxis2();
	struct object *prof, *p[2], *cprof, *polyprof(), *straighten();
	struct polygondomain *poly, *uspaxis8();
	struct chromcands *chcp,chc[MAXOBJS];
	struct cand *cndp,*cndp1,cnd[MAXCANDS];
	struct pframe f;
	int i, j, dist, number, numberused, nlab, ao, amin, imin;
	int QUANTUM, PROF_CANDS, BOUND_CANDS, totalcands;
	char *cell,s[256],*timestring,*outfilename;
	FILE *infile, *pf, *outfile, *fin;
	long daytime;
	int DISPLAY, PRINT, BPOLMOD, STRAIGHT, TYPE3;
	struct polygondomain *boundary;
	struct histogramdomain *hist, *hist_diff, *bound_curv();

	TYPE3 = DISPLAY = PRINT = BPOLMOD = PROF_CANDS = BOUND_CANDS = STRAIGHT = 0;
	totalcands = 0;
#ifdef VAX
	daytime = time();
	timestring = ctime(&daytime);
	timestring[24] = ' ';
#endif VAX
	outfile = NULL;
	outfilename = "<no output>";
	pf = stdout;

	/*
	 * defaults for QUANTUM
	 */
	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 '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 'S':
			STRAIGHT++;
			break;
		case 'Y':
			TYPE3++;
			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;
	}
	if (PROF_CANDS == 0 && BOUND_CANDS == 0)
		PROF_CANDS = BOUND_CANDS = 1;
	if (DISPLAY != 0) {
		dinit();
		lwidth(3);
		user(40,0);
		colourlut();
		fin = fopen("/dev/tty","r");
	}
	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 = 0;
	chcp = chc;
	cndp = cnd;
	while ((obj = readobj(infile)) != NULL) {
		if (obj->plist != NULL && (obj->plist->number <= 0 || obj->plist->number > MAXOBJS))
			swabplist(obj->plist);
		/*
		 * 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 totcands %d\n",number,obj->plist->area,obj->plist->mass,(int)(cndp-cnd));
			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);
			if (STRAIGHT) {
				spobj = straighten(spobj, mpol);
				freeobj(mpol);
				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);
#ifdef FIPTEST
printhist(stderr,hist);
#endif FIPTEST
			
			/*
			 * identify composite/overlap objects
			 */
		/*	bound_comp(obj, spobj, boundary, hist, hist_diff);	*/
			
			/*
			 * fill up top of chromcands structure
			 */
			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);
			}
			if (TYPE3)
				type3_cand(prof,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 */
			}
			/*
			 * display candidates if required
			 */
			if (DISPLAY >= 1)
 			    display1(BOUND_CANDS,PROF_CANDS,cell,number,spobj,mpol,prof,boundary,&f,chcp,timestring,SCALE);

			/*--------------------------------------------------------------*/
			/*
			 * 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:
			    case 3:
			        cpol = centpoly(mpol,cndp->profpos);
			        break;
			    case 2:
			        cpol = boundpoly(boundary, cndp, spobj);
				if (cpol == NULL) 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:
			    case 3:
				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 crossing line and profile if required
			 */
			    if (DISPLAY >= 1)
			    	display2(cndp,cpol,cprof, mpol,&f,i,SCALE);
			    freeobj(cpol);
			    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++;
		}
		/*----------------------------------------------------------------------*/
		freeobj(obj);
		if (DISPLAY) {
			update();
			fprintf(stderr,"Type return to continue\n");
			getc(fin);
		}
	}

	/*
	 * 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);
	if (DISPLAY >= 1)
		finish();
}




struct object *straighten(obj, poly)
struct object *obj, *poly;
{
	struct object *strobj, *makemain(), *threshold();
	struct intervaldomain *idom,*makedomain();
	struct valuetable *vtbl,*makevaluetb();
	struct polygondomain *p, *ply, *makepolydmn();
	struct iwspace iwsp;
	struct gwspace gwsp;
	int strwidth,halfwidth,sline,mlines,i,*g;
	float x,y,sintheta,costheta;

	p = (struct polygondomain *) poly->idom;
	if (poly->type != 10 || p->type != 2) {
		fprintf(stderr,"not a type 2 polygon\n");
		exit(91);
	}
	mlines = p->nvertices;

	/*
	 * approximate mean width of object
	 * this is used as half-width of straightened object
	 */
	halfwidth = area(obj)/mlines;
	if (halfwidth > 10)
		halfwidth = 10;
	strwidth = 1 + 2*halfwidth;

	/*
	 * make a rectangular object with an empty grey table
	 */
	idom = makedomain(2,1,mlines,1,strwidth);
	vtbl = makevaluetb(1,1,mlines,0,NULL);
	strobj = makemain(1,idom,vtbl,NULL,obj);
	vtbl->original = strobj;
/* have to make grey-table one bigger than might be thought necessary
as somewhere there is a loss of sync !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
	g = (int *) Calloc(mlines*strwidth+1,sizeof(int));
	for (i=1; i<=mlines; i++) {
		makevalueline(vtbl,i,1,strwidth,g);
		g += strwidth;
	}

	/*
	 * start transformation
	 */
	initgreyscan(strobj,&iwsp,&gwsp);
	for (sline=1; sline<=mlines; sline++) {
		/*
	 	* get next column, line (x,y) in unit-spaced steps
	 	* along smoothed mid-points polygon, and local orientation.
	 	* sline is line number in straightened object
	 	*/
		setnextpoint(p,&x,&y,&sintheta,&costheta,sline);

		/*
		 *  rotate and fill rotated interval centred at x,y.
		 */
		nextgreyinterval(&iwsp);
		g = gwsp.grintptr;
		rotate(obj,x,y,sintheta,costheta,halfwidth,g);
	}
	/*
	 * threshold to remove unset points
	 */
	strobj = threshold(strobj,1);
	return(strobj);
}

#ifdef POLYTELL
polytell(p)
struct polygondomain *p;
{
	int i;
	struct ivertex *vtx = p->vtx;
	printf("polytell: %d vertices (max %d)\n",p->nvertices,p->maxvertices);
	for (i=0; i<p->nvertices; i++,vtx++)
	printf("%5d %5d\n",vtx->vtX,vtx->vtY);
	fflush(stdout);
}
#endif POLYTELL
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);
}



#ifdef FIPTEST

printhist(pf,hist)
FILE *pf;
struct histogramdomain *hist;
{
	int i,j,k,*h;
	if (hist->type != 1) {
		fprintf(pf,"Bound_curv histogram not type 1\n");
		return;
	}
	h = hist->hv;
	for (i=0; i<hist->npoints; i++,h++) {
		fprintf(pf,"%4d ",*h);
		if (i>0 && (i+1)%8 == 0)
			fprintf(pf,"\n");
	}
	fprintf(pf,"\n");
}

#endif FIPTEST