/*
 * nxpclass.c
 *
 * classify x-probed chromosomes with features
 *
 * nxpclass [-m] [-p] classparams [inputdata]
 *
 * Flags:
 *	-m	known to be of male sex.  Otherwise assumed female.
 *	-p	use pooled covariance matrices
 *	-d	display chromosomes classified as X
 *	-D	display X chromosomes in cell context
 */
#include <stdio.h>
#include <wstruct.h>
#include <chromanal.h>
#include <otypes.h>
#include "xprobe.h"
#include <vcv.h>

#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;

/* logs of prior probabilities */
static float priorp[MAXCLASS] = {
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 
1.0, 1.0, 1.0, 1.0};
static int nclass,ndim,nsfv,sfv[MAXDIM];

extern int POOL;

main(argc,argv)
char **argv;
{
	FILE *fclass;
	int c,i,j,n,nc,trueclass,class,vclass,ivclass,xrank;
	char celln[100], cellwas[100];
	struct xfeats xf[MAXRANK];
	struct chromosome *obj, *objl[MAXRANK], *readobj();
	struct chromosome *cellobj[MAXOBJS];
	struct chromplist *plist;
	struct pframe f;
	FILE *in, *cin, *tin;
	int SEX = SEX_FEMALE;
	int DISPLAY = 0, PAUSE = 0;

	
	setbuf(stderr,NULL);

	while (argc > 1 && argv[1][0] == '-') {
		switch (argv[1][1]) {
			case 'p':
				POOL = 1;
				break;
			case 'm':
				SEX = SEX_MALE;
				break;
			case 'd':
				DISPLAY = 1;
				break;
			case 'D':
				DISPLAY = 2;
				break;
			case 'P':
				PAUSE = 1;
				tin = fopen("/dev/tty","r");
				setbuf(tin,NULL);
				break;
			default:
				fprintf(stderr,"%s not understood\n",argv[1]);
				exit(1);
		}
		argc--;
		argv++;
	}
	/* read class mean vectors and variance-covariance matrix */
	fclass = fopen(argv[1],"r");
	if (fclass == NULL) {
		fprintf(stderr,"Cannot open classifier parameter file\n");
		exit(2);
	}
	readclassparams(fclass,&nclass,&ndim,sfv,&nsfv);

	if (argc > 2) {
		in = fopen(argv[2],"r");
		if (in == NULL) {
			fprintf(stderr,"Cannot open file %s\n",argv[2]);
			exit(2);
		}
	}
	else
		in = stdin;

	if (DISPLAY) {
		setout("stosun -p270,0 -d1 -r870,870",'p');
		dinit();
		defaultlut();
		lut_overlay();
		sleep(4);	/* give stosun time to get up */
		lwidth(1);
	}
	/*
	 * Read chromosome data file.  Divide into cells
	 * according to cell-name stored in WDD slots of plist.
	 * Classify by multivariate Gaussian maximum likelihood,
	 * choose the best 1 or 2 in each cell and output.
	 * Display if desired.
	 */
	obj = readobj(in);
	if (obj != NULL)
		strcpy(cellwas,&obj->plist->wdd[0]);
	n = 0;
	do {
		if (obj != NULL) {
			strcpy(celln,&obj->plist->wdd[0]);
		}
		if (obj == NULL || strcmp(celln,cellwas)) {
			if (DISPLAY == 1) {
				colour(255);
				clear();
				defaultlut();
				lut_overlay();
				xpdispinit(objl,n,&f);
			}
			else if (DISPLAY == 2) {
				cin = fopen(cellwas,"r");
				if (cin != NULL) {
					nc = 0;
					while ((cellobj[nc]=readobj(cin)) != NULL)
						nc++;
					fclose(cin);
					colour(255);
					clear();
					xpdispinit(cellobj,nc,&f);
					for (i=0; i<nc; i++) {
						picframe(cellobj[i],&f);
						freeobj(cellobj[i]);
					}
				}
				else {
					xpdispinit(objl,n,&f);
				}
			}
			xclassify(objl,xf,n,SEX);
			for (i=0; i<n; i++) {
				if (DISPLAY == 1 || (DISPLAY == 2 && cin == NULL)) {
					colour(248);
					picframe(objl[i],&f);
				}
				if (xf[i].vcvclass == 1) {
					c = 4;
					writeobj(stdout,objl[i]);
					fprintf(stderr,"%s :%d",&(objl[i]->plist->wdd[0]),objl[i]->plist->number);
					if (xfragile(objl[i])) {
						objl[i]->plist->btype |= FRAXA;
						c = 2;
						fprintf(stderr," - fragile ?");
						if (DISPLAY == 2 && cin != NULL)
							cornerdisp(objl[i]);
					}
					if (DISPLAY)
						xpdisp(objl[i],&f,c,1);
					fprintf(stderr,"\n");
				}
				else if (DISPLAY)
					xpdisp(objl[i],&f,1,1);
			}
			if (DISPLAY) {
				update();
				if (PAUSE) {
					fprintf(stderr,"type return to continue\n");
					while(getc(tin) != '\n')
						;
				}
				else
					sleep(3);
			}
			strcpy(cellwas,celln);
			n = 0;
		}
		else {
			objl[n++] = obj;
			obj = readobj(in);
		}
	} while (obj != NULL);
	if (DISPLAY)
		finish();
}


xclassify(ol,xf,n,s)
struct chromosome **ol;
struct xfeats *xf;
int n,s;
{
	int i;
	for (i=0; i<n; i++)
		vcvclassxobj(ol[i],xf+i);
}


vcvclassxobj(o,xf)
register struct chromosome *o;
register struct xfeats *xf;
{
	float fvec[MAXDIM],svec[MAXDIM];
	struct chromcl chrcl;
	int i;

	getxpfeats(o,xf);
	fxfeat(fvec,xf);
	for (i=0; i<nsfv; i++)
		svec[i] = fvec[sfv[i]];
	/* likelihood classifier */
	xf->vcvclass = classdist(svec,priorp,chrcl.cllik,nclass,nsfv);
	xf->vcvscore = chrcl.cllik[0].lval;
}


int xfragile(obj)
struct chromosome *obj;
{
	struct object *spobj, *mpol, *cvh, *rect;
	struct object *fipconvhull(), *spinsqueeze(), *minwrect();
	struct object *poly, *centpoly(), *chromaxis2();
	struct object *prof, *p[5], *nprofile();
	struct histogramdomain *h;
	int cl,r1,r2,ar,np,maxpos,peak;
	double os9frig;

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

	/*
	 * find axis
	 */
	mpol = chromaxis2(spobj);

	/*
	 * 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);

	prof = p[4];
	histonormalise(prof);
	h = (struct histogramdomain *)prof->idom;
	np = h->npoints;
	r1 = tipratio(prof,maxpos);
	prof = p[GPROF];
	histonormalise(prof);
	r2 = tipratio(prof,maxpos);
	ar = 10 * np * np / area(obj);
	if (ar >= 25 && r1*r2 <= 5000)
		cl = 1;
	else
		cl = 0;
	freeobj(p[0]);
	freeobj(p[1]);
	freeobj(p[2]);
	freeobj(p[3]);
	freeobj(p[4]);
	freeobj(cvh);
	freeobj(mpol);
	freeobj(rect);
	freeobj(spobj);
	return(cl);
}


cornerdisp(obj)
struct object *obj;
{
	struct pframe f;
	struct intervaldomain *idom = obj->idom;
	f.dx = 3500;
	f.dy = 3500;
	f.ix = 4;
	f.iy = 3;
	f.scale = 3;
	f.ox = (idom->kol1 + idom->lastkl)/2;
	f.oy = (idom->line1 + idom->lastln)/2;
	colour(248);
	picframe(obj,&f);
	xpdisp(obj,&f,2,1);
}