/*
 * dcfeatures.c		Jim Piper	08-09-86
 *
 * Dicentric finding system - measure features of a candidate centromere.
 *
 * The input parameters are (i) the "chromcands" structure for the current
 * object, (ii) the "cand" structure for the current candidate, (iii) the
 * "crossing profile" for the current candidate, (iv) the current object's
 * density profile.
 *
 * MODIFICATIONS
 * -------------
 * jimp 07-06-90	Add cprofcurv() function (for second time).
 *			Need to re-implement kurtosis() yet!
 * jimp	08-04-88	Added mom1ratio feature
 * jimp 02-03-88	Add profcurv feature
 * jimp	13-05-87	detect zero crossing profile, deal appropriately ?
 */
#include <stdio.h>
#include <wstruct.h>
#include <dicen.h>


dcfeatures(chc,cnd,crossprof,prof)
struct chromcands *chc;
struct cand *cnd;
struct object *crossprof, *prof;
{
	struct object *cvhh, *cvhhisto();
	/*
	 * original features
	 */
	cnd->u.sy.profdepth = chc->profmax - cnd->u.sy.profval;
	cnd->u.sy.width = truewidth(crossprof);
	if (cnd->u.sy.width > 1) {
		cnd->u.sy.modality = modality(crossprof,CQUANTUM);
		cnd->u.sy.skewness = skewness(crossprof); 
		cnd->u.sy.cprofpeak = peakval(crossprof);
		cnd->u.sy.cprofsum = cprofsum(crossprof);
		cvhh = cvhhisto(crossprof);
		cnd->u.sy.cvhdeficiency = cvhharea(cvhh) - cnd->u.sy.cprofsum;
		cnd->u.sy.cvhhwdiff = cvhhwdiff(cvhh,crossprof);
		cnd->u.sy.profcurv = secondderiv(prof,cnd->profpos,4,chc->profmax);
		cnd->u.sy.mom1ratio = mom1ratio(crossprof);
		cnd->u.sy.cprofcurv = cprofcurv(crossprof);
		freeobj(cvhh);
	}
	else if (cnd->u.sy.width == 1) {
		cnd->u.sy.modality = 1;
		cnd->u.sy.skewness = 0;
		cnd->u.sy.cprofpeak = peakval(crossprof);
		cnd->u.sy.cprofsum = cprofsum(crossprof);
		cnd->u.sy.cvhdeficiency = 0;
		cnd->u.sy.cvhhwdiff = 0;
		cnd->u.sy.profcurv = 1000;
		cnd->u.sy.mom1ratio = 0;
		cnd->u.sy.cprofcurv = 0;
	}
	else {
		cnd->u.sy.modality = 1;
		cnd->u.sy.skewness = 0; 
		cnd->u.sy.cprofpeak = 0;
		cnd->u.sy.cprofsum = 1;		/* prevent division-by-zero */
		cnd->u.sy.cvhdeficiency = 0;
		cnd->u.sy.cvhhwdiff = 0;
		cnd->u.sy.profcurv = 1000;
		cnd->u.sy.mom1ratio = 0;
		cnd->u.sy.cprofcurv = 0;
	}
	/*
	 * derived features
	 */
	cnd->u.sy.nprofval = 100 * cnd->u.sy.profval / chc->profmax;
	cnd->u.sy.ncprofpeak = 100 * cnd->u.sy.cprofpeak / chc->density;
	cnd->u.sy.ncvhdeficiency = 100 * cnd->u.sy.cvhdeficiency / cnd->u.sy.cprofsum;
}

/*
 * measure skewness as 100*|centroid-center| of crossing profile
 */
skewness(p)
struct object *p;
{
	register struct histogramdomain *hdom;
	register int c,*h,i,m;
	int l;

	hdom = (struct histogramdomain *)p->idom;
	h = hdom->hv;
	l = hdom->npoints;
	c = m = 0;
	for (i=0; i<l; i++) {
		m += *h;
		c += i * *h++;
	}
	if (m == 0) {
		fprintf(stderr,"zero crossing profile sum\n");
		return(0);
	}
	else {
		c -= l*m/2;
		c = 100*c/m;
		if (c < 0)
			c = -c;
		return(c);
	}
}

/*
 * measure curvature at centroid
 */
cprofcurv(p)
struct object *p;
{
	register struct histogramdomain *hdom;
	register int c,*h,i,m;
	int l,j,s1,s2;

	hdom = (struct histogramdomain *)p->idom;
	h = hdom->hv;
	l = hdom->npoints;
	c = m = 0;
	for (i=0; i<l; i++) {
		m += *h;
		c += i * *h++;
	}
	if (m == 0)
		return(0);
	c /= m;
	h = hdom->hv;
	s1=s2=0;
	for (i=c-1; i<=c+1; i++) {
		if (i >= 0 && i < l)
			s1 += h[i];
		j = i-3;
		if (j >= 0 && j < l)
			s2 += h[j];
		j = i+3;
		if (j >= 0 && j < l)
			s2 += h[j];
	}
	s1 *= 2;
	if (s1 == 0)
		return(0);
	else
		return(100*(s1-s2)/s1);
}

/*
 * Compute true width of chromosome at crossing profile
 * by deleting zero points at either end.  Perhaps a better
 * algorithm is called for here.....
 */
truewidth(p)
struct object *p;
{
	register struct histogramdomain *hdom;
	register int i,w,*h;
	hdom = (struct histogramdomain *)p->idom;
	w = hdom->npoints;
	h = hdom->hv + w - 1;
	while (*h-- == 0 && w>0)
		w--;
	if (w == 0)
		return(0);
	h = hdom->hv;
	while (*h++ == 0 && w>0)
		w--;
	return(w);
}

/*
 * compute modality of crossing profile, using a quantum.
 */
modality(prof,quantum)
struct object *prof;
{
	register struct histogramdomain *hdom;
	register int i,l,*h,max;
	int m,curmax,curmin,p,nmax;
	hdom = (struct histogramdomain *)prof->idom;
	h = hdom->hv;
	l = hdom->npoints;
	max = 0;
	for (i=0; i<l; i++) {
		if (*h > max)
			max = *h;
		h++;
	}
	quantum = max*quantum/100;
	/*
	 * search profile for extrema
	 */
	nmax = 0;
	m = 1;
	curmax = 0;
	h = hdom->hv;
	for (i=0; i<hdom->npoints; i++) {
		switch (m) {
		case 1:			/* looking for maximum */
			if (*h > curmax) {
				curmax = *h;
			} else if (*h < curmax - quantum) {
				nmax++;
				m = -m;
				curmin = *h;
			}
			break;
		case -1:		/* looking for minimum */
			if (*h < curmin) {
				curmin = *h;
			} else if (*h > curmin + quantum) {
				m = -m;
				curmax = *h;
			}
			break;
		}
		h++;
	}
	/*
	 * If still searching for a maximum at end of profile,
	 * then see if one is found if we imagine a zero-value
	 * point stuck on the end.
	 */
	if (m == 1 && curmax > quantum)
		nmax++;
	return(nmax);
}

/*
 * Find the maximum value of profile p
 */
peakval(p)
struct object *p;
{
	register struct histogramdomain *hdom;
	register int i,l,max,*h;
	hdom = (struct histogramdomain *)p->idom;
	h = hdom->hv;
	l = hdom->npoints;
	max = 0;
	for (i=0; i<l; i++,h++)
		if (*h > max)
			max = *h;
	return(max);
}


/*
 * Find the integral of profile p
 */
cprofsum(p)
struct object *p;
{
	register struct histogramdomain *hdom;
	register int i,l,sum,*h;
	hdom = (struct histogramdomain *)p->idom;
	h = hdom->hv;
	l = hdom->npoints;
	sum = 0;
	for (i=0; i<l; i++)
		sum += *h++;
	return(sum);
}

int secondderiv(prof,pos,wid,max)
struct object *prof;
int pos,wid;
{
	int diff, *hv;
	struct histogramdomain *hdom = (struct histogramdomain *)prof->idom;
	
	if (pos-wid < 0 || pos+wid > hdom->npoints-1)
		return(0);
	hv = hdom->hv;
	diff = 2*hv[pos] - hv[pos-wid] - hv[pos+wid];
	diff = 1000 + 1000*diff/max;
	return(diff);
}



/*
 * Find the ratio of 1st to 0th central moment of histogram (profile) p
 */
mom1ratio(p)
struct object *p;
{
	register struct histogramdomain *hdom;
	register int l,sum,sumd,*h;
	float centroid,x,mom1;
	hdom = (struct histogramdomain *)p->idom;
	h = hdom->hv;
	sum = 0;
	sumd = 0;
	for (l=hdom->npoints; l>0; l--) {
		sumd += l * *h;
		sum += *h++;
	}
	centroid = (float) sumd / sum;
	h = hdom->hv;
	l = hdom->npoints;
	mom1 = 0.0;
	for (l=hdom->npoints; l>0; l--) {
		x = l - centroid;
		if (x < 0)
			x = -x;
		mom1 += x * *h++;
	}
	mom1 /= sum;
	return((int)(100.0 * mom1));
}