gnuastro-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[gnuastro-commits] master 498e1c3: MakeCatalog: correct treating of NaN


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master 498e1c3: MakeCatalog: correct treating of NaN values pixels
Date: Sat, 23 Nov 2019 11:53:02 -0500 (EST)

branch: master
commit 498e1c38b26722c15e6a0d32c363d7fb4a2b4793
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    MakeCatalog: correct treating of NaN values pixels
    
    Until now, MakeCatalog's `p->hasblank' was mistakenly only set based on the
    labels image, not the values image. Also, if there were NaN pixels in the
    values image (and not the labels image), then the brightness and magnitude
    error columns and the S/N column would become NaN.
    
    With this commit, we now check the Sky value and the Sky standard deviation
    value to be NaN before actually using it. It was thus necessary to keep a
    count of how many pixels go into each. Also, when estimating `CCOL_SUM_VAR'
    and `OCOL_SUM_VAR' (the sum of the variances including the values dataset)
    only non-NaN value and STD pixels are used.
    
    This bug was reported by Joseph Putko.
    
    This fixes bug #57293.
---
 NEWS                    |  1 +
 bin/mkcatalog/columns.c |  8 ++---
 bin/mkcatalog/main.h    | 12 ++++---
 bin/mkcatalog/parse.c   | 85 ++++++++++++++++++++++++++++++++++---------------
 bin/mkcatalog/ui.c      |  2 +-
 5 files changed, 73 insertions(+), 35 deletions(-)

diff --git a/NEWS b/NEWS
index b8c275e..f7f3c65 100644
--- a/NEWS
+++ b/NEWS
@@ -132,6 +132,7 @@ See the end of the file for license conditions.
   bug #57164: MakeCatalog crashes when a label isn't in the dataset.
   bug #57180: MakeCatalog reporting infinity S/N when --instd isn't an image.
   bug #57200: Generated pkgconfig must request wcslib, not wcs.
+  bug #57293: NaN value for brightness-related columns when values have NaN.
 
 
 
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index a2f2443..666f1c3 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -2074,12 +2074,12 @@ columns_fill(struct mkcatalog_passparams *pp)
           break;
 
         case UI_KEY_SKY:
-          ((float *)colarr)[oind] = MKC_RATIO(oi[OCOL_SUMSKY], oi[OCOL_NUM]);
+          ((float *)colarr)[oind] = MKC_RATIO(oi[OCOL_SUMSKY], 
oi[OCOL_NUMSKY]);
           break;
 
         case UI_KEY_STD:
           ((float *)colarr)[oind] = sqrt( MKC_RATIO( oi[ OCOL_SUMVAR ],
-                                                     oi[ OCOL_NUM    ]) );
+                                                     oi[ OCOL_NUMVAR ]) );
           break;
 
         case UI_KEY_SEMIMAJOR:
@@ -2323,12 +2323,12 @@ columns_fill(struct mkcatalog_passparams *pp)
 
           case UI_KEY_SKY:
             ((float *)colarr)[cind] = MKC_RATIO( ci[ CCOL_SUMSKY],
-                                                 ci[ CCOL_NUM] );
+                                                 ci[ CCOL_NUMSKY] );
             break;
 
           case UI_KEY_STD:
             ((float *)colarr)[cind] = sqrt( MKC_RATIO( ci[ CCOL_SUMVAR ],
-                                                       ci[ CCOL_NUM    ] ));
+                                                       ci[ CCOL_NUMVAR ] ));
             break;
 
           case UI_KEY_SEMIMAJOR:
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 0aa3f72..e281449 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -75,7 +75,7 @@ enum objectcols
     OCOL_NUM,            /* Number of values used in this object.     */
     OCOL_NUMXY,          /* Number of values in the first two dims.   */
     OCOL_SUM,            /* Sum of (value-sky) in object.             */
-    OCOL_SUM_VAR,        /* Varience of sum (for brightness error).   */
+    OCOL_SUM_VAR,        /* Variance including values (not just sky). */
     OCOL_MEDIAN,         /* Median of value in object.                */
     OCOL_VX,             /* Sum of (value-sky) * x.                   */
     OCOL_VY,             /* Sum of (value-sky) * y.                   */
@@ -84,7 +84,9 @@ enum objectcols
     OCOL_VYY,            /* Sum of (value-sky) * y * y.               */
     OCOL_VXY,            /* Sum of (value-sky) * x * y.               */
     OCOL_SUMSKY,         /* Sum of sky value on this object.          */
+    OCOL_NUMSKY,         /* Number of sky value on this object.       */
     OCOL_SUMVAR,         /* Sum of sky variance value on this object. */
+    OCOL_NUMVAR,         /* Number of sky value on this object.       */
     OCOL_SUMWHT,         /* Sum of positive image pixels.             */
     OCOL_NUMWHT,         /* Number of positive pixels used for wht.   */
     OCOL_GX,             /* Geometric center of object in X.          */
@@ -119,7 +121,7 @@ enum clumpcols
     CCOL_NUM,            /* Number of values used in clump.           */
     CCOL_NUMXY,          /* Number of values only in first two dims.  */
     CCOL_SUM,            /* River subtracted brightness.              */
-    CCOL_SUM_VAR,        /* Variance of sum (for brightness error).   */
+    CCOL_SUM_VAR,        /* Variance including values (not just sky). */
     CCOL_MEDIAN,         /* Median of values in clump.                */
     CCOL_RIV_NUM,        /* Num river pixels around this clump.       */
     CCOL_RIV_SUM,        /* Sum of rivers around clump.               */
@@ -130,8 +132,10 @@ enum clumpcols
     CCOL_VXX,            /* Sum of flux*x*x of this clump.            */
     CCOL_VYY,            /* Sum of flux*y*y of this clump.            */
     CCOL_VXY,            /* Sum of flux*x*y of this clump.            */
-    CCOL_SUMSKY,         /* Sum of sky value on this object.          */
-    CCOL_SUMVAR,         /* Sum of sky variance value on this object. */
+    CCOL_SUMSKY,         /* Sum of sky value on this clump.           */
+    CCOL_NUMSKY,         /* Number of sky value on this clump.        */
+    CCOL_SUMVAR,         /* Sum of sky variance value on this clump.  */
+    CCOL_NUMVAR,         /* Number of sky variance value on this clump.*/
     CCOL_SUMWHT,         /* Sum of positive image pixels for wht.     */
     CCOL_NUMWHT,         /* Num of positive image pixels for wht.     */
     CCOL_GX,             /* Geometric center of clump in X.           */
diff --git a/bin/mkcatalog/parse.c b/bin/mkcatalog/parse.c
index 1229de9..8ab9831 100644
--- a/bin/mkcatalog/parse.c
+++ b/bin/mkcatalog/parse.c
@@ -441,10 +441,10 @@ parse_objects(struct mkcatalog_passparams *pp)
   double *oi=pp->oi;
   gal_data_t *xybin=NULL;
   size_t *tsize=pp->tile->dsize;
-  uint8_t *xybinarr=NULL, *u, *uf;
-  float var, sval, *V=NULL, *SK=NULL, *ST=NULL;
+  uint8_t *u, *uf, goodvalue, *xybinarr=NULL;
   size_t d, pind=0, increment=0, num_increment=1;
   int32_t *O, *OO, *C=NULL, *objarr=p->objects->array;
+  float var, sval, varval, skyval, *V=NULL, *SK=NULL, *ST=NULL;
   float *std=p->std?p->std->array:NULL, *sky=p->sky?p->sky->array:NULL;
 
   /* If tile processing isn't necessary, set `tid' to a blank value. */
@@ -563,8 +563,12 @@ parse_objects(struct mkcatalog_passparams *pp)
 
 
               /* Value related measurements. */
+              goodvalue=0;
               if( p->values && !( p->hasblank && isnan(*V) ) )
                 {
+                  /* For the standard-deviation measurements later. */
+                  goodvalue=1;
+
                   /* General flux summations. */
                   if(xybin) xybinarr[ pind ]=2;
                   if(oif[ OCOL_NUM    ]) oi[ OCOL_NUM     ]++;
@@ -608,13 +612,19 @@ parse_objects(struct mkcatalog_passparams *pp)
 
 
               /* Sky value based measurements. */
-              if(p->sky)
-                if(oif[ OCOL_SUMSKY ])
-                  oi[ OCOL_SUMSKY  ] += ( pp->st_sky
-                                          ? *SK             /* Full array   */
-                                          : ( p->sky->size>1
-                                              ? sky[tid]    /* Tile         */
-                                              : sky[0] ) ); /* Single value */
+              if(p->sky && oif[ OCOL_SUMSKY ])
+                {
+                  skyval = ( pp->st_sky
+                             ? (isnan(*SK)?0:*SK)               /* Full array  
*/
+                             : ( p->sky->size>1
+                                 ? (isnan(sky[tid])?0:sky[tid]) /* Tile        
*/
+                                 : sky[0] ) );                  /* Single 
value*/
+                  if(!isnan(skyval))
+                    {
+                      oi[ OCOL_NUMSKY  ]++;
+                      oi[ OCOL_SUMSKY  ] += skyval;
+                    }
+                }
 
 
               /* Sky standard deviation based measurements.*/
@@ -622,7 +632,11 @@ parse_objects(struct mkcatalog_passparams *pp)
                 {
                   sval = pp->st_std ? *ST : (p->std->size>1?std[tid]:std[0]);
                   var = p->variance ? sval : sval*sval;
-                  if(oif[ OCOL_SUMVAR ]) oi[ OCOL_SUMVAR  ] += var;
+                  if(oif[ OCOL_SUMVAR ] && (!isnan(var)))
+                    {
+                      oi[ OCOL_NUMVAR  ]++;
+                      oi[ OCOL_SUMVAR  ] += var;
+                    }
                   /* For each pixel, we have a sky contribution to the
                      counts and the signal's contribution. The standard
                      deviation in the sky is simply `sval', but the
@@ -633,9 +647,11 @@ parse_objects(struct mkcatalog_passparams *pp)
                      absolute value, because especially as the signal gets
                      noisy there will be negative values, and we don't want
                      them to decrease the variance. */
-                  if(oif[ OCOL_SUM_VAR ])
-                    oi[ OCOL_SUM_VAR ] += ( (p->variance ? var : sval)
-                                            + fabs(*V) );
+                  if(oif[ OCOL_SUM_VAR ] && goodvalue)
+                    {
+                      varval=p->variance ? var : sval;
+                      if(!isnan(varval)) oi[ OCOL_SUM_VAR ] += varval + 
fabs(*V);
+                    }
                 }
             }
 
@@ -718,10 +734,10 @@ parse_clumps(struct mkcatalog_passparams *pp)
   gal_data_t *xybin=NULL;
   size_t *tsize=pp->tile->dsize;
   int32_t *O, *OO, *C=NULL, nlab;
-  uint8_t *u, *uf, *cif=p->ciflag;
-  float var, sval, *V=NULL, *SK=NULL, *ST=NULL;
+  uint8_t *u, *uf, goodvalue, *cif=p->ciflag;
   size_t nngb=gal_dimension_num_neighbors(ndim);
   size_t i, ii, d, pind=0, increment=0, num_increment=1;
+  float var, sval, varval, skyval, *V=NULL, *SK=NULL, *ST=NULL;
   int32_t *objects=p->objects->array, *clumps=p->clumps->array;
   float *std=p->std?p->std->array:NULL, *sky=p->sky?p->sky->array:NULL;
 
@@ -855,8 +871,12 @@ parse_clumps(struct mkcatalog_passparams *pp)
 
                   /* Value related measurements, see `parse_objects' for
                      comments. */
+                  goodvalue=0;
                   if( p->values && !( p->hasblank && isnan(*V) ) )
                     {
+                      /* For the standard-deviation measurement. */
+                      goodvalue=1;
+
                       /* Fill in the necessary information. */
                       if(cif[ CCOL_NUM   ]) ci[ CCOL_NUM ]++;
                       if(cif[ CCOL_SUM   ]) ci[ CCOL_SUM ] += *V;
@@ -882,13 +902,19 @@ parse_clumps(struct mkcatalog_passparams *pp)
                     }
 
                   /* Sky based measurements. */
-                  if(p->sky)
-                    if(cif[ CCOL_SUMSKY ])
-                      ci[ CCOL_SUMSKY  ] += ( pp->st_sky
-                                              ? *SK             /* Full */
-                                              : ( p->sky->size>1
-                                                  ? sky[tid]    /* Tile */
-                                                  : sky[0] ) ); /* 1 value */
+                  if(p->sky && cif[ CCOL_SUMSKY ])
+                    {
+                      skyval = ( pp->st_sky
+                                 ? *SK             /* Full. */
+                                 : ( p->sky->size>1
+                                     ? sky[tid]    /* Tile. */
+                                     : sky[0] ) ); /* 1 value. */
+                      if(!isnan(skyval))
+                        {
+                          ci[ CCOL_NUMSKY  ]++;
+                          ci[ CCOL_SUMSKY  ] += skyval;
+                        }
+                    }
 
                   /* Sky Standard deviation based measurements, see
                      `parse_objects' for comments. */
@@ -898,10 +924,17 @@ parse_clumps(struct mkcatalog_passparams *pp)
                                ? *ST
                                : (p->std->size>1 ? std[tid] : std[0]) );
                       var = p->variance ? sval : sval*sval;
-                      if(cif[ CCOL_SUMVAR  ]) ci[ CCOL_SUMVAR ] += var;
-                      if(cif[ CCOL_SUM_VAR ])
-                        ci[ CCOL_SUM_VAR ] += ( (p->variance ? var : sval)
-                                                + fabs(*V) );
+                      if(cif[ CCOL_SUMVAR  ] && (!isnan(var)))
+                        {
+                          ci[ CCOL_NUMVAR ]++;
+                          ci[ CCOL_SUMVAR ] += var;
+                        }
+                      if(cif[ CCOL_SUM_VAR ] && goodvalue)
+                        {
+                          varval=p->variance ? var : sval;
+                          if(!isnan(varval))
+                            ci[ CCOL_SUM_VAR ] += varval + fabs(*V);
+                        }
                     }
                 }
 
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index adcfe46..c6c1a55 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -1159,7 +1159,7 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
       /* Initially, `p->hasblank' was set based on the objects image, but
          it may happen that the objects image only has zero values for
          blank pixels, so we'll also do a check on the input image. */
-      p->hasblank = gal_blank_present(p->objects, 1);
+      p->hasblank = gal_blank_present(p->values, 1);
 
       /* Reset the units of the value-based columns if the input dataset
          has defined units. */



reply via email to

[Prev in Thread] Current Thread [Next in Thread]