src/fits2jpeg.c

   1 /***************************************************************************
   2  *   Copyright (C) 2012 by CADS/UV Software Team,                          *
   3  *                         Indian Institute of Astrophysics                *
   4  *                         Bangalore 560034                                *
   5  *                         cads_AT_iiap.res.in                             *
   6  *                                                                         *
   7  *   This program is free software; you can redistribute it and/or modify  *
   8  *   it under the terms of the GNU General Public License as published by  *
   9  *   the Free Software Foundation; either version 2 of the License, or     *
  10  *   (at your option) any later version.                                   *
  11  *                                                                         *
  12  *   This program is distributed in the hope that it will be useful,       *
  13  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
  14  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
  15  *   GNU General Public License for more details.                          *
  16  *                                                                         *
  17  *   You should have received a copy of the GNU General Public License     *
  18  *   along with this program; if not, write to the                         *
  19  *   Free Software Foundation, Inc.,                                       *
  20  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
  21  ***************************************************************************/
  22 /*---------------------------------------------------------------------------
  23   This program will write out FITS data into a jpeg file.
  24   Assume that data are in first extension and are two dimensional
  25
  26   Libraries needed are CFITSIO, JPEG
  27   Original author: Jayant Murthy
  28   Last Modification: 11/05/02
  29
  30   Reconstruction and additional features: Reks
  31  *---------------------------------------------------------------------------*/
  32 /*---------------------------------------------------------------------------
  33   Modification History
  34
  35   Reks, 5th October 2007: v0.01
  36   1. configure checks for jpeglib.
  37      optional extra arguments --with-jpeglib=<location>
  38                               --with-cfitsio=<location>
  39   2. Makes blank images! datamin and max are not set..
  40
  41   Reks, 6th October 2007: v0.10
  42   1. Added 2 functions to fetch max and min from image
  43   2. user can choose various scales like log/square/square-root
  44   3. Increased verbosity, added banner, help message, etc
  45   4. Inputs through getopt. will Ask if mandatory input (fits file) is not
  46      provided.
  47   5. fitsio2
  48
  49   Reks, 10th October 2007: v0.9.0 (close to 1.0)
  50   1. Found our niche :-) Can handle wild cards at command line,
  51      batch processing.
  52   2. Now users cannot specify the output file name though.. Its fixed to
  53      <input_file_root>.jpg.
  54   3. No getopt for input file.. getopt is limited to various image scaling
  55      options.
  56
  57  Reks, 17th October 2007, v0.9.5
  58   1. Merged all image scaling options into one input param: -s
  59      optarg (from getopt) will decide scaling function. Now there is
  60      enough room for additional flags (in future).
  61
  62  Reks, 22 June 2010 v1.0.0
  63   1. Removed datamax() and datamin(). That was required only once!
  64   2. option histoeq shifted to image enhancement option
  65   3. Added "normalize" option
  66      Now users can do scaling + enhance (if at all someone feels like it)
  67   4. All options tested. First stable release :)
  68
  69   Reks 28 June 2012
  70   1. rebranded to CADS/UVS
  71  *-------------------------------------------------------------------------*/
  72
  73 #ifdef HAVE_CONFIG_H
  74 #include <config.h>
  75 #endif
  76
  77 /*Header Definitions*/
  78 #include <stdio.h>
  79 #include <stdlib.h>
  80 #include <string.h>
  81 #include <math.h>
  82 #include <libgen.h>
  83 #include <dirent.h>
  84 #include <fitsio2.h>
  85 #include <jpeglib.h>
  86
  87 #define PROGRAM "fits2jpeg"
  88 #define r2d (90./atan2(1,0))
  89 #define MAX_TEXT 150
  90 int my_getopt(int, char * const *, const char *);
  91 void signals_handler(int);
  92 void set_signals(void);
  93
  94 char *optarg;
  95 int optindex;
  96 /*---------------------------------------------------------------------------*/
  97 void PRINT_BANNER()
  98 {
  99     printf("\n");
 100     printf(" %s %s\n", PROGRAM, VERSION);
 101     printf(" Converts fits image files to jpeg\n");
 102     printf("-----------------------------------------------------------\n");
 103 }
 104 /*---------------------------------------------------------------------------*/
 105 void usage()
 106 {
 107     printf ("\n Usage: fits2jpeg [-s <options>] <fits_file>             \n");
 108     printf ("  Options are:                                             \n");
 109     printf ("    -h help                                                \n");
 110     printf ("    -s <scale_type>                                        \n");
 111     printf ("       scale for output image, where <scale_type> can be:  \n");
 112     printf ("         linear         Linear scale, default            \n\n");
 113     printf ("         sqroot         for square root scale              \n");
 114     printf ("         square         for quadratic scale                \n");
 115     printf ("         cubic          for cubic scale                    \n");
 116     printf ("         log            for log scale                      \n");
 117     printf ("    -e <operation>                                         \n");
 118     printf ("       Image enhancements, where <operation> can be:       \n");
 119     printf ("         normalize      for linear histogram stretch       \n");
 120     printf ("         equalize       for histogram equalization         \n");
 121     printf ("  Output will be written to <fits_file_root>.jpg. Wild card\n");
 122     printf ("  entries allowed in <fits_file>; For eg: *.fits, m31*.fits\n");
 123     printf ("  ngc???.fits etc.                                         \n");
 124     printf ("-----------------------------------------------------------\n");
 125     printf ("Copyright (c) 2012 The CADS Team, IIAp. [GPL v3.0 or later]\n");
 126     printf ("Report Bugs to: http://cads.iiap.res.in/bugzilla           \n");
 127     printf ("Documentation : http://cads.iiap.res.in/software         \n\n");
 128     exit(1);
 129 }
 130 /*---------------------------------------------------------------------------*/
 131
 132 /*--------------------------- Begin Main Progam -----------------------------*/
 133 int main(int argc, char *argv[], char *envp[])
 134 {
 135     fitsfile *fptr;
 136
 137     int scale = 0, process = 0, status = 0, nfound, anynull;
 138     unsigned int i = 0, j = 0, npixels = 0;
 139     long naxes[2], row_stride;
 140     float datamin = 0.0, datamax = 0.0;
 141     float nullval = 0.0, scl_data;
 142     float *data;
 143     float tmp = 0.0;
 144     float hist[256] = {0.0}, cumhist[256] = {0.0};
 145     char jpeg_file_name[MAX_TEXT] = "", fits_file_name[MAX_TEXT] = "";
 146
 147     int opt;                                                   /* For getopt */
 148     extern char *optarg;
 149     extern int optindex;
 150
 151     FILE *jpeg_file;
 152
 153     struct jpeg_error_mgr jerr;
 154     struct jpeg_compress_struct cinfo;
 155     int JMAXVAL = 255;                    /* Max value for greyscale in jpeg */
 156     JSAMPROW row_pointer[1];
 157     JSAMPLE *image_buffer;
 158     /*---------------- End of variable initialization -----------------------*/
 159
 160     PRINT_BANNER();        /* Header: Talks about program, version & purpose */
 161     set_signals();                                  /* Trap un-natural exits */
 162
 163     /*-------------------- Parse command line inputs ------------------------*/
 164
 165     if (argc < 2) usage();  /* People who don't give any arguments need help */
 166     while ( ( opt = my_getopt ( argc, argv, "s:e:h" ) ) != EOF )
 167         switch ( opt )
 168         {
 169         case 's':
 170             if (strcmp(optarg, "sqroot")   ==0) scale = 1;
 171             if (strcmp(optarg, "square")   ==0) scale = 2;
 172             if (strcmp(optarg, "cube")     ==0) scale = 3;
 173             if (strcmp(optarg, "log")      ==0) scale = 4;
 174             break;
 175
 176         case 'e':
 177             if (strcmp(optarg, "normalize")==0) process = 1;
 178             if (strcmp(optarg, "equalize") ==0) process = 2;
 179             break;
 180
 181         case 'h':
 182             usage();
 183             break;
 184
 185         case '*':
 186             usage();
 187             break;
 188         }
 189
 190     if ((argc - optindex) == 0) usage();       /* Somebody gave only options */
 191
 192     /*---------------------- Begin processing files -------------------------*/
 193
 194     /* Getopt would have processed all argument with a "-" in front. optindex
 195        is the counter, which keeps the number of options supplied. Whatever
 196        remains in the list of command line arguments are the fits files.
 197        NOTE: POSIX systems will do the regex matching and expand wild card
 198              entries for us. Even windoze XP Pro is known to do that...      */
 199     for (j = optindex; j < argc; j++)
 200     {
 201
 202         strcpy(fits_file_name, argv[j]);
 203
 204         strcpy(jpeg_file_name, fits_file_name);
 205         strtok(jpeg_file_name, ".");
 206         strcat(jpeg_file_name, ".jpg");
 207
 208         fits_open_file(&fptr, fits_file_name, READONLY, &status);
 209         fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxes, &nfound, &status);
 210         if (status)
 211         {
 212             fprintf(stderr, "ERROR  : Could not open file: %s\n",
 213                     fits_file_name);
 214             exit(EXIT_FAILURE);
 215         }/*Endif*/
 216
 217         /* Read in data */
 218         npixels = naxes[0] * naxes[1];
 219         data = (float *) malloc(sizeof(float) * npixels);
 220
 221         nullval = 0;
 222         if (fits_read_img(fptr, TFLOAT, 1, npixels, &nullval, data, &anynull,
 223                           &status))
 224         {
 225             fprintf(stderr, "ERROR  : No valid fits image data in %s\n",
 226                     fits_file_name);
 227             exit(EXIT_FAILURE);
 228         }/*Endif*/
 229
 230         fits_close_file(fptr, &status);
 231         /*Close data file*/
 232         printf("INFO   : data input from %s\n", fits_file_name);
 233
 234         /* find min & max in image */
 235         datamax = -1.0e9;
 236         datamin = +1.0e9;
 237         for (i = 0; i < npixels; ++i)
 238         {
 239             if (data[i] > datamax) datamax = data[i];
 240             if (data[i] < datamin) datamin = data[i];
 241         } /*endfor*/
 242
 243         /* Convert data into bytscaled values for jpeg file                 */
 244         /* the dynamic range is reduced to 255 for jpeg                     */
 245         scl_data = (datamax - datamin)/(float)JMAXVAL;
 246         for (i = 0; i < npixels; ++i)
 247             data[i] = (data[i] - datamin)/scl_data;
 248
 249         /* All data is now squeezed into the range 0 - 255                   */
 250         /* NOTE: At this point onwards min & max is 0 and 255 respectively   */
 251
 252
 253         /* Now we need to look at user options..                             */
 254         /* 1. scale it as per user's requirement.                            */
 255         /* 2. image enhancing operations                                     */
 256                                                     /* Allocate image buffer */
 257         image_buffer = (unsigned char *) malloc(sizeof(char) * npixels);
 258         scl_data = 1.0;
 259         switch (scale)
 260         {
 261         case 1 :                                              /* Square root */
 262             printf("INFO   : Using square-root scale\n");
 263             scl_data = sqrt((float)JMAXVAL)/(float)JMAXVAL;
 264             for (i = 0; i < npixels; ++i)
 265                 image_buffer[i] = (int)(sqrt(data[i])/scl_data);
 266             break;
 267
 268         case 2 :                                                   /* Square */
 269             printf("INFO   : Using quadratic scale\n");
 270             scl_data = pow((float)JMAXVAL,2)/(float)JMAXVAL;
 271             for (i = 0; i < npixels; ++i)
 272                 image_buffer[i] = (int)abs((pow(data[i],2) - 1.0)/scl_data);
 273             break;
 274
 275         case 3 :                                                    /* Cubic */
 276             printf("INFO   : Using cubic scale\n");
 277             scl_data = pow((float)JMAXVAL,3)/(float)JMAXVAL;
 278             for (i = 0; i < npixels; ++i)
 279                 image_buffer[i] = (int)abs((pow(data[i],3) - 1.0)/scl_data);
 280             break;
 281
 282         case 4 :                                                      /* log */
 283             printf("INFO   : Using log scale\n");
 284             scl_data = log(1.0 + (float)JMAXVAL)/(float)JMAXVAL;
 285             for (i = 0; i < npixels; ++i)
 286                 image_buffer[i] = (int)((log(abs(data[i]) + 1.0))/scl_data);
 287             break;
 288
 289         default:                                   /* Linear scale (min-max) */
 290             for (i = 0; i < npixels; ++i)
 291                 image_buffer[i] = (int)(data[i]);
 292             break;
 293         }
 294
 295         /* initialize image histogram. ensure all are zeroes in hist[]       */
 296         for (i = 0; i <= JMAXVAL; ++i) hist[i] = 0;
 297
 298         /* construct the image histogram */
 299         tmp = 1.0/(float)npixels;
 300         for (i = 0; i <= npixels; ++i)
 301             hist[(int)floor(data[i])] += tmp;
 302
 303         /* And the cumulative histogram */
 304         cumhist[0] = hist[0];
 305         for (i = 1; i <= JMAXVAL; ++i)
 306             cumhist[i] += cumhist[i - 1] + hist[i];
 307
 308         switch (process)
 309         {
 310         case 1 :                                         /* contrast stretch */
 311             printf("INFO   : Performing contrast stretch (normalization) \n");
 312
 313             datamax = (float)JMAXVAL;
 314             datamin = 0.0;
 315
 316             /* We need to go through the cumulative histogram to pick the
 317                appropriate values for datamin and datamax                    */
 318             i = 0;
 319             while (i < JMAXVAL)
 320             {
 321                 if (cumhist[i] >= 0.01)
 322                 {
 323                     datamin = (float) i;
 324                     break;
 325                 }
 326                 i++;
 327             }
 328             i = JMAXVAL;
 329             while (i > 0)
 330             {
 331                 if (cumhist[i] <= 0.99)
 332                 {
 333                     datamax = (float) i;
 334                     break;
 335                 }
 336                 i--;
 337             }
 338             scl_data = (datamax - datamin)/(float)JMAXVAL;
 339             for (i = 0; i < npixels; ++i)
 340             {
 341                 if (image_buffer[i] >= datamax)
 342                     image_buffer[i] = JMAXVAL;
 343                 else if (image_buffer[i] <= datamin)
 344                     image_buffer[i] = 0;
 345                 else
 346                     image_buffer[i]=(int)abs((image_buffer[i]-datamin)/scl_data);
 347             }
 348             break;
 349
 350
 351         case 2 :                                   /* histogram equalization */
 352             printf("INFO   : Performing Histogram Equalization\n");
 353             for (i = 0; i <  npixels; ++i)
 354                 image_buffer[i] = cumhist[image_buffer[i]]*255;
 355             break;
 356         }
 357
 358         /*Write out data into JPEG file*/
 359         cinfo.err = jpeg_std_error(&jerr);             /* JPEG error handler */
 360         jpeg_create_compress(&cinfo);                 /* JPEG initialization */
 361         jpeg_file = fopen(jpeg_file_name, "wb");/* Open JPEG file for writing*/
 362         jpeg_stdio_dest(&cinfo, jpeg_file); /* Send compressed data to stdio */
 363         cinfo.image_width = naxes[0];                         /* Image width */
 364         cinfo.image_height = naxes[1];                       /* Image height */
 365         cinfo.input_components = 1;            /* Number of colors per pixel */
 366         cinfo.in_color_space = JCS_GRAYSCALE;   /* colorspace of input image */
 367         jpeg_set_defaults(&cinfo);                      /* Set JPEG defaults */
 368         jpeg_start_compress(&cinfo, TRUE);       /* default data compression */
 369         row_stride = naxes[0];            /* JSAMPLEs per row inimage buffer */
 370
 371         /* Now we have to turn the data upside down                          */
 372         while (cinfo.next_scanline < cinfo.image_height)
 373         {
 374             row_pointer[0] = &image_buffer[(cinfo.image_height -
 375                                             cinfo.next_scanline) *
 376                                             row_stride];
 377             (void) jpeg_write_scanlines(&cinfo, row_pointer, 1);
 378         }/*Loop through file writing one line at a time*/
 379
 380         jpeg_finish_compress(&cinfo);                  /* Finish compression */
 381         fclose(jpeg_file);                                     /* Close file */
 382         jpeg_destroy_compress(&cinfo);                     /* Release memory */
 383
 384         free(data);
 385         free(image_buffer);
 386         printf("INFO   : wrote output to %s\n", jpeg_file_name);
 387     }
 388     exit(EXIT_SUCCESS);
 389 } /*End of program*/