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