From b13ab2fa9923c6d8dce84e4533105a2aeae1c838 Mon Sep 17 00:00:00 2001 From: Ricardo Zanmar Date: Sun, 5 Jun 2011 22:00:25 +0200 Subject: [PATCH] histogram based min max values. environment vals The y axis was inverted becase the way the fits files are saved. This has been fixed. I now use DATAMIN DATAMAX environmental variables to pass values to the loader. At the moment this only works if both are defined but it might be useful to split this. If these variables don't exist the we do histogram and base the min and max on this, very similar to ds9's 99.5 scaling. --- Makefile | 6 ++- loader_fits.c | 144 ++++++++++++++++++++++++++++++++++++++++---------- 2 files changed, 121 insertions(+), 29 deletions(-) diff --git a/Makefile b/Makefile index 705fb85..2979010 100644 --- a/Makefile +++ b/Makefile @@ -1,10 +1,12 @@ +IMLIB_PATH=../imlib2 + all: fits.so fits.so: loader_fits.o gcc -shared $< -lcfitsio -lSM -lICE -Wl,-soname -Wl,fits.so -o $@ -loader_fits.o: loader_fits.c loader_common.h config.h - gcc -Wall -g -DHAVE_CONFIG_H -c -fPIC -DPIC -o $@ -I. -I../imlib2 -I../imlib2/src -I../imlib2/src/lib $< +loader_fits.o: loader_fits.c config.h + gcc -Wall -g -DHAVE_CONFIG_H -c -fPIC -DPIC -o $@ -I. -I$(IMLIB_PATH) -I$(IMLIB_PATH)/src -I$(IMLIB_PATH)/src/lib -I$(IMLIB_PATH)/src/modules/loaders $< config.h: touch config.h diff --git a/loader_fits.c b/loader_fits.c index 3641a3d..be76102 100644 --- a/loader_fits.c +++ b/loader_fits.c @@ -1,5 +1,11 @@ -#include "loader_common.h" -#include "fitsio.h" +#include +#include + +//TODO: add a color table. +//invert gray scale, Check if feh or qiv can do this. +// add option to check for datamin and datamax on the header +// make it more robust. datamax > datamin +// write the fits writer! #define SWAP32(x) (x) = \ ((((x) & 0x000000ff ) << 24) |\ @@ -11,17 +17,19 @@ char load(ImlibImage * im, ImlibProgressFunction progress, char progress_granularity, char immediate_load) { + fprintf(stderr, "fits load function\n"); + int w = 0, h = 0, alpha = 0; //FILE *f; fitsfile *f; - int status; + int status = 0; if (im->data) return 0; //f = fopen(im->real_file, "rb"); - if ( fits_open_file( &f, im->real_file, READONLY, &status) ) + if ( fits_open_file( &f, im->real_file, READONLY, &status) ) { return 0; - + } //if (!f) // return 0; @@ -78,17 +86,18 @@ load(ImlibImage * im, ImlibProgressFunction progress, DATA8 gray; int y, x, pl = 0; char pper = 0; - float *buffer, datamin, datamax, nullval = 0; + double *buffer, datamin, datamax, nullval = 0; int anynull, nbuffer; long fpixel = 1; - - //gray = ceil(floatvalue * 256) - //*ptr = ( ( ( ( (0xff << 8) | gray) << 8) | gray) << 8) | gray + long *bin, nmax = 0, ndataok, sumpix = 0; + int nbin, k, kmin, mostfreq = 0; + double binw; + char *envdmin, *envdmax; /* must set the im->data member before callign progress function */ ptr = im->data = malloc(w * h * sizeof(DATA32)); nbuffer = w; - buffer = malloc( nbuffer * sizeof( float ) ); + buffer = malloc( nbuffer * sizeof( double ) ); if (!im->data) { @@ -96,30 +105,110 @@ load(ImlibImage * im, ImlibProgressFunction progress, //fclose(f); return 0; } - datamin = 1E8; - datamax = -1E8; + + envdmin = getenv( "DATAMIN" ); + envdmax = getenv( "DATAMAX" ); + if( envdmin != NULL && envdmax != NULL ){ + datamin = atof( envdmin ); + datamax = atof( envdmax ); + fprintf(stderr, "ENV:data min %g, max %g\n", datamin, datamax);} + else{ + + datamin = HUGE_VALF; + datamax = -HUGE_VALF; + + //Find min max + for (y = 0; y < h; y++) + { + fits_read_img( f, TDOUBLE, fpixel, nbuffer, &nullval, + buffer, &anynull, &status); + for( x = 0; x < w; x++ ){ + if( buffer[ x ] > datamax ) + datamax = buffer[ x ]; + if( buffer[ x ] < datamin ) + datamin = buffer[ x ]; + } + fpixel += nbuffer; + } + fprintf(stderr, "data min %g, max %g\n", datamin, datamax); + + nbin = 4 * sqrt( datamax - datamin ); + if( nbin > 5000 ) nbin = 1000; + binw = ( datamax - datamin ) / nbin; + fprintf(stderr, "nbin and bin width:%d %g\n", nbin, binw); + bin = calloc( nbin, sizeof( long ) ); + + //Find histogram + fpixel = 1; + for (y = 0; y < h; y++) + { + fits_read_img( f, TDOUBLE, fpixel, nbuffer, &nullval, + buffer, &anynull, &status); + for( x = 0; x < w; x++ ){ + for( k = 0; k < nbin; k++ ){ + if( buffer[ x ] >= k * binw && + buffer[ x ] < ( k + 1 ) * binw ){ + bin[ k ] += 1; + break; + } + } + } + fpixel += nbuffer; + } + //Find most frequent pixel value + for( k = 0; k < nbin; k++ ){ + if( bin[ k ] > nmax ){ + nmax = bin[ k ]; + mostfreq = k; + } + fprintf(stderr, "%ld ", bin[ k ] ); + } + fprintf(stderr, "mostfreq and nmax = %d %ld ", mostfreq, nmax ); + + //Set new DATAMIN. look for FWHM on a log scale + nmax = log( nmax ) / 2; + datamin = ( mostfreq + 0.5 ) * binw; + for( k = mostfreq - 1; k >= 0; k-- ){ + if( log( bin[ k ] ) < nmax ){ + kmin = k + 1; + datamin = ( (kmin) + 0.5 ) * binw; + break; + } + } + //Set new DATAMAX when we have 99.5% of the pixels accounted for + ndataok = h * w; + for( k = 0; k < kmin; k++) + ndataok -= bin[ k ]; + //ndataok = h * w; //To get values closer to ds9's uncomment this + for( k = kmin; k < nbin; k++){ + sumpix += bin[ k ]; + if( 1.0 * sumpix / ndataok > 0.997 ) + break; + } + datamax = ( ( k + 1 ) + 0.5 ) * binw; + } + + // use the new datamin datamax values to rescale and create imlib im + ptr += w * h - 1; + fpixel = 1; for (y = 0; y < h; y++) { - //fread(ptr, im->w, 4, f); - fits_read_img( f, TFLOAT, fpixel, nbuffer, &nullval, + fits_read_img( f, TDOUBLE, fpixel, nbuffer, &nullval, buffer, &anynull, &status); - for( x = 0; x < w; x++ ){ - if( buffer[ x ] > datamax ) - datamax = buffer[ x ]; - if( buffer[ x ] < datamin ) - datamin = buffer[ x ]; - } - for( x = 0; x < w; x++ ){ - gray = ceil( ( buffer[ x ] - datamin ) / ( datamax - datamin ) * 256 ); - *ptr = ( ( ( ( (0xff << 8) | gray) << 8) | gray) << 8) | gray ; - ptr++; + for( x = w - 1; x >= 0; x-- ){ + if( buffer[ x ] <= datamin ) + gray = 0; + else + if( buffer[ x ] >= datamax ) + gray = 255; + else + gray = ceil(( buffer[ x ] - datamin ) / ( datamax - datamin ) * 255 ); + *ptr = (((( (0xff << 8) | gray) << 8) | gray) << 8) | gray ; + ptr--; } - fpixel += nbuffer; - //*ptr = y; - //ptr += im->w; if (progress) { char per; @@ -141,6 +230,7 @@ load(ImlibImage * im, ImlibProgressFunction progress, } } } + fprintf(stderr, "data min %g, max %g\n", datamin, datamax); } //fclose(f); fits_close_file(f, &status); -- 2.32.0.93.g670b81a890