diff --git a/src/image.c b/src/image.c
new file mode 100644 (file)
index 0000000..7a068e9
--- /dev/null
@@ -0,0 +1,158 @@
+/***************************************************************************
+ * This file is a part of CADS/UVS fits2jpeg conversion software           *
+ *   Copyright (C) 2012 by CADS/UV Software Team,                          *
+ *                         Indian Institute of Astrophysics                *
+ *                         Bangalore 560034                                *
+ *                         cads_AT_iiap.res.in                             *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ *   This program is distributed in the hope that it will be useful,       *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU General Public License for more details.                          *
+ *                                                                         *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+
+/*Header Definitions*/
+#include "fits2jpeg.h"
+
+void scale_image(int scale, int npixels,
+                 float *data, JSAMPLE *image_buffer)
+{
+    unsigned int i = 0;
+    int JMAXVAL = 255;
+    float datamax = 0.0, datamin = 0.0, tmp = 0.0;
+    float hist[256] = {0.0}, cumhist[256] = {0.0};
+    float scl_data = 0.0;
+
+
+    /* first find min & max in data */
+    datamax = -1.0 * FLT_MAX;
+    datamin = FLT_MAX;
+    for (i = 0; i < npixels; ++i)
+    {
+        if (data[i] > datamax) datamax = data[i];
+        if (data[i] < datamin) datamin = data[i];
+    } /*endfor*/
+
+    /* Convert data into bytscaled values for jpeg file                 */
+    /* the dynamic range is reduced to 255 for jpeg                     */
+    scl_data = (datamax - datamin)/(float)JMAXVAL;
+    for (i = 0; i < npixels; ++i)
+        data[i] = (data[i] - datamin)/scl_data;
+
+    /* All data is now squeezed into the range 0 - 255                   */
+    /* NOTE: At this point onwards min & max is 0 and 255 respectively   */
+    datamax = (float)JMAXVAL;
+    datamin = 0.0;
+
+    /* initialize image histogram. ensure all are zeroes in hist[]       */
+    /*-------------------------------------------------------------------*/
+    for (i = 0; i <= JMAXVAL; ++i) hist[i] = 0;
+
+    /* construct the image histogram */
+    tmp = 1.0/(float)npixels;
+    for (i = 0; i <= npixels; ++i)
+        hist[(int)floor(data[i])] += tmp;
+
+    /* And the cumulative histogram */
+    cumhist[0] = hist[0];
+    for (i = 1; i <= JMAXVAL; ++i)
+        cumhist[i] += cumhist[i - 1] + hist[i];
+
+    /* Linear scale (min-max) : This is the default scaling
+     * if we dont generate image_buffer here, histo-eq will fail */
+    for (i = 0; i < npixels; ++i)
+        image_buffer[i] = (int)(data[i]);
+
+    /*-------------------------------------------------------------------*/
+
+
+    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);
+            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);
+            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);
+            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);
+            break;
+
+        case 5 :
+            /* contrast stretch */
+            printinfo("Performing histogram stretch (normalization)");
+
+            /* We need to go through the cumulative histogram to pick the
+             *  appropriate values for datamin and datamax               */
+            i = 0;
+            while (i < JMAXVAL)
+            {
+                if (cumhist[i] >= 0.01)
+                {
+                    datamin = (float) i;
+                    break;
+                }
+                i++;
+            }
+            i = JMAXVAL;
+            while (i > 0)
+            {
+                if (cumhist[i] <= 0.99)
+                {
+                    datamax = (float) i;
+                    break;
+                }
+                i--;
+            }
+            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;
+                else
+                    image_buffer[i] = (int) abs((image_buffer[i] - datamin)/scl_data);
+            }
+            break;
+
+        case 6 :
+            /* histogram equalization */
+            printinfo("Performing Histogram Equalization");
+            for (i = 0; i <  npixels; ++i)
+                image_buffer[i] = cumhist[image_buffer[i]] * 255;
+            break;
+        default :
+            printinfo("Using linear scale");
+            break;
+    }
+
+}