/********************************/
/*	areas	valdes 8 30 82	*/
/*				*/
/*	Subroutines for areas	*/
/********************************/
#include	"focas1.h"


struct areas    area[ARMAX];
unsigned char  *silh;
int     xwmax, ywmax;
int     xfltsz, yfltsz;
float   xwavg, ywavg;

static int swap = 0;

/*Find areas in a directory search and return file descriptor*/
FILE *
fndar (arfl, mode)
char   *arfl;
int     mode;
{
    FILE    *fd;
    short   head, temp;


    switch (mode) {
    case 0:
    case 2:
        if ((fd = fopen (arfl, FOPEN_RW)) == NULL)
	    focaserr (1, "fndar: Can't open area file ", arfl);
	fseek (fd, 0L, 0);
	if (fread (&head, 2, 1, fd) != 1)
	    focaserr (1, "fndar: Can't read area file ", arfl);
	if (head == 0)
	    fseek (fd, 0L, 0);
	else if (head != 1) {
	    swap = 1;
	    swaps (&head, &temp);
	    if (temp != 1)
		focaserr (2, "fndar: Bad area file header in ", arfl);
	}
	break;
    case 1:
        if ((fd = fopen (arfl, FOPEN_WO)) == NULL) {
	    c_delete (arfl);
	    if ((fd = fopen (arfl, FOPEN_WO)) == NULL) {
		focaserr (3, "fndar: Can't create ", arfl);
		return (NULL);
	    }
	}
	if (swap) {
	    temp = 1;
	    swaps (&temp, &head);
	} else
	    head = 1;
	fwrite (&head, 2, 1, fd);
	break;
    }

    return (fd);
}

/*Make detection area within rct from ap*/
long
detect (ap, rct, a)
struct smblk   *ap;
short   rct[4];
struct areas   *a;
{
    long    ar;
    short   xa, xb, xc, xd, y, nrlc;
    short   ytop, ybot;
    struct smblk   *lp;

    a -> aflgs &= ~SIZE;
    a -> xmin = MAXBUF;
    a -> xmax = 0;
    nrlc = 0;
    ar = 0;
    for (lp = ap -> smb2.lnp -> smb1.b; lp != 0; lp = lp -> smb1.b) {
	xa = lp -> smb2.xl;
	xb = lp -> smb3.xr;
	y = lp -> smb4.y;
	if ((y < rct[1]) || (y > rct[3]))
	    continue;
	if ((xb < rct[0]) || (xa > rct[2]))
	    continue;
	if (nrlc == 0)
	    a -> d.ychst = ytop = y;
	a -> d.xla[nrlc] = xc = max (rct[0], xa);
	a -> d.xra[nrlc] = xd = min (rct[2], xb);
	a -> d.nx[ytop - y] = nrlc + 1;
	ar += xd - xc + 1;
	if (xc < a -> xmin)
	    a -> xmin = xc;	/* min x of entity */
	if (xd > a -> xmax)
	    a -> xmax = xd;	/* max x of entity */
	ybot = y;
	nrlc++;
	if (nrlc == MXRLC)
	    break;
    }
    a -> rct[0] = max (a -> xmin - sp.buf, xfltsz);
    a -> rct[1] = max (ybot - sp.buf, yfltsz);
    a -> rct[2] = min (a -> xmax + sp.buf, sp.kpts - 1 - xfltsz);
    a -> rct[3] = min (ytop + sp.buf, sp.kscan - 1 - yfltsz);
    a -> rct[0] = max (a -> rct[0], rct[0]);
    a -> rct[1] = max (a -> rct[1], rct[1]);
    a -> rct[2] = min (a -> rct[2], rct[2]);
    a -> rct[3] = min (a -> rct[3], rct[3]);
    a -> d.nlines = ytop - ybot + 1;
    if (nrlc == MXRLC)
	a -> aflgs |= SIZE;
    if ((ytop > (sp.kscan - yfltsz - 2)) ||
	    (ybot < (yfltsz + 1)) ||
	    (a -> xmin < (xfltsz + 3)) ||
	    (a -> xmax > (sp.kpts - xfltsz - 4))
	)
	a -> aflgs |= EDGE;
    if (a -> d.nlines < 1)
	focaserr (4, "detect: Error in detection area", "");
    return (ar);
}

/*Find isophotal boundary*/
isophote (a)
struct areas   *a;
{
    int     i, j, k, l;
    int     dx, dy, x0, xl, xr, begin, ar;
    unsigned    size;
    struct arblk   *d, *is;
    char   *fvalloc ();
    unsigned char  *cp;

    d = &a -> d;
    is = &a -> i;
    x0 = a -> xmin - 1;
    dx = a -> xmax + 2 - x0;
    dy = d -> nlines + 2;
    size = dx * dy;
    if (size > MXSZ) {
	a -> aflgs |= SIZE;
	fprintf (stderr, "isophote: Image too large\n");
	fprintf (stderr, "xmin = %d xmax = %d ychst = %d nlines = %d\n",
		a -> xmin, a -> xmax, d -> ychst, d -> nlines);
	return (1);
    }
    if (silh == 0) {
	silh = (unsigned char *) fvalloc (MXSZ, sizeof (char));
	if (silh == 0)
	    focaserr (4, "isophote: Not enough memory", "");
    }
 /* set up silhouette */
    ar = 0;
    cp = silh;
    for (i = 0; i < size; i++)
	cp[i] = 0;
    for (i = 0, l = 0; i < d -> nlines; i++) {
	cp += dx;
	for (; l < d -> nx[i]; l++) {
	    xl = d -> xla[l] - x0;
	    xr = d -> xra[l] - x0;
	    for (k = xr; k >= xl; k--) {
		cp[k] = 1;
		ar++;
	    }
	}
    }

    begin = d -> xra[0] - x0 + dx;
    silh[begin] = 2;
    if (ar > 1) {
	for (i = begin, j = 0, k = 0; j != begin; i = j, k = (k + 6) % 8) {
	    for (;; k = (++k) % 8) {
		switch (k) {
		    case 0: 
			j = i - dx;
			break;
		    case 1: 
			j = i - dx + 1;
			break;
		    case 2: 
			j = i + 1;
			break;
		    case 3: 
			j = i + dx + 1;
			break;
		    case 4: 
			j = i + dx;
			break;
		    case 5: 
			j = i + dx - 1;
			break;
		    case 6: 
			j = i - 1;
			break;
		    case 7: 
			j = i - dx - 1;
			break;
		}
		if (silh[j] != 0)
		    break;
	    }
	    silh[j] = 2;
	}
    }
    is -> nlines = d -> nlines;
    is -> ychst = d -> ychst;
    cp = silh;
    for (i = 0, l = 0; i < is -> nlines; i++) {
	cp += dx;
	for (j = dx - 2; j > 0; j--) {
	    if (cp[j] == 2) {
		is -> xra[l] = j + x0;
		for (; cp[j] == 2; j--);
		is -> xla[l] = j + 1 + x0;
		is -> nx[i] = ++l;
		if (l == MXRLC) {
		    a -> aflgs |= SIZE;
		    fprintf (stderr, "isophote:  Warning too many rlcs\n");
		    return (1);
		}
	    }
	}
    }
    return (0);
}

imisophotes (cfd, afd, imiso, grow)
FILE	*cfd;
FILE	*afd;
struct imiso *imiso;
float	grow;
{
    struct objrec   ob;
    struct areas    a;
    struct arblk    *iso;
    struct liso *liso;
    char	*calloc(), *malloc(), *realloc();
    short	i, j, k, l, n, x0, y0, nl, nc, nbytes;
    long   totarea();

    x0 = imiso -> x0;
    y0 = imiso -> y0;
    nc = imiso -> nx;
    nl = imiso -> ny;

    imiso -> liso = (struct liso *) calloc (nl, sizeof(struct liso));
    if (imiso -> liso == 0)
	focaserr (1, "Memory allocation error", ""); 

    fseek (cfd, szcathdr(), 0);
    for (;;) {
	if (rdcatob (cfd, 0L, &ob))
	    break;
	if (ffilter (&ob) != 1)
	    continue;
	if (rdarea (afd, ob.entnum, ob.subent, ob.arpos, &a))
	    continue;

	if (grow <= 1.)
	    iso = &a.d;
	else {
	    (void) totarea (&a, grow);
	    iso = &a.t;
	}

        for (i = 0, j = iso->ychst, k = 0; i < iso->nlines; i++, j--) {
	    l = j - y0;
	    if (l < 0 || l >= nl) {
		k += iso->nx[i];
		continue;
	    }
	    liso = &imiso -> liso[l];
	    if (liso -> nalloc < liso -> nused + iso->nx[i]) {
		liso -> nalloc = liso -> nused + iso->nx[i] + 20;
		nbytes = liso -> nalloc * sizeof (short);
		if (liso -> xla == 0) {
		    liso -> xla = (short *) malloc (nbytes);
		    liso -> xra = (short *) malloc (nbytes);
		} else {
		    liso -> xla = (short *) realloc (liso -> xla, nbytes);
		    liso -> xra = (short *) realloc (liso -> xra, nbytes);
		}
		if (liso -> xla == 0 || liso -> xra == 0)
		    focaserr (1, "Memory allocation error", ""); 
	    }
	    n = liso -> nused;
	    for (; k < iso->nx[i]; k++) {
		if (iso->xra[k] < x0 || iso->xla[k] >= x0 + nc)
		    continue;
		liso -> xla[n] = iso->xla[k];
		liso -> xra[n] = iso->xra[k];
		n++;
	    }
	    liso -> nused = n;
	}
    }
}

long
totarea (a, grow)
struct areas   *a;
float   grow;
{
    int     i, j, k, l, xmin, xmax;
    long    ar, ar0;
    struct arblk   *d, *t;

    d = &a -> d;
    t = &a -> t;
    xmin = sp.kpts + 1;
    xmax = -1;
    j = k = 0;

    /* Remove concavities */
    t -> ychst = d -> ychst;
    t -> nlines = d -> nlines;
    for (i = 0; i < t -> nlines; i++) {
	 t -> nx[i] = i + 1;
	 t -> xla[i] = d -> xla[d -> nx[i] - 1];
	 t -> xra[i] = (i == 0 ? d -> xra[i] : d -> xra[d -> nx[i-1]]);
	if (t  -> xla[i] < xmin)
	    xmin = t -> xla[j = i];
	if (t -> xra[i] > xmax)
	    xmax = t -> xra[k = i];
    }
    do {
	for (i = j - 1, l = 0; i > 0; i--)
	    if (t -> xla[i] > t -> xla[i-1]) {
		t -> xla[i] = t -> xla[i-1];
		l++;
	    }
    } while (l != 0);
    do {
	for (i = j + 1, l = 0; i < t -> nlines - 1; i++)
	    if (t -> xla[i] > t -> xla[i+1]) {
		t -> xla[i] = t -> xla[i+1];
		l++;
	    }
    } while (l != 0);
    do {
	for (i = k - 1, l = 0; i > 0; i--)
	    if (t -> xra[i] < t -> xra[i-1]) {
		t -> xra[i] = t -> xra[i-1];
		l++;
	    }
    } while (l != 0);
    do {
	for (i = k + 1, l = 0; i < t -> nlines - 1; i++)
	    if (t -> xra[i] < t -> xra[i+1]) {
		t -> xra[i] = t -> xra[i+1];
		l++;
	    }
    } while (l != 0);

    /* Evaluate area */
    for (ar0 = 0, i = 0; i < t -> nlines; i++)
	ar0 += t -> xra[i] - t -> xla[i] + 1;

    /* Grow boundary with at least one additional ring */
    do {
	for (i = t -> nlines - 1; i >= 0; i--) {
	    t -> xla[i+1] = t -> xla[i] - 1;
	    t -> xra[i+1] = t -> xra[i] + 1;
	}
	t -> xla[0] = t -> xla[1] + 1;
	t -> xra[0] = t -> xra[1] - 1;
	t -> xla[t -> nlines+1] = t -> xla[t -> nlines] + 1;
	t -> xra[t -> nlines+1] = t -> xra[t -> nlines] - 1;
	j++;
	k++;
	t -> ychst++;
	t -> nlines += 2;
	t -> nx[t -> nlines-2] = t -> nlines - 1;
	t -> nx[t -> nlines-1] = t -> nlines;
	for (i = j - 1; i > 0; i--)
	    if (t -> xla[i] > t -> xla[i+1] + 1)
		t -> xla[i] = t -> xla[i+1] + 1;
	for (i = j + 1; i < t -> nlines - 1; i++)
	    if (t -> xla[i] > t -> xla[i-1] + 1)
		t -> xla[i] = t -> xla[i-1] + 1;
	for (i = k - 1; i > 0; i--)
	    if (t -> xra[i] < t -> xra[i+1] - 1)
		t -> xra[i] = t -> xra[i+1] - 1;
	for (i = k + 1; i < t -> nlines - 1; i++)
	    if (t -> xra[i] < t -> xra[i-1] - 1)
		t -> xra[i] = t -> xra[i-1] - 1;
	for (ar = 0, i = 0; i < t -> nlines; i++)
	    ar += t -> xra[i] - t -> xla[i] + 1;
    } while (ar < grow * ar0);

    /* If shrinking remove rings but require at list one pixel be left */
    if (grow <= 1. && t -> nlines > 2) {
	do {
	    for (i = 0; i < t -> nlines - 2; i++) {
		if (t -> xra[i] - t -> xla[i] + 1 > 2) {
		    t -> xla[i] = t -> xla[i+1] + 1;
		    t -> xra[i] = t -> xra[i+1] - 1;
		}
	    }
	    j--;
	    k--;
	    t -> ychst--;
	    t -> nlines -= 2;
	    for (ar = 0, i = 0; i < t -> nlines; i++)
		ar += t -> xra[i] - t -> xla[i] + 1;
	} while (ar > grow * ar0 && t -> nlines > 2);
	/* Add one ring back */
	do {
	    for (i = t -> nlines - 1; i >= 0; i--) {
		t -> xla[i+1] = t -> xla[i] - 1;
		t -> xra[i+1] = t -> xra[i] + 1;
	    }
	    t -> xla[0] = t -> xla[1] + 1;
	    t -> xra[0] = t -> xra[1] - 1;
	    t -> xla[t -> nlines+1] = t -> xla[t -> nlines] + 1;
	    t -> xra[t -> nlines+1] = t -> xra[t -> nlines] - 1;
	    j++;
	    k++;
	    t -> ychst++;
	    t -> nlines += 2;
	    t -> nx[t -> nlines-2] = t -> nlines - 1;
	    t -> nx[t -> nlines-1] = t -> nlines;
	    for (ar = 0, i = 0; i < t -> nlines; i++)
		ar += t -> xra[i] - t -> xla[i] + 1;
	} while (ar < grow * ar0);
    }
    return (ar);
}

/*writes area to file descriptor fd and returns the offset*/
/*to beginning of the write*/
long
wrarea (fd, ob, a)
FILE    *fd;
struct objrec  *ob;
struct areas   *a;
{
    short   i, stemp;
    long    afpos, ltemp;
    struct arblk   *d;
    register short *sptr1, *sptr2;

    d = &a -> d;

    arorg (a, sp.xs, sp.ys);
    fseek (fd, 0L, 2);
    afpos = ftell (fd);

    if (swap) {
	sptr1 = &stemp;
	swaps (&ob -> entnum, sptr1);
	fwrite (sptr1, 2, 1, fd);
	swapl (&ob -> subent, &ltemp);
	fwrite (&ltemp, 4, 1, fd);

	sptr2 = a -> rct;
	swaps (sptr2++, sptr1);
	fwrite (sptr1, 2, 1, fd);
	swaps (sptr2++, sptr1);
	fwrite (sptr1, 2, 1, fd);
	swaps (sptr2++, sptr1);
	fwrite (sptr1, 2, 1, fd);
	swaps (sptr2++, sptr1);
	fwrite (sptr1, 2, 1, fd);
	swaps (&a -> xmin, sptr1);
	fwrite (sptr1, 2, 1, fd);
	swaps (&a -> xmax, sptr1);
	fwrite (sptr1, 2, 1, fd);

	swaps (&d -> ychst, sptr1);
	fwrite (sptr1, 2, 1, fd);
	swaps (&d -> nlines, sptr1);
	fwrite (sptr1, 2, 1, fd);

	sptr2 = d -> nx;
	for (i = 0; i < d -> nlines; i++) {
	    swaps (sptr2++, sptr1);
	    fwrite (sptr1, 2, 1, fd);
	}
	sptr2 = d -> xla;
	for (i = 0; i < d -> nx[d -> nlines-1]; i++) {
	    swaps (sptr2++, sptr1);
	    fwrite (sptr1, 2, 1, fd);
	}
	sptr2 = d -> xra;
	for (i = 0; i < d -> nx[d -> nlines-1]; i++) {
	    swaps (sptr2++, sptr1);
	    fwrite (sptr1, 2, 1, fd);
	}
    } else {
        fwrite (&ob -> entnum, 2, 1, fd);
        fwrite (&ob -> subent, 4, 1, fd);

        fwrite (a -> rct, 2, 4, fd);
        fwrite (&a -> xmin, 2, 1, fd);
        fwrite (&a -> xmax, 2, 1, fd);

        fwrite (&d -> ychst, 2, 1, fd);
        fwrite (&d -> nlines, 2, 1, fd);
        fwrite (&d -> nx[0], 2, d -> nlines, fd);
        fwrite (&d -> xla[0], 2, d -> nx[d -> nlines - 1], fd);
        fwrite (&d -> xra[0], 2, d -> nx[d -> nlines - 1], fd);
    }

    return (afpos);
}

/* Change the area descriptor origin */
arorg (a, xs, ys)
struct areas   *a;
int     xs, ys;
{
    int     i;
    struct arblk   *d;

    d = &a -> d;

    a -> rct[0] += xs;
    a -> rct[1] += ys;
    a -> rct[2] += xs;
    a -> rct[3] += ys;
    a -> xmin += xs;
    a -> xmax += xs;
    d -> ychst += ys;
    for (i = 0; i < d -> nx[d -> nlines - 1]; i++) {
	d -> xla[i] += xs;
	d -> xra[i] += xs;
    }
}

/*reads area from file descriptor fd starting at offset afpos*/
rdarea (fd, ent, sbent, arpos, a)
FILE    *fd;
short   ent;
long    sbent;
long    arpos;
struct areas   *a;
{
    int	    subent, ltemp;
    struct  arblk *d;
    short   i, entnum, oldsubent, stemp;
    register short *sptr1, *sptr2;

    d = &a -> d;

    if (fseek (fd, arpos, 0) == -1) {
	fprintf (stderr, "rdarea: Error on seek\n");
	return (1);
    }

    if (swap) {
	sptr1 = &stemp;
        fread (sptr1, 2, 1, fd);
	swaps (sptr1, &entnum);
        fread (&ltemp, 4, 1, fd);
	swapl (&ltemp, &subent);
        if ((entnum != ent) || (subent != sbent)) {
	    fseek (fd, -4, 1);
	    fread (sptr1, 2, 1, fd);
	    swaps (sptr1, &oldsubent);
	    if ((entnum != ent) || (oldsubent != sbent)) {
	        fprintf (stderr, "rdarea: Error on read\n");
	        fprintf (stderr, "  catalog: %d.%-d   areas: %d.%-d\n",
		        ent, sbent, entnum, subent);
	        return (1);
	    }
        }

	sptr2 = a -> rct;
        fread (sptr1, 2, 1, fd);
	swaps (sptr1, sptr2++);
        fread (sptr1, 2, 1, fd);
	swaps (sptr1, sptr2++);
        fread (sptr1, 2, 1, fd);
	swaps (sptr1, sptr2++);
        fread (sptr1, 2, 1, fd);
	swaps (sptr1, sptr2++);
        fread (sptr1, 2, 1, fd);
	swaps (sptr1, &a -> xmin);
        fread (sptr1, 2, 1, fd);
	swaps (sptr1, &a -> xmax);

        fread (sptr1, 2, 1, fd);
	swaps (sptr1, &d -> ychst);
        fread (sptr1, 2, 1, fd);
	swaps (sptr1, &d -> nlines);

	sptr2 = d -> nx;
	for (i=0; i < d -> nlines; i++) {
            fread (sptr1, 2, 1, fd);
	    swaps (sptr1, sptr2++);
	}
	sptr2 = d -> xla;
	for (i=0; i < d -> nx[d -> nlines-1]; i++) {
            fread (sptr1, 2, 1, fd);
	    swaps (sptr1, sptr2++);
	}
	sptr2 = d -> xra;
	for (i=0; i < d -> nx[d -> nlines-1]; i++) {
            fread (sptr1, 2, 1, fd);
	    swaps (sptr1, sptr2++);
	}
    } else {
        fread (&entnum, 2, 1, fd);
        fread (&subent, 4, 1, fd);
        if ((entnum != ent) || (subent != sbent)) {
	    fseek (fd, -4, 1);
	    fread (&oldsubent, 2, 1, fd);
	    if ((entnum != ent) || (oldsubent != sbent)) {
	        fprintf (stderr, "rdarea: Error on read\n");
	        fprintf (stderr, "  catalog: %d.%-d   areas: %d.%-d\n",
		        ent, sbent, entnum, subent);
	        return (1);
	    }
        }

        fread (a -> rct, 2, 4, fd);
        fread (&a -> xmin, 2, 1, fd);
        fread (&a -> xmax, 2, 1, fd);

        fread (&d -> ychst, 2, 1, fd);
        fread (&d -> nlines, 2, 1, fd);
        fread (&d -> nx[0], 2, d -> nlines, fd);
        fread (&d -> xla[0], 2, d -> nx[d -> nlines - 1], fd);
        fread (&d -> xra[0], 2, d -> nx[d -> nlines - 1], fd);
    }

    arorg (a, -sp.xs, -sp.ys);
    (void) isophote (a);
    (void) totarea (a, sp.grow);
    return (0);
}

/*Checks area id from file descriptor fd starting at offset afpos*/
ckarea (fd, ob)
FILE    *fd;
struct objrec  *ob;
{
    int	    subent, ltemp;
    short   i, entnum, oldsubent, stemp;
    register short *sptr1;


    if (fseek (fd, ob -> arpos, 0) == -1) {
	fprintf (stderr, "ckarea: Error on seek\n");
	return (1);
    }

    if (swap) {
	sptr1 = &stemp;
        fread (sptr1, 2, 1, fd);
	swaps (sptr1, &entnum);
        fread (&ltemp, 4, 1, fd);
	swapl (&ltemp, &subent);
        if ((entnum != ob -> entnum) || (subent != ob -> subent)) {
	    fseek (fd, -4, 1);
	    fread (sptr1, 2, 1, fd);
	    swaps (sptr1, &oldsubent);
	    if ((entnum != ob -> entnum) || (oldsubent != ob -> subent)) {
	        fprintf (stderr, "ckarea: Error on read\n");
	        fprintf (stderr, "  catalog: %d.%-d   areas: %d.%-d\n",
		        ob -> entnum, ob -> subent, entnum, subent);
	        return (1);
	    }
        }
    } else {
        fread (&entnum, 2, 1, fd);
        fread (&subent, 4, 1, fd);
        if ((entnum != ob -> entnum) || (subent != ob -> subent)) {
	    fseek (fd, -4, 1);
	    fread (&oldsubent, 2, 1, fd);
	    if ((entnum != ob -> entnum) || (oldsubent != ob -> subent)) {
	        fprintf (stderr, "ckarea: Error on read\n");
	        fprintf (stderr, "  catalog: %d.%-d   areas: %d.%-d\n",
		        ob -> entnum, ob -> subent, entnum, subent);
	        return (1);
	    }
        }
    }
    return (0);
}

dmparea (a)
struct areas   *a;
{
    int     i, j;
    struct arblk   *d, *is, *t;


    d = &a -> d;
    is = &a -> i;
    t = &a -> t;
    fprintf (stderr, "  rct = %4d %4d %4d %4d\n",
	    a -> rct[0], a -> rct[1], a -> rct[2], a -> rct[3]);
    fprintf (stderr, "  min = %4d max = %4d\n", a -> xmin, a -> xmax);
    fprintf (stderr, "  d: ychst = %4d   nlines = %4d\n", d -> ychst, d -> nlines);
    for (i = 0, j = 0; i < d -> nlines; i++) {
	fprintf (stderr, "      i = %2d  xl = %4d   xr = %4d\n", i, d -> xla[j], d -> xra[j]);
	for (j++; j < d -> nx[i]; j++)
	    fprintf (stderr, "             xl = %4d   xr = %4d\n", d -> xla[j], d -> xra[j]);
    }
    fprintf (stderr, "  i: ychst = %4d   nlines = %4d\n", is -> ychst, is -> nlines);
    for (i = 0, j = 0; i < is -> nlines; i++) {
	fprintf (stderr, "      i = %2d  xl = %4d   xr = %4d\n", i, is -> xla[j], is -> xra[j]);
	for (j++; j < is -> nx[i]; j++)
	    fprintf (stderr, "             xl = %4d   xr = %4d\n", is -> xla[j], is -> xra[j]);
    }
    fprintf (stderr, "  t: ychst = %4d   nlines = %4d\n", t -> ychst, t -> nlines);
    for (i = 0, j = 0; i < t -> nlines; i++) {
	fprintf (stderr, "      i = %2d  xl = %4d   xr = %4d\n", i, t -> xla[j], t -> xra[j]);
	for (j++; j < t -> nx[i]; j++)
	    fprintf (stderr, "             xl = %4d   xr = %4d\n", t -> xla[j], t -> xra[j]);
    }
}

min (r, s)
int     r, s;
{
    return (r < s ? r : s);
}

max (r, s)
int     r, s;
{
    return (r > s ? r : s);
}