From b9b24b1dd9abdf9801f186f1136739858539f779 Mon Sep 17 00:00:00 2001
From: Ricardo Zanmar <rzanmar@gmail.com>
Date: Mon, 27 Jun 2011 16:37:24 +0200
Subject: [PATCH] Fixed bug on the histogram calculation.

I was forgetting to add datamin to the datamin and datamax values drawn
from the histogram. This was evident on images with a very small range
of values (bias images).
---
 loader_fits.c | 46 +++++++++++++++-------------------------------
 1 file changed, 15 insertions(+), 31 deletions(-)

diff --git a/loader_fits.c b/loader_fits.c
index d090c56..ce925f5 100644
--- a/loader_fits.c
+++ b/loader_fits.c
@@ -88,6 +88,7 @@ load(ImlibImage * im, ImlibProgressFunction progress,
         int                y, x, pl = 0;
         char               pper = 0;
         double             *buffer, datamin, datamax, nullval = 0, *subim;
+        double             zscalemin, zscalemax;
         int                anynull, nbuffer;
         long               fpixel = 1;
         long               *bin, nmax = 0, ndataok, sumpix = 0;
@@ -150,7 +151,7 @@ load(ImlibImage * im, ImlibProgressFunction progress,
             fprintf(stderr, "subimage has %d pixels out of %d... %d %d\n", nsub, w*h, x/stepx, y/stepy);
 
             nbin = 4 * sqrt( datamax - datamin );
-            if( nbin > 5000 ) nbin = 1000;
+            if( nbin > 1000 ) nbin = 1000;
             binw = ( datamax - datamin ) / nbin;
             fprintf(stderr, "nbin and bin width:%d %g\n", nbin, binw);
             bin = calloc( nbin, sizeof( long ) );
@@ -158,31 +159,13 @@ load(ImlibImage * im, ImlibProgressFunction progress,
             //Find histogram
             for( y = 0; y < nsub; y++ ){
                  for( k = 0; k < nbin; k++ ){
-                     if( subim[ y ] >= k*binw && subim[ y ] < ( k+1 )*binw ){
+            //         if( y==nsub/2 )fprintf(stderr,"%g %g %g\n",k*binw + datamin, ( k+1 )*binw + datamin, subim[ y ]);
+                     if( subim[ y ] >= k*binw + datamin && subim[ y ] < ( k+1 )*binw + datamin ){
                          bin[ k ] += 1;
                          break;
                      }
                  }
             }
-            //free( subim );
-/*
-            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++ ){
@@ -196,22 +179,23 @@ load(ImlibImage * im, ImlibProgressFunction progress,
 
             //Set new DATAMIN. look for FWHM on a log scale 
             nmax = log( nmax ) / 2;
-            datamin = ( mostfreq + 0.5 ) * binw;
+            //datamin = ( mostfreq + 0.5 ) * binw;
+            zscalemin = ( mostfreq + 0.5 ) * binw + datamin;
             for( k = mostfreq - 1; k >= 0; k-- ){
                 if( bin[ k ] > 0 )
                     if( log( bin[ k ] ) < nmax ){
                         kmin = k + 1;
-                        datamin = ( (kmin) + 0.5 ) * binw;
+                        zscalemin = ( (kmin) + 0.5 ) * binw + datamin;
                         break;
                     }
             }
-            fprintf(stderr, "\n after fwfm log datamin, kmin: %g %d %d ", datamin,kmin, k );
+            fprintf(stderr, "\n after fwfm log datamin, kmin: %g %d %d ", zscalemin,kmin, k );
             if( k <= 0 ){   //case were we have one single peak surrounded by 0
                kmin = mostfreq; 
-               datamin = ( (kmin) + 0.5 ) * binw;
+               zscalemin = ( (kmin) + 0.5 ) * binw + datamin;
             }
 
-            fprintf(stderr, "\n after fwfm log datamin, kmin: %g %d ", datamin,kmin );
+            fprintf(stderr, "\n after fwfm log datamin, kmin: %g %d ", zscalemin,kmin );
 
             //Set new DATAMAX when we have 99.5% of the pixels accounted for
             //ndataok = h * w;
@@ -219,18 +203,18 @@ load(ImlibImage * im, ImlibProgressFunction progress,
             for( k = 0; k < kmin; k++)
                 ndataok -= bin[ k ];
             //ndataok = h * w; //To get values closer to ds9's uncomment this
-            ndataok = 0.997 * ndataok;
+            ndataok = 0.995 * ndataok;
             for( k = kmin; k < nbin; k++){
                 sumpix += bin[ k ];
                 if( sumpix > ndataok )
                     break;
             }
-            datamax = ( ( k + 1 ) + 0.5 ) * binw;
+            zscalemax = ( ( k + 1 ) + 0.5 ) * binw + zscalemin;
         }
 
         //free( bin );
 
-        // use the new datamin datamax values to rescale and create imlib im
+        // use the new zscale min/max values to rescale and create imlib im
         ptr += w * h - 1;
         fpixel = 1;
         color = 1;
@@ -241,7 +225,7 @@ load(ImlibImage * im, ImlibProgressFunction progress,
                           buffer, &anynull, &status);
 
              for( x = w - 1; x >= 0; x-- ){
-                 getRGB( buffer[ x ], datamin, datamax, color, inv, &r, &g, &b);
+                 getRGB( buffer[ x ], zscalemin, zscalemax, color, inv, &r, &g, &b);
                  if( x == w / 100 && y == h / 100 ) fprintf( stderr, "\n%d %d %d %g\n",r,g,b, buffer[ x ]);
                  *ptr = (((( (0xff << 8) | r) << 8) | g) << 8) | b ;
                  ptr--;
@@ -269,7 +253,7 @@ load(ImlibImage * im, ImlibProgressFunction progress,
                     }
                }
           }
-          fprintf(stderr, "\ndata min %g, max %g\n", datamin, datamax);
+          fprintf(stderr, "\ndata min %g, max %g\n", zscalemin, zscalemax);
           //free( buffer );
           //free( ptr );
      }
-- 
2.32.0.93.g670b81a890