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 /*---------------------------------------------------------------------------*
  28  * SCALE_PIXELS: Changes the pixel scale to linear/log/sqroot/etc..
  29  *---------------------------------------------------------------------------*/
  30 void scale_pixels(int scale, unsigned int npixels, float *data,
  31                  JSAMPLE ** image_buffer)
  32 {
  33     unsigned int i = 0;
  34     int JMAXVAL = 255;
  35     float datamax = 0.0, datamin = 0.0, tmp = 0.0;
  36     float hist[256] = {0.0}, cumhist[256] = {0.0};
  37     float scl_data = 0.0;
  38
  39
  40     /* first find min & max in data                                          */
  41     datamax = -1.0 * FLT_MAX;
  42     datamin = FLT_MAX;
  43     for (i = 0; i < npixels; ++i)
  44     {
  45         if (data[i] > datamax) datamax = data[i];
  46         if (data[i] < datamin) datamin = data[i];
  47     } /*endfor*/
  48
  49
  50     /* Convert data into bytscaled values for jpeg file                     */
  51     /* the dynamic range is reduced to 255 for jpeg                         */
  52     scl_data = (datamax - datamin)/(float)JMAXVAL;
  53
  54     for (i = 0; i < npixels; ++i)
  55         data[i] = (data[i] - datamin)/scl_data;
  56
  57
  58     /* All data is now squeezed into the range 0 - 255                       */
  59     /* NOTE: At this point onwards min & max is 0 and 255 respectively       */
  60     datamax = (float)JMAXVAL;
  61     datamin = 0.0;
  62
  63     /* initialize image histogram. ensure all are zeroes in hist[]           */
  64     /*-----------------------------------------------------------------------*/
  65     for (i = 0; i <= JMAXVAL; ++i) hist[i] = 0;
  66
  67     /* construct the image histogram */
  68     tmp = 1.0/(float)npixels;
  69     for (i = 0; i <= npixels; ++i)
  70         hist[(int)floor(data[i])] += tmp;
  71
  72     /* And the cumulative histogram */
  73     cumhist[0] = hist[0];
  74     for (i = 1; i <= JMAXVAL; ++i)
  75         cumhist[i] += cumhist[i - 1] + hist[i];
  76
  77     /* Allocate image buffer */
  78     (*image_buffer) = malloc(sizeof(unsigned char) * npixels);
  79
  80
  81     /* Linear scale (min-max) : This is the default scaling
  82      * histo-eq will fail if we dont generate image_buffer here              */
  83     for (i = 0; i < npixels; ++i)
  84         (*image_buffer)[i] = (int)(data[i]);
  85
  86     /*-----------------------------------------------------------------------*/
  87
  88
  89     switch (scale)
  90     {
  91         case 1 :                                              /* Square root */
  92             printinfo("Using square-root scale");
  93             scl_data = sqrt((float)JMAXVAL)/(float)JMAXVAL;
  94             for (i = 0; i < npixels; ++i)
  95                 (*image_buffer)[i] = (int)(sqrt(data[i])/scl_data);
  96             break;
  97
  98         case 2 :                                                   /* Square */
  99             printinfo("Using quadratic scale");
 100             scl_data = pow((float)JMAXVAL,2)/(float)JMAXVAL;
 101             for (i = 0; i < npixels; ++i)
 102                 (*image_buffer)[i] = (int)abs((pow(data[i],2) - 1.0)/scl_data);
 103             break;
 104
 105         case 3 :                                                    /* Cubic */
 106             printinfo("Using cubic scale");
 107             scl_data = pow((float)JMAXVAL,3)/(float)JMAXVAL;
 108             for (i = 0; i < npixels; ++i)
 109                 (*image_buffer)[i] = (int)abs((pow(data[i],3) - 1.0)/scl_data);
 110             break;
 111
 112         case 4 :                                                      /* log */
 113             printinfo("Using log scale");
 114             scl_data = log(1.0 + (float)JMAXVAL)/(float)JMAXVAL;
 115             for (i = 0; i < npixels; ++i)
 116                 (*image_buffer)[i] = (int)((log(abs(data[i]) + 1.0))/scl_data);
 117             break;
 118
 119         case 5 :
 120             /* contrast stretch */
 121             printinfo("Performing histogram stretch (normalization)");
 122
 123             /* We need to go through the cumulative histogram to pick the
 124              *  appropriate values for datamin and datamax               */
 125             i = 0;
 126             while (i < JMAXVAL)
 127             {
 128                 if (cumhist[i] >= 0.01)
 129                 {
 130                     datamin = (float) i;
 131                     break;
 132                 }
 133                 i++;
 134             }
 135             i = JMAXVAL;
 136             while (i > 0)
 137             {
 138                 if (cumhist[i] <= 0.99)
 139                 {
 140                     datamax = (float) i;
 141                     break;
 142                 }
 143                 i--;
 144             }
 145             scl_data = (datamax - datamin)/(float)JMAXVAL;
 146             for (i = 0; i < npixels; ++i)
 147             {
 148                 if ((*image_buffer)[i] >= datamax)
 149                     (*image_buffer)[i] = JMAXVAL;
 150                 else if ((*image_buffer)[i] <= datamin)
 151                     (*image_buffer)[i] = 0;
 152                 else
 153                     (*image_buffer)[i] = (int) abs(((*image_buffer)[i]
 154                                     - datamin)/scl_data);
 155             }
 156             break;
 157
 158         case 6 :
 159             /* histogram equalization */
 160             printinfo("Performing Histogram Equalization");
 161             for (i = 0; i <  npixels; ++i)
 162                 (*image_buffer)[i] = cumhist[(*image_buffer)[i]] * JMAXVAL;
 163             break;
 164         default :
 165             printinfo("Using linear scale");
 166             break;
 167     }
 168 }
 169
 170 /*---------------------------------------------------------------------------*
 171  * RESIZE_IMAGE: Scales down/up the image_buffer using bilinear scaling
 172  * Based on an article by "John" at
 173  * http://tech-algorithm.com/articles/bilinear-image-scaling/
 174  *---------------------------------------------------------------------------*/
 175 void resize_image(long *xdim, long *ydim, float zoomfact,
 176                   JSAMPLE ** image_buffer)
 177 {
 178     int offset = 0, index = 0;
 179     int A, B, C, D, x, y, gray;
 180     JSAMPLE *buff;
 181     unsigned int i = 0, j = 0;
 182     unsigned long npixels = 0;
 183     long w = *xdim, h = *ydim;
 184     long zxdim = 0, zydim = 0;
 185     float xdiff, ydiff, xratio, yratio;
 186
 187     zxdim  = (int)(w * zoomfact);
 188     zydim  = (int)(h * zoomfact);
 189
 190     npixels= zxdim * zydim;
 191
 192     xratio = ((float)(w - 1))/zxdim;
 193     yratio = ((float)(h - 1))/zydim;
 194
 195                             /* allocate space for *buff */
 196     buff   = malloc(sizeof(unsigned char) * zxdim * zydim);
 197
 198     index  = 0;
 199     offset = 0;
 200     for (i = 0; i < zydim; i++)
 201     {
 202         y     = (int)(yratio * i);
 203         ydiff = (yratio * i) - y;
 204
 205         for (j = 0; j < zxdim; j++)
 206         {
 207             x = (int)(xratio * j);
 208
 209             xdiff = (xratio * j) - x;
 210             index = y * w + x;
 211
 212             A = (*image_buffer)[index]         & 0xff;
 213             B = (*image_buffer)[index + 1]     & 0xff;
 214             C = (*image_buffer)[index + w]     & 0xff;
 215             D = (*image_buffer)[index + w + 1] & 0xff;
 216
 217             gray = (int)(A * (1 - xdiff) * (1 - ydiff)
 218                  +       B * (xdiff)     * (1 - ydiff)
 219                  +       C * (ydiff)     * (1 - xdiff)
 220                  +       D * (xdiff)     * (ydiff)
 221                     );
 222             buff[offset++] = gray;
 223         }
 224     }
 225     *xdim = zxdim;
 226     *ydim = zydim;
 227     (*image_buffer) = realloc((*image_buffer), sizeof(unsigned char) * npixels);
 228     for (i = 0; i <  npixels; ++i)
 229         (*image_buffer)[i] = buff[i];
 230     free(buff);
 231 }