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     for (i = 0; i < npixels; ++i)
  83         data[i] = (data[i] - datamin)/scl_data;
  84
  85
  86     /* All data is now squeezed into the range 0 - 255                       */
  87     /* NOTE: At this point onwards min & max is 0 and 255 respectively       */
  88     datamax = (float)JMAXVAL;
  89     datamin = 0.0;
  90
  91     /* initialize image histogram. ensure all are zeroes in hist[]           */
  92     /*-----------------------------------------------------------------------*/
  93     for (i = 0; i <= JMAXVAL; ++i) hist[i] = 0;
  94
  95     /* construct the image histogram */
  96     tmp = 1.0/(float)npixels;
  97     for (i = 0; i <= npixels; ++i)
  98         hist[(int)floor(data[i])] += tmp;
  99
 100     /* And the cumulative histogram */
 101     cumhist[0] = hist[0];
 102     for (i = 1; i <= JMAXVAL; ++i)
 103         cumhist[i] += cumhist[i - 1] + hist[i];
 104
 105     /* Allocate image buffer */
 106     (*image_buffer) = malloc(sizeof(unsigned char) * npixels);
 107
 108
 109     /* Linear scale (min-max) : This is the default scaling
 110      * histo-eq will fail if we dont generate image_buffer here              */
 111     for (i = 0; i < npixels; ++i)
 112         (*image_buffer)[i] = (int)(data[i]);
 113
 114     /*-----------------------------------------------------------------------*/
 115
 116
 117     switch (scale)
 118     {
 119         case 1 :                                              /* Square root */
 120             scl_data = sqrt((float)JMAXVAL)/(float)JMAXVAL;
 121             for (i = 0; i < npixels; ++i)
 122                 (*image_buffer)[i] = (int)(sqrt(data[i])/scl_data);
 123             break;
 124
 125         case 2 :                                                   /* Square */
 126             scl_data = pow((float)JMAXVAL,2)/(float)JMAXVAL;
 127             for (i = 0; i < npixels; ++i)
 128                 (*image_buffer)[i] = (int)abs((pow(data[i],2) - 1.0)/scl_data);
 129             break;
 130
 131         case 3 :                                                    /* Cubic */
 132             scl_data = pow((float)JMAXVAL,3)/(float)JMAXVAL;
 133             for (i = 0; i < npixels; ++i)
 134                 (*image_buffer)[i] = (int)abs((pow(data[i],3) - 1.0)/scl_data);
 135             break;
 136
 137         case 4 :                                                      /* log */
 138             scl_data = log(1.0 + (float)JMAXVAL)/(float)JMAXVAL;
 139             for (i = 0; i < npixels; ++i)
 140                 (*image_buffer)[i] = (int)((log(abs(data[i]) + 1.0))/scl_data);
 141             break;
 142
 143         case 5 :                                         /* contrast stretch */
 144             /* We need to go through the cumulative histogram to pick the
 145              *  appropriate values for datamin and datamax                   */
 146             i = 0;
 147             while (i < JMAXVAL)
 148             {
 149                 if (cumhist[i] >= 0.01)
 150                 {
 151                     datamin = (float) i;
 152                     break;
 153                 }
 154                 i++;
 155             }
 156             i = JMAXVAL;
 157             while (i > 0)
 158             {
 159                 if (cumhist[i] <= 0.99)
 160                 {
 161                     datamax = (float) i;
 162                     break;
 163                 }
 164                 i--;
 165             }
 166             scl_data = (datamax - datamin)/(float)JMAXVAL;
 167             for (i = 0; i < npixels; ++i)
 168             {
 169                 if ((*image_buffer)[i] >= datamax)
 170                     (*image_buffer)[i] = JMAXVAL;
 171                 else if ((*image_buffer)[i] <= datamin)
 172                     (*image_buffer)[i] = 0;
 173                 else
 174                     (*image_buffer)[i] = (int) abs(((*image_buffer)[i]
 175                                     - datamin)/scl_data);
 176             }
 177             break;
 178
 179         case 6 :                                   /* histogram equalization */
 180             for (i = 0; i <  npixels; ++i)
 181                 (*image_buffer)[i] = cumhist[(*image_buffer)[i]] * JMAXVAL;
 182             break;
 183         default :
 184             break;
 185     }
 186 }
 187
 188 /*---------------------------------------------------------------------------*
 189  * RESIZE_IMAGE: Scales down/up the image_buffer using bilinear scaling
 190  * Based on an article by "John" at
 191  * http://tech-algorithm.com/articles/bilinear-image-scaling/
 192  *---------------------------------------------------------------------------*/
 193 void resize_image(long *xdim, long *ydim, float zoomfact,
 194                   JSAMPLE ** image_buffer)
 195 {
 196     int offset = 0, index = 0;
 197     int A, B, C, D, x, y, gray;
 198     JSAMPLE *buff;
 199     unsigned int i = 0, j = 0;
 200     unsigned long npixels = 0;
 201     long w = *xdim, h = *ydim;
 202     long zxdim = 0, zydim = 0;
 203     float xdiff, ydiff, xratio, yratio;
 204
 205     zxdim  = (int)(w * zoomfact);
 206     zydim  = (int)(h * zoomfact);
 207
 208     npixels= zxdim * zydim;
 209
 210     xratio = ((float)(w - 1))/zxdim;
 211     yratio = ((float)(h - 1))/zydim;
 212
 213                             /* allocate space for *buff */
 214     buff   = malloc(sizeof(unsigned char) * zxdim * zydim);
 215
 216     index  = 0;
 217     offset = 0;
 218     for (i = 0; i < zydim; i++)
 219     {
 220         y     = (int)(yratio * i);
 221         ydiff = (yratio * i) - y;
 222
 223         for (j = 0; j < zxdim; j++)
 224         {
 225             x = (int)(xratio * j);
 226
 227             xdiff = (xratio * j) - x;
 228             index = y * w + x;
 229
 230             A = (*image_buffer)[index]         & 0xff;
 231             B = (*image_buffer)[index + 1]     & 0xff;
 232             C = (*image_buffer)[index + w]     & 0xff;
 233             D = (*image_buffer)[index + w + 1] & 0xff;
 234
 235             gray = (int)(A * (1 - xdiff) * (1 - ydiff)
 236                  +       B * (xdiff)     * (1 - ydiff)
 237                  +       C * (ydiff)     * (1 - xdiff)
 238                  +       D * (xdiff)     * (ydiff)
 239                     );
 240             buff[offset++] = gray;
 241         }
 242     }
 243     *xdim = zxdim;
 244     *ydim = zydim;
 245     (*image_buffer) = realloc((*image_buffer), sizeof(unsigned char) * npixels);
 246     for (i = 0; i <  npixels; ++i)
 247         (*image_buffer)[i] = buff[i];
 248     free(buff);
 249 }