/*************************************************************************** * This file is a part of CADS/UVS fits2jpeg conversion software * * Copyright (C) 2012 by CADS/UV Software Team, * * Indian Institute of Astrophysics * * Bangalore 560034 * * cads_AT_iiap.res.in * * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the * * Free Software Foundation, Inc., * * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * ***************************************************************************/ /*Header Definitions*/ #include "fits2jpeg.h" /*---------------------------------------------------------------------------* * READ_FITS: To reads data from a fits image file.. (isn't that obvious?) *---------------------------------------------------------------------------*/ void read_fits(char * fits_file_name, long * xdim, long * ydim, float ** data) { fitsfile *fptr; int status = 0, nfound, anynull; long naxes[2]; long npixels; float nullval = 0.0; fits_open_file(&fptr, fits_file_name, READONLY, &status); fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxes, &nfound, &status); if (status) printerro(strcat(fits_file_name, " <-- Failed to open the file")); /* Read in data */ npixels = naxes[0] * naxes[1]; (*data) = malloc(sizeof(float) * npixels); nullval = 0; if (fits_read_img(fptr, TFLOAT, 1, npixels, &nullval, (*data), &anynull, &status)) printerro(strcat(fits_file_name, " has no valid fits image data")); *xdim = naxes[0]; *ydim = naxes[1]; fits_close_file(fptr, &status); } /*---------------------------------------------------------------------------* * SCALE_PIXELS: Changes the pixel scale to linear/log/sqroot/etc.. *---------------------------------------------------------------------------*/ void scale_pixels(int scale, unsigned int npixels, float *data, JSAMPLE ** image_buffer) { unsigned int i = 0; int JMAXVAL = 255; float datamax = 0.0, datamin = 0.0, tmp = 0.0; float hist[256] = {0.0}, cumhist[256] = {0.0}; float scl_data = 0.0; /* first find min & max in data */ datamax = -1.0 * FLT_MAX; datamin = FLT_MAX; for (i = 0; i < npixels; ++i) { if (data[i] > datamax) datamax = data[i]; if (data[i] < datamin) datamin = data[i]; } /*endfor*/ /* Convert data into bytscaled values for jpeg file */ /* the dynamic range is reduced to 255 for jpeg */ scl_data = (datamax - datamin)/(float)JMAXVAL; /* we will end up with segfaults if scl_data = 0 */ if (scl_data == 0) scl_data = 1; for (i = 0; i < npixels; ++i) data[i] = (data[i] - datamin)/scl_data; /* All data is now squeezed into the range 0 - 255 */ /* NOTE: At this point onwards min & max is 0 and 255 respectively */ datamax = (float)JMAXVAL; datamin = 0.0; /* initialize image histogram. ensure all are zeroes in hist[] */ /*-----------------------------------------------------------------------*/ for (i = 0; i <= JMAXVAL; ++i) hist[i] = 0; /* construct the image histogram */ tmp = 1.0/(float)npixels; for (i = 0; i <= npixels; ++i) hist[(int)floor(data[i])] += tmp; /* And the cumulative histogram */ cumhist[0] = hist[0]; for (i = 1; i <= JMAXVAL; ++i) cumhist[i] += cumhist[i - 1] + hist[i]; /* Allocate image buffer */ (*image_buffer) = malloc(sizeof(unsigned char) * npixels); /* Linear scale (min-max) : This is the default scaling * histo-eq will fail if we dont generate image_buffer here */ for (i = 0; i < npixels; ++i) (*image_buffer)[i] = (int)(data[i]); /*-----------------------------------------------------------------------*/ switch (scale) { case 1 : /* Square root */ scl_data = sqrt((float)JMAXVAL)/(float)JMAXVAL; for (i = 0; i < npixels; ++i) (*image_buffer)[i] = (int)(sqrt(data[i])/scl_data); break; case 2 : /* Square */ scl_data = pow((float)JMAXVAL,2)/(float)JMAXVAL; for (i = 0; i < npixels; ++i) (*image_buffer)[i] = (int)abs((pow(data[i],2) - 1.0)/scl_data); break; case 3 : /* Cubic */ scl_data = pow((float)JMAXVAL,3)/(float)JMAXVAL; for (i = 0; i < npixels; ++i) (*image_buffer)[i] = (int)abs((pow(data[i],3) - 1.0)/scl_data); break; case 4 : /* log */ scl_data = log(1.0 + (float)JMAXVAL)/(float)JMAXVAL; for (i = 0; i < npixels; ++i) (*image_buffer)[i] = (int)((log(abs(data[i]) + 1.0))/scl_data); break; case 5 : /* contrast stretch */ /* We need to go through the cumulative histogram to pick the * appropriate values for datamin and datamax */ i = 0; while (i < JMAXVAL) { if (cumhist[i] >= 0.01) { datamin = (float) i; break; } i++; } i = JMAXVAL; while (i > 0) { if (cumhist[i] <= 0.99) { datamax = (float) i; break; } i--; } scl_data = (datamax - datamin)/(float)JMAXVAL; for (i = 0; i < npixels; ++i) { if ((*image_buffer)[i] >= datamax) (*image_buffer)[i] = JMAXVAL; else if ((*image_buffer)[i] <= datamin) (*image_buffer)[i] = 0; else (*image_buffer)[i] = (int) abs(((*image_buffer)[i] - datamin)/scl_data); } break; case 6 : /* histogram equalization */ for (i = 0; i < npixels; ++i) (*image_buffer)[i] = cumhist[(*image_buffer)[i]] * JMAXVAL; break; default : break; } } /*---------------------------------------------------------------------------* * RESIZE_IMAGE: Scales down/up the image_buffer using bilinear scaling * Based on an article by "John" at * http://tech-algorithm.com/articles/bilinear-image-scaling/ *---------------------------------------------------------------------------*/ void resize_image(long *xdim, long *ydim, float zoomfact, JSAMPLE ** image_buffer) { int offset = 0, index = 0; int A, B, C, D, x, y, gray; JSAMPLE *buff; unsigned int i = 0, j = 0; unsigned long npixels = 0; long w = *xdim, h = *ydim; long zxdim = 0, zydim = 0; float xdiff, ydiff, xratio, yratio; zxdim = (int)(w * zoomfact); zydim = (int)(h * zoomfact); npixels= zxdim * zydim; xratio = ((float)(w - 1))/zxdim; yratio = ((float)(h - 1))/zydim; /* allocate space for *buff */ buff = malloc(sizeof(unsigned char) * zxdim * zydim); index = 0; offset = 0; for (i = 0; i < zydim; i++) { y = (int)(yratio * i); ydiff = (yratio * i) - y; for (j = 0; j < zxdim; j++) { x = (int)(xratio * j); xdiff = (xratio * j) - x; index = y * w + x; A = (*image_buffer)[index] & 0xff; B = (*image_buffer)[index + 1] & 0xff; C = (*image_buffer)[index + w] & 0xff; D = (*image_buffer)[index + w + 1] & 0xff; gray = (int)(A * (1 - xdiff) * (1 - ydiff) + B * (xdiff) * (1 - ydiff) + C * (ydiff) * (1 - xdiff) + D * (xdiff) * (ydiff) ); buff[offset++] = gray; } } *xdim = zxdim; *ydim = zydim; (*image_buffer) = realloc((*image_buffer), sizeof(unsigned char) * npixels); for (i = 0; i < npixels; ++i) (*image_buffer)[i] = buff[i]; free(buff); }