index d1c22f0..819158b 100644 (file)
 #include "fits2jpeg.h"
 
 /*---------------------------------------------------------------------------*
+ * READ_FITS: To reads data from a fits image file.. (isn't that obvious?)
+ *---------------------------------------------------------------------------*/
+void read_fits(char * fits_file_name, long * xdim, long * ydim, float ** data)
+{
+    fitsfile *fptr;
+    int status = 0, nfound, anynull;
+    long naxes[2];
+    long npixels;
+    float nullval = 0.0;
+
+    fits_open_file(&fptr, fits_file_name, READONLY, &status);
+    fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxes, &nfound, &status);
+    if (status)
+        printerro(strcat(fits_file_name, " <-- Failed to open the file"));
+
+    /* Read in data */
+    npixels = naxes[0] * naxes[1];
+    (*data) = malloc(sizeof(float) * npixels);
+
+    nullval = 0;
+    if (fits_read_img(fptr, TFLOAT, 1, npixels, &nullval, (*data), &anynull,
+        &status))
+        printerro(strcat(fits_file_name, " has no valid fits image data"));
+
+    *xdim = naxes[0];
+    *ydim = naxes[1];
+
+    fits_close_file(fptr, &status);
+}
+
+/*---------------------------------------------------------------------------*
  * SCALE_PIXELS: Changes the pixel scale to linear/log/sqroot/etc..
  *---------------------------------------------------------------------------*/
-void scale_pixels(int scale, unsigned int npixels,
-                 float *data, JSAMPLE *image_buffer)
+void scale_pixels(int scale, unsigned int npixels, float *data,
+                 JSAMPLE ** image_buffer)
 {
     unsigned int i = 0;
     int JMAXVAL = 255;
@@ -51,6 +82,9 @@ void scale_pixels(int scale, unsigned int npixels,
     /* the dynamic range is reduced to 255 for jpeg                         */
     scl_data = (datamax - datamin)/(float)JMAXVAL;
 
+    /* we will end up with segfaults if scl_data = 0                        */
+    if (scl_data == 0) scl_data = 1;
+
     for (i = 0; i < npixels; ++i)
         data[i] = (data[i] - datamin)/scl_data;
 
@@ -74,10 +108,14 @@ void scale_pixels(int scale, unsigned int npixels,
     for (i = 1; i <= JMAXVAL; ++i)
         cumhist[i] += cumhist[i - 1] + hist[i];
 
+    /* Allocate image buffer */
+    (*image_buffer) = malloc(sizeof(unsigned char) * npixels);
+
+
     /* Linear scale (min-max) : This is the default scaling
      * histo-eq will fail if we dont generate image_buffer here              */
     for (i = 0; i < npixels; ++i)
-        image_buffer[i] = (int)(data[i]);
+        (*image_buffer)[i] = (int)(data[i]);
 
     /*-----------------------------------------------------------------------*/
 
@@ -85,39 +123,32 @@ void scale_pixels(int scale, unsigned int npixels,
     switch (scale)
     {
         case 1 :                                              /* Square root */
-            printinfo("Using square-root scale");
             scl_data = sqrt((float)JMAXVAL)/(float)JMAXVAL;
             for (i = 0; i < npixels; ++i)
-                image_buffer[i] = (int)(sqrt(data[i])/scl_data);
+                (*image_buffer)[i] = (int)(sqrt(data[i])/scl_data);
             break;
 
         case 2 :                                                   /* Square */
-            printinfo("Using quadratic scale");
             scl_data = pow((float)JMAXVAL,2)/(float)JMAXVAL;
             for (i = 0; i < npixels; ++i)
-                image_buffer[i] = (int)abs((pow(data[i],2) - 1.0)/scl_data);
+                (*image_buffer)[i] = (int)abs((pow(data[i],2) - 1.0)/scl_data);
             break;
 
         case 3 :                                                    /* Cubic */
-            printinfo("Using cubic scale");
             scl_data = pow((float)JMAXVAL,3)/(float)JMAXVAL;
             for (i = 0; i < npixels; ++i)
-                image_buffer[i] = (int)abs((pow(data[i],3) - 1.0)/scl_data);
+                (*image_buffer)[i] = (int)abs((pow(data[i],3) - 1.0)/scl_data);
             break;
 
         case 4 :                                                      /* log */
-            printinfo("Using log scale");
             scl_data = log(1.0 + (float)JMAXVAL)/(float)JMAXVAL;
             for (i = 0; i < npixels; ++i)
-                image_buffer[i] = (int)((log(abs(data[i]) + 1.0))/scl_data);
+                (*image_buffer)[i] = (int)((log(abs(data[i]) + 1.0))/scl_data);
             break;
 
-        case 5 :
-            /* contrast stretch */
-            printinfo("Performing histogram stretch (normalization)");
-
+        case 5 :                                         /* contrast stretch */
             /* We need to go through the cumulative histogram to pick the
-             *  appropriate values for datamin and datamax               */
+             *  appropriate values for datamin and datamax                   */
             i = 0;
             while (i < JMAXVAL)
             {
@@ -141,33 +172,32 @@ void scale_pixels(int scale, unsigned int npixels,
             scl_data = (datamax - datamin)/(float)JMAXVAL;
             for (i = 0; i < npixels; ++i)
             {
-                if (image_buffer[i] >= datamax)
-                    image_buffer[i] = JMAXVAL;
-                else if (image_buffer[i] <= datamin)
-                    image_buffer[i] = 0;
+                if ((*image_buffer)[i] >= datamax)
+                    (*image_buffer)[i] = JMAXVAL;
+                else if ((*image_buffer)[i] <= datamin)
+                    (*image_buffer)[i] = 0;
                 else
-                    image_buffer[i] = (int) abs((image_buffer[i]
+                    (*image_buffer)[i] = (int) abs(((*image_buffer)[i]
                                     - datamin)/scl_data);
             }
             break;
 
-        case 6 :
-            /* histogram equalization */
-            printinfo("Performing Histogram Equalization");
+        case 6 :                                   /* histogram equalization */
             for (i = 0; i <  npixels; ++i)
-                image_buffer[i] = cumhist[image_buffer[i]] * JMAXVAL;
+                (*image_buffer)[i] = cumhist[(*image_buffer)[i]] * JMAXVAL;
             break;
         default :
-            printinfo("Using linear scale");
             break;
     }
-
 }
 
 /*---------------------------------------------------------------------------*
- * RESIZE_IMAGE: Scales down/up the image_buffer
+ * RESIZE_IMAGE: Scales down/up the image_buffer using bilinear scaling
+ * Based on an article by "John" at
+ * http://tech-algorithm.com/articles/bilinear-image-scaling/
  *---------------------------------------------------------------------------*/
-void resize_image(long *xdim, long *ydim, float zoomfact, JSAMPLE *image_buffer)
+void resize_image(long *xdim, long *ydim, float zoomfact,
+                  JSAMPLE ** image_buffer)
 {
     int offset = 0, index = 0;
     int A, B, C, D, x, y, gray;
@@ -186,8 +216,8 @@ void resize_image(long *xdim, long *ydim, float zoomfact, JSAMPLE *image_buffer)
     xratio = ((float)(w - 1))/zxdim;
     yratio = ((float)(h - 1))/zydim;
 
-                                     /* allocate space for *buff */
-    buff   = (unsigned char *) malloc(sizeof(char) * zxdim * zydim);
+                            /* allocate space for *buff */
+    buff   = malloc(sizeof(unsigned char) * zxdim * zydim);
 
     index  = 0;
     offset = 0;
@@ -203,10 +233,10 @@ void resize_image(long *xdim, long *ydim, float zoomfact, JSAMPLE *image_buffer)
             xdiff = (xratio * j) - x;
             index = y * w + x;
 
-            A = image_buffer[index]         & 0xff;
-            B = image_buffer[index + 1]     & 0xff;
-            C = image_buffer[index + w]     & 0xff;
-            D = image_buffer[index + w + 1] & 0xff;
+            A = (*image_buffer)[index]         & 0xff;
+            B = (*image_buffer)[index + 1]     & 0xff;
+            C = (*image_buffer)[index + w]     & 0xff;
+            D = (*image_buffer)[index + w + 1] & 0xff;
 
             gray = (int)(A * (1 - xdiff) * (1 - ydiff)
                  +       B * (xdiff)     * (1 - ydiff)
@@ -218,10 +248,8 @@ void resize_image(long *xdim, long *ydim, float zoomfact, JSAMPLE *image_buffer)
     }
     *xdim = zxdim;
     *ydim = zydim;
-    image_buffer = realloc(image_buffer, sizeof(char) * npixels);
-    if (!image_buffer)
-        printerro("Failed to allocate memory");
-
+    (*image_buffer) = realloc((*image_buffer), sizeof(unsigned char) * npixels);
     for (i = 0; i <  npixels; ++i)
-        image_buffer[i] = buff[i];
+        (*image_buffer)[i] = buff[i];
+    free(buff);
 }