/*
 * nxp.c
 *
 *******************************************************************
 * WARNING - Uses chromosome property list in non-standard fashion *
 *******************************************************************
 * Test bed for X-centromere labelling detection.
 * Measures features and runs box classifier.
 * Selected objects are labelled with cell-name
 * and output; feature values are also output
 * in "vcv" format for training the quadratic
 * classifier.
 *
 * xprobe [-d] [-m] [-f] [-o selobjfile] [-a featfile] [file]
 *
 * -d display on a Sun workstation screen.
 * -m known to be male
 * -f known to be female
 * -o append selected objects to "selobjfile" (only
 *    if named file provided for input).
 * -a append "vcv" feature values for each object to "featfile"
 */

#include <stdio.h>
#include <wstruct.h>
#include <chromanal.h>
#include <otypes.h>
#include "xprobe.h"

int wh[2] = {-2,2};	/* max density, second moment */

int MODE = 1;
float xscale = 1.0;
float yscale = 0.75;
static struct chromosome *parl[1000] = {0};

main(argc,argv)
char **argv;
{
	struct object *spobj, *mpol, *cvh, *rect;
	struct chromosome *obj, *objl[1000], *readobj();
	struct object *fipconvhull(), *spinsqueeze(), *minwrect();
	struct object *poly, *centpoly(), *chromaxis2();
	struct object *prof, *p[3], *nprofile();
	struct object *polyprof(), *polyp;
	struct histogramdomain *h;
	struct polygondomain *pdom, *uspaxis8();
	struct pframe f;
	struct xfeats xfl[200],*xflp;
	int i, j, l, numarea, number, cpos, maxpos, ca, cl, ydiff, peak;
	int c,deriv2,xcentromere,cia,pia,mediana,maxobjs,areao,difpos;
	int areal[200];
	double os9frig;
	char s[200], *celln;
	unsigned char uc;
	FILE *in,*fout,*oout, *tin;
	int intsort();
	int FRESULTS = 0, ORESULTS = 0, DISPLAY = 0;
	int SEX = SEX_UNKNOWN;

	setbuf(stderr,NULL);

	while (argc > 1 && argv[1][0] == '-') {
		switch (argv[1][1]) {
			case 'o':
				ORESULTS++;
				oout = fopen(argv[2],"a");
				if (oout == NULL) {
					fprintf(stderr,"Cannot open %s for appending box-selected objects\n",argv[2]);
					exit(1);
				}
				argc--;
				argv++;
				break;
			case 'a':
				FRESULTS++;
				fout = fopen(argv[2],"a");
				if (fout == NULL) {
					fprintf(stderr,"Cannot open %s for appending box-selected data\n",argv[2]);
					exit(1);
				}
				argc--;
				argv++;
				break;
			case 'd':
				DISPLAY++;
				tin = fopen("/dev/tty","r");
				setbuf(tin,NULL);
				break;
			case 'm':
				SEX = SEX_MALE;
				break;
			case 'f':
				SEX = SEX_FEMALE;
				break;
			default:
				fprintf(stderr,"%s not understood\n",argv[1]);
				exit(1);
		}
		argc--;
		argv++;
	}
	if (argc > 1) {
		in = fopen(argv[1],"r");
		celln = argv[1];
	}
	else {
		in = stdin;
		celln = "stdin";
		if (ORESULTS) {
			fprintf(stderr,"no named input file, sono output objects\n");
			ORESULTS = 0;
		}
	}
	numarea = number = 0;

	/*
	 * read objects, reject if unsuitable, find areas
	 */
	while ((obj = readobj(in)) != NULL) {
		swabplist(obj->plist);
		parl[obj->plist->number] = obj;
		/* reject non-chromosome */
		if (obj->plist->otype > 1)
			continue;
		areal[numarea++] = obj->plist->area = area(obj);
		/* reject manual overlap split */
		if (obj->plist->pnumber != 0 && (obj->plist->history[0]&2))
			continue;
		/* reject auto-split if overlap */
		if (obj->plist->pnumber != 0
		  && parl[obj->plist->pnumber]
		  && parl[obj->plist->pnumber]->plist->otype >= 3)
			continue;
		objl[number] = obj;
		number++;
	}
	if (DISPLAY) {
		setout("stosun -p270,0 -d1 -r870,870",'p');
		dinit();
		defaultlut();
		lut_overlay();
		lwidth(3);
		xpdispinit(objl,number,&f);
		for (i=0; i<number; i++)
			picframe(objl[i],&f);
		/*
		 * Allow display to catch up !
		 */
		update();
		sleep(4);
	}

	/*
	 * rank area list by size order and find median for normalisation
	 */
	qsort(areal,numarea,sizeof(int),intsort);
	mediana = areal[numarea/2];
	maxobjs = number;

	/*
	 * measure features, fill up the "xfeats" structure
	 */
	for (number=0; number<maxobjs; number++) {
		obj = objl[number];
		xflp = xfl+number;
		xflp->number = obj->plist->number;
		xflp->Cbtype = obj->plist->Cbtype;
		xflp->relarea = areao = 100*area(obj)/mediana;
		if (areao < 40 || areao > 180)
			continue;
		if (DISPLAY)
			xpdisp(obj,&f,1,1);
		cvh = fipconvhull(obj);
		rect = minwrect(cvh);

		/*
		 * rotate to vertical
		 */
		os9frig = ((struct frect *)rect->idom)->rangle;
		spobj = spinsqueeze(obj,os9frig,xscale,yscale);

		mpol = chromaxis2(spobj);
	
		/*
		 * moment profile computation and smoothing
		 */
		multiprofiles(spobj,mpol,p,wh,2);
		shorten(mpol,p[0],p[1]);
		whistosmooth(p[1],p[0],3); 

		/*
		 * locate centromere in un-mirrored profile with
		 * quantum parameter set to be 10% of peak value.
		 */
		prof = p[1];
		cpos = cmerelocation(prof,10,0);

		/*
		 * analyse density profile for maximum
		 * and other properties
		 */
		prof = p[0];
		maxpos = maxlocation(prof,&peak);
		if (cpos >= 0) {
			difpos = maxpos-cpos;
			cia = ((struct fvertex *)((struct polygondomain *)(mpol->idom))->vtx)[cpos].vtY;
			cia = cindexa(spobj,cia);
		}
		else {
			difpos = 0;
			cia = -1;
		}
		pia = ((struct fvertex *)((struct polygondomain *)(mpol->idom))->vtx)[maxpos].vtY;
		pia = cindexa(spobj,pia);
		if (cia > 50) {
			cia = 100 - cia;
			pia = 100 - pia;
			difpos = -difpos;
		}
		else if (cia == -1 && pia > 50)
			pia = 100 - pia;
		deriv2 = secondderiv(prof,maxpos,3);

		xflp->cindexa = cia;
		xflp->pindexa = pia;
		xflp->difpos = difpos;
		xflp->profpeak = peak;
		xflp->deriv2 = deriv2;
		greyhistofeats(obj,xflp);

		freeobj(p[0]);
		freeobj(p[1]);
		freeobj(cvh);
		freeobj(mpol);
		freeobj(rect);
		freeobj(spobj);
		if (DISPLAY)
			xpdisp(obj,&f,1,0);
	}

	/*
	 * Box classify
	 */
	xpboxclass(xfl,maxobjs);

	/*
	 * Display selected chromosomes,
	 * write features if required
	 */

	if (DISPLAY) {
		tsize(80,240);
		intens(1);
		colour(248);
		moveto(0,3960);
		pixbegin(8,1,1,4095,135);
		pixel(80);
		pixend();
		colour(2);
		dispstring("Believed to be X",100,4000);
		colour(1);
		dispstring("Believed NOT to be X",1500,4000);
		colour(4);
		dispstring("Possibly X",3200,4000);
	}
	xflp = xfl;
	for (i=0; i<number; i++,xflp++) {
		if (xflp->btype == 1) {
			obj = parl[xflp->number];
			if (DISPLAY) {
				if (xflp->xrank == 1)
					c = 2;
				else if (xflp->xrank >= 3)
					c = 1;
				else
					switch (SEX) {
					case SEX_UNKNOWN:
						c = 4;
						break;
					case SEX_MALE:
						c = 1;
						break;
					case SEX_FEMALE:
						c = 2;
						break;
					}
				xpdisp(obj,&f,c,1);
			}
			if (FRESULTS)
				writexpvcvfeat(fout,xflp);
			if (ORESULTS) {
				strcpy(&(obj->plist->wdd[0]),argv[1]);
				putxpfeats(obj,xflp);
				writeobj(oout,obj);
			}
		}
	}
	if (DISPLAY) {
		fprintf(stderr,"type return to continue\n");
		while(getc(tin) != '\n')
			;
		finish();
	}
}

xpboxclass(xf,number)
struct xfeats *xf;
int number;
{
	int i,j;
	struct xfeats *xfp;
	int xfpcomp(),vcvcomp(),rank,max2d,xrank;
	/*
	 * massage features as necessary
	 */
	xfp = xf;
	max2d = 0;
	for (i=0; i<number; i++,xfp++)
		if (max2d < xfp->deriv2)
			max2d = xfp->deriv2;
	xfp = xf;
	for (i=0; i<number; i++,xfp++)
		xfp->deriv2 = 100 * xfp->deriv2 / max2d;

	/*
	 * sort by density profile peak, code rank.
	 */
	qsort(xf,number,sizeof(struct xfeats),xfpcomp);

	/*
	 * Classification:
	 * (1) First MAXRANK only of those which satisfy wide pindexa bounds
	 * (2) relarea limits
	 * (3) cindexa limits, or pindexa limits if cindexa == -1
	 * (4) difpos limit if cindexa != -1
	 */
	xfp = xf;
	for (i=1, rank=1, xrank=1; xrank<=XRANK && i<=number && rank<=MAXRANK; i++,rank++,xfp++) {
		xfp->btype = 2;
		if (xfp->pindexa < PILB || xfp->pindexa > PIUB) {
			rank--;
			continue;
		}
		if (xfp->relarea < AREAMIN || xfp->relarea > AREAMAX ||
		    xfp->deriv2 < D2MIN)
			continue;
		if (xfp->cindexa != -1) {
			if (xfp->cindexa < CIMIN)
				continue;
			if (xfp->cindexa > CIMAX) {
				if (xfp->pindexa > PIMAX || xfp->difpos < 0 ||
				    xfp->difpos > DIFF)
					continue;
			}
			if (xfp->difpos < -DIFF || xfp->difpos > DIFF)
				continue;
		}
		else if (xfp->pindexa < PIMIN || xfp->pindexa > PIMAX)
			continue;
		/*
		 * Potential probed X found
		 */
		xfp->xrank = xrank++;
		xfp->btype = 1;
		xfp->rank = rank;
	}
}

int xfpcomp(x1,x2)
register struct xfeats *x1,*x2;
{
	return(x2->profpeak - x1->profpeak);
}