/*
 * profile.c		Jim Piper		20/1/84
 *
 * Extract profiles from banded chromosomes and similar objects.
 * Input is object and mid-points polygon.
 *
 * Also extract mid-points polygon from object know to be rotated
 * so that original guesstimate of axis direction is vertical
 *
 * MODIFICATIONS
 * -------------
 * 13-02-91 jimp	Added profiles -4, -3 to multiprofiles() --
 *			-4 : left/right density difference
 *			-3 : ratio of -4 to integrated density
 */
#include <stdio.h>
#include <math.h>
#include <wstruct.h>
#include <chromanal.h>


struct object *midptspoly(obj)
struct object *obj;
{
	struct object *p, *makemain();
	struct polygondomain *pdom, *makepolydmn();
	register struct ivertex *vtx;
	struct iwspace iwsp;
	register int kl, n;
	n = 1 + obj->idom->lastln - obj->idom->line1;
	pdom = makepolydmn(1,NULL,0,n,1);
	pdom->nvertices = n;
	vtx = pdom->vtx;
	p = makemain(10,pdom,NULL,NULL,NULL);
	initrasterscan(obj,&iwsp,0);
	while (nextinterval(&iwsp) == 0) {
		if (iwsp.nwlpos >= 1) {
			kl = iwsp.lftpos;
			/*
			 * if there is a gap in the object then must
			 * reduce the number of vertices by the number
			 * of lines skipped
			 */
			pdom->nvertices -= (iwsp.nwlpos - 1);
		}
		if (iwsp.intrmn == 0) {
			vtx->vtX = (kl + iwsp.rgtpos)/2;
			vtx->vtY = iwsp.linpos;
			vtx++;
		}
	}
	return(p);
}


modifytips(p)
struct object *p;
{
	register int i,n,m,x;
	register struct ivertex *vtx;
	struct polygondomain *pdom;

	pdom = (struct polygondomain *)p->idom;
	vtx = pdom->vtx;
	n = pdom->nvertices;
	m = n/8;
	if (m<2)
		m = 2;
	x = 2 * vtx[m].vtX;
	for (i=0; i<m; i++)
		vtx[i].vtX = x - vtx[2*m-i].vtX;
	x = 2 * vtx[n-m-1].vtX;
	for (i=0; i<m; i++)
		vtx[n-1-i].vtX = x - vtx[n-1-2*m+i].vtX;
}


removeblips(p)
struct object *p;
{
	register int i,n;
	register struct ivertex *vtx;
	struct polygondomain *pdom;

	pdom = (struct polygondomain *)p->idom;
	vtx = pdom->vtx;
	n = pdom->nvertices;
	/* alter position of atypical, isolated points */
	for (i=0; i<n-2; i++) {
		/* look at x coordinates
		/* if x coordinates of next but adjoining points the same, force
		/* x coordinate of middle point to be same
		/* (cheaper than further checks) */
		if (vtx[i].vtX == vtx[i+2].vtX)
			vtx[i+1].vtX = vtx[i].vtX;
	}
}



/*
 * return a "histogram" type object as the profile.
 *
 * "mpol" MUST be a mid-points polygon derived from the object, i.e.
 * it must have one vertex per line of the object.
 *
 * "moment" selects profile type :
 *      -4	density difference left/right
 *      -3	ratio of density difference left/right to integrated density
 *	-2	max density
 *	-1	mean density (width defined as sum of non-zero points)
 *	 0	integrated density
 *	 1	first moment (sigma |md| / sigma m) (cf. centroid)
 *	 2	second moment (sigma mdd / sigma m) (radius of gyration)
 */
struct object *profile(obj, mpol, moment)
struct object *obj, *mpol;
{
	struct ivertex *vtx;
	register GREY gv, *g, maxg;
	register int l1, k, i, x, sum, sumd, d, mom;
	int n, *hv, sump;
	struct iwspace iwsp;
	struct gwspace gwsp;
	struct object *prof, *makemain();
	struct histogramdomain *hdom;
	
	vtx = ((struct polygondomain *)mpol->idom)->vtx;
	l1 = obj->idom->line1;
	n = 1 + obj->idom->lastln - l1;
	hdom = (struct histogramdomain *) Malloc(sizeof(struct histogramdomain) + n * sizeof(int));
	hdom->type = 1;
	hdom->linkcount = 0;
	hdom->r = 0;	/* display with vertical base */
	hdom->l = l1;
	hdom->k = obj->idom->kol1;
	hdom->npoints= n;
	hv = hdom->hv = (int *) (hdom + 1);
	prof = makemain(13,hdom,NULL,NULL,NULL);
	initgreyscan(obj,&iwsp,&gwsp);
	while (nextgreyinterval(&iwsp) == 0) {
		g = gwsp.grintptr;
		if (iwsp.nwlpos != 0) {
			maxg = 0;
			sum = 0;
			sump = 0;
			sumd = 0;
			x = vtx->vtX;
		}
		for (k=iwsp.lftpos; k<=iwsp.rgtpos; k++) {
			gv = *g++;
			if (gv > maxg)
				maxg = gv;
			sum += gv;
			if (gv > 0)
				sump++;
			if (moment > 0) {
				mom = gv;
				d = k-x;
				if (d < 0)
					d = -d;
				for (i=0; i<moment; i++)
					mom *= d;
				sumd += mom;
			}
		}
		if (iwsp.intrmn == 0) {
			if (moment == -1 && sump > 0)
				sum /= sump;
			*hv = (moment == -2? maxg: (moment <= 0? sum: (sum==0? 0: 10*sumd/sum)));
			hv++;
			vtx++;
		}
	}
	return(prof);
}




/*
 * Make profile from object and input mid-points polygon.
 * Proceed as follows :
 *	"mpol" is a unit-spaced polygon.
 *	Step along it in unit steps, and find perpendicular direction.
 *	At each step, compute desired profile value.
 * Format same as profile routine in "profile.c", but mpol can extend
 * beyond object boundaries without catastrophic effects.
 */
struct object *nprofile(obj, mpol, moment)
struct chromosome *obj;
struct object *mpol;
{
	struct object *prof;
	multiprofiles(obj,mpol,&prof,&moment,1);
	return(prof);
}


/*
 * 02-06-86	jimp&AUC	2nd moment about slice centroid
 *
 * make multiple profiles as specified in array "moments", which should be
 * monotonic increasing, minimum value >= -1.
 * This saves considerable time in the extraction of slices perpendicular
 * to medial axis.  The maximum number of simultaneous profiles is MAXSPROF.
 */
#define MAXSPROF	5
multiprofiles(obj, mpol, prof, moments, nmoments)
struct chromosome *obj;
struct object *mpol, **prof;
int *moments;
{
	register int d,g,i,sum,maxg;
	int centroid,j, k, mlines, sline, sumd[MAXSPROF], sump, sumq, mom;
	int sumleft, sumright;
	int areao, halfwidth, *hv[MAXSPROF], rbuf[2*MAXHALFWIDTH+1];
	float x,y,sintheta,costheta;
	struct object *makemain();
	struct histogramdomain *hdom;
	struct polygondomain *p;

	p = (struct polygondomain *) mpol->idom;
	if (mpol->type != 10 || p->type != 2) {
		fprintf(stderr,"not a type 2 polygon\n");
		exit(91);
	}

	/*
	 * approximate 1.5 mean width of object
	 * to provide a limiting distance from axis for profile extraction
	 */
	mlines = p->nvertices;
	areao = (obj->plist==NULL || obj->plist->area<=0)? area(obj): obj->plist->area;
	halfwidth = 2*areao/(3*mlines);
	if (halfwidth > MAXHALFWIDTH)
		halfwidth = MAXHALFWIDTH;

	for (i=0; i<nmoments; i++) {
		hdom = (struct histogramdomain *) Calloc(sizeof(struct histogramdomain) + mlines , sizeof(int));
		hdom->type = 1;
		hdom->linkcount = 0;
		hdom->r = 0;	/* display with vertical base */
		hdom->l = obj->idom->line1;
		hdom->k = obj->idom->kol1;
		hdom->npoints= mlines;
		hv[i] = hdom->hv = (int *) (hdom + 1);
		prof[i] = makemain(13,hdom,NULL,NULL,NULL);
	}

	for (sline=1; sline<=mlines; sline++) {
		/*
	 	* get next column, line (x,y) in unit-spaced steps
	 	* along smoothed mid-points polygon, and local orientation.
	 	* sline is line number in straightened object
	 	*/
		setnextpoint(p,&x,&y,&sintheta,&costheta,sline);
		/*
		 * rotate into buffer.
		 */
		rotate(obj,x,y,sintheta,costheta,halfwidth,rbuf);
		/*
		 * compute density profile value
		 */
		maxg = sum = sump = centroid = sumleft = sumright = 0;
		for (i=0; i<nmoments; i++)
			sumd[i] = 0;
		for (k=0; k<=halfwidth*2; k++) {
			g = rbuf[k];
			if (g > maxg)
				maxg = g;
			sum += g;
			if (k <= halfwidth)
				sumleft += g;
			else
				sumright += g;
			if (g > 0)
				sump++;
			mom = g;
			d = halfwidth - k;
			centroid += d*g;
			if (d < 0)
				d = -d;
			i = 0;
			for (j=0; j<nmoments; j++) {
				for (; i<moments[j]; i++)
					mom *= d;
				sumd[j] += mom;
			}
		}
		sumq = sum;
		/*
		 * modify 2nd moment profile to be about slice centroid
		 */
		if (sum != 0) 
			for (j=0; j<nmoments; j++) {
				if (moments[j] == -1 && sump > 0)
					sum /= sump;
				if (moments[j] == 2)
					sumd[j] -= centroid*centroid/sum;
		}
		for (i=0; i<nmoments; i++) {
			switch (moments[i]) {
				case -4:
					*hv[i] = (sumleft - sumright);
					break;
				case -3:
					*hv[i] = sumq==0? 0: 100*(sumleft - sumright)/sumq;
					break;
				case -2:
					*hv[i] = maxg;
					break;
				case -1:
				case 0:
					*hv[i] = sum;
					break;
				case 1:
				default:
					*hv[i] = sumq==0? 0: 10*sumd[i]/sumq;
					break;
			}
			hv[i]++;
		}
	}
	return;
}



setnextpoint(p,xx,yy,sini,cosi,sln)
struct polygondomain *p;
float *xx,*yy,*sini,*cosi;
{
	struct fvertex *fv;
	int int1,int2;
	fv = (struct fvertex *) p->vtx + sln - 1;
	*xx=fv->vtX;
	*yy=fv->vtY;
	int1=sln-HALF_SLOPE_NHD;
	int2=sln+HALF_SLOPE_NHD;
	if (int1 < 1) {
		int1=1;
		int2++;
	}
	if (int2 > p->nvertices) {
		int2=p->nvertices;
		if (int1 > 1)
			int1--;
	}
	slope(p,int1,int2,sini,cosi);
}


rotate(obj,x2,y2,sin2,cos2,halfwidth,g)
struct object *obj;
float x2,y2,sin2,cos2;
int *g;
{
	register i;
	int xxx,yyy,xinc,yinc,xbar,ybar,wd1,wd2;
	int isin2,icos2,ix3,iy3,ixfrac,iyfrac;
	int ixx,iyy,jxx,jyy;
	GREY g4[4];
	icos2=181*cos2+0.48;
	isin2=181*sin2+0.48;
	ix3=x2;
	iy3=y2;
	/*
	 * fractional remainder, out of 181
	 */
	ixfrac=181*(x2-ix3);
	iyfrac=181*(y2-iy3);
	/*
	 * everything will now be done centered round ix3,iy3.
	 */
	for (i= -halfwidth; i<=halfwidth; i++) {
		/*
		 * coordinates, minus ix3, iy3 respectively, *181
		 */
		xxx=ixfrac+i*icos2;
		yyy=iyfrac+i*isin2;
		ixx=xxx/181;
		iyy=yyy/181;
		xinc=xxx-ixx*181;
		if (xinc < 0) {
			xinc += 181;
			ixx--;
		}
		xbar=181-xinc;
		yinc=yyy-iyy*181;
		if (yinc < 0) {
			yinc += 181;
			iyy--;
		}
		ybar=181-yinc;
		/*
		 * xinc, xbar, yinc, ybar should now all be non-negative
		 */
		jxx=ix3+ixx;
		jyy=iy3+iyy;
		/*
		 * interpolated grey value
		 */
		grey4val(obj,jyy,jxx,g4);
		*g++ = (g4[0] * ybar*xbar + g4[1] * ybar*xinc
		     + g4[2] * yinc*xbar + g4[3] * yinc*xinc
		     + 16383) / 32767;
	}
}

/*
 * compute estimate of slope of polygon vertices between lf and lr
 * using a gradient-finding convolution independently in X and Y.
 */
static
slope(p,lf,lr,sinsl,cossl)
struct polygondomain *p;
float *sinsl,*cossl;
{
	register struct fvertex *fv, *gv;
	register float s,c,sscc;
	register int i;

	s = 0.0;
	c = 0.0;
	fv = (struct fvertex *) p->vtx + lf - 1;
	gv = (struct fvertex *) p->vtx + lr - 1;
	lf = (lr-lf+1)/2;
	for (i=0; i<lf; i++) {
		c += gv->vtY - fv->vtY;
		s += fv->vtX - gv->vtX;
		gv--;
		fv++;
	}
	sscc = sqrt(s*s + c*c);
	*sinsl = s/sscc;
	*cossl = c/sscc;
}