[gnuastro-commits] master e0c5f91: Library (arithmetic.h): NaN used in multioperands with float output

[Mohammad Akhlaghi]

branch: master
commit e0c5f91dcc9cf52e5c619efe378288a10dbb1275
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Library (arithmetic.h): NaN used in multioperands with float output
    
    Until now, when the inputs to the 'mean' operator were integer types, and
    some of the output pixels weren't covered by any of the inputs, instead of
    a NaN value, those pixels contained large negative values (for example
    -2.14748e+09 when the inputs are 32-bit integers).
    
    When the inputs to the 'mean' operator don't contain any data in a pixel,
    the output pixel should have a NaN value (the output of 'mean' is always a
    floating point numerical data type, irrespective of the input's data
    type). The same applies to the 'sum', 'std', 'median' and all 'sigclip-*'
    operators.
    
    The issue was caused because we were simply writing the input's blank value
    in these pixels of the output (hence a value of ~2 billion for a signed
    integer). But this was wrong, because the output of these operators doesn't
    have the same type as the input.
    
    With this commit, for these operators instead of using the input's blank
    value, we are simply writing a NaN when there is no data in that pixel.
    
    This bug was reported by Zohre Ghaffari.
    
    This fixes bug #61698.
---
 NEWS             |  5 +++++
 lib/arithmetic.c | 21 +++++++++++----------
 2 files changed, 16 insertions(+), 10 deletions(-)

diff --git a/NEWS b/NEWS
index 921ec29..2c5a1c0 100644
--- a/NEWS
+++ b/NEWS
@@ -144,6 +144,11 @@ See the end of the file for license conditions.
               extension while it should have been deleted.
   bug #61598: ds9-region script not using precise WCS values, found and
               fixed by Raúl Infante-Sainz.
+  bug #61698: Mean operator of Arithmetic (for stacking many images) not
+              returning NaN for blank regions when input is integer. This
+              fix also resolves the same problem in the 'sum', 'std',
+              'median' and all 'sigclip-*' operators). This bug was
+              reported by Zohre Ghaffari.
 
 
 
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 5616799..bff7456 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -1131,8 +1131,8 @@ struct multioperandparams
             /* Use in sum if necessary. */                              \
             if(use) { sum += a[i][j]; ++n; }                            \
           }                                                             \
-        o[j] = n ? sum : b;                                             \
-      }                                                                 \
+        o[j] = n ? sum : NAN; /* Not using 'b', because input type */   \
+      }           /* may be integer, while output is always float. */   \
   }
 
 
@@ -1165,8 +1165,8 @@ struct multioperandparams
             /* Calculate the mean if necessary. */                      \
             if(use) { sum += a[i][j]; ++n; }                            \
           }                                                             \
-        o[j] = n ? sum/n : b;                                           \
-      }                                                                 \
+        o[j] = n ? sum/n : NAN; /* Not using 'b', because input type */ \
+      }             /* may be integer, while output is always float. */ \
   }
 
 
@@ -1204,8 +1204,8 @@ struct multioperandparams
                 ++n;                                                    \
               }                                                         \
           }                                                             \
-        o[j] = n ? sqrt( (sum2-sum*sum/n)/n ) : b;                      \
-      }                                                                 \
+        o[j] = n ? sqrt( (sum2-sum*sum/n)/n ) : NAN; /* Not using 'b' */ \
+      } /* because input may be integer, but output is always float. */ \
   }
 
 
@@ -1247,8 +1247,8 @@ struct multioperandparams
             o[j] = n%2 ? pixs[n/2] : (pixs[n/2] + pixs[n/2-1])/2 ;      \
           }                                                             \
         else                                                            \
-          o[j]=b;                                                       \
-      }                                                                 \
+          o[j]=NAN; /* Not using 'b' because input may be integer */    \
+      }                             /* but output is always float.*/    \
                                                                         \
     /* Clean up. */                                                     \
     free(pixs);                                                         \
@@ -1348,8 +1348,8 @@ struct multioperandparams
             cont->size=cont->dsize[0]=p->dnum;                          \
           }                                                             \
         else                                                            \
-          o[j]=b;                                                       \
-      }                                                                 \
+          o[j]=NAN; /* Not using 'b' because input can be inter, but */ \
+      }                                   /* output is always float. */ \
                                                                         \
     /* Clean up (note that 'pixs' is inside of 'cont'). */              \
     gal_data_free(cont);                                                \
@@ -1363,6 +1363,7 @@ struct multioperandparams
     TYPE b, **a;                                                        \
     gal_data_t *tmp;                                                    \
     size_t i=0, tind;                                                   \
+                                                                        \
     /* Allocate space to keep the pointers to the arrays of each. */    \
     /* Input data structure. The operators will increment these */      \
     /* pointers while parsing them. */                                  \