/*
 * detect - FOCAS detection
 */
#include	"focas1.h"
#define		SKY	8

PIXEL  *tmpbfr;
extern	PIXEL * rcdptr[];
float  *dsky, th0, th1, wts, wts2;
int     yp, inclns = LFWRD + LBACK + 1;
int    nent, nsb;
int     xfltsz, yfltsz;
PIXEL  *lnptrs[NYFLT];
float  *smthln;

dodetect (catalog, first, last, minsky, maxsky, skfl)
char   *catalog, *skfl;
int    *first, *last;
float  *minsky, *maxsky;
{
    char   str[80];
    PIXEL  *cp;
    float  a0, a1, b;
    int    i, j, nln, rlc0[4], rlc1[4], x, nsky = 0;
    float  nn, ns;
    int    xlen, ylen, skyline;
    float  skyot;
    static struct imgdef    skphdr;

 /* program initialization */
    cfd = catopen (catalog, 3);
    pfd = fndfld (sp.ptfl, 0);
    afd = fndar (sp.arfl, 1);
    sp.xs = pthdr.crpix1;
    sp.ys = pthdr.crpix2;
    sp.kpts = pthdr.naxis1;
    sp.kscan = pthdr.naxis2;
    memory (sp.kpts);

    if (*first < 1)
	*first = 1;
    if (*last < 1)
	*last = sp.kscan;
    if (*minsky < pthdr.minpxl)
	*minsky = pthdr.minpxl;
    if (*maxsky > pthdr.maxpxl)
	*maxsky = pthdr.maxpxl;
    if (*minsky > *maxsky)
	focaserr (1, "Error in sky limits", "");

    for (i = 0, wts = 0., wts2 = 0.; i < sp.fy; i++)
	for (j = 0; j < sp.fx; j++) {
	    wts += dfltr[i][j];
	    wts2 += dfltr[i][j] * dfltr[i][j];
	}
    xfltsz = sp.fx / 2;
    yfltsz = sp.fy / 2;
    xlen = sp.kpts - xfltsz;
    ylen = *last - yfltsz;
    skyline = 0;
    if (sp.thrln0 < 0.)
	th0 = -sp.thrln0 * wts;
    else
	th0 = sp.thrln0 * sp.skyhw * sqrt (wts2);
    if (sp.thrln1 < 0.)
	th1 = -sp.thrln1 * wts;
    else
	th1 = sp.thrln1 * sp.skyhw * sqrt (wts2);
    nent = 0;
    nsb = 0;

    if ((*first-1+LFWRD >= sp.kscan) || (*first-1+yfltsz >= ylen))
	focaserr (1,
	    "Insufficient number of lines in detection image region", "");

 /* set up initial values of sky array */
    for (nln = *first - 1; nln <= *first - 1 + LFWRD; nln++)
	rdimage (pfd, 0, nln, pthdr.naxis1, &pthdr, rcdptr[nln % inclns]);
    for (i = 0; i < sp.fy; i++)
	lnptrs[i] = rcdptr[(*first - 1 + i) % inclns];
    cmmnt ("detect");
    intsky (*minsky, *maxsky, SKY);
    wtcathdr (cfd, 0);

 /* Set sky header */
    sfd = 0;
    if (skfl[0] != 0) {
	skphdr = pthdr;
	sprintf (str, "Sky for %s", sp.pltnm);
	strncpy (skphdr.object, str, 31);
	skphdr.cdelt1 *= SKY;
	skphdr.cdelt2 *= SKY;
	skphdr.naxis1 = (skphdr.naxis1 + SKY - 1) / SKY;
	skphdr.naxis2 = (skphdr.naxis2 + SKY - 1) / SKY;
	skphdr.maxpxl = -MAXPXL;
	skphdr.minpxl = MAXPXL;
	sfd = imopen (skfl, &skphdr, 1);
    }

 /* main loop, segmentation program */
    for (yp = *first - 1 + yfltsz; yp < ylen; yp++) {
	rlc0[0] = rlc0[1] = -1;
	rlc0[2] = yp;
	rlc1[0] = rlc1[1] = -1;
	rlc1[2] = yp;
	if (nln < sp.kscan) {
	    rdimage (pfd, 0, nln, pthdr.naxis1, &pthdr,
		    rcdptr[nln % inclns]);
	    nln++;
	}
	for (i = 0; i < sp.fy; i++)
	    lnptrs[i] = rcdptr[(yp - yfltsz + i) % inclns];
	pixel (lnptrs);

    /* segmentation proceeds on x traversing loop */
	for (x = xfltsz, nn = ns = 0.; x < xlen; x++) {
	    j = x / SKY;
	    if (((x - xfltsz) % SKY) == 0) {
		if (sp.thrln0 < 0.)
		    a0 = th0;
		else
		    a0 = dsky[j] + th0;
		if (sp.thrln1 < 0.)
		    a1 = th1;
		else
		    a1 = dsky[j] - th1;
	    }

	/* Detection filter value */
	    b = smthln[x];

	    if (b > a0)		/* point is star - galaxy */
		if (rlc0[0] < 0)
		    rlc0[0] = rlc0[1] = x;
		else
		    rlc0[1] = x;/* right x value updated */

	    else
		if (b < a1)	/* point is dark */
		    if (rlc1[0] < 0)
			rlc1[0] = rlc1[1] = x;
		    else
			rlc1[1] = x;/* right x value updated */

		else {		/* point is sky */
		/* if previous point was in entity output it */
		    if (rlc0[0] >= 0) {
			add (rlc0, 0);
			rlc0[0] = -1;
		    }
		    if (rlc1[0] >= 0) {
			add (rlc1, 1);
			rlc1[0] = -1;
		    }
		/* update sky value */
		    nn++;
		    ns += b;
		}
	    if (((x - xfltsz) % SKY) == 0) {
		j = j > 0 ? j : 1;
		if (nn > sp.nupdte)
		    dsky[j] = sp.sb * (dsky[j + 1] + dsky[j - 1]) +
			sp.sa * ns / nn + (1.0 - sp.sa - sp.sb - sp.sb) *
			dsky[j];
		else
		    dsky[j] = sp.sb * (dsky[j + 1] + dsky[j - 1]) +
			(1.0 - sp.sb - sp.sb) * dsky[j];
		ns = nn = 0.;
	    }
	}
	dsky[j + 1] = dsky[j];
	dsky[0] = dsky[1];
	if (rlc0[0] >= 0)
	    add (rlc0, 0);
	if (rlc1[0] >= 0)
	    add (rlc1, 1);
	linend ();
	if (sfd > 0) {
	    cp = tmpbfr;
	    if (nsky == 0)
		for (i = 0; i < skphdr.naxis1; i++)
		    *cp++ = dsky[i];
	    else
		for (i = 0; i < skphdr.naxis1; i++)
		    *cp++ += dsky[i];
	    nsky++;
	    if ((yp % SKY) == SKY - 1) {
		cp = tmpbfr;
		skyot = 1 / wts / nsky;
		for (i = 0; i < skphdr.naxis1; i++, cp++)
		    *cp = *cp * skyot;
		wtimage (sfd, 0, skyline++, skphdr.naxis1, &skphdr, tmpbfr);
		pixrange (&skphdr, skphdr.naxis1, tmpbfr);
		nsky = 0;
	    }
	}
    }
    yp += sp.fy;
    linend ();

    if (sfd > 0) {
	if (nsky > 0) {
	    cp = tmpbfr;
	    skyot = 1 / wts / nsky;
	    for (i = 0; i < skphdr.naxis1; i++, cp++)
		*cp = *cp * skyot;
	    wtimage (sfd, 0, skyline++, skphdr.naxis1, &skphdr, tmpbfr);
	    pixrange (&skphdr, skphdr.naxis1, tmpbfr);
	}
	for (; skyline < skphdr.naxis2;)
	    wtimage (sfd, 0, skyline++, skphdr.naxis1, &skphdr, tmpbfr);
	skphdr.naxis2 = 0;
	skphdr.date[0] = 0;
	setimhdr (sfd, &skphdr);
	imclose (sfd);
    }
    imclose (pfd);
    fclose (afd);
    fclose (cfd);
}

pixel (lines)
PIXEL * lines[];
{
    register    PIXEL * ptr;
    register float *smthptr;
    register int    sum1, sum2, sum3, sum4;
    register int    i, j, k;
    int     imax;

    smthptr = smthln + xfltsz;
    imax = sp.kpts - 2 * xfltsz;
    if (sp.pflags & FLTR) {
	lines[0]++;
	lines[4]++;
	for (i = 0; i < imax; i++) {
	    ptr = lines[0]++;
	    sum1 = *ptr++;
	    sum2 = *ptr++;
	    sum1 += *ptr;
	    ptr = lines[1]++;
	    sum1 += *ptr++;
	    sum2 += *ptr++;
	    sum3 = *ptr++;
	    sum2 += *ptr++;
	    sum1 += *ptr;
	    ptr = lines[2]++;
	    sum2 += *ptr++;
	    sum3 += *ptr++;
	    sum4 = *ptr++;
	    sum3 += *ptr++;
	    sum2 += *ptr;
	    ptr = lines[3]++;
	    sum1 += *ptr++;
	    sum2 += *ptr++;
	    sum3 += *ptr++;
	    sum2 += *ptr++;
	    sum1 += *ptr;
	    ptr = lines[4]++;
	    sum1 += *ptr++;
	    sum2 += *ptr++;
	    sum1 += *ptr;
	    *smthptr++ = sum1 + (sum2 << 1) + 3 * sum3 + (sum4 << 2);
	}
	return;
    }
    for (i = 0; i < imax; i++) {
	*smthptr = 0.;
	for (j = 0; j < sp.fy; j++) {
	    ptr = lines[j]++;
	    for (k = 0; k < sp.fx; k++) {
		*smthptr += dfltr[j][k] * *ptr++;
	    }
	}
	smthptr++;
    }
}