/*
 * xprobe.c
 *
 * Test bed for X-centromere labelling detection.
 *
 * xprobe [-v] [-S] [-D] [-d] [-o featfile] [file]
 *
 * -v - don't rotate objects to vertical first
 * -S selects objects according to criterion in "xdecide()"
 * -D displays each object with feature values
 * -d displays selected objects with feature values
 * -o write 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 */


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 deriv2,xcentromere,cia,pia,mediana,maxobjs,areao,difpos;
	int areal[200];
	double os9frig;
	char s[200], *celln;
	FILE *in,*out;
	int intsort();
	int VERT = 0;
	int SELECT = 0;
	int DISPLAY = 0;
	int RESULTS = 0;

	setbuf(stderr,NULL);

	while (argc > 1 && argv[1][0] == '-') {
		switch (argv[1][1]) {
			case 'v':
				VERT = 1;
				break;
			case 'S':
				SELECT = 1;
				DISPLAY = 0;
				break;
			case 'D':
				DISPLAY = 1;
				SELECT = 0;
				break;
			case 'd':
				DISPLAY = 2;
				SELECT = 0;
				break;
			case 'o':
				RESULTS++;
				out = fopen(argv[2],"w");
				if (out == NULL) {
					fprintf(stderr,"Cannot open %s for results\n",argv[2]);
					exit(1);
				}
				argc--;
				argv++;
				break;
		}
		argc--;
		argv++;
	}
	if (argc > 1) {
		in = fopen(argv[1],"r");
		celln = argv[1];
	}
	else {
		in = stdin;
		celln = "stdin";
	}
	if (DISPLAY) {
		dinit();
		user(40,0);
		tsize(50,50);
	}
	numarea = number = 0;
	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
/*		  && (obj->plist->history[0]&020)		*/
		  && parl[obj->plist->pnumber]
		  && parl[obj->plist->pnumber]->plist->otype >= 3)
			continue;
		objl[number] = obj;
		number++;
	}
	qsort(areal,numarea,sizeof(int),intsort);
	mediana = areal[numarea/2];
	maxobjs = number;
	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) {
		f.type = 60;
		f.scale = 4;
		f.ix = 8;
		f.iy = 8;
		f.dx = 2560;
		f.dy = 1024;
		if (VERT)
			spobj = (struct object *) obj;
		else {
			cvh = fipconvhull(obj);
			rect = minwrect(cvh);

			/*
			 * rotate to vertical
			 */
			os9frig = ((struct frect *)rect->idom)->rangle;
			spobj = spinsqueeze(obj,os9frig,xscale,yscale);
		}
		f.ox = spobj->idom->kol1;
		f.oy = spobj->idom->line1;
		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);
		xcentromere = xdecide(xflp);
		/*
		 * write features if required
		 */
		if (RESULTS)
			writexfeat(out,xflp);

		/*
		 * display things if object selected
		 */
		if (xcentromere || DISPLAY == 2) {
		if (SELECT)
			writeobj(stdout,spobj);
		else if (DISPLAY >= 1) {
		histonormalise(prof);
		clear();
		f.dx -= 768;
		h = (struct histogramdomain *) prof->idom;
		ydiff = h->npoints - (spobj->idom->lastln-spobj->idom->line1+1);
		ydiff *= f.scale * f.iy / 2;
		f.dy -= ydiff;
		moveto(f.dx,f.dy);
		moveby(0,-100);
#ifdef TEXT
		text("mean");
#endif TEXT
		moveby(0,-70);
#ifdef TEXT
		text("density");
#endif TEXT
		moveby(0,-70);
#ifdef TEXT
		text("profile");*/
#endif TEXT
		moveby(0,240);
		lineby(0,h->npoints*f.scale*f.iy);
		f.ix = 2;
		picframe(prof,&f);

		prof = p[1];
		f.dx -= 769;
		cpos++;
		moveto(64,64);
#ifdef TEXT
		text(s);
#endif TEXT
		moveto(f.dx,f.dy+cpos*f.scale*f.iy);
		moveby(-600,0);
#ifdef TEXT
		text("centromere");
#endif TEXT
		moveby(600,0);
		lineby(1450,0);

		histonormalise(prof);
		h = (struct histogramdomain *) prof->idom;
		moveto(f.dx,f.dy);
		moveby(0,-100);
#ifdef TEXT
		text("shape");
#endif TEXT
		moveby(0,-70);
#ifdef TEXT
		text("profile");
#endif TEXT
		moveby(0,170);
		lineby(0,h->npoints*f.scale*f.iy);
		f.ix = 2;
		picframe(prof,&f);


		f.dy += ydiff;
		f.dx += 1536;
		f.ix = 8;
		picframe(spobj,&f);
		picframe(mpol,&f);

		if (cpos >= 0) {
			f.ix = f.iy = 1;
			f.ox *= 8;
			f.oy *= 8;
			f.dx += 4*f.scale;
			poly = centpoly(mpol,cpos);
			pdom = (struct polygondomain *) poly->idom;
			poly->idom = (struct intervaldomain *) uspaxis8(pdom);
			picframe(poly,&f);
			free(pdom);
		}
		sleep(2);
	}
		}
		freeobj(p[0]);
		freeobj(p[1]);
		freeobj(cvh);
		freeobj(mpol);
		freeobj(rect);
		freeobj(spobj);
		}
		freeobj(obj);
	}
	if (DISPLAY)
		finish();
}