path: root/vendor/voclient/voapps/lib/voFITS.c

#include <string.h>
#include <stdio.h>
#include <math.h>
#include <unistd.h>
#include "voApps.h"
#include "fitsio.h"


#define  MAX_IMAGES      20480           /* max images to process        */

#define  dabs(x)         ((x<0.0?-x:x))


/**
 *  Private methods.
 */
static void  vot_printFrameInfo (FILE *fd, frameInfo *im);
static  int  vot_getFrameWcs (fitsfile *fptr, frameInfo *info);

extern  int  vot_fileType (char *name);




/**
 *  VOT_IMAGEINFO -- Get information about a FITS file structure and WCS.
 *  
 *  @fn      info = vot_imageInfo (char *name, int do_all)
 *
 *  @brief            Get information about a FITS file structure and WCS.
 *  @param   name     name of FITS file to open
 *  @param   do_all   do all extensions in an MEF file?
 *  @return  	      pointer to ImInfo structure
 */
ImInfo *
vot_imageInfo (char *name, int do_all)
{
    fitsfile *fptr;
    ImInfo   *info = (ImInfo *) NULL;
    long      naxes[3] = {0, 0, 0}, nrows=0;
    int	      nextns=0, naxis=0, bitpix=0, extnum=0;
    int       hdupos=0, hdutype=0, ncols=0, i=0, status=0;
    double    cxsum=0.0, cysum=0.0, rxsum=0.0, rysum=0.0;


    /*  Check for file existence.
     */
    if (access (name, F_OK) != 0) {
	fprintf (stderr, "Error: cannot open image '%s'\n", name);
	return ((ImInfo *) NULL);
    }
    if (vot_fileType (name) != VOT_FITS) {
	fprintf (stderr, "Error: file '%s' is not a FITS image\n", name);
	return ((ImInfo *) NULL);
    }

    info = (ImInfo *) calloc (1, sizeof (ImInfo));
    if (fits_open_file (&fptr, name, READONLY, &status) == 0) {
	fits_get_num_hdus (fptr, &nextns, &status);
        fits_get_hdu_num (fptr, &hdupos);    /* get the current HDU position */

	strncpy (info->imname, name, strlen(name));
	info->nextend = (nextns - 1);
	info->extns = (frameInfo *) calloc (nextns, sizeof (frameInfo));

        for (; !status; hdupos++) {    /* Main loop for each HDU */
            fits_get_hdu_type (fptr, &hdutype, &status); /* Get the HDU type */

            if (hdutype == IMAGE_HDU) {    /* primary array or image HDU */
                fits_get_img_param (fptr, 10, &bitpix, &naxis, naxes, &status);

		if (hdupos == 0) {				/* PHU	*/
		    info->frame.is_image = 1;
		    info->frame.is_table = 0;
		    info->frame.naxis    = naxis;
		    info->frame.bitpix   = bitpix;
		    for (i=0; i < naxis; i++)
			info->frame.naxes[i] = naxes[i];

		} else {					/* EHU	*/
		    extnum = hdupos - 1;
		    info->extns[extnum].is_image = 1;
		    info->extns[extnum].is_table = 0;
		    info->extns[extnum].extnum   = extnum;
		    info->extns[extnum].naxis    = naxis;
		    info->extns[extnum].bitpix   = bitpix;
		    for (i=0; i < naxis; i++)
			info->extns[extnum].naxes[i] = naxes[i];

		    if (vot_getFrameWcs (fptr, &info->extns[extnum]) == 0)
		        info->extns[extnum].has_wcs = 1;
		}

            } else {					/* a table HDU */
		if (hdupos > 0) {				/* EHU	*/
                    fits_get_num_rows (fptr, &nrows, &status);
                    fits_get_num_cols (fptr, &ncols, &status);

		    extnum = hdupos - 1;
		    info->extns[extnum].is_image = 0;
		    info->extns[extnum].is_table = 1;
		    info->extns[extnum].naxis    = 2;
		    info->extns[extnum].bitpix   = 0;
		    info->extns[extnum].naxes[0] = ncols;
		    info->extns[extnum].naxes[1] = nrows;
		}

#ifdef GET_TBLINFO
                printf ("%s Table:  %d columns x %ld rows\n", 
		    ((hdutype==ASCII_TBL) ? "ASCII" : "Binary"); ncols, nrows);
                printf(" COL NAME             FORMAT\n");
                for (i = 1; i <= ncols; i++) {
                    fits_make_keyn ("TTYPE", i, keyname,&status);
                    fits_read_key (fptr, TSTRING, keyname,colname,NULL,&status);
                    fits_make_keyn ("TFORM", i, keyname,&status);
                    fits_read_key (fptr, TSTRING, keyname,coltype,NULL,&status);
                    printf(" %3d %-16s %-16s\n", i, colname, coltype);
                }
#endif
            }

	    /*  Move to next extension.
	     */
            fits_movrel_hdu (fptr, 1, NULL, &status); 
        }

        if (status == END_OF_FILE)
            status = 0;    			/* reset normal error 	   */

    } else if (status)    			/* print any error message */
        fits_report_error (stderr, status);


    /*  Compute the values for the entire frame.
     */
    info->frame.lx = info->frame.ly =  360.0;
    info->frame.ux = info->frame.uy = -360.0;

    info->frame.cx     = info->extns[0].cx;
    info->frame.cy     = info->extns[0].cy;
    info->frame.lx     = info->extns[0].lx;
    info->frame.ly     = info->extns[0].ly;
    info->frame.ux     = info->extns[0].ux;
    info->frame.uy     = info->extns[0].uy;
    info->frame.rotang = info->extns[0].rotang;
    info->frame.xrval  = info->extns[0].xrval;
    info->frame.yrval  = info->extns[0].yrval;
    info->frame.xrpix  = info->extns[0].xrpix;
    info->frame.yrpix  = info->extns[0].yrpix;
    info->frame.radius = info->extns[0].radius;
		    
    if (nextns == 1) {
        for (i=0; i < naxis; i++)
	    info->frame.naxes[i] = info->extns[0].naxes[i];
    }
    memcpy (&info->frame.xc[0], &info->extns[0].xc[0], (sizeof(double)*4));
    memcpy (&info->frame.yc[0], &info->extns[0].yc[0], (sizeof(double)*4));

    for (i=1; i < nextns; i++) {
	if (info->extns[i].lx < info->frame.lx) 
	    info->frame.lx = info->extns[i].lx;
	if (info->extns[i].ly < info->frame.ly) 
	    info->frame.ly = info->extns[i].ly;

	if (info->extns[i].ux > info->frame.ux) 
	    info->frame.ux = info->extns[i].ux;
	if (info->extns[i].uy > info->frame.uy) 
	    info->frame.uy = info->extns[i].uy;
	cxsum += info->extns[i].cx;
	cysum += info->extns[i].cy;

	rxsum += info->extns[i].naxes[0];
	rysum += info->extns[i].naxes[1];

	info->frame.scale  = info->extns[i].scale;
	info->frame.rotang = info->extns[i].rotang;
	strcpy (info->frame.ctype, info->extns[i].ctype);
    }

    if (nextns > 1) {
        info->frame.xrval = info->frame.cx = (cxsum / (double) nextns);
        info->frame.yrval = info->frame.cy = (cysum / (double) nextns);
        info->frame.xrpix = info->frame.radius / 
		(info->frame.scale / 3600.) / 2.0;
        info->frame.yrpix = info->frame.radius / 
		(info->frame.scale / 3600.) / 2.0;
    } else {
        info->frame.xrval = info->frame.cx = info->extns[0].cx;
        info->frame.yrval = info->frame.cy = info->extns[0].cy;
        info->frame.xrpix = info->extns[0].xrpix;
        info->frame.yrpix = info->extns[0].yrpix;
	strcpy (info->frame.ctype, info->extns[0].ctype);
    }

    info->frame.width  = ((info->frame.ux+360.) - (info->frame.lx+360.));
    info->frame.height = ((info->frame.uy+ 90.) - (info->frame.ly+ 90.));
    info->frame.radius = sqrt (
	(info->frame.cx - info->frame.lx) * 
	(info->frame.cx - info->frame.lx) + 
	(info->frame.cy - info->frame.ly) * 
	(info->frame.cy - info->frame.ly) );
    info->frame.xc[0]  = info->frame.lx;   info->frame.yc[0] = info->frame.ly;
    info->frame.xc[1]  = info->frame.lx;   info->frame.yc[1] = info->frame.uy;
    info->frame.xc[2]  = info->frame.ux;   info->frame.yc[2] = info->frame.uy;
    info->frame.xc[3]  = info->frame.ux;   info->frame.yc[3] = info->frame.ly;


    fits_close_file (fptr, &status);
    return ( (ImInfo *) info);
}


/**
 *  VOT_IMAGENEXTNS -- Get the number of extensions in an MEF file.
 *
 *  @fn      nextn = vot_imageNExtns (char *name)
 *
 *  @brief            Get the number of extensions in an MEF file.
 *  @param   name     name of FITS file to open
 *  @return           number of extensions in an MEF file
 */
int
vot_imageNExtns (char *image)
{
    fitsfile *fptr;
    int   nextns = -1, status = 0;


    if (fits_open_file (&fptr, image, READONLY, &status) == 0) {
	fits_get_num_hdus (fptr, &nextns, &status);
        fits_close_file (fptr, &status);
    }

    return (nextns);
}


/**
 *  VOT_PRINTIMAGEINFO -- Print the image information.
 *
 *  @fn      vot_printImageInfo (FILE *fd, ImInfo *im)
 *
 *  @brief            Print the image information.
 *  @param   fd       output file descriptor (or stdout/stderr)
 *  @param   im       image information structure
 *  @return           nothing
 */
void
vot_printImageInfo (FILE *fd, ImInfo *im)
{
    register int i;

    fprintf (fd, "Name: %s    nextns = %d\n", im->imname, im->nextend);

    fprintf (fd, "Frame:\n");   im->frame.extnum = -1;
        vot_printFrameInfo (fd, &im->frame);

    for (i=1; i < im->nextend; i++)
	vot_printFrameInfo (fd, &im->extns[i]);
}


/**
 *  VOC_FREEIMAGEINFO -- Free the image information structure.
 *
 *  @fn      vot_freeImageInfo (ImInfo *im)
 *
 *  @brief            Free the image information structure.
 *  @param   im       image information structure
 *  @return           nothing
 */
void
vot_freeImageInfo (ImInfo *im)
{
    free ((void *) im->extns);            /* extension structs    */
    free ((void *) im);                   /* main image structs   */
}



/****************************************************************************
 ***   				Private procedures
 ****************************************************************************/

/**
 *  VOT_PRINTFRAMEINFO -- Print information about a specific frame.
 */
static void
vot_printFrameInfo (FILE *fd, frameInfo *im)
{
    if (im->extnum >= 0) {
        fprintf (fd, "\nExt: %d  dims[%d] = %d %d %d\t",
            im->extnum, im->naxis, im->naxes[0], im->naxes[1], im->naxes[2]);
        fprintf (fd, "  image: %d   table: %d   has_wcs: %d   flip: %d\n",
            im->is_image, im->is_table, im->has_wcs, im->axflip);
    }

    fprintf (fd, 
	"  center:  %8.4f %8.4f   ll: %8.4f %8.4f    ur: %8.4f %8.4f\n",
        im->cx, im->cy, im->lx, im->ly, im->ux, im->uy);
    fprintf (fd, "  corners: %8.4f %8.4f       %8.4f %8.4f\n",
	im->xc[1], im->yc[1], im->xc[2], im->yc[2]);
    fprintf (fd, "\t   %8.4f %8.4f       %8.4f %8.4f\n",
        im->xc[0], im->yc[0], im->xc[3], im->yc[3]);
    fprintf (fd, 
	"  crval: %8.4f %8.4f   crpix: %8.4f %8.4f  ctype: '%s'\n",
        im->xrval, im->yrval, im->xrpix, im->yrpix, im->ctype);
    fprintf (fd, 
	"  w/h: %8.4f %8.4f   radius: %8.4f  rot: %8.4f  scale: %8.4f\n",
        dabs(im->width), dabs(im->height), im->radius, im->rotang, im->scale); 
}


/**
 *  VOT_GETFRAMEWCS -- Get the WCS information for a given frame.
 */
static int  
vot_getFrameWcs (fitsfile *fptr, frameInfo *info)
{
    double   xrval=0.0, yrval=0.0, xrpix=0.0, yrpix=0.0, xpix=0.0, ypix=0.0;
    double   xinc=0.0, yinc=0.0, rot=0.0, scale=0.0, xrot=0.0, yrot=0.0;
    double   cx=0.0, cy=0.0, lx=0.0, ly=0.0, ux=0.0, uy=0.0;
    double   cd11=0.0, cd12=0.0, cd21=0.0, cd22=0.0, cdelt1=0.0, cdelt2=0.0;
    int      i, axflip=0, status = 0;
    char     str[32], ctype[5], comment[80];


    /*  Get the header WCS keywords.
     */
    fits_read_img_coord (fptr, &xrval, &yrval, &xrpix,
               &yrpix, &xinc, &yinc, &rot, ctype, &status);

    info->xrval    = xrval;
    info->yrval    = yrval;
    info->xrpix    = xrpix;
    info->yrpix    = yrpix;
    info->has_wcs  = 1;
    info->is_image = 1;
    strcpy (info->ctype, ctype);

    xpix = (double) 0.5;                        /*   Lower-left         */
    ypix = (double) 0.5;
    status = 0;
    fits_pix_to_world (xpix, ypix, xrval, yrval, xrpix, yrpix, xinc, yinc, 
	rot, ctype, &info->xc[0], &info->yc[0], &status);

    xpix = (double) 0.5;                        /*   Upper-left         */
    ypix = (double) info->naxes[1] - 0.5;
    status = 0;
    fits_pix_to_world (xpix, ypix, xrval, yrval, xrpix, yrpix, xinc, yinc, 
	rot, ctype, &info->xc[1], &info->yc[1], &status);

    xpix = (double) info->naxes[0] - 0.5;       /*   Upper-right        */
    ypix = (double) info->naxes[1] - 0.5;
    status = 0;
    fits_pix_to_world (xpix, ypix, xrval, yrval, xrpix, yrpix, xinc, yinc, 
	rot, ctype, &info->xc[2], &info->yc[2], &status);

    xpix = (double) info->naxes[0] - 0.5;       /*   Lower-right        */
    ypix = (double) 0.5;
    status = 0;
    fits_pix_to_world (xpix, ypix, xrval, yrval, xrpix, yrpix, xinc, yinc, 
	rot, ctype, &info->xc[3], &info->yc[3], &status);

    xpix = (double) info->naxes[0] / 2. - 0.5;	/*  Center		*/
    ypix = (double) info->naxes[1] / 2. - 0.5;
    status = 0;
    fits_pix_to_world (xpix, ypix, xrval, yrval, xrpix, yrpix, xinc, yinc, 
	rot, ctype, &info->cx, &info->cy, &status);
    info->cx = xrval;
    info->cy = yrval;

    /*  Get center and cone radius.
     */
    lx = info->lx = info->xc[0];		/*  Lower-left		*/
    ly = info->ly = info->yc[0];
    ux = info->ux = info->xc[2];		/*  Upper-right		*/
    uy = info->uy = info->yc[2];
    cx = info->cx;
    cy = info->cy;
    status = 0;
    if ((i = fits_read_key_dbl(fptr, "CD1_1", &cd11, comment, &status))==0) {
        fits_read_key_dbl (fptr, "CD1_2", &cd12, comment, &status);
        fits_read_key_dbl (fptr, "CD2_1", &cd21, comment, &status);
        fits_read_key_dbl (fptr, "CD2_2", &cd22, comment, &status);

        scale = 3600.0 * sqrt ((cd11*cd11+cd21*cd21+cd12*cd12+cd22*cd22) / 2.);
	xrot  = dabs (atan2 ( cd21, cd11));
	yrot  = dabs (atan2 (-cd12, cd22));
	rot   = (xrot + yrot) / 2.0;
    } else {
	/*  Old-style keywords.
	 */
	status = 0;
        if (!fits_read_key_dbl (fptr, "CDELT1", &cdelt1, comment, &status)) {
             fits_read_key_dbl (fptr, "CDELT2", &cdelt2, comment, &status);

	    scale = 3600.0 * sqrt ((cdelt1*cdelt1 + cdelt2*cdelt2) / 2.);

            if (!fits_read_key_dbl (fptr, "CROTA1", &xrot, comment, &status)) {
                 fits_read_key_dbl (fptr, "CROTA2", &yrot, comment, &status);
		     rot   = (xrot + yrot) / 2.0;
            }
        } else
	    info->has_wcs  = 0;
    }

    status = 0;
    memset (str, 0, 32);
    if ((i = fits_read_key_str(fptr, "CTYPE1", str, comment, &status))==0) {
	if (strncasecmp (str,"DEC",3) == 0 || strncasecmp (str,"LAT",3) == 0)
	    axflip = 1;
    }

    /*  For a bad/approximate WCS, compute in rough coords.
     */
    if ( (lx == cx) && (ly == cy) ) {
	double s = (scale / 3600.0);

        lx = info->lx = cx - (info->naxes[0]/2 * s); /* Lower-left	*/
        ly = info->ly = cy - (info->naxes[1]/2 * s);
        ux = info->ux = cx + (info->naxes[0]/2 * s); /* Upper-right	*/
        uy = info->uy = cy + (info->naxes[1]/2 * s);

        info->xc[0] = cx - (info->naxes[0]/2 * s); /* Lower-left	*/
        info->yc[0] = cy - (info->naxes[1]/2 * s);

        info->xc[1] = cx - (info->naxes[0]/2 * s); /* Upper-left	*/
        info->yc[1] = cy + (info->naxes[1]/2 * s);

        info->xc[2] = cx + (info->naxes[0]/2 * s); /* Upper-right	*/
        info->yc[2] = cy + (info->naxes[1]/2 * s);

        info->xc[3] = cx + (info->naxes[0]/2 * s); /* Lower-right	*/
        info->yc[3] = cy - (info->naxes[1]/2 * s);
    }

    /*  FIXME -- Doens't handle rotation properly. */
    info->width  = dabs((ux - lx));		/* in degrees		*/
    info->height = dabs((uy - ly));
    info->radius = sqrt ((cx - lx) * (cx - lx) + (cy - ly) * (cy - ly));
    info->rotang = rot;
    info->scale  = scale;
    info->axflip = axflip;

    return (info->has_wcs);
}


/*****************************************************************************
 *  Program main.
 ****************************************************************************/
#ifdef UNIT_TEST
int 
main (int argc, char *argv[])
{
    ImInfo *img = (ImInfo *) NULL;


    if ( (img = vot_imageInfo (argv[1], 1)) ) {
	vot_printImageInfo (stdout, img);
	vot_freeImageInfo (img);
    }

    return (0);
}
#endif