/********************************/
/*	rdfield	valdes 1 7 83	*/
/*				*/
/*	Read field image data	*/
/********************************/
#include	"focas1.h"
int     inclns;
PIXEL * rcdptr[LFWRD + LBACK + 1];

/*Find field in a directory search and return file descriptor*/
fndfld (ptfl, mode)
char   *ptfl;
int     mode;
{
    int     fd;
    char    imfile[40];

    if ((fd = imopen (ptfl, &pthdr, mode)) > 0)
	return (fd);
    /*
    if (access (ptfl, 0) == 0) {
	if ((fd = imopen (ptfl, &pthdr, mode)) > 0)
	    return (fd);
	focaserr (1, "fndfld: Can't access image file ", ptfl);
    }
    strcpy (imfile, ptfl);
    strcat (imfile, ".imh");
    if (access (imfile, 0) == 0) {
	if ((fd = imopen (imfile, &pthdr, mode)) > 0)
	    return (fd);
	focaserr (2, "fndfld: Can't access image file ", imfile);
    }
    strcpy (imfile, ptfl);
    strcat (imfile, ".fits");
    if (access (imfile, 0) == 0) {
	if ((fd = imopen (imfile, &pthdr, mode)) > 0)
	    return (fd);
	focaserr (2, "fndfld: Can't access image file ", imfile);
    }
    */
    fprintf (stderr, "fndfld: Can't access image file %s\n", ptfl);
    fflush (stderr);
    for (;;) {
	fprintf (stderr, "Enter image filename (<cr> to exit): ");
	fflush (stderr);
	if (strlen (gets (imfile))) {
	    if ((fd = imopen (imfile, &pthdr, mode)) > 0)
		return (fd);
	    else
		fprintf (stderr, "fndfld: Can't access image file %s\n", imfile);
	} else
	    focaserr (1, "No image", "");
    }
}

/*Read an image within rectangle rct and return pointer to image structure*/
/*Return: - 1 = error, 0 = new image, 1 = image unchanged*/
rdimg (fd, rct, p)
int     fd;
short   rct[4];
struct image   *p;
{
    register    PIXEL * ptr;
    int     i, j, k;
    char   *fvalloc ();

 /* Check desired rectangle */
    if ((rct[0] < 0) || (rct[1] < 0) || (rct[2] < 0) || (rct[3] < 0))
	return (-1);
    if ((rct[0] >= pthdr.naxis1) || (rct[1] >= pthdr.naxis2) ||
	    (rct[2] >= pthdr.naxis1) || (rct[3] >= pthdr.naxis2))
	return (-1);
    if ((rct[0] > rct[2]) || (rct[1] > rct[3]))
	return (-1);

 /* Check if last image is the same */
    if ((p->x0 == rct[0]) && (p->y0 == rct[1]) &&
	    (p->dx == rct[2] - rct[0] + 1) && (p->dy == rct[3] - rct[1] + 1))
	return (1);

    p->x0 = rct[0];
    p->dx = rct[2] - p->x0 + 1;
    p->y0 = rct[1];
    p->dy = rct[3] - p->y0 + 1;
    if (p->pic == 0) {
	p->mxsz = pthdr.naxis1 * pthdr.naxis2;
	p->mxsz = MXSZ < p->mxsz ? MXSZ : p->mxsz;
	for (i = p->mxsz; p->pic == 0 && i > 0; i *= .1) {
	    p->mxsz = i;
	    p->pic = (PIXEL *) fvalloc (p->mxsz, sizeof (PIXEL));
	}
    }
    if (p->dx * p->dy > p->mxsz) {
	fprintf (stderr, "rdimg: Image too large\n");
	return (-1);
    }
    p->min = MAXPXL;
    p->max = -MAXPXL;
    for (i = p->y0, j = 0; i < p->y0 + p->dy; i++, j += p->dx) {
	ptr = p->pic + j;
	rdimage (fd, p->x0, i, p->dx, &pthdr, ptr);
	for (k = 0; k < p->dx; k++) {
	    p->min = p->min < *ptr ? p->min : *ptr;
	    p->max = p->max > *ptr ? p->max : *ptr;
	    ptr++;
	}
    }
    return (0);
}

/*Return minimum density within isophote*/
interior (a, p, blank)
struct areas   *a;
struct image   *p;
int     blank;
{
    int     i, j, k, nx, yi, xl, xr;
    PIXEL * cp, dmin;

    dmin = sp.satr;
    for (i = p->dy - 1, j = 0; i >= 0; i--) {
	cp = p->pic + i * p->dx;
	yi = a->d.ychst - p->y0 - i;
	if ((yi >= 0) && (yi < a->d.nlines)) {
	    nx = yi == 0 ? a->d.nx[yi] : a->d.nx[yi] - a->d.nx[yi - 1];
	    for (k = p->dx - 1; j < a->d.nx[yi]; j++) {
		xl = a->d.xla[j] - p->x0;
		xr = a->d.xra[j] - p->x0;
		for (; k > xr; k--)
		    cp[k] = blank;
		for (; k >= xl; k--) {
		    if ((cp[k] < dmin) && (cp[k] > blank))
			dmin = cp[k];
		}
		if (nx == 1)
		    for (; k >= 0; k--)
			cp[k] = blank;
		nx--;
	    }
	}
	else
	    for (k = 0; k < p->dx; k++)
		cp[k] = blank;
    }
    return ((int) dmin);
}

/*replace an image with sky*/
imgsub (a, p, grow)
struct areas   *a;
struct image   *p;
float   grow;
{
    int     i, j, k, nx, yi, xl, xr, peak, max, dmin, dmax;
    struct arblk   *ar;
    int     nhstgrm, hstgrm[512];
    PIXEL * cp, start, step, mean1, mean2, mean;
    float   bin;

    if (meanskyd (p, &a->d, &mean1) < 1)
	return (1);

    if (meanisphtd (p, &a->i, &mean2) < 1)
	return (1);;

    if (mean1 < mean2)
	mean = mean1;
    else
	mean = mean2;

    step = sp.skybw * sp.skyhw;
    if (step < 1)
	step = 1;
    nhstgrm = 512;
    start = mean - (nhstgrm / 2) * step;

 /* Determine the density histogram */

    if (skyhst (&a->d, p, hstgrm, nhstgrm, step, start) < 1)
	return (1);

    for (i = 0, max = 0; i < 512; i++)
	if (hstgrm[i] > max)
	    max = hstgrm[peak = i];
    dmin = peak - 2 * sp.skyhw / step;
    dmax = peak + 2 * sp.skyhw / step;
    if (dmin < 0)
	dmin = 0;
    if (dmax < 0)
	dmax = 0;
    if (dmin >= nhstgrm)
	dmin = nhstgrm - 1;
    if (dmax >= nhstgrm)
	dmax = nhstgrm - 1;
    if (dmax < dmin) {
	peak = dmin;
	dmin = dmax;
	dmax = peak;
    }

    totarea (a, grow);
    ar = &a->t;
    for (i = p->dy - 1, j = 0; i >= 0; i--) {
	cp = p->pic + i * p->dx;
	yi = ar->ychst - p->y0 - i;
	if ((yi >= 0) && (yi < ar->nlines)) {
	    nx = yi == 0 ? ar->nx[yi] : ar->nx[yi] - ar->nx[yi - 1];
	    for (k = p->dx - 1; j < ar->nx[yi]; j++) {
		xl = ar->xla[j] - p->x0;
		xr = ar->xra[j] - p->x0;
		for (; k > xr; k--);
		for (; k >= xl; k--) {
		    cp[k] = start +
			step * ransky (dmin, dmax, max, hstgrm);
		}
		if (nx == 1)
		    for (; k >= 0; k--);
		nx--;
	    }
	}
    }
    return (0.);
}

/*Put zeros outside aperture of radius r centered at (x, y)*/
/*and return luminosity interior*/
float   aperture (x, y, r, ob, p)
int     x, y;
int     r;
struct objrec  *ob;
struct image   *p;
{
    int     a, i, j, xc, yc;
    float   lum, r2, i2;
    PIXEL * cp;

    xc = x - p->x0;
    yc = y - p->y0;
    r2 = r * r;
    lum = 0.;
    for (i = 0; i < p->dy; i++) {
	cp = p->pic + i * p->dx;
	if ((i2 = (i - yc) * (i - yc)) > r2)
	    a = -1;
	else
	    a = sqrt (r2 - i2);
	for (j = 0; j < xc - a; j++)
	    cp[j] = 0;
	for (; (j <= xc + a) && (j < p->dx); j++)
	    lum += film (cp[j]) - ob->sbr;
	for (; j < p->dx; j++)
	    cp[j] = 0;
    }
    return (lum);
}

/*write standard file with image
wtimg(s, o, p)
struct sgparm	*s;
struct objrec	*o;
struct image	*p;
{
	int	i;
	char	file[20];
	struct imgdef	phdr;
	int	fd;
	PIXEL	*cp;

	sprintf(file, "img%-d.%-d.%-d", s->subcat, o->entnum, o->subent);
	phdr = pthdr;
	phdr.maxpxl = -32000;
	phdr.minpxl = 32000;
	phdr.naxis1 = p->dx;
	phdr.naxis2 = p->dy;
	updimhdr( &phdr);
	strncpy(phdr.object, file, 32);
	sprintf(phdr.comment, "subfield at %d %d, res = 1\nclass %s mag %f", 
	p->x0, p->y0, o->class, o->mag);
	fd = imopen(file, &phdr, 1);
	for (i = 0; i < p->dy; i++) {
		cp = p->pic + i * p->dx;
		wtimage(fd, 0, i, p->dx, &phdr, cp);
		pixrange( &phdr, p->dx, cp);
	}
	wtimhdr( &phdr, fd);
	imclose(fd);
}
*/

/*read standard file with image
rdimg1(s, ent, sbent, p)
struct sgparm	*s;
short	ent;
long	sbent;
struct image	*p;
{
	int	i, fd;
	char	file[20];
	struct imgdef	phdr;

	sprintf(file, "img%-d.%-d.%-d", s->subcat, ent, sbent);
	fd = imopen(file, &phdr, 0);
	p->x0 = phdr.crpix1;
	p->y0 = phdr.crpix2;
	p->dx = phdr.naxis1;
	p->dy = phdr.naxis1;
        if (p->pic == 0) {
	    p->mxsz = pthdr.naxis1 * pthdr.naxis2;
	    p->mxsz = MXSZ < p->mxsz ? MXSZ : p->mxsz;
	    for (i = p->mxsz; p->pic == 0 && i > 0; i *= .1) {
	        p->mxsz = i;
	        p->pic = (PIXEL *) fvalloc (p->mxsz, sizeof (PIXEL));
	    }
        }
	if (p->dx * p->dy > p->mxsz)
	    focaserr (1, "rdimg1: Image too large", "");
	rdimage(fd, 0, 0, p->dx * p->dy, &phdr, p->pic);
	imclose(fd);
	return(0);
}
*/

/*Add an image within rectangle rct and return pointer to image structure*/
addimg (fd, q)
int     fd;
struct image   *q;
{
    int     i, j;
    static struct image padd;
    char   *fvalloc ();

    if (padd.pic == 0) {
	padd.pic = (PIXEL *) fvalloc (MXSZ, sizeof (PIXEL));
	if (padd.pic == 0)
	    focaserr (2, "addimg: Not enough memory");
    }
    for (i = 0, j = 0; i < q->dy; i++, j += q->dx)
	rdimage (fd, q->x0, i, q->dx, &pthdr, padd.pic);
    for (i = 0; i < q->dy * q->dx; i++)
	padd.pic[i] += q->pic[i];
    for (i = 0, j = 0; i < q->dy; i++, j += q->dx)
	wtimage (fd, q->x0, i, q->dx, &pthdr, padd.pic);
}

/*Substitute an image rct and return pointer to image structure*/
subimg (fd, p)
int     fd;
struct image   *p;
{
    int     i, j;

    for (i = 0, j = 0; i < p->dy; i++, j += p->dx)
	wtimage (fd, p->x0, i + p->y0, p->dx, &pthdr, p->pic + j);
}

prntp (fd, p)
FILE * fd;
struct image   *p;
{
    int     i, j;
    PIXEL * cp;

    cp = p->pic;
    for (i = 0; i < p->dy; i++) {
	for (j = 0; j < p->dx; j++)
	    fprintf (fd, "%3d ", *cp++);
	fprintf (fd, "\n");
    }
}