/*
 * xprofs.c
 *
 * Test bed for centromere-labelled fragile X detection
 *
 * xprofs [-D] [file]
 *
 * -D displays each object with feature values
 */

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

/*
 * some display control definitions
 */
#define SCALE	4	/* true DIGS scale divided by 4 */
#define XORIG	1200
#define YORIG	2800
#define XN	3
#define YN	3
#define XINC	1300
#define YINC	1300

#define GPROF 2		/* mean grey */
/* max grey, mean grey, total grey, 1st mom, 2nd mom */
static int wh[] = {-2,-1,0,1,2};

float xscale = 1.0;
float yscale = 0.75;

main(argc,argv)
char **argv;
{
	struct object *spobj, *mpol, *cvh, *rect;
	struct chromosome *obj, *readobj();
	struct object *fipconvhull(), *spinsqueeze(), *minwrect();
	struct object *poly, *centpoly(), *chromaxis2();
	struct object *prof, *p[5], *nprofile();
	struct object *polyprof(), *polyp;
	struct object *dilation(), *erosion(), *d1,*d2,*e1,*e2;
	struct object *skeleton(), *skobj;
	struct histogramdomain *h;
	struct polygondomain *pdom, *uspaxis8();
	struct pframe f;
	struct valuetable *newvaluetb();
	int i, j, l, fiy, numarea, number, cpos, invert, maxpos, ca, cl, ydiff, peak;
	int deriv2,xcentromere,cia,pia,mediana,maxobjs,areao,difpos,dcount;
	int end1,end2;
	double os9frig;
	char s[200], *celln;
	FILE *in,*out;
	int intsort();
	int DISPLAY = 1;

	while (argc > 1 && argv[1][0] == '-') {
		switch (argv[1][1]) {
			case 'D':
				DISPLAY = 0;
				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(16*SCALE,16*SCALE);
	}
	numarea = number = 0;
	while ((obj = readobj(in)) != NULL) {
		swabplist(obj->plist);
		/* reject non-chromosome */
		if (obj->plist->otype > 1)
			continue;
#ifdef NOOVS
		/* 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]&8))
			continue;
#endif NOOVS

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

		/*
		 * find axis
		 */
#ifdef CLEAN
		d1 = dilation(spobj);
		e1 = erosion(d1);
		e2 = erosion(e1);
		d2 = dilation(e2);
		mpol = chromaxis2(d2);
#else CLEAN
		mpol = chromaxis2(spobj);
#endif CLEAN	
		/*
		 * moment profile computation and smoothing
		 */
		multiprofiles(spobj,mpol,p,wh,5);
		multishorten(mpol,p[2],p,5);
		histosmooth(p[2],2);
		whistosmooth(p[3],p[2],3); 
		whistosmooth(p[4],p[2],3); 
		maxpos = maxlocation(p[0],&peak);

		if (DISPLAY >= 1) {
			if (dcount % (XN*YN) == 0) {
				if (dcount > 0) {
					sleep(5);
					clear();
				}
				moveto(8,4040);
				text("2nd moment, max grey");
				moveto(3800,4040);
				text("page");
				dispnum(1+dcount/(XN*YN),4000,4040);
			}
			f.type = 60;
			f.scale = SCALE;
			f.ix = 4;
			fiy = f.iy = 4;
			f.dx = XORIG + XINC*(dcount % XN);
			f.dy = YORIG - YINC*((dcount/XN) % YN);
			dcount++;
			f.ox = spobj->idom->kol1;
			f.oy = spobj->idom->line1;

			h = (struct histogramdomain *) p[2]->idom;
			invert = -1;
			if (maxpos*2 < h->npoints) {
				invert = 1;
				f.iy = -4;
				f.oy = spobj->idom->lastln;
			}

			f.dx -= 48*SCALE;
			ydiff = h->npoints - (spobj->idom->lastln-spobj->idom->line1+1);
			ydiff *= f.scale * fiy / 2;
			f.dy -= ydiff;
			moveto(f.dx,f.dy);
			if (obj->plist->Cbtype&FRAXA) {
				moveby(0,-24*SCALE);
				text("fra(X)");
				moveby(0,24*SCALE);
			}
			else if (obj->plist->Cbtype&XPROBE) {
				moveby(0,-24*SCALE);
				text("X");
				moveby(0,24*SCALE);
			}

#ifdef ALLPROFS
			for (i=0; i<5; i++) {
				f.dx -= 50*SCALE;
				prof = p[i];
				histonormalise(prof);
				f.ix = 1;
				picframe(prof,&f);
			}
			f.dx += 298*SCALE;
#endif ALLPROFS

			f.dx -= 160*SCALE;
			moveto(f.dx,f.dy);
			lineby(0,h->npoints*fiy*f.scale);
			moveto(f.dx,f.dy);
			moveby(8,((h->npoints+9)/10) * fiy*f.scale);
			lineby(-16,0);
			prof = p[4];
			histonormalise(prof);
			f.ix = 1;
			picframe(prof,&f);
#ifndef NOMIRROR
			submirrorprof(prof);
			picframe(prof,&f);
#endif NOMIRROR
			end1 = invert*profend(prof);
			dispnum(end1,f.dx,f.dy-24*SCALE);

			f.dx -= 60*SCALE;
			moveto(f.dx,f.dy);
			lineby(0,h->npoints*fiy*f.scale);
			moveto(f.dx,f.dy);
			moveby(8,((h->npoints+9)/10) * fiy*f.scale);
			lineby(-16,0);
			prof = p[GPROF];
			histonormalise(prof);
#ifndef NOMIRROR
			submirrorprof(prof);
#endif NOMIRROR
			f.ix = 1;
			picframe(prof,&f);
			end2 = invert*profend(prof);
			dispnum(end2,f.dx,f.dy-24*SCALE);


			f.dx += 60*SCALE;
			f.dx += 160*SCALE;
			dispclass(obj,f.dx-40*SCALE,f.dy-24*SCALE,end1,end2);

			f.dy += ydiff;
			f.ix = 4;
			picframe(spobj,&f);
/*
			skobj = skeleton(spobj,2);
			skobj->vdom = newvaluetb(skobj,1,0);
			setval(skobj,255);
			picframe(skobj,&f);
			freeobj(skobj);
*/
			picframe(mpol,&f);
		}
		freeobj(p[0]);
		freeobj(p[1]);
		freeobj(p[2]);
		freeobj(p[3]);
		freeobj(p[4]);
		freeobj(cvh);
		freeobj(mpol);
		freeobj(rect);
		freeobj(spobj);
		freeobj(obj);
#ifdef CLEAN	
		freeobj(d1); freeobj(d2); freeobj(e1); freeobj(e2);
#endif CLEAN	
	}
	if (DISPLAY)
		finish();
}


/*
 * compute profile of average values of profile
 * subtracted from its mirror image
 */
submirrorprof(p)
struct object *p;
{
	struct histogramdomain *hdom;
	int *ev,*hv,n,v;
	if (p->type != 13) {
		fprintf(stderr,"SUBMIRRORPROF: not a histogram\n");
		exit(1);
	}
	hdom = (struct histogramdomain *)p->idom;
	hv = hdom->hv;
	n = hdom->npoints;
	ev = hv + n - 1;
	n = (n+1)/2;
	for(; n>0; n--,hv++,ev--) {
		v = (*hv - *ev);
		*hv = v;
		*ev = -v;
	}
}


/*
 * compute mean of first 10% of profile values
 */
profend(p)
struct object *p;
{
	struct histogramdomain *hdom;
	int *hv,i,n,v;
	if (p->type != 13) {
		fprintf(stderr,"PROFEND: not a histogram\n");
		exit(1);
	}
	hdom = (struct histogramdomain *)p->idom;
	hv = hdom->hv;
	n = hdom->npoints;
	v = 0;
	i = n = (n+9)/10;
	for (; i>0; i--,hv++)
		v += *hv;
	return(v/n);
}


dispclass(obj,x,y,e1,e2)
struct chromosome *obj;
{
	static int c=1;
	fprintf(stderr,"%2d ",c++);
	if ((e2 >= 5 && e1 >= 10 && e1 >= e2) ||
	    (e1 >= e2 && e1+e2 >= 25) ||
	    (e2 >= 20 && e1 >= 20)) {
		fprintf(stderr,"****");
		dispstring("****",x,y);
	}
	else
		fprintf(stderr,"-   ");
	if (obj->plist->Cbtype&FRAXA)
		fprintf(stderr," fra(X)");
	else if (obj->plist->Cbtype&XPROBE)
		fprintf(stderr," X");
	fprintf(stderr,"\n");
}