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 ***************************************************************************/
24 /*Header Definitions*/
25 #include "fits2jpeg.h"
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;
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"));
43 /* Read in data */
44 npixels = naxes[0] * naxes[1];
45 (*data) = malloc(sizeof(float) * npixels);
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"));
52 *xdim = naxes[0];
53 *ydim = naxes[1];
55 fits_close_file(fptr, &status);
56 }
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;
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*/
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;
85 /* we will end up with segfaults if scl_data = 0 */
86 if (scl_data == 0) scl_data = 1;
88 for (i = 0; i < npixels; ++i)
89 data[i] = (data[i] - datamin)/scl_data;
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;
97 /* initialize image histogram. ensure all are zeroes in hist[] */
98 /*-----------------------------------------------------------------------*/
99 for (i = 0; i <= JMAXVAL; ++i) hist[i] = 0;
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;
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];
111 /* Allocate image buffer */
112 (*image_buffer) = malloc(sizeof(unsigned char) * npixels);
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]);
120 /*-----------------------------------------------------------------------*/
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;
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;
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;
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;
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;
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 }
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;
211 zxdim = (int)(w * zoomfact);
212 zydim = (int)(h * zoomfact);
214 npixels= zxdim * zydim;
216 xratio = ((float)(w - 1))/zxdim;
217 yratio = ((float)(h - 1))/zydim;
219 /* allocate space for *buff */
220 buff = malloc(sizeof(unsigned char) * zxdim * zydim);
222 index = 0;
223 offset = 0;
224 for (i = 0; i < zydim; i++)
225 {
226 y = (int)(yratio * i);
227 ydiff = (yratio * i) - y;
229 for (j = 0; j < zxdim; j++)
230 {
231 x = (int)(xratio * j);
233 xdiff = (xratio * j) - x;
234 index = y * w + x;
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;
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 }