src/image.c

   1 /***************************************************************************
   2  * This file is a part of CADS/UVS fits2jpeg conversion software           *
   3  *   Copyright (C) 2012 by CADS/UV Software Team,                          *
   4  *                         Indian Institute of Astrophysics                *
   5  *                         Bangalore 560034                                *
   6  *                         cads_AT_iiap.res.in                             *
   7  *                                                                         *
   8  *   This program is free software; you can redistribute it and/or modify  *
   9  *   it under the terms of the GNU General Public License as published by  *
  10  *   the Free Software Foundation; either version 2 of the License, or     *
  11  *   (at your option) any later version.                                   *
  12  *                                                                         *
  13  *   This program is distributed in the hope that it will be useful,       *
  14  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
  15  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
  16  *   GNU General Public License for more details.                          *
  17  *                                                                         *
  18  *   You should have received a copy of the GNU General Public License     *
  19  *   along with this program; if not, write to the                         *
  20  *   Free Software Foundation, Inc.,                                       *
  21  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
  22  ***************************************************************************/
  23
  24 /*Header Definitions*/
  25 #include "fits2jpeg.h"
  26
  27 void scale_image(int scale, int npixels,
  28                  float *data, JSAMPLE *image_buffer)
  29 {
  30     unsigned int i = 0;
  31     int JMAXVAL = 255;
  32     float datamax = 0.0, datamin = 0.0, tmp = 0.0;
  33     float hist[256] = {0.0}, cumhist[256] = {0.0};
  34     float scl_data = 0.0;
  35
  36
  37     /* first find min & max in data */
  38     datamax = -1.0 * FLT_MAX;
  39     datamin = FLT_MAX;
  40     for (i = 0; i < npixels; ++i)
  41     {
  42         if (data[i] > datamax) datamax = data[i];
  43         if (data[i] < datamin) datamin = data[i];
  44     } /*endfor*/
  45
  46     /* Convert data into bytscaled values for jpeg file                 */
  47     /* the dynamic range is reduced to 255 for jpeg                     */
  48     scl_data = (datamax - datamin)/(float)JMAXVAL;
  49     for (i = 0; i < npixels; ++i)
  50         data[i] = (data[i] - datamin)/scl_data;
  51
  52     /* All data is now squeezed into the range 0 - 255                   */
  53     /* NOTE: At this point onwards min & max is 0 and 255 respectively   */
  54     datamax = (float)JMAXVAL;
  55     datamin = 0.0;
  56
  57     /* initialize image histogram. ensure all are zeroes in hist[]       */
  58     /*-------------------------------------------------------------------*/
  59     for (i = 0; i <= JMAXVAL; ++i) hist[i] = 0;
  60
  61     /* construct the image histogram */
  62     tmp = 1.0/(float)npixels;
  63     for (i = 0; i <= npixels; ++i)
  64         hist[(int)floor(data[i])] += tmp;
  65
  66     /* And the cumulative histogram */
  67     cumhist[0] = hist[0];
  68     for (i = 1; i <= JMAXVAL; ++i)
  69         cumhist[i] += cumhist[i - 1] + hist[i];
  70
  71     /* Linear scale (min-max) : This is the default scaling
  72      * if we dont generate image_buffer here, histo-eq will fail */
  73     for (i = 0; i < npixels; ++i)
  74         image_buffer[i] = (int)(data[i]);
  75
  76     /*-------------------------------------------------------------------*/
  77
  78
  79     switch (scale)
  80     {
  81         case 1 :                                              /* Square root */
  82             printinfo("Using square-root scale");
  83             scl_data = sqrt((float)JMAXVAL)/(float)JMAXVAL;
  84             for (i = 0; i < npixels; ++i)
  85                 image_buffer[i] = (int)(sqrt(data[i])/scl_data);
  86             break;
  87
  88         case 2 :                                                   /* Square */
  89             printinfo("Using quadratic scale");
  90             scl_data = pow((float)JMAXVAL,2)/(float)JMAXVAL;
  91             for (i = 0; i < npixels; ++i)
  92                 image_buffer[i] = (int)abs((pow(data[i],2) - 1.0)/scl_data);
  93             break;
  94
  95         case 3 :                                                    /* Cubic */
  96             printinfo("Using cubic scale");
  97             scl_data = pow((float)JMAXVAL,3)/(float)JMAXVAL;
  98             for (i = 0; i < npixels; ++i)
  99                 image_buffer[i] = (int)abs((pow(data[i],3) - 1.0)/scl_data);
 100             break;
 101
 102         case 4 :                                                      /* log */
 103             printinfo("Using log scale");
 104             scl_data = log(1.0 + (float)JMAXVAL)/(float)JMAXVAL;
 105             for (i = 0; i < npixels; ++i)
 106                 image_buffer[i] = (int)((log(abs(data[i]) + 1.0))/scl_data);
 107             break;
 108
 109         case 5 :
 110             /* contrast stretch */
 111             printinfo("Performing histogram stretch (normalization)");
 112
 113             /* We need to go through the cumulative histogram to pick the
 114              *  appropriate values for datamin and datamax               */
 115             i = 0;
 116             while (i < JMAXVAL)
 117             {
 118                 if (cumhist[i] >= 0.01)
 119                 {
 120                     datamin = (float) i;
 121                     break;
 122                 }
 123                 i++;
 124             }
 125             i = JMAXVAL;
 126             while (i > 0)
 127             {
 128                 if (cumhist[i] <= 0.99)
 129                 {
 130                     datamax = (float) i;
 131                     break;
 132                 }
 133                 i--;
 134             }
 135             scl_data = (datamax - datamin)/(float)JMAXVAL;
 136             for (i = 0; i < npixels; ++i)
 137             {
 138                 if (image_buffer[i] >= datamax)
 139                     image_buffer[i] = JMAXVAL;
 140                 else if (image_buffer[i] <= datamin)
 141                     image_buffer[i] = 0;
 142                 else
 143                     image_buffer[i] = (int) abs((image_buffer[i] - datamin)/scl_data);
 144             }
 145             break;
 146
 147         case 6 :
 148             /* histogram equalization */
 149             printinfo("Performing Histogram Equalization");
 150             for (i = 0; i <  npixels; ++i)
 151                 image_buffer[i] = cumhist[image_buffer[i]] * 255;
 152             break;
 153         default :
 154             printinfo("Using linear scale");
 155             break;
 156     }
 157
 158 }