58800fbbd6b5b946f88862be75e6356153d77cff
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
26   Libraries needed are CFITSIO, JPEG
27   Original author: Jayant Murthy
28   Last Modification: 11/05/02
30   Reconstruction and additional features: Reks
31  *---------------------------------------------------------------------------*/
32 /*---------------------------------------------------------------------------
33   Modification History
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..
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
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.
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).
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 :)
69   Reks 28 June 2012
70   1. rebranded to CADS/UVS
72   Reks 05 July 2012 v2.0 plans
73   1. Merged scale & operations to single flag
74   TODO
75       1. New option to specify min-max for image
76       2. specify output directory
77       3. output filename for single input file
78       4. Image resizing options
80  *-------------------------------------------------------------------------*/
82 #ifdef HAVE_CONFIG_H
83 #include <config.h>
84 #endif
86 /*Header Definitions*/
87 #include <stdio.h>
88 #include <stdlib.h>
89 #include <string.h>
90 #include <math.h>
91 #include <libgen.h>
92 #include <dirent.h>
93 #include <fitsio2.h>
94 #include <jpeglib.h>
96 #define PROGRAM "fits2jpeg"
97 #define r2d (90./atan2(1,0))
98 #define MAX_TEXT 150
99 int my_getopt(int, char * const *, const char *);
100 void signals_handler(int);
101 void set_signals(void);
103 char *optarg;
104 int optindex;
105 /*---------------------------------------------------------------------------*/
106 void PRINT_BANNER()
108     printf("\n");
109     printf(" %s %s\n", PROGRAM, VERSION);
110     printf(" Converts fits image files to jpeg\n");
111     printf("-----------------------------------------------------------\n");
113 /*---------------------------------------------------------------------------*/
114 void usage()
116     printf ("\n Usage: fits2jpeg [-s <options>] <fits_file>             \n");
117     printf ("  Options are:                                             \n");
118     printf ("    -h help                                                \n");
119     printf ("    -s <scale_type>                                        \n");
120     printf ("       scale for output image, where <scale_type> can be:  \n");
121     printf ("         linear         Linear scale, default              \n");
122     printf ("         sqroot         for square root scale              \n");
123     printf ("         square         for quadratic scale                \n");
124     printf ("         cubic          for cubic scale                    \n");
125     printf ("         log            for log scale                      \n");
126     printf ("         normalize      for linear histogram stretch       \n");
127     printf ("         equalize       for histogram equalization         \n");
128     printf ("  Output will be written to <fits_file_root>.jpg. Wild card\n");
129     printf ("  entries allowed in <fits_file>; For eg: *.fits, m31*.fits\n");
130     printf ("  ngc???.fits etc.                                         \n");
131     printf ("-----------------------------------------------------------\n");
132     printf ("Copyright (c) 2012 The CADS Team, IIAp. [GPL v3.0 or later]\n");
133     printf ("Report Bugs to: http://cads.iiap.res.in/bugzilla           \n");
134     printf ("Documentation : http://cads.iiap.res.in/software         \n\n");
135     exit(1);
137 /*---------------------------------------------------------------------------*/
139 /*--------------------------- Begin Main Progam -----------------------------*/
140 int main(int argc, char *argv[], char *envp[])
142     fitsfile *fptr;
144     int scale = 0, process = 0, status = 0, nfound, anynull;
145     unsigned int i = 0, j = 0, npixels = 0;
146     long naxes[2], row_stride;
147     float datamin = 0.0, datamax = 0.0;
148     float nullval = 0.0, scl_data;
149     float *data;
150     float tmp = 0.0;
151     float hist[256] = {0.0}, cumhist[256] = {0.0};
152     char jpeg_file_name[MAX_TEXT] = "", fits_file_name[MAX_TEXT] = "";
154     int opt;                                                   /* For getopt */
155     extern char *optarg;
156     extern int optindex;
158     FILE *jpeg_file;
160     struct jpeg_error_mgr jerr;
161     struct jpeg_compress_struct cinfo;
162     int JMAXVAL = 255;                    /* Max value for greyscale in jpeg */
163     JSAMPROW row_pointer[1];
164     JSAMPLE *image_buffer;
165     /*---------------- End of variable initialization -----------------------*/
167     PRINT_BANNER();        /* Header: Talks about program, version & purpose */
168     set_signals();                                  /* Trap un-natural exits */
170     /*-------------------- Parse command line inputs ------------------------*/
172     if (argc < 2) usage();  /* People who don't give any arguments need help */
173     while ( ( opt = my_getopt ( argc, argv, "s:e:r:h" ) ) != EOF )
174         switch ( opt )
175         {
176         case 's':
177             if      (strcmp(optarg, "sqroot")   ==0) scale = 1;
178             else if (strcmp(optarg, "square")   ==0) scale = 2;
179             else if (strcmp(optarg, "cube")     ==0) scale = 3;
180             else if (strcmp(optarg, "log")      ==0) scale = 4;
181             else if (strcmp(optarg, "normalize")==0) scale = 5;
182             else if (strcmp(optarg, "equalize") ==0) scale = 6;
183             else
184             {
185                 fprintf(stderr, "ERROR  : Unrecognized option : %s\n", optarg);
186                 exit(EXIT_FAILURE);
187             }
188             break;
190         case 'r':
191             /*
192              * read a string like <min>:<max>
193              * set a flag. If this flag is on, do a 'reverse map'
194              * instead of search for datamin & datamax
195              */
196             break;
198         case 'h':
199             usage();
200             break;
202         case '*':
203             usage();
204             break;
205         }
207     if ((argc - optindex) == 0) usage();       /* Somebody gave only options */
209     /*---------------------- Begin processing files -------------------------*/
211     /* Getopt would have processed all argument with a "-" in front. optindex
212        is the counter, which keeps the number of options supplied. Whatever
213        remains in the list of command line arguments are the fits files.
214        NOTE: POSIX systems will do the regex matching and expand wild card
215              entries for us. Even windoze XP Pro is known to do that...      */
216     for (j = optindex; j < argc; j++)
217     {
219         strcpy(fits_file_name, argv[j]);
221         strcpy(jpeg_file_name, fits_file_name);
222         strtok(jpeg_file_name, ".");
223         strcat(jpeg_file_name, ".jpg");
225         fits_open_file(&fptr, fits_file_name, READONLY, &status);
226         fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxes, &nfound, &status);
227         if (status)
228         {
229             fprintf(stderr, "ERROR  : Could not open file: %s\n",
230                     fits_file_name);
231             exit(EXIT_FAILURE);
232         }/*Endif*/
234         /* Read in data */
235         npixels = naxes[0] * naxes[1];
236         data = (float *) malloc(sizeof(float) * npixels);
238         nullval = 0;
239         if (fits_read_img(fptr, TFLOAT, 1, npixels, &nullval, data, &anynull,
240                           &status))
241         {
242             fprintf(stderr, "ERROR  : No valid fits image data in %s\n",
243                     fits_file_name);
244             exit(EXIT_FAILURE);
245         }/*Endif*/
247         fits_close_file(fptr, &status);
248         /*Close data file*/
249         printf("INFO   : data input from %s\n", fits_file_name);
251         /* IF no range is provided, find min & max in image */
252         datamax = -1.0e9;
253         datamin = +1.0e9;
254         for (i = 0; i < npixels; ++i)
255         {
256             if (data[i] > datamax) datamax = data[i];
257             if (data[i] < datamin) datamin = data[i];
258         } /*endfor*/
261         /* Convert data into bytscaled values for jpeg file                 */
262         /* the dynamic range is reduced to 255 for jpeg                     */
263         scl_data = (datamax - datamin)/(float)JMAXVAL;
264         for (i = 0; i < npixels; ++i)
265             data[i] = (data[i] - datamin)/scl_data;
267         /* All data is now squeezed into the range 0 - 255                   */
268         /* NOTE: At this point onwards min & max is 0 and 255 respectively   */
271         /* Now we need to look at user options..                             */
272         /* 1. scale it as per user's requirement.                            */
273         /* 2. image enhancing operations                                     */
274                                                     /* Allocate image buffer */
275         image_buffer = (unsigned char *) malloc(sizeof(char) * npixels);
276         scl_data = 1.0;
279         /* initialize image histogram. ensure all are zeroes in hist[]       */
280         /*-------------------------------------------------------------------*/
281         for (i = 0; i <= JMAXVAL; ++i) hist[i] = 0;
283         /* construct the image histogram */
284         tmp = 1.0/(float)npixels;
285         for (i = 0; i <= npixels; ++i)
286             hist[(int)floor(data[i])] += tmp;
288         /* And the cumulative histogram */
289         cumhist[0] = hist[0];
290         for (i = 1; i <= JMAXVAL; ++i)
291             cumhist[i] += cumhist[i - 1] + hist[i];
292         /*-------------------------------------------------------------------*/
295         switch (scale)
296         {
297         case 1 :                                              /* Square root */
298             printf("INFO   : Using square-root scale\n");
299             scl_data = sqrt((float)JMAXVAL)/(float)JMAXVAL;
300             for (i = 0; i < npixels; ++i)
301                 image_buffer[i] = (int)(sqrt(data[i])/scl_data);
302             break;
304         case 2 :                                                   /* Square */
305             printf("INFO   : Using quadratic scale\n");
306             scl_data = pow((float)JMAXVAL,2)/(float)JMAXVAL;
307             for (i = 0; i < npixels; ++i)
308                 image_buffer[i] = (int)abs((pow(data[i],2) - 1.0)/scl_data);
309             break;
311         case 3 :                                                    /* Cubic */
312             printf("INFO   : Using cubic scale\n");
313             scl_data = pow((float)JMAXVAL,3)/(float)JMAXVAL;
314             for (i = 0; i < npixels; ++i)
315                 image_buffer[i] = (int)abs((pow(data[i],3) - 1.0)/scl_data);
316             break;
318         case 4 :                                                      /* log */
319             printf("INFO   : Using log scale\n");
320             scl_data = log(1.0 + (float)JMAXVAL)/(float)JMAXVAL;
321             for (i = 0; i < npixels; ++i)
322                 image_buffer[i] = (int)((log(abs(data[i]) + 1.0))/scl_data);
323             break;
325         case 5 :
326             /* contrast stretch */
327             printf("INFO   : Performing histogram stretch (normalization)\n");
329             datamax = (float)JMAXVAL;
330             datamin = 0.0;
332             /* We need to go through the cumulative histogram to pick the
333              *  appropriate values for datamin and datamax               */
334             i = 0;
335             while (i < JMAXVAL)
336             {
337                 if (cumhist[i] >= 0.01)
338                 {
339                     datamin = (float) i;
340                     break;
341                 }
342                 i++;
343             }
344             i = JMAXVAL;
345             while (i > 0)
346             {
347                 if (cumhist[i] <= 0.99)
348                 {
349                     datamax = (float) i;
350                     break;
351                 }
352                 i--;
353             }
354             scl_data = (datamax - datamin)/(float)JMAXVAL;
355             for (i = 0; i < npixels; ++i)
356             {
357                 if (image_buffer[i] >= datamax)
358                     image_buffer[i] = JMAXVAL;
359                 else if (image_buffer[i] <= datamin)
360                     image_buffer[i] = 0;
361                 else
362                     image_buffer[i] = (int) abs((image_buffer[i]
363                                     - datamin)/scl_data);
364             }
365             break;
367         case 6 :
368            /* histogram equalization */
369             printf("INFO   : Performing Histogram Equalization\n");
370             for (i = 0; i <  npixels; ++i)
371                 image_buffer[i] = cumhist[image_buffer[i]]*255;
372             break;
374         default:                                   /* Linear scale (min-max) */
375             for (i = 0; i < npixels; ++i)
376                 image_buffer[i] = (int)(data[i]);
377             break;
378         }
381         /*Write out data into JPEG file*/
382         cinfo.err = jpeg_std_error(&jerr);             /* JPEG error handler */
383         jpeg_create_compress(&cinfo);                 /* JPEG initialization */
384         jpeg_file = fopen(jpeg_file_name, "wb");/* Open JPEG file for writing*/
385         jpeg_stdio_dest(&cinfo, jpeg_file); /* Send compressed data to stdio */
386         cinfo.image_width = naxes[0];                         /* Image width */
387         cinfo.image_height = naxes[1];                       /* Image height */
388         cinfo.input_components = 1;            /* Number of colors per pixel */
389         cinfo.in_color_space = JCS_GRAYSCALE;   /* colorspace of input image */
390         jpeg_set_defaults(&cinfo);                      /* Set JPEG defaults */
391         jpeg_start_compress(&cinfo, TRUE);       /* default data compression */
392         row_stride = naxes[0];            /* JSAMPLEs per row inimage buffer */
394         /* Now we have to turn the data upside down                          */
395         while (cinfo.next_scanline < cinfo.image_height)
396         {
397             row_pointer[0] = &image_buffer[(cinfo.image_height -
398                                             cinfo.next_scanline) *
399                                             row_stride];
400             (void) jpeg_write_scanlines(&cinfo, row_pointer, 1);
401         }/*Loop through file writing one line at a time*/
403         jpeg_finish_compress(&cinfo);                  /* Finish compression */
404         fclose(jpeg_file);                                     /* Close file */
405         jpeg_destroy_compress(&cinfo);                     /* Release memory */
407         free(data);
408         free(image_buffer);
409         printf("INFO   : wrote output to %s\n", jpeg_file_name);
410     }
411     exit(EXIT_SUCCESS);
412 } /*End of program*/