aboutsummaryrefslogtreecommitdiff
path: root/vendor/voclient/voapps/lib/voFITS.c
diff options
context:
space:
mode:
Diffstat (limited to 'vendor/voclient/voapps/lib/voFITS.c')
-rw-r--r--vendor/voclient/voapps/lib/voFITS.c467
1 files changed, 467 insertions, 0 deletions
diff --git a/vendor/voclient/voapps/lib/voFITS.c b/vendor/voclient/voapps/lib/voFITS.c
new file mode 100644
index 00000000..b6d9c48e
--- /dev/null
+++ b/vendor/voclient/voapps/lib/voFITS.c
@@ -0,0 +1,467 @@
+#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
generated by cgit (git 2.34.1) at 2025-09-20 16:39:04 -0400