/*
 * bound_cands.c	Simon Towers	09-02-87
 *
 * finds candidate centromeres by boundary analysis
 *
 * 01-06-89	ji/jimp	Reverse curvature to suit output of "curvat()"
 *			(in "bound_curv()" and "curvhist()").
 * 17-05-88	simon	bug fix in curvhist()
 * 14-04-88	simon	bug fix in centromere()
 * 15-03-88	simon	new boundary
 * 23-02-88	simon	Freeing of curve structs after use.
 * 02-06-87	jimp	Free-ing of tmp->assoc removed (removed "obj" !!).
 */
 
#include    <stdio.h>
#include    <wstruct.h>
#include    <dicen.h>
#include    <chromanal.h>
#define     LONG    5
#define     SHORT   3

struct curve {
    int diff;
    int index;
    int x;
    int y;
    int dir;
    int chg;
    char type;      /*  l: left centromere point
                        r: right     "       "
                        a: arm link
                        j: joint                    */
    char type2;
    struct curve *next;
    };
int length;


boundary_cands(obj, boundary, hist, hist_diff, chcand)
struct object *obj;
struct polygondomain *boundary;
struct histogramdomain *hist, *hist_diff;
struct chromcands *chcand;
{
    struct curve *cur1, *cur, *curvhist();	


/* construct boundary and histograms */
    cur1 = curvhist(boundary, hist, hist_diff);
	
/* identify centromere candidates */
    centromere(obj, cur1, chcand, boundary);

    while (cur1 != NULL) {
        cur = cur1;
    	cur1 = cur1->next;
    	Free(cur);
    }
}


struct curve *curvhist(pdom, hdom, hdomdiff)
struct histogramdomain *hdomdiff, *hdom;
struct polygondomain *pdom;
{
	struct curve *cur, *cur1;
	int *hvd, *hv, i, peak, position, direction, nvtx, p2;

	nvtx = hdomdiff->npoints;
	hvd = hdomdiff->hv;
	hv = hdom->hv;
	
	cur1 = cur = NULL;
	for (i=0 ; i<nvtx ; i++) {
		if (hvd[i] > 0) {
			peak = hvd[i];
			direction = hv[i];
			position = i + 1;
			while (hvd[++i] > 0) {
				if (hvd[i] > peak) {
					peak = hvd[i];
					direction = hv[i];
					position = i + 1;
				}
			}
			if (cur == NULL)
				cur1 = cur = (struct curve *) Calloc(1, sizeof(struct curve));
			else {
				cur->next = (struct curve *) Calloc(1, sizeof(struct curve ));
				cur = cur->next;
			}
			cur->diff = -peak;
			cur->index = position - 1;
			cur->dir = direction;
			if (position > pdom->nvertices) position = 1;
			cur->x = (pdom->vtx + position - 1)->vtX;
			cur->y = (pdom->vtx + position - 1)->vtY;
			p2 = position + length;
			if (p2 > nvtx) p2 -= nvtx;
			cur->chg = hv[p2-1];
			p2 = position - length;
			if (p2 < 1) p2 += nvtx;
			cur->chg -= hv[p2-1];
			if (cur->chg >= 180) cur->chg -= 360;
			if (cur->chg <= -180) cur->chg += 360;
			typeit(cur);
		}
	}

	return(cur1);
}


angle(x, y)
{
    double dx, dy, atan2();

if (x < -9 || x > 9 || y < -9 || y > 9)
fprintf(stderr,"ANGLE: out of bounds (%d,%d)\n",x,y);
    dx = (double) x;
    dy = (double) y;
    if (x == 0) {
        if (y < 0)
            dx = -90.;
        else
            dx = 90.;
        }
    else
        dx = atan2(dy, dx) * 180. / 3.14159265;
    if (y < 0) dx += 360.;
#ifdef FIPTEST
fprintf(stderr,"%4d %4d %4d\n",x,y,(int)dx);
#endif FIPTEST

    return((int) dx);
}

typeit(cur)
struct curve *cur;
{
    /* determine type */
    if (cur->diff <= -60 && cur->chg <= -60
        && cur->chg + cur->diff <= -150
        && length == LONG
        && ((cur->dir>=35 && cur->dir<=145)
        || (cur->dir>=215 && cur->dir<=325)))
        cur->type = 'j';
    else if (cur->dir>=35 && cur->dir<=145)
        cur->type = 'r';
    else if (cur->dir>=215 && cur->dir<=325)
        cur->type = 'l';
    else
        cur->type = 'a';
    cur->type2 = ' ';
}

centromere(obj, cur1, chcand, pdom)
struct object *obj;
struct curve *cur1;
struct chromcands *chcand;
struct polygondomain *pdom;
{
    struct cand *cands;
    struct curve *cur, *cur2, *curmax;
    int centroflag, dist, a, i, j, yl, yr, yl2, yr2;
    char lcflag;

    /* locate centromere candidates */
    for (centroflag=0 ; centroflag<=2 ; centroflag++) {
        /* find first relevant angle */
        curmax = NULL;
        cur = cur1;
        while (cur != NULL) {
            if (cur->type2==' ' && (cur->type=='l' || cur->type=='r')) {
                curmax = cur;
                cur = NULL;
            }
            else
                cur = cur->next;
        }
        /* find largest angle */
        if (curmax != NULL)
            cur = curmax->next;
        else
            break;
        while (cur != NULL) {
            if (cur->diff<curmax->diff && cur->type2==' ' && (cur->type=='l' || cur->type=='r'))
                curmax = cur;
            cur = cur->next;
        }
        lcflag = 'l';
        if (curmax->type == 'l') lcflag='r';
        sprintf(&(curmax->type2), "%1d", centroflag+1);

        /* write info about candidate centromere */
	cands = (chcand->cands) + (chcand->ncands)++;
        cands->candtype = 2;
        cands->candnum = chcand->ncands;
        if (lcflag == 'r') {
            cands->lboundpos = curmax->index;
            cands->boundmaxcurv = 1;
        }
        else {
            cands->rboundpos = curmax->index;
            cands->boundmaxcurv = 2;
        }
        cands->u.sy.curvmax = cands->u.sy.curvsum = curmax->diff;
	cands->u.sy.curvmin = 0;

        /* search for matching angle on other side */
        cur = cur1;
        cur2 = NULL;
        dist = 500;
        while (cur != NULL) {
            if (cur->type2==' ' && cur->type == lcflag) {
                if ((a=abs(cur->y - curmax->y)) < dist) {
                    dist = a;
                    cur2 = cur;
                }
                else if (a == dist && cur->diff < cur2->diff)
                    cur2 = cur;
            }
            cur = cur->next;
        }
        /* check centromere ok, and produce vector */
        if (dist <= (length*2 - 1)) {
            sprintf(&(cur2->type2), "%1d", centroflag+1);
            if (lcflag == 'r')
                cands->rboundpos = cur2->index;
            else
                cands->lboundpos = cur2->index;
            cands->u.sy.curvmin = cur2->diff;
	    cands->u.sy.curvsum += cur2->diff;
        }
        else {
            if (lcflag == 'r')
                cands->rboundpos = -1;
            else
                cands->lboundpos = -1;
            if (length == SHORT) break;
        }
    }
    
    if (!centroflag) return;
/* check for one-sided centromeres near ends */
    for (i=centroflag-1 ; i>=0 ; i--) {
    	if ((cands-i)->candtype != 2) continue;
    	if ((cands-i)->lboundpos == -1 || (cands-i)->rboundpos == -1) {
    	    if ((cands-i)->boundmaxcurv == 1) 
    	        yl = (cands-i)->lboundpos;
    	    else
    	        yl = (cands-i)->rboundpos;
    	    if (abs((pdom->vtx+yl)->vtY - obj->idom->line1)<6
    	        || abs((pdom->vtx+yl)->vtY - obj->idom->lastln)<6)
    	        (cands-i)->candtype = -1;
    	}
    }

/* check for close centromeres */
    if (centroflag < 2) return;
    for (i=centroflag-1 ; i>0 ; i--) {
        if ((cands-i)->candtype != 2) continue;
        if ((yl=(cands-i)->lboundpos) == -1)
            yl = (cands-i)->rboundpos;
        yl = (pdom->vtx + yl)->vtY;
        if ((yr=(cands-i)->rboundpos) == -1)
            yr = (cands-i)->lboundpos;
        yr = (pdom->vtx + yr)->vtY;
        j = i - 1;
        while (j >= 0) {
            if ((yl2=(cands-j)->lboundpos) == -1)
                yl2 = (cands-j)->rboundpos;
            yl2 = (pdom->vtx + yl2)->vtY;
            if ((yr2=(cands-j)->rboundpos) == -1)
                yr2 = (cands-j)->lboundpos;
            yr2 = (pdom->vtx + yr2)->vtY;
            if ((cands-j)->candtype==2 && (abs(yl-yl2)<length*2 || abs(yr-yr2)<length*2)) {
                a = (cands-i)->u.sy.curvmax + (cands-i)->u.sy.curvmin
                    - (cands-j)->u.sy.curvmax - (cands-j)->u.sy.curvmin;
                if (a >= 0) {
                    (cands-i)->candtype = -1;
                    j--;
                }
                else
                    (cands-j)->candtype = -1;
                if (--centroflag < 2) return;
            }
            j--;
        }
    }
}

struct histogramdomain *bound_curv(obj, pdom, hist_diff)
struct object *obj;
struct polygondomain **pdom;
struct histogramdomain **hist_diff;
{
	struct histogramdomain *hdom, *makehistodmn();
	struct ivertex *vtx1, *vtx2;
	int a, b, c, d, i, nvtx, *hv, *hvd;
	struct object *tmp, *tmp2, *ibound(), *unitbound();
	struct object *curvat();
	
	tmp = ibound(obj, 1);
	tmp2 = unitbound(tmp);
	freeobj(tmp);
	tmp = curvat(tmp2);
	curvsmooth(tmp, 2);
	*hist_diff = (struct histogramdomain *) tmp->idom;
	*pdom = (struct polygondomain *) tmp2->idom;

	length = LONG;
	nvtx = (*pdom)->nvertices;
	if (nvtx/length <= 20) length = SHORT;
	hdom = makehistodmn(1, nvtx);
	hv = hdom->hv;
	hdom->r = 1;
	hdom->k = 0;
	hdom->l = 0;

	hvd = (*hist_diff)->hv;
	vtx1 = (*pdom)->vtx;
	for (i=1 ; i <= nvtx ; i++, vtx1++) {
/*		*hvd++ *= -1;		removed 01-06-89	*/
		vtx2 = vtx1 - length;
		if (i <= length) vtx2 += nvtx-1;
	        a = vtx1->vtX - vtx2->vtX;
        	b = vtx1->vtY - vtx2->vtY;
		a = angle(a, b);
        	vtx2 = vtx1 + length;
        	if (i+length > nvtx) vtx2 -= nvtx-1;
	        c = vtx2->vtX - vtx1->vtX;
        	d = vtx2->vtY - vtx1->vtY;
		c = angle(c, d);
		if ((a<90 && c>270) || (c<90 && a>270)) {
			a = (a + c + 360) / 2;
			if (a > 359) a -= 360;
			*hv++ = a;
		}
		else
			*hv++ = (a + c) / 2;
	}
	return(hdom);
}