/*
 * yprofs.c
 *
 * Test bed for centromere-labelled fragile X detection
 *
 * yprofs [file]
 *
 * outputs feature vector (preceded by class: 1 fragile, 2 not fragile)
 */

#include <stdio.h>
#include <wstruct.h>
#include <chromanal.h>
#include <otypes.h>
#include "xprobe.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;

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 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 r1,r2,ar,np;
	double os9frig;
	char s[200], *celln;
	FILE *in,*out;
	int intsort();
	int DISPLAY = 1;

	if (argc > 1) {
		in = fopen(argv[1],"r");
		celln = argv[1];
	}
	else {
		in = stdin;
		celln = "stdin";
	}

	numarea = number = 0;
	while ((obj = readobj(in)) != NULL) {
		swabplist(obj->plist);
		/* reject non-chromosome */
		if (obj->plist->otype > 1)
			continue;

		/*
		 * 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);
		cl = 2;
		if (obj->plist->Cbtype & FRAXA)
			cl = 1;
		ar = 10 * np * np / area(obj);
/*		printf("%1d %3d %3d\n",cl,r1*r2/10,ar);	*/
		printf("%1d ",cl);
		if (ar >= 25 && r1*r2 <= 5000) {
			if (cl == 1)
				printf("TP");
			else
				printf("fp");
		}
		else {
			if (cl == 1)
				printf("fn");
		}
		printf("\n");

		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);
	}
}