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