[gnuastro-commits] master eae75ee3: Library (statistics.h): std-dev calculation puts value in double

[Mohammad Akhlaghi]

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

    Library (statistics.h): std-dev calculation puts value in double
    
    Until now, when calculating the standard deviation, to find the sum of the
    squares, we would simply multiply the values in their native data type,
    then put the value in a 'double' type. However, this would cause integer
    overflow when the values were larger than the square root of the maximum
    value! Therefore a next sanity check would produce a zero standard
    deviation.
    
    With this commit, to avoid such situations, a temporary variable of type
    'double' has been introduced in the two functions that do standard
    deviation parsing. The raw value (in any type) is first written into this
    variable, and then used in the sum and sum of squares. This fixes the
    problem.
    
    This fixes bug #63013.
---
 NEWS             |  2 ++
 lib/statistics.c | 20 ++++++++++++++------
 2 files changed, 16 insertions(+), 6 deletions(-)

diff --git a/NEWS b/NEWS
index 8d687764..68e1eb3a 100644
--- a/NEWS
+++ b/NEWS
@@ -178,6 +178,8 @@ See the end of the file for license conditions.
               Found and fixed by Sepideh Eskandarlou.
   bug #62944: No warning when the option value isn't immediately after the
               equal sign in long format. Found by Faezeh Bijarchian.
+  bug #63013: Sigma clip segmentation fault when input has an integer type
+              with values close to saturation-level.
 
 
 
diff --git a/lib/statistics.c b/lib/statistics.c
index 977394eb..910040de 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -227,7 +227,7 @@ gal_data_t *
 gal_statistics_std(gal_data_t *input)
 {
   size_t dsize=1, n=0;
-  double s=0.0f, s2=0.0f, *o;
+  double v, *o, s=0.0f, s2=0.0f;
   gal_data_t *out=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &dsize,
                                  NULL, 1, -1, 1, NULL, NULL, NULL);
 
@@ -246,8 +246,12 @@ gal_statistics_std(gal_data_t *input)
     /* More than one element. */
     default:
 
-      /* Find the 's' and 's2' by parsing the data. */
-      GAL_TILE_PARSE_OPERATE(input, out, 0, 1, {++n; s += *i; s2 += *i * *i;});
+      /* Parse the data to measure 's' and 's2'. Its important to put each
+         value into a 'double' type variable ('v') before multiplying (for
+         's2') because the multiplication of integer types close to their
+         limits will cause overflow and thus an unreasonable output). */
+      GAL_TILE_PARSE_OPERATE(input, out, 0, 1,
+                             {++n; v=*i; s+=v; s2+=v*v;});
 
       /* Write the standard deviation. */
       o[0] = gal_statistics_std_from_sums(s, s2, n);
@@ -269,7 +273,7 @@ gal_data_t *
 gal_statistics_mean_std(gal_data_t *input)
 {
   size_t dsize=2, n=0;
-  double s=0.0f, s2=0.0f, *o;
+  double v, *o, s=0.0f, s2=0.0f;
   gal_data_t *out=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &dsize,
                                  NULL, 1, -1, 1, NULL, NULL, NULL);
 
@@ -291,8 +295,12 @@ gal_statistics_mean_std(gal_data_t *input)
     /* More than one element. */
     default:
 
-      /* Parse the data. */
-      GAL_TILE_PARSE_OPERATE(input, out, 0, 1, {++n; s += *i; s2 += *i * *i;});
+      /* Parse the data. Its important to put each value into a 'double'
+         type variable ('v') before multiplying (for 's2') because the
+         multiplication of integer types close to their limits will cause
+         overflow and thus an unreasonable output). */
+      GAL_TILE_PARSE_OPERATE(input, out, 0, 1,
+                             {++n; v=*i; s+=v; s2+=v*v;});
 
       /* Write the mean */
       o[0]=s/n;