/****************************************/
/*	imageio                         */
/*					*/
/*	Image io subroutines		*/
/****************************************/
#include	"image.h"


char   *funits[] = {		/* List of possible units	 */
    "",				/* Unitless			 */
    "U",			/* Intensity in U bandpass	 */
    "B",			/* Intensity in B bandpass	 */
    "J",			/* Intensity in J bandpass	 */
    "V",			/* Intensity in V bandpass	 */
    "F",			/* Intensity in F bandpass	 */
    "PDS",			/* PDS density units		 */
    "RA",			/* Right Ascension in degrees	 */
    "DEC",			/* Declination in degrees	 */
    "ARC SEC",			/* Angular arc second units	 */
    0
};


/* Open image file and return fd */
imopen (file, imhdr, mode)
char   *file;
struct imgdef  *imhdr;
int     mode;
/* If mode = 0  Read image header */
/* If mode = 1  Write image header */
/* If mode = 2  If file exists read header otherwise write header */
{
    int     i, fdim, acmode, hdrlen=0;

    switch (mode) {
	case 0: 
	    acmode = READ_ONLY;
	    iferr (fdim = immap (c_sppstr (file), &acmode, &hdrlen)) {
		focaserr (1, "imopen: Can't read ", file);
		return (NULL);
	    }
	    rdimhdr (fdim, imhdr);
	    break;
	case 1: 
	    acmode = NEW_FILE;
	    iferr (fdim = immap (c_sppstr (file), &acmode, &hdrlen)) {
		iferr (imdelete (c_sppstr (file)))
		    ;
		iferr (fdim = immap (c_sppstr (file), &acmode, &hdrlen)) {
		    focaserr (2, "imopen: Can't write ", file);
		    return (NULL);
		}
	    }
	    i = 2; imaddi (&fdim, "naxis", &i);
    	    i = imhdr->naxis1; imaddi (&fdim, "naxis1", &i);
    	    i = imhdr->naxis2; imaddi (&fdim, "naxis2", &i);
	    if (imhdr->bitpix == 8 || imhdr->bitpix == 16) {
		imhdr->bytepix = 2;
		imhdr->bitpix = 16;
		i = TY_SHORT;
	    } else {
		imhdr->bytepix = 4;
		imhdr->bitpix = 32;
		i = TY_REAL;
	    }
	    imaddi (&fdim, "pixtype",  &i);

	    wtimhdr (fdim, imhdr);
	    break;
	case 2: 
	    acmode = READ_WRITE;
	    iferr (fdim = immap (c_sppstr (file), &acmode, &hdrlen)) {
		acmode = NEW_FILE;
		iferr (fdim = immap (c_sppstr (file), &acmode, &hdrlen)) {
		    focaserr (3, "imopen: Can't write image ", file);
		    return (NULL);
		}
		i = 2; imaddi (&fdim, "naxis", &i);
		i = imhdr->naxis1; imaddi (&fdim, "naxis1",  &i);
		i = imhdr->naxis2; imaddi (&fdim, "naxis2",  &i);
		if (imhdr->bitpix == 8 || imhdr->bitpix == 16) {
		    imhdr->bytepix = 2;
		    imhdr->bitpix = 16;
		    i = TY_SHORT;
		} else {
		    imhdr->bytepix = 4;
		    imhdr->bitpix = 32;
		    i = TY_REAL;
		}
		imaddi (&fdim, "pixtype",  &i);

		wtimhdr (fdim, imhdr);
	    }
	    rdimhdr (fdim, imhdr);
	    break;
    }
    return (fdim);
}

/* Close image file */
imclose (fdim)
int    fdim;
{
    imunmap (&fdim);
}


/*Read image header and return read status*/
rdimhdr (fdim, imhdr)
int     fdim;
struct imgdef  *imhdr;
{
    int     i;
    float   r;
    char    sval[80];

    imhdr->simple[0] = 'T';
    imhdr->blksiz = 2880;
    imgeti (&fdim, "naxis", &i); imhdr->naxis = i;
    imgeti (&fdim, "naxis1", &i); imhdr->naxis1 = i;
    imgeti (&fdim, "naxis2", &i); imhdr->naxis2 = i;
    imgeti (&fdim, "pixtype", &i);
	if (i == 3) {
	    imhdr->bitpix = 16;
	    imhdr->bytepix = 2;
	} else {
	    imhdr->bitpix = 32;
	    imhdr->bytepix = 4;
	}
    /*
    iferr (imgeti (&fdim, "dflags", &i))
	imhdr->dflags = 0;
    else
	imhdr->dflags = i;
    iferr (imgeti (&fdim, "maptyp", &i))
	imhdr->maptyp = 0;
    else
	imhdr->maptyp = i;
    */
    iferr (imgetr (&fdim, "maxpixval", &r))
	imhdr->maxpxl = MAXPXL;
    else
	imhdr->maxpxl = r;
    iferr (imgetr (&fdim, "minpixval", &r))
	imhdr->minpxl = -MAXPXL;
    else
	imhdr->minpxl = r;
    if (imhdr->maxpxl == 0 && imhdr->minpxl == 0) {
	imhdr->maxpxl = MAXPXL;
	imhdr->minpxl = -MAXPXL;
    }
    iferr (imgetr (&fdim, "bscale", &r))
	imhdr->bscale = 1.;
    else
	imhdr->bscale = r;
    iferr (imgetr (&fdim, "bzero", &r))
	imhdr->bzero =  0.;
    else
	imhdr->bzero =  r;
    iferr (imgetr (&fdim, "ltm1_1", &r))
	r = 1.;
    if (r == 0.)
	focaserr (1,
	    "rdcathdr: ltm1_1 is zero.  Possibly need to reset WCS", "");
    imhdr->cdelt1 = 1 / r;
    iferr (imgetr (&fdim, "ltv1", &r))
	r = 0.;
    imhdr->crpix1 = imhdr->cdelt1 * (1 - r);
    /*
    iferr (imgetr (&fdim, "crval1", &r))
	imhdr->crval1 = 0.;
    else
	imhdr->crval1 = r;
    */
    iferr (imgetr (&fdim, "ltm2_2", &r))
	r = 1.;
    if (r == 0.)
	focaserr (1,
	    "rdcathdr: ltm2_2 is zero.  Possibly need to reset WCS", "");
    imhdr->cdelt2 = 1 / r;
    iferr (imgetr (&fdim, "ltv2", &r))
	r = 0.;
    imhdr->crpix2 = imhdr->cdelt2 * (1 - r);
    /*
    iferr (imgetr (&fdim, "crval2", &r))
	imhdr->crval2 = 0.;
    else
	imhdr->crval2 = r;
    */
    iferr (imgeti (&fdim, "blank",  &i))
	imhdr->blank = -(MAXPXL-1);
    else
	imhdr->blank = i;
    iferr (imgetr (&fdim, "dscale", &r))
	imhdr->dscale = 1.;
    else
	imhdr->dscale = r;
    iferr (imgetr (&fdim, "dzero",  &r))
	imhdr->dzero = 0.;
    else
	imhdr->dzero = r;
    iferr (imgstr (&fdim, "title", sval))
	imhdr->object[0] = 0;
    else
	strncpy (imhdr->object, sval, 32);
    iferr (imgstr (&fdim, "origin", sval))
	imhdr->origin[0] = 0;
    else
	strncpy (imhdr->origin, sval, 32);
    iferr (imgstr (&fdim, "date-obs", sval))
	imhdr->date[0] = 0;
    else
	strncpy (imhdr->date, sval, 32);
    iferr (imgstr (&fdim, "observer", sval))
	imhdr->observer[0] = 0;
    else
	strncpy (imhdr->observer, sval, 32);
    imhdr->comment[0] = 0;
}

/*Write image header and return write status*/
wtimhdr (fdim, imhdr)
int     fdim;
struct imgdef  *imhdr;
{
    int    i;
    float  r;

    r = imhdr->maxpxl; imaddr (&fdim, "maxpixval", &r);
    r = imhdr->minpxl; imaddr (&fdim, "minpixval", &r);
    r = 1-imhdr->crpix1/imhdr->cdelt1; imaddr (&fdim, "ltv1",  &r);
    r = 1/imhdr->cdelt1; imaddr (&fdim, "ltm1_1",  &r);
    r = 1-imhdr->crpix2/imhdr->cdelt2; imaddr (&fdim, "ltv2",  &r);
    r = 1/imhdr->cdelt2; imaddr (&fdim, "ltm2_2",  &r);
    r = imhdr->dscale; imaddr (&fdim, "dscale",  &r);
    r = imhdr->dzero;  imaddr (&fdim, "dzero",   &r);
    /*
    i = imhdr->dflags; imaddi (&fdim, "dflags",  &i);
    i = imhdr->maptyp; imaddi (&fdim, "maptyp",  &i);
    r = imhdr->bscale; imaddr (&fdim, "bscale",  &r);
    r = imhdr->bzero;  imaddr (&fdim, "bzero",   &r);
    r = imhdr->crval1; imaddr (&fdim, "crval1",  &r);
    r = imhdr->crval2; imaddr (&fdim, "crval2",  &r);
    i = imhdr->blank;  imaddi (&fdim, "blank",   &i);
    imastr (&fdim, "title", imhdr->object);
    imastr (&fdim, "origin", imhdr->origin);
    imastr (&fdim, "date-obs", imhdr->date);
    imastr (&fdim, "observer", imhdr->observer);
    */
}

/*Read image file*/
rdimage (fdim, x, y, pixels, imhdr, buf)
int     fdim;
int     x, y, pixels;
struct imgdef  *imhdr;
PIXEL * buf;
{
    XPOINTER rbuf;
    XREAL *ptr;
    int x1, x2, y1, y2;

    x1 = x + 1;
    x2 = x1 + pixels - 1;
    y1 = y + 1;
    y2 = y1;
    iferr (rbuf = imgs2r (&fdim, &x1, &x2, &y1, &y2)) {
	focaserr (4, "rdimage: Error reading image", "");
	return (1);
    }
    ptr = &Memr[rbuf];
    for (x1 = 0; x1 < pixels; x1++)
	buf[x1] = *ptr++;
    return (0);
}

/*Write image file*/
wtimage (fdim, x, y, pixels, imhdr, buf)
int     fdim;
int     x, y, pixels;
struct imgdef  *imhdr;
PIXEL  *buf;
{
    XPOINTER rbuf;
    float *ptr;
    int x1, x2, y1, y2;

    x1 = x + 1;
    x2 = x1 + pixels - 1;
    y1 = y + 1;
    y2 = y1;
    iferr (rbuf = imps2r (&fdim, &x1, &x2, &y1, &y2)) {
	focaserr (5, "wtimage: Error writing image", "");
	return (1);
    }
    ptr = &Memr[rbuf];
    for (x1 = 0; x1 < pixels; x1++)
	*ptr++ = buf[x1];
    return (0);
}


/*Set initial header*/
setimhdr (fd, imhdr)
int     fd;
struct imgdef  *imhdr;
{
    if (imhdr->origin[0] == 0)
	sprintf (imhdr->origin, "KPNO");
    if (imhdr->date[0] == 0)
	updimhdr (imhdr);
    if (imhdr->cdelt1 == 0)
	imhdr->cdelt1 = 1.;
    if (imhdr->cdelt2 == 0)
	imhdr->cdelt2 = 1.;
    if ((imhdr->maxpxl == 0) && (imhdr->minpxl == 0)) {
	imhdr->maxpxl = -MAXPXL;
	imhdr->minpxl = MAXPXL;
    }
    if (imhdr->bscale == 0.) {
	imhdr->bscale = 1.;
	imhdr->bzero = 0.;
    }
    if (imhdr->dscale == 0.) {
	imhdr->dscale = imhdr->bscale;
	imhdr->dzero = imhdr->bzero;
    }
    wtimhdr (fd, imhdr);
}

/*Update header*/
updimhdr (imhdr)
struct imgdef  *imhdr;
{
    c_cnvtime (c_clktime(0), imhdr->date, 31);
}


dmpimhdr (imhdr)
struct imgdef  *imhdr;
{

    /*
    printf ("Object:   %s\n", imhdr->object);
    printf ("Observer: %s\n", imhdr->observer);
    printf ("Date:     %s\n", imhdr->date);
    printf ("Origin:   %s\n", imhdr->origin);
    printf ("\n");
    */
    printf ("Bitpix:   %d\n", imhdr->bitpix);
    printf ("Bytepix:  %d\n", imhdr->bytepix);
    printf ("Naxis1:   %d\n", imhdr->naxis1);
    printf ("Naxis2:   %d\n", imhdr->naxis2);
    printf ("Maxpxl:   %f\n", imhdr->maxpxl);
    printf ("Minpxl:   %f\n", imhdr->minpxl);
    /*
    printf ("Bscale:   %f\n", imhdr->bscale);
    printf ("Bzero:    %f\n", imhdr->bzero);
    */
    printf ("Dscale:   %f\n", imhdr->dscale);
    printf ("Dzero:    %f\n", imhdr->dzero);
    /* printf ("Crval1:   %f\n", imhdr->crval1); */
    printf ("Crpix1:   %f\n", imhdr->crpix1);
    printf ("Cdelt1:   %f\n", imhdr->cdelt1);
    /* printf ("Crval2:   %f\n", imhdr->crval2); */
    printf ("Crpix2:   %f\n", imhdr->crpix2);
    printf ("Cdelt2:   %f\n", imhdr->cdelt2);
    /*
    printf ("Blank:    %d\n", imhdr->blank);
    printf ("Dflags:   %d\n", imhdr->dflags);
    printf ("Maptyp:   %d\n", imhdr->maptyp);
    */
    prcmmnt (imhdr);
}

/*Set range of pixel values in line*/
pixrange (imhdr, n, buf)
struct imgdef  *imhdr;
int     n;
PIXEL * buf;
{
    int     i;
    long    j;

    for (i = 0; i < n; i++) {
	j = buf[i];
	if (j == imhdr->blank)
	    continue;
	imhdr->maxpxl = j > imhdr->maxpxl ? j : imhdr->maxpxl;
	imhdr->minpxl = j < imhdr->minpxl ? j : imhdr->minpxl;
    }
}

/*Modify an image header*/
modimhdr (imhdr)
struct imgdef  *imhdr;
{
    char    str[81];

    /*
    printf ("Object   ( %32s ) = ", imhdr->object);
    if (strlen (gets (str)))
	strncpy (imhdr->object, str, 32);
    printf ("Observer ( %32s ) = ", imhdr->observer);
    if (strlen (gets (str)))
	strncpy (imhdr->observer, str, 32);
    printf ("Date     ( %32s ) = ", imhdr->date);
    if (strlen (gets (str)))
	strncpy (imhdr->date, str, 32);
    printf ("Origin   ( %32s ) = ", imhdr->origin);
    if (strlen (gets (str)))
	strncpy (imhdr->origin, str, 32);
    */
    printf ("Bitpix   ( %32d )\n", imhdr->bitpix);
    printf ("Bytepix  ( %32d )\n", imhdr->bytepix);
    printf ("Naxis1   ( %32d )\n", imhdr->naxis1);
    printf ("Naxis2   ( %32d )\n", imhdr->naxis2);
    printf ("Minpxl   ( %32f )\n", imhdr->minpxl);
    printf ("Maxpxl   ( %32f )\n", imhdr->maxpxl);
    /*
    printf ("Bscale   ( %32f ) = ", imhdr->bscale);
    if (strlen (gets (str)))
	sscanf (str, "%hf", &imhdr->bscale);
    printf ("Bzero    ( %32f ) = ", imhdr->bzero);
    if (strlen (gets (str)))
	sscanf (str, "%hf", &imhdr->bzero);
    */
    printf ("Dscale   ( %32f ) = ", imhdr->dscale);
    if (strlen (gets (str)))
	sscanf (str, "%hf", &imhdr->dscale);
    printf ("Dzero    ( %32f ) = ", imhdr->dzero);
    if (strlen (gets (str)))
	sscanf (str, "%hf", &imhdr->dzero);
    /*
    printf ("Crval1   ( %32f ) = ", imhdr->crval1);
    if (strlen (gets (str)))
	sscanf (str, "%hf", &imhdr->crval1);
    */
    printf ("Crpix1   ( %32f ) = ", imhdr->crpix1);
    if (strlen (gets (str)))
	sscanf (str, "%hf", &imhdr->crpix1);
    printf ("Cdelt1   ( %32f ) = ", imhdr->cdelt1);
    if (strlen (gets (str)))
	sscanf (str, "%hf", &imhdr->cdelt1);
    /*
    printf ("Crval2   ( %32f ) = ", imhdr->crval2);
    if (strlen (gets (str)))
	sscanf (str, "%hf", &imhdr->crval2);
    */
    printf ("Crpix2   ( %32f ) = ", imhdr->crpix2);
    if (strlen (gets (str)))
	sscanf (str, "%hf", &imhdr->crpix2);
    printf ("Cdelt2   ( %32f ) = ", imhdr->cdelt2);
    if (strlen (gets (str)))
	sscanf (str, "%hf", &imhdr->cdelt2);
    /*
    printf ("Blank    ( %32d ) = ", imhdr->blank);
    if (strlen (gets (str)))
	sscanf (str, "%d", &imhdr->blank);
    if (strlen (gets (str)))
	sscanf (str, "%hd", &imhdr->dflags);
    printf ("Maptyp   ( %32d ) = ", imhdr->maptyp);
    if (strlen (gets (str)))
	sscanf (str, "%hd", &imhdr->maptyp);
    */
}

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

 /* Check desired rectangle */
    if ((rct[0] < 0) || (rct[1] < 0) || (rct[2] < 0) || (rct[3] < 0))
	return (-1);
    if ((rct[0] >= imhdr->naxis1) || (rct[1] >= imhdr->naxis2) ||
	    (rct[2] >= imhdr->naxis1) || (rct[3] >= imhdr->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->dx * p->dy > MXSZ) {
	fprintf (stderr, "rdrct: Image too large\n");
	return (-1);
    }
    if (p->pic == 0) {
	p->pic = (PIXEL *) calloc (MXSZ, sizeof (PIXEL));
	if (p->pic == 0) {
	    fprintf (stderr, "rdrct: Not enough memory\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, imhdr, 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);
}

/* Add comment */
addcmmnt (imhdr, string)
struct imgdef  *imhdr;
char   *string;
{
    int     i, j;

    i = strlen (imhdr->comment);
    j = strlen (string) + 1;
    if (i + j > 255)
	return;
    if (imhdr->comment[0] != 0)
	strcat (imhdr->comment, "\n");
    strcat (imhdr->comment, string);
}

/* Print comments */
prcmmnt (imhdr)
struct imgdef  *imhdr;
{
    printf ("Comments:\n%s\n", imhdr->comment);
}