/*
 * xtra_feats.c		Jim Piper		10 February 1987
 *
 * Compute the features missing from boundary candidates and
 * profile candidates respectively.
 *
 * MODIFICATIONS
 * -------------
 * 10-04-91	jimp	Negated curvatures in xtra_b_feats(), to complete
 *			fix of 01-06-89 -- which arises in "bound_cands.c"
 * 01-06-89	ji/jimp	Reverse curvature to suit output of "curvat()".
 * 02-03-88	jimp	added chcp->u.sy.profcurv to xtra_p_feats,
 *			need extra parameter chc as a result.
 */
#include <stdio.h>
#include <wstruct.h>
#include <chromanal.h>
#include <dicen.h>


/*
 * Find the locations where the crossing profile crosses
 * the piece-wise linear boundary.
 * Insert these locations and the corresponding curvature
 * features in the candidate feature structure.
 * This routine expects the following conventions:
 *	cprof to be integer histogram of crossing densities of cpol, used to
 *		find the end-points.
 *	cpol to be integer-scaled-by-8 polygon, unit-spaced
 *	bound to be an integer polygon domain of the piece-wise-linear boundary
 *	hist to be in integer histogramdomain of the curvature of bound
 * The algorithm is to find the end-points of cpol as the first from either
 * end with non-zero cprof values, then to scan bound for the nearest points.
 */
xtra_b_feats(cndp,cprof,cpol,bound,hist)
register struct cand *cndp;
struct object *cprof,*cpol;
struct polygondomain *bound;
struct histogramdomain *hist;
{
	register int c,*hv;
	struct polygondomain *cpdom;
	struct histogramdomain *chdom;
	register struct ivertex *cvtx;

	cpdom = (struct polygondomain *)cpol->idom;
	chdom = (struct histogramdomain *)cprof->idom;

	/*
	 * if not a valid candidate, forget it
	 */
	if (cndp->candtype < 0)
		return;

	/*
	 * find first end-point coordinate pair
	 * check that there IS a non-zero point in crossing profile,
	 * otherwise abandon this candidate (and print out why??)
	 */
	hv = chdom->hv;
	cvtx = cpdom->vtx;
	c = 0;
	while (*hv == 0) {
		hv++;
		cvtx++;
		if (++c == chdom->npoints) {
			fprintf(stderr,"Abandoning type %d candidate - no crossing profile\n",cndp->candtype);
			cndp->candtype = -1;
			return;
		}
	}
	/*
	 * Now find nearest point in boundary polygon to coordinate pair
	 */
	cndp->lboundpos = near_intersect(cvtx->vtX,cvtx->vtY,bound);

	/*
	 * repeat with the second endpoint
	 */
	hv = chdom->hv + chdom->npoints - 1;
	cvtx = cpdom->vtx + cpdom->nvertices - 1;
	while (*hv == 0) {
		hv--;
		cvtx--;
	}
	cndp->rboundpos = near_intersect(cvtx->vtX,cvtx->vtY,bound);

	/*
	 * find the boundary curvatures, find the larger (negative = concave)
	 * curvature, fill in the data structure.
	 */
	cndp->boundmaxcurv = 1;
	cndp->u.sy.curvmax = - hist->hv[cndp->lboundpos];
	cndp->u.sy.curvmin = - hist->hv[cndp->rboundpos];
	if (cndp->u.sy.curvmax == cndp->u.sy.curvmin)
		cndp->boundmaxcurv = 3;
	else if (cndp->u.sy.curvmax > cndp->u.sy.curvmin) {
		c = cndp->u.sy.curvmax;
		cndp->u.sy.curvmax = cndp->u.sy.curvmin;
		cndp->u.sy.curvmin = c;
		cndp->boundmaxcurv = 2;
	}
	cndp->u.sy.curvsum = cndp->u.sy.curvmax + cndp->u.sy.curvmin;
}

/*
 * find serial location of nearest point in boundary domain "bound"
 * (integer-scaled-by-1) to point (x,y) (integer-scaled-by-8).
 */
static near_intersect(x,y,bound)
register int x,y;
register struct polygondomain *bound;
{
	register int i,dx,dy;
	int imin,dist,mindist;
	register struct ivertex *vtx = bound->vtx;
	mindist = 10000000;
	imin = 0;
	for (i=0; i<bound->nvertices; i++,vtx++) {
		dx = 8*vtx->vtX - x;
		dy = 8*vtx->vtY - y;
		dist = isqrt(dx*dx + dy*dy);
		if (dist < mindist) {
			mindist = dist;
			imin = i;
		}
	}
	return(imin);
}

/*
 * Find the location where the crossing profile crosses
 * the chromosome medial axis.  Fill in the location
 * and the longitudinal profile value.
 * This routine expects the following conventions:
 *	cpol to be integer-scaled-by-8 polygon
 *	mpol to be real polygon
 *	prof to be integer histogram
 */
xtra_p_feats(chc,cndp,cpol,mpol,prof)
struct chromcands *chc;
struct cand *cndp;
struct object *cpol, *mpol, *prof;
{
	struct ivertex *ivtx;
	struct fvertex *fvtx;
	struct polygondomain *cpdom,*mpdom;
	struct histogramdomain *hdom;
	int i,j,mini,minj;
	float xdiff,ydiff,diff,diffmin;
	cpdom = (struct polygondomain *)cpol->idom;
	mpdom = (struct polygondomain *)mpol->idom;
	hdom = (struct histogramdomain *)prof->idom;
	if (mpdom->type != 2 || cpdom->type != 1) {
		fprintf(stderr,"Xtra_p_feats: unexpected polygon types %d %d\n",mpdom->type,cpdom->type);
		return;
	}
	ivtx = cpdom->vtx;
	fvtx = (struct fvertex *)mpdom->vtx;
	diffmin = 1000000000.0;
	/*
	 * This brute-force search for the minimum distance between
	 * a point in either polygon needs to be improved - and can be !!
	 */
	for (i=0; i<cpdom->nvertices; i++)
		for (j=0; j<mpdom->nvertices; j++) {
			xdiff = ivtx[i].vtX - 8.0 * fvtx[j].vtX;
			ydiff = ivtx[i].vtY - 8.0 * fvtx[j].vtY;
			diff = xdiff*xdiff + ydiff*ydiff;
			if (diff < diffmin) {
				diffmin = diff;
				mini = i;
				minj = j;
			}
		}
	cndp->profpos = minj;
	cndp->u.sy.profval = hdom->hv[minj];
	cndp->u.sy.profcurv = secondderiv(prof,cndp->profpos,4,chc->profmax);
}