/*
 * dcdemo.c		16-07-87		Jim Piper
 */

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

#define RADIANS 57.296

int Cbtype[MAXOBJS];
struct pframe f;

main(argc,argv)
char **argv;
{
	int 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;

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

	scale = -1;
	nonum = 1;
	setout(argv[1],'p');
	if (argc > 2) {
		cell = argv[2];
		infile = fopen(cell,"r");
		if (infile == NULL) {
			fprintf(stderr,"Can't open %s\n",cell);
			exit(1);
		}
	} else {
		cell = "<stdin";
		infile = stdin;
	}
	while((obj = readobj(infile)) != NULL) {
		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);
			}
		}
	}
fprintf(stderr,"objects read\n");
	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) {
		dinit();
		lwidth(1);
		colourlut();
#ifdef LOGO
		user(80,0);	/* for laser printer */
		clear();
		logo();
		sprintf(s,"%s: file %s",timestring,cell);
		colour(248);
		dispstring(s,24,24);
#endif LOGO
		colour(248);
		tsize(40,40);
		for (objnum=1; objnum<=number; objnum++)
			picframe(objlist[objnum-1],&f);
		if (!nonum)
			for (objnum=1; objnum<=number; objnum++)
				frameobjnum(objnum,objlist[objnum-1],&f);
		tsize(70,80);
		colour(2);
		dispstring("Candidate metacentric centromeres",32,4000);
		numberused = dcmain(objlist,number,chc,cnd);
		update();
		fprintf(stderr,"pausing for 5 seconds - if you wish to WAIT\n");
		sleep(5);
		colour(1);
		dispstring("Selected centromeres",2600,4000);
		classmain(chc,cnd,numberused);
		finish();
	}
}


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


min(i,j)
{
	return ( i<j? i: j);
}
/*
 * dccf.c		Jim Piper		21-01-87
 *
 * find and display dicentric centromere candidates
 *
 * dccf [-d <stoxxxx>] [-p] [-q <QUANTUM>] [-P] [-B] [-o outputfile] [file]
 *
 * MODIFICATIONS
 * -------------
 * 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)
 */


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

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


dcmain(objl,totalnumber,chc,cnd)
struct chromcands *chc;
struct cand *cnd;
struct chromosome **objl;
{
	struct object *spobj, *mpol;
	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;
	int i, j, number, numberused, sernum;
	int QUANTUM, PROF_CANDS, BOUND_CANDS;
	char *cell,s[256],*timestring,*outfilename;
	FILE *infile, *pf, *outfile;
	long daytime;
	int DISPLAY = 0;
	int PRINT = 0;
	struct polygondomain *boundary;
	struct histogramdomain *hist, *hist_diff, *bound_curv();

	PROF_CANDS = BOUND_CANDS = 0;
	pf = stdout;

	/*
	 * defaults for QUANTUM
	 */
	QUANTUM = 12;
	PROF_CANDS++;
	DISPLAY++;
	if (DISPLAY != 0) {
		colourlut();
	}
	numberused = number = sernum = 0;
	chcp = chc;
	cndp = cnd;
	while (sernum < totalnumber) {
fprintf(stderr,"obj serial %d\n",sernum);
		obj = objl[sernum++];
		/*
		 * give it a propertylist and serial number if necessary
		 */
		if (obj->plist == NULL) {
			obj->plist = makechromplist();
			++number;
		}
		else
#ifdef NUMBWORKS
			number = obj->plist->number;
#else
			{
			++number;
		}
#endif
		/*
		 * 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 analyses objects believed to be chromosome (i.e.
		 * NOT composite, noise, overlap, etc..) by any previous pass.
		 */
		else if (obj->plist->otype <= 1) {
			Cbtype[numberused] = obj->plist->Cbtype;
			numberused++;

			/*
			 * rotate to vertical, find axis,
			 * find density and 1st moment profile
			 */
			spobj = dcvertical(obj,FIP);
			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, 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;
			chcp->area = plist->area;
			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);
			}
			/*
			 * display candidates if required
			if (DISPLAY >= 1)
 			    display1(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:
			        cpol = centpoly(mpol,cndp->profpos);
			        break;
			    case 2:
			        cpol = boundpoly(boundary, cndp, spobj);
				if (cpol == NULL) continue;
			        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(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);
			 */
			    cndp->u.af[15] = (struct object *) cpol;
			    cndp->u.af[16] = (struct object *) obj;
			    dispcand(obj,cpol,cndp,&f,2,1);
			    freeobj(cprof);
			}

			chcp++;

			freeobj(spobj);
			freeobj(mpol);
			freeobj(p[0]);
			freeobj(p[1]);
			free(boundary);
			free(hist);
			free(hist_diff);
		}
	}

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

/*
 * dcclass.c		Jim Piper		21-01-87
 *
 * perform box classification on dicentric candidate/feature data
 *
 * dcclass [-f n] [-p] [-m <PC1> <PC2>] [file]
 *
 * -f n		selects whether classification on type 1 candidates (n=1),
 *		type 2 candidates (n=2), or some other candidate selection
 *		method (to be determined).
 */




/* which features (numbering from one, as in dicen.h) */
int nfeat = 9;
int selfeat[MAXFEATURES] = {1,3,5,6,8,9,13,14,15};
int selllim[MAXFEATURES] = {1,1,1,0,0,0, 1, 1, 0};


classmain(chc,cnd,numberused)
struct chromcands *chc;
struct cand *cnd;
{
	struct object *spobj, *mpol, *cvh, *rect;
	struct chromosome *obj;
	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;
	struct cand *cndp;
	int i, j, number;
	int PC1, PC2;
	char *cell,s[256];
	FILE *infile, *pf;
	int DISPLAY = 0;
	int PRINT = 0;
	int CANDSEL = 1;

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

	pf = stdout;
	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");
	}

	/*
	 * normalise features
	 */
	dcnormalise(chc,numberused);

	/*
	 * do the classification
	 */
	dcboxclass(1,0,chc,numberused,nfeat,selfeat,selllim,PC1,PC2,PRINT,pf);

	/*
	 * display the found true and false positives
	 * and the false negatives
	 */
	dcselcanddisp(chc,numberused,&f);
}

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->ctype == 1) {
				nc++;
				dispcand(cndp->u.af[16],cndp->u.af[15],cndp,f,2,0);
				dispcand(cndp->u.af[16],cndp->u.af[15],cndp,f,1,1);
			}
		}
		if (nc > 1) {
			cndp--;
			picframe(convhull(cndp->u.af[16]),f);
		}
	}
}