/********************************/
/*    evlobj			*/
/*    Evaluate object		*/
/********************************/
#include    "focas1.h"
#define	BIN    10		/* Bin width is sigma from skyhst */
unsigned char  *silh;		/* Allocated in isophote() */

int     skytype;		/* 0 = mean of histogram, 1 = mode of histogram
				   */
static struct image pevl;

evlobj (ob, area, diag, fd)
struct objrec  *ob;
struct areas   *area;
int     diag;
FILE * fd;
{
    PIXEL * cp;
    struct arblk   *ar;
    int     i, j, k, l, n, xi, yi, xl, xr, xci, yci;
    int     dx, dy, xw, yw, nxw, nyw;
    int     hstgrm[512];
    PIXEL step, start;
    float   x, y, xc, yc, a, b, c, r;
    float   xv[3], yv[3], am[3][3], bm[3][3];
    float   s1, s2;
    float   fastr ();


    isophote (area);

    ob->eflgs &= ~DARK;
    ob->eflgs &= ~PEAK;
    ob->eflgs &= ~EVAL;

    if (ob->nsbr < 10) {
	if (diag)
	    fprintf (fd, "evlobj: Failed because nsbr = %d < 10\n", ob->nsbr);
	return (1);
    }

    if (rdimg (pfd, area->rct, &pevl) < 0)
	return (0);

    ar = &area->d;
    xv[0] = yv[0] = 1.0;
    s1 = s2 = 0.;
    for (i = 0; i < 3; i++)
	for (j = 0; j < (3 - i); j++)
	    am[i][j] = bm[i][j] = 0.0;

 /* Compute moments of object intensity within detection area */

    for (i = pevl.dy - 1, l = 0; i >= 0; i--) {
	yi = ar->ychst - pevl.y0 - i;
	if ((yi >= 0) && (yi < ar->nlines)) {
	    for (; l < ar->nx[yi]; l++) {
		xl = ar->xla[l] - pevl.x0;
		xr = ar->xra[l] - pevl.x0;
		cp = pevl.pic + i * pevl.dx + xl;
		y = i;
		for (j = 1; j < 3; j++)
		    yv[j] = yv[j - 1] * y;
		for (xi = xl; xi <= xr; xi++) {
		    x = xi;
		    j = *cp++;
		    if (j > sp.satr)
			ob->eflgs |= PEAK;
		    a = film (j);
		    s1 += a;
		    s2 += a * a;
		    for (j = 1; j < 3; j++)
			xv[j] = xv[j - 1] * x;
		    for (j = 0; j < 3; j++) {
			y = yv[j];
			for (k = 0; k < (3 - j); k++) {
			    b = y * xv[k];
			    am[j][k] += a * b;
			    bm[j][k] += b;
			}
		    }
		}
	    }
	}
    }
    if (bm[0][0] > 0.) {
	s1 /= bm[0][0];
	s2 = s2 / bm[0][0] - s1 * s1;
	if (s2 > 0.)
	    s2 = sqrt (s2);
	else
	    s2 = 0.;
	ob->sLi = s2;
    }

    for (j = 0; j < 3; j++)
	for (k = 0; k < (3 - j); k++)
	    am[j][k] -= ob->sbr * bm[j][k];
    if (am[0][0] <= 0.0)
	ob->eflgs |= DARK;

 /* Peak position and peak luminosity */
    c = 0.;
    xc = pevl.dx /2.;
    yc = pevl.dy /2.;
    for (yci = pevl.dy - 1, l = 0; yci >= 0; yci--) {
	yi = ar->ychst - pevl.y0 - yci;
	if ((yi >= 0) && (yi < ar->nlines)) {
	    for (; l < ar->nx[yi]; l++) {
		xl = ar->xla[l] - pevl.x0;
		xr = ar->xra[l] - pevl.x0;
		cp = pevl.pic + yci * pevl.dx + xl;
		for (xci = xl; xci <= xr; xci++) {
		    b = film (*(cp - pevl.dx - 1));
		    b += film (*(cp - pevl.dx));
		    b += film (*(cp - pevl.dx + 1));
		    b += film (*(cp - 1));
		    b += film (*(cp));
		    b += film (*(cp + 1));
		    b += film (*(cp + pevl.dx - 1));
		    b += film (*(cp + pevl.dx));
		    b += film (*(cp + pevl.dx + 1));
		    b -= 9 * ob->sbr;
		    if (b > c) {
			c = b;
			xc = xci;
			yc = yci;
		    }
		    cp++;
		}
	    }
	}
    }

    xci = xc;
    yci = yc;
    ob->Lc = c;
    ob->xc = lpcoords (xc + pevl.x0, 1, 1);
    ob->yc = lpcoords (yc + pevl.y0, 2, 1);
    ob->ra = ob->xc * sp.xfrm[0][0] + ob->yc * sp.xfrm[0][1] + sp.xfrm[0][2];
    ob->dec = ob->xc * sp.xfrm[1][0] + ob->yc * sp.xfrm[1][1] + sp.xfrm[1][2];

 /* Central moments */
    if (bm[0][0] != 0.) {
	ob->cx = lpcoords (bm[0][1] / bm[0][0] + pevl.x0, 1, 1);
	ob->cy = lpcoords (bm[1][0] / bm[0][0] + pevl.y0, 2, 1);
	ob->yy = (bm[2][0] + yc * yc * bm[0][0] - 2.0 * yc * bm[1][0]) /
	    bm[0][0];
	ob->xx = (bm[0][2] + xc * xc * bm[0][0] - 2.0 * xc * bm[0][1]) /
	    bm[0][0];
	ob->xy = (bm[1][1] + xc * yc * bm[0][0] - xc * bm[1][0] -
		yc * bm[0][1]) / bm[0][0];
    }
    if (am[0][0] != 0.) {
	ob->icx = lpcoords (am[0][1] / am[0][0] + pevl.x0, 1, 1);
	ob->icy = lpcoords (am[1][0] / am[0][0] + pevl.y0, 2, 1);
	ob->iyy = (am[2][0] + yc * yc * am[0][0] - 2.0 * yc * am[1][0]) /
	    am[0][0];
	ob->ixx = (am[0][2] + xc * xc * am[0][0] - 2.0 * xc * am[0][1]) /
	    am[0][0];
	ob->ixy = (am[1][1] + xc * yc * am[0][0] - xc * am[1][0] -
	    yc * am[0][1]) / am[0][0];
    }

 /* test for noise object

    if ((ob->ixx <= 0.) || (ob->iyy <= 0.)) {
	if (diag)
	    fprintf (fd,
		    "evlobj: Failed at x = %6.1f y = %6.1f because ixx = %g or iyy = %g\n",
		    ob->xc, ob->yc, ob->ixx, ob->iyy);
	return (1);
    }
  */

 /* Total luminosity */
    if (ob->subent == 0) {
	ar = &area->t;
	ob->Ltotal = 0.;
	for (i = pevl.dy - 1, l = 0; i >= 0; i--) {
	    yi = ar->ychst - pevl.y0 - i;
	    if ((yi >= 0) && (yi < ar->nlines)) {
		xl = ar->xla[yi] - pevl.x0;
		xr = ar->xra[yi] - pevl.x0;
		cp = pevl.pic + i * pevl.dx + xl;
		k = xr - xl;
		for (j = 0; j <= k; j++)
		    ob->Ltotal += film (*cp++) - ob->sbr;
	    }
	}
    }

    ob->Li = am[0][0];

    if (ob->Li <= 0.) {
	if (ob->Li < 0.)
	    ob->mag = sp.magfst - log (-ob->Li) / 0.92103;
	else
	    ob->mag = 100.;
    }
    else
	ob->mag = sp.magfst - log (ob->Li) / 0.92103;

 /* ob->r1 = 0.; */
    ob->ir1 = ob->ir3 = ob->ir4 = 0.;

 /* Radial moments */
    ar = &area->d;
    for (i = pevl.dy - 1, l = 0; i >= 0; i--) {
	yi = ar->ychst - pevl.y0 - i;
	if ((yi >= 0) && (yi < ar->nlines)) {
	    for (; l < ar->nx[yi]; l++) {
		xl = ar->xla[l] - pevl.x0;
		xr = ar->xra[l] - pevl.x0;
		cp = pevl.pic + i * pevl.dx + xl;
		for (k = xl; k <= xr; k++) {
		    r = fastr (k - xci, i - yci);
		    c = film (*cp++) - ob->sbr;
		/* ob->r1 += r; */
		    ob->ir1 += r * c;
		    ob->ir3 += r * r * r * c;
		    ob->ir4 += r * r * r * r * c;
		}
	    }
	}
    }

 /* ob->r1 /= bm[0][0]; */
    if (am[0][0] != 0.) {
	ob->ir1 /= am[0][0];
	ob->ir3 /= am[0][0];
	ob->ir4 /= am[0][0];
    }
  /*
    if ((ob->ir1 <= 0.) || (ob->ir3 <= 0.) || (ob->ir4 <= 0.)) {
	if (diag)
	    fprintf (fd, "evlobj: Failed at x = %6.1f y = %6.1f because ir1 = %g or ir3 = %g or ir4 = %g\n",
		    ob->xc, ob->yc, ob->ir1 / ob->Li, ob->ir3 / ob->Li,
		    ob->ir4 / ob->Li);
	return (1);
    }
  */

  /* Fixed aperature luminosity */
    ob->Lfca = 0.;
    for (i = pevl.dy - 1; i >= 0; i--) {
	cp = pevl.pic + i * pevl.dx;
	for (j = 0; j < pevl.dx; j++) {
	    r = fastr (j - xci, i - yci);
	    c = film (*cp++) - ob->sbr;
	    if (r < sp.rfca)
		ob->Lfca += c;
	}
    }

 /* Isophotal parameters */
    ar = &area->i;
    n = 0;
    ob->ispht = 0.;
    for (i = pevl.dy - 1, l = 0; i >= 0; i--) {
	yi = ar->ychst - pevl.y0 - i;
	if ((yi >= 0) && (yi < ar->nlines)) {
	    for (; l < ar->nx[yi]; l++) {
		xl = ar->xla[l] - pevl.x0;
		xr = ar->xra[l] - pevl.x0;
		cp = pevl.pic + i * pevl.dx + xl;
		for (j = xl; j <= xr; j++) {
		    ob->ispht += film (*cp++) - ob->sbr;
		    n++;
		}
	    }
	}
    }
    ob->ispht /= n;

 /* Do reality test */
    /*
    if (ob->ssbr <= 0.)
	ob->prob = -100.;
    else
        ob->prob = (ob->Lc / 9 - ob->ispht) / (ob->ssbr / 3);
    */
    ob->prob = (ob->Lc / 9 - ob->ispht) / (sp.skyhw / 3);

    if (ob->prob < sp.sig) {
	if (diag) fprintf(fd,
	    "evlobj: Failed sig. test at x = %6.1f y = %6.1f\n",
	    ob->xc, ob->yc);
        return(1);
    }

 /* Determine average x and y widths */
    ar = &area->d;
    dx = area->xmax - area->xmin + 3;
    dy = ar->nlines + 2;
    ob->xavg = ob->yavg = xw = yw = nxw = nyw = 0;
    for (i = 0; i < dy; i++) {
	for (j = 0, xw = 0; j < dx; j++) {
	    if (silh[i * dx + j] > 0)
		xw++;
	    else {
		if (xw > 0) {
		    ob->xavg += xw;
		    nxw++;
		}
		xw = 0;
	    }
	}
    }
    for (j = 0; j < dx; j++) {
	for (i = 0, yw = 0; i < dy; i++) {
	    if (silh[i * dx + j] > 0)
		yw++;
	    else {
		if (yw > 0) {
		    ob->yavg += yw;
		    nyw++;
		}
		yw = 0;
	    }
	}
    }
    ob->xavg /= nxw;
    ob->yavg /= nyw;
    ob->eflgs |= EVAL;
    return (0);
}