1 #include <loader_common.h>
6 // - read the image only once and recycle the ptr array if the image is not
8 // - biseearch for the lower limit. Would converge faster in flat images where the max is towards the end of the histogram. this images take a second more.
9 // add option to check for datamin and datamax on the header
10 // make it more robust. datamax > datamin
11 // write the fits writer!
13 #define SWAP32(x) (x) = \
14 ((((x) & 0x000000ff ) << 24) |\
15 (((x) & 0x0000ff00 ) << 8) |\
16 (((x) & 0x00ff0000 ) >> 8) |\
17 (((x) & 0xff000000 ) >> 24))
19 /* return the difference in time, in milliseconds */
20 double timediff_ms(const struct timeval *start, const struct timeval *stop)
22 return (stop->tv_sec - start->tv_sec)*1000.0 + (stop->tv_usec - start->tv_usec)/1000.0;
25 void getRGB( double, double, double, short, short, DATA8 *, DATA8 *, DATA8 * );
28 load(ImlibImage * im, ImlibProgressFunction progress,
29 char progress_granularity, char immediate_load)
32 int w = 0, h = 0, alpha = 0;
39 //f = fopen(im->real_file, "rb");
40 if ( fits_open_file( &f, im->real_file, READONLY, &status) ) {
51 //char buf[256], buf2[256];
53 /* read the NAXIS1 and NAXIS2 keyword to get image size */
54 if ( fits_read_keys_lng( f, "NAXIS", 1, 2, naxes, &nfound, &status ) ){
55 fits_close_file(f, &status);
63 if (!fgets(buf, 255, f))
69 sscanf(buf, "%s %i %i %i", buf2, &w, &h, &alpha);
70 if (strcmp(buf2, "ARGB"))
76 if (!IMAGE_DIMENSIONS_OK(w, h))
78 fits_close_file(f, &status);
87 SET_FLAG(im->flags, F_HAS_ALPHA);
89 UNSET_FLAG(im->flags, F_HAS_ALPHA);
90 im->format = strdup("fits");
93 if (((!im->data) && (im->loader)) || (immediate_load) || (progress))
99 double *buffer, datamin, datamax, nullval = 0, *subim;
100 double zscalemin, zscalemax;
101 int anynull, nbuffer;
103 long *bin, nmax = 0, ndataok, sumpix = 0;
104 int nbin, k, kmin, mostfreq = 0, nsub;
106 char *envdmin, *envdmax;
107 short color, inv, stepx, stepy;
109 struct timeval start, end;
112 gettimeofday(&start, NULL);
114 fprintf(stderr, "%s\n", im->real_file );
115 /* must set the im->data member before callign progress function */
116 ptr = im->data = malloc(w * h * sizeof(DATA32));
118 buffer = malloc( nbuffer * sizeof( double ) );
122 fits_close_file(f, &status);
127 envdmin = getenv( "DATAMIN" );
128 envdmax = getenv( "DATAMAX" );
129 if( envdmin != NULL && envdmax != NULL ){
130 datamin = atof( envdmin );
131 datamax = atof( envdmax );
132 fprintf(stderr, "ENV:data min %g, max %g\n", datamin, datamax);}
136 datamax = -HUGE_VALF;
140 if( stepx == 0 ) stepx = 1;
141 if( stepy == 0 ) stepy = 1;
143 fprintf(stderr, "subsample image: %d %d\n", stepx, stepy );
145 //Find min max and save a subsample image
148 subim = malloc( (w / stepx + 1)*(h / stepy +1)*sizeof( double ) );
149 fprintf(stderr, "image: %d %d %d \n", (w / stepx + 1) * (h / stepy + 1), w/stepx+1,h/stepy+1);
150 for (y = 0; y < h; y += stepy )
152 fits_read_img( f, TDOUBLE, fpixel, nbuffer, &nullval,
153 buffer, &anynull, &status);
154 for( x = 0; x < w; x += stepx ){
155 subim[ nsub ] = buffer[ x ];
157 if( buffer[ x ] > datamax )
158 datamax = buffer[ x ];
159 if( buffer[ x ] < datamin )
160 datamin = buffer[ x ];
162 fpixel += nbuffer * stepy;
164 fprintf(stderr, "data min %g, max %g\n", datamin, datamax);
165 fprintf(stderr, "subimage has %d pixels out of %d... %d %d\n", nsub, w*h, x/stepx, y/stepy);
167 nbin = 4 * sqrt( datamax - datamin );
168 if( nbin > 1000 ) nbin = 1000;
169 binw = ( datamax - datamin ) / nbin;
170 fprintf(stderr, "nbin and bin width:%d %g\n", nbin, binw);
171 bin = calloc( nbin, sizeof( long ) );
174 for( y = 0; y < nsub; y++ ){
175 for( k = 0; k < nbin; k++ ){
176 // if( y==nsub/2 )fprintf(stderr,"%g %g %g\n",k*binw + datamin, ( k+1 )*binw + datamin, subim[ y ]);
177 if( subim[ y ] >= k*binw + datamin && subim[ y ] < ( k+1 )*binw + datamin ){
184 //Find most frequent pixel value
185 for( k = 0; k < nbin; k++ ){
186 if( bin[ k ] > nmax ){
190 fprintf(stderr, "%ld ", bin[ k ] );
192 fprintf(stderr, "\nmostfreq and nmax = %d %ld ", mostfreq, nmax );
194 //Set new DATAMIN. look for FWHM on a log scale
195 nmax = log( nmax ) / 2;
196 //datamin = ( mostfreq + 0.5 ) * binw;
197 zscalemin = ( mostfreq + 0.5 ) * binw + datamin;
198 for( k = mostfreq - 1; k >= 0; k-- ){
200 if( log( bin[ k ] ) < nmax ){
202 zscalemin = ( (kmin) + 0.5 ) * binw + datamin;
206 fprintf(stderr, "\n after fwfm log datamin, kmin: %g %d %d ", zscalemin,kmin, k );
207 if( k <= 0 ){ //case were we have one single peak surrounded by 0
209 zscalemin = ( (kmin) + 0.5 ) * binw + datamin;
212 fprintf(stderr, "\n after fwfm log datamin, kmin: %g %d ", zscalemin,kmin );
214 //Set new DATAMAX when we have 99.5% of the pixels accounted for
217 for( k = 0; k < kmin; k++)
219 //ndataok = h * w; //To get values closer to ds9's uncomment this
220 ndataok = 0.995 * ndataok;
221 for( k = kmin; k < nbin; k++){
223 if( sumpix > ndataok )
226 zscalemax = ( ( k + 1 ) + 0.5 ) * binw + datamin;
231 // use the new zscale min/max values to rescale and create imlib im
236 for (y = 0; y < h; y++)
238 fits_read_img( f, TDOUBLE, fpixel, nbuffer, &nullval,
239 buffer, &anynull, &status);
241 for( x = w - 1; x >= 0; x-- ){
242 getRGB( buffer[ x ], zscalemin, zscalemax, color, inv, &r, &g, &b);
243 if( x == w / 100 && y == h / 100 ) fprintf( stderr, "\n%d %d %d %g\n",r,g,b, buffer[ x ]);
244 *ptr = (((( (0xff << 8) | r) << 8) | g) << 8) | b ;
254 per = (char)((100 * y) / im->h);
255 if (((per - pper) >= progress_granularity) ||
259 if (!progress(im, per, 0, (y - l), im->w, l))
261 fits_close_file(f, &status);
270 fprintf(stderr, "\ndata min %g, max %g\n", zscalemin, zscalemax);
274 gettimeofday(&end, NULL);
275 delta_t = timediff_ms( &start, &end );
276 fprintf(stderr,"\n --> Used %.2f msec.\n", delta_t);
279 fits_close_file(f, &status);
284 save(ImlibImage * im, ImlibProgressFunction progress, char progress_granularity)
288 int y, pl = 0, alpha = 0;
292 /* no image data? abort */
295 f = fopen(im->real_file, "wb");
298 if (im->flags & F_HAS_ALPHA)
300 fprintf(f, "ARGB %i %i %i\n", im->w, im->h, alpha);
302 for (y = 0; y < im->h; y++)
304 fwrite(ptr, im->w, 4, f);
311 per = (char)((100 * y) / im->h);
312 if (((per - pper) >= progress_granularity) || (y == (im->h - 1)))
315 if (!progress(im, per, 0, (y - l), im->w, l))
331 formats(ImlibLoader * l)
333 char *list_formats[] = { "fits", "fts", "fit" };
338 l->num_formats = (sizeof(list_formats) / sizeof(char *));
339 l->formats = malloc(sizeof(char *) * l->num_formats);
340 for (i = 0; i < l->num_formats; i++)
341 l->formats[i] = strdup(list_formats[i]);
345 void getRGB( double val, double min, double max, short color, short inv,
346 DATA8 *r, DATA8 *g, DATA8 *b ){
349 const float slsR[200]={0, 0.043442, 0.086883, 0.130325, 0.173767, 0.217208,
350 0.260650, 0.304092, 0.347533, 0.390975, 0.434417, 0.477858, 0.521300,
351 0.506742, 0.492183, 0.477625, 0.463067, 0.448508, 0.433950, 0.419392,
352 0.404833, 0.390275, 0.375717, 0.361158, 0.346600, 0.317717, 0.288833,
353 0.259950, 0.231067, 0.202183, 0.173300, 0.144417, 0.115533, 0.086650,
354 0.057767, 0.028883, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
355 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
356 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.033167,
357 0.066333, 0.099500, 0.132667, 0.165833, 0.199000, 0.232167, 0.265333,
358 0.298500, 0.331667, 0.364833, 0.398000, 0.430950, 0.463900, 0.496850,
359 0.529800, 0.562750, 0.595700, 0.628650, 0.661600, 0.694550, 0.727500,
360 0.760450, 0.793400, 0.810617, 0.827833, 0.845050, 0.862267, 0.879483,
361 0.896700, 0.913917, 0.931133, 0.948350, 0.965567, 0.982783, 1, 0.995725,
362 0.991450, 0.987175, 0.982900, 0.978625, 0.974350, 0.970075, 0.965800,
363 0.961525, 0.957250, 0.952975, 0.948700, 0.952975, 0.957250, 0.961525,
364 0.965800, 0.970075, 0.974350, 0.978625, 0.982900, 0.987175, 0.991450,
365 0.995725, 1, 0.998342, 0.996683, 0.995025, 0.993367, 0.991708, 0.990050,
366 0.988392, 0.986733, 0.985075, 0.983417, 0.981758, 0.980100, 0.955925,
367 0.931750, 0.907575, 0.883400, 0.859225, 0.835050, 0.810875, 0.786700,
368 0.762525, 0.738350, 0.714175, 0.690000, 0.715833, 0.741667, 0.767500,
369 0.793333, 0.819167, 0.845000, 0.870833, 0.896667, 0.922500, 0.948333,
370 0.974167, 1, 1, 1, 1, 1, 1, 1, 1};
372 const float slsG[200]={0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
373 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.019817,
374 0.039633, 0.059450, 0.079267, 0.099083, 0.118900, 0.138717, 0.158533,
375 0.178350, 0.198167, 0.217983, 0.237800, 0.268533, 0.299267, 0.330000,
376 0.360733, 0.391467, 0.422200, 0.452933, 0.483667, 0.514400, 0.545133,
377 0.575867, 0.606600, 0.631733, 0.656867, 0.682000, 0.707133, 0.732267,
378 0.757400, 0.782533, 0.807667, 0.832800, 0.857933, 0.883067, 0.908200,
379 0.901908, 0.895617, 0.889325, 0.883033, 0.876742, 0.870450, 0.864158,
380 0.857867, 0.851575, 0.845283, 0.838992, 0.832700, 0.832308, 0.831917,
381 0.831525, 0.831133, 0.830742, 0.830350, 0.829958, 0.829567, 0.829175,
382 0.828783, 0.828392, 0.828000, 0.834167, 0.840333, 0.846500, 0.852667,
383 0.858833, 0.865000, 0.871167, 0.877333, 0.883500, 0.889667, 0.895833,
384 0.902000, 0.902000, 0.902000, 0.902000, 0.902000, 0.902000, 0.902000,
385 0.902000, 0.902000, 0.902000, 0.902000, 0.902000, 0.902000, 0.897133,
386 0.892267, 0.887400, 0.882533, 0.877667, 0.872800, 0.867933, 0.863067,
387 0.858200, 0.853333, 0.848467, 0.843600, 0.824892, 0.806183, 0.787475,
388 0.768767, 0.750058, 0.731350, 0.712642, 0.693933, 0.675225, 0.656517,
389 0.637808, 0.619100, 0.600408, 0.581717, 0.563025, 0.544333, 0.525642,
390 0.506950, 0.488258, 0.469567, 0.450875, 0.432183, 0.413492, 0.394800,
391 0.361900, 0.329000, 0.296100, 0.263200, 0.230300, 0.197400, 0.164500,
392 0.131600, 0.098700, 0.065800, 0.032900, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
393 0, 0, 0, 0.083333, 0.166667, 0.250000, 0.333333, 0.416667, 0.500000,
394 0.583333, 0.666667, 0.750000, 0.833333, 0.916667, 1, 1, 1, 1, 1, 1, 1,
397 const float slsB[200]={0, 0.052883, 0.105767, 0.158650, 0.211533, 0.264417,
398 0.317300, 0.370183, 0.423067, 0.475950, 0.528833, 0.581717, 0.634600,
399 0.640217, 0.645833, 0.651450, 0.657067, 0.662683, 0.668300, 0.673917,
400 0.679533, 0.685150, 0.690767, 0.696383, 0.702000, 0.712192, 0.722383,
401 0.732575, 0.742767, 0.752958, 0.763150, 0.773342, 0.783533, 0.793725,
402 0.803917, 0.814108, 0.824300, 0.838942, 0.853583, 0.868225, 0.882867,
403 0.897508, 0.912150, 0.926792, 0.941433, 0.956075, 0.970717, 0.985358,
404 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.975300, 0.950600, 0.925900,
405 0.901200, 0.876500, 0.851800, 0.827100, 0.802400, 0.777700, 0.753000,
406 0.728300, 0.703600, 0.676675, 0.649750, 0.622825, 0.595900, 0.568975,
407 0.542050, 0.515125, 0.488200, 0.461275, 0.434350, 0.407425, 0.380500,
408 0.354858, 0.329217, 0.303575, 0.277933, 0.252292, 0.226650, 0.201008,
409 0.175367, 0.149725, 0.124083, 0.098442, 0.072800, 0.066733, 0.060667,
410 0.054600, 0.048533, 0.042467, 0.036400, 0.030333, 0.024267, 0.018200,
411 0.012133, 0.006067, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.003983,
412 0.007967, 0.011950, 0.015933, 0.019917, 0.023900, 0.027883, 0.031867,
413 0.035850, 0.039833, 0.043817, 0.047800, 0.051600, 0.055400, 0.059200,
414 0.063000, 0.066800, 0.070600, 0.074400, 0.078200, 0.082000, 0.085800,
415 0.089600, 0.093400, 0.085617, 0.077833, 0.070050, 0.062267, 0.054483,
416 0.046700, 0.038917, 0.031133, 0.023350, 0.015567, 0.007783, 0, 0, 0,
417 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
418 0.083333, 0.166667, 0.250000, 0.333333, 0.416667, 0.500000, 0.583333,
419 0.666667, 0.750000, 0.833333, 0.916667, 1, 1, 1, 1, 1, 1, 1, 1};
427 fac = ( val - min ) / ( max - min );
429 if( inv ) fac = 1 - fac;
431 if( color == 0 ){ // gray
432 *r = ceil ( fac * 255 );
437 if( color == 1 ){ // color SLS
438 ind = ceil( fac * 199 ); // color table has 200 elements
439 *r = ceil( slsR[ ind ] * 255 );
440 *g = ceil( slsG[ ind ] * 255 );
441 *b = ceil( slsB[ ind ] * 255 );