/********************************/
/*	msplit	valdes 8 30 82	*/
/*				*/
/*	Split images by		*/
/*	increasing isophotes	*/
/********************************/
#include	"focas1.h"
/*
#include	"iisfv.h"
*/

struct objrec   ob[ARMAX];
struct areas    area[ARMAX];
int     nent, nsb;
int     yp;

static struct image p;
static int  blank;

msplit (o, a, nlevels, diag)
struct objrec  *o;
struct areas   *a;
int     nlevels;
int     diag;
{
    int     i, n;
    int     x, y, rlc[4], d, dmin, dmax;
    float   sum, Ltotal;
    short   rct[4];
    PIXEL * cp, step;

 /* Test for objects not to be split */

    if (o -> eflgs & DARK)
	return (0);
    if ((o -> eflgs & EVAL) == 0)
	return (0);
    if (o -> eflgs & SIZE)
	return (0);

 /* Find the smallest rectangle containing the isophote and read the image */

    rct[0] = a -> xmin;
    rct[1] = a -> d.ychst - (a -> d.nlines - 1);
    rct[2] = a -> xmax;
    rct[3] = a -> d.ychst;
    if (rdimg (pfd, rct, &p) < 0)
	return (0);

 /* Set the exterior pixels to blank and find minimum density */

    blank = p.min - 1;
    nent = o -> entnum;
    nsb = 1;

 /* Determine the density step and limits */

    step = sp.sptbw * sp.skyhw;
    if (step < 1)
	step = 1;
    dmin = interior (a, &p, blank);
    dmax = dmin + nlevels * sp.skyhw;

 /* Step through the increasing thresholds */

    for (d = dmin; d <= dmax; d += step) {
    /* Run detector */
	for (y = p.y0; y < p.y0 + p.dy; y++) {
	    yp = y;
	    rlc[0] = rlc[1] = -1;
	    rlc[2] = yp;
	    cp = p.pic + (y - p.y0) * p.dx;
	    for (x = p.x0; x < p.x0 + p.dx; x++) {
		if (cp[x - p.x0] > d)
		    if (rlc[0] < 0)
			rlc[0] = rlc[1] = x;
		    else
			rlc[1] = x;
		else {
		    if (rlc[0] >= 0) {
			add (rlc, 0);
			rlc[0] = -1;
		    }
		    cp[x - p.x0] = blank;
		}
	    }
	    if (rlc[0] >= 0)
		add (rlc, 0);
	    linend ();
	}
	yp += 5;
	linend ();

	/*
	if (diag)
	    dspspt (&p);
	*/

    /* Object has disappeared */
	if (nsb == 1)
	    break;

    /* Successful split */
	if (nsb > 3)
	    break;

    /* Possible successful split */
	if (nsb == 3) {
	    if (evlobj (&ob[1], &area[1], 0, stdout) == 0) {
		strcpy (o -> class, "m");
		break;
	    }
	}

    /* Only original object detected */
	nsb = 1;
    }

    if (nsb == 1)
	o -> eflgs |= SNGL;
    else {
	strcpy (o -> class, "m");
	sum = 0.;
	for (i = 1; i < nsb; i++) {
	    ob[i].entnum = ob[0].entnum;
	    ob[i].subent = 10 * ob[0].subent + i;
	    strcpy (ob[i].class, "u");
	    ob[i].scale = ob[i].frac = 0.;
	    sum += ob[i].Li;
	}
    /* Luminosity estimate */
	if (o -> subent == 0)
	    /* Ltotal = o -> Li; */
	    Ltotal = o->Ltotal;
	else
	    Ltotal = o -> Ltotal;
	for (i = 1; i < nsb; i++) {
	    ob[i].Ltotal = Ltotal * ob[i].Li / sum;
	}
    }
    n = nsb - 1 > 1 ? nsb - 1 : 0;
    return (n);
}

sptevl (diag)
int     diag;
{
    int     j, k, l, yi, xl, xr;
    PIXEL * cp;

    if (ob[nsb].area < sp.minarea)
	goto remove;
    if (nsb == 1) {		/* Just check area */
	return (0);
    } else if (nsb == 2) {		/* Evaluate smallest area object */
	if (ob[1].area < ob[2].area) {
	    /*
	    ob[3] = ob[2];
	    area[3] = area[2];
	    ob[2] = ob[1];
	    area[2] = area[1];
	    ob[1] = ob[3];
	    area[1] = area[3];
	    */
	    bytecopy (&ob[2], &ob[3], sizeof (struct objrec));
	    bytecopy (&ob[1], &ob[2], sizeof (struct objrec));
	    bytecopy (&ob[3], &ob[1], sizeof (struct objrec));
	    bytecopy (&area[2], &area[3], sizeof (struct areas));
	    bytecopy (&area[1], &area[2], sizeof (struct areas));
	    bytecopy (&area[3], &area[1], sizeof (struct areas));
	}
	if (evlobj (&ob[2], &area[2], diag, stdout) == 0)
	    return (0);
	else
	    goto remove;
    } else {			/* Evaluate all */
	if ((ob[1].eflgs & EVAL) == 0) {/* Evaluate smallest area */
	    if (ob[1].area < ob[nsb].area) {
		/*
		ob[4] = ob[3];
		area[4] = area[3];
		ob[3] = ob[1];
		area[3] = area[1];
		ob[1] = ob[4];
		area[1] = area[4];
		*/
	        bytecopy (&ob[3], &ob[4], sizeof (struct objrec));
	        bytecopy (&ob[1], &ob[3], sizeof (struct objrec));
	        bytecopy (&ob[4], &ob[1], sizeof (struct objrec));
	        bytecopy (&area[3], &area[4], sizeof (struct areas));
	        bytecopy (&area[1], &area[3], sizeof (struct areas));
	        bytecopy (&area[4], &area[1], sizeof (struct areas));
	    }
	    if (evlobj (&ob[3], &area[3], diag, stdout) == 1)
		goto remove;
	    /*
	    ob[4] = ob[3];
	    area[4] = area[3];
	    ob[3] = ob[1];
	    area[3] = area[1];
	    ob[1] = ob[4];
	    area[1] = area[4];
	    */
	    bytecopy (&ob[3], &ob[4], sizeof (struct objrec));
	    bytecopy (&ob[1], &ob[3], sizeof (struct objrec));
	    bytecopy (&ob[4], &ob[1], sizeof (struct objrec));
	    bytecopy (&area[3], &area[4], sizeof (struct areas));
	    bytecopy (&area[1], &area[3], sizeof (struct areas));
	    bytecopy (&area[4], &area[1], sizeof (struct areas));
	}
	if (evlobj (&ob[nsb], &area[nsb], diag, stdout) == 0)
	    return (0);
    }

remove: 
    for (j = p.dy - 1, k = 0; j >= 0; j--) {
	yi = area[nsb].d.ychst - p.y0 - j;
	if ((yi >= 0) && (yi < area[nsb].d.nlines)) {
	    cp = p.pic + j * p.dx;
	    for (; k < area[nsb].d.nx[yi]; k++) {
		xl = area[nsb].d.xla[k] - p.x0;
		xr = area[nsb].d.xra[k] - p.x0;
		for (l = xl; l <= xr; l++)
		    cp[l] = blank;
	    }
	}
    }
    return (1);
}

bytecopy (a, b, len)
char	*a, *b;
int	len;
{
    register char *ap, *bp;
    register int i;

    ap = a;
    bp = b;
    for (i=0; i<len; i++)
	*bp++ = *ap++;
}