[gnuastro-commits] master 8276762e: Library (statistics.h): clipping accounts for a spread of zero

[Mohammad Akhlaghi]

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

    Library (statistics.h): clipping accounts for a spread of zero
    
    Until now, when the spread (standard deviation in sigma-clipping or MAD in
    MAD-clipping) became zero, the clipping function would fall into an
    infinite loop, resulting in a final error about the size of the inputs!
    
    With this commit, when checking if we have passed the tolerance, we also
    check if the spread has become zero or not. If so, then the clipping stops.
    
    While investigating the cause of the problem above (and checking if there
    is any memory leaks), I noticed that when writing the value of linked-list
    options, we weren't allocating the space for the value. This has also been
    corrected in 'options.c' in this commit.
---
 lib/arithmetic.c |  8 ++++++--
 lib/options.c    | 10 ++++++----
 lib/statistics.c | 28 +++++++++++++++-------------
 3 files changed, 27 insertions(+), 19 deletions(-)

diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 4eedefd1..c76d983e 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -2037,7 +2037,7 @@ struct multioperandparams
         n=0;                                                            \
         j=tprm->indexs[tind];                                           \
                                                                         \
-        /* Read the necessay values from each input. */                 \
+        /* Put the values from each input into a single array. */       \
         for(i=0;i<p->dnum;++i) pixs[n++]=a[i][j];                       \
                                                                         \
         /* If there are any usable elements, do the measurement. */     \
@@ -2085,7 +2085,7 @@ struct multioperandparams
                       p->operator);                                     \
               }                                                         \
                                                                         \
-            /* If we are on clip-dilate operators, keep the values. */  \
+            /* If we are on clip-fill operators, keep the values. */    \
             if(center)                                                  \
               {                                                         \
                 center[j]=carr[GAL_STATISTICS_CLIP_OUTCOL_MEDIAN];      \
@@ -2537,6 +2537,10 @@ arithmetic_multioperand_clip_fill(struct 
multioperandparams *p,
   gal_data_free(fp.upper);    /* We are freeing these here to avoid*/
   gal_data_free(fp.lower);    /* keeping extra RAM. */
 
+  /* Free the center and spreads (no longer necessary). */
+  gal_data_free(p->center); p->center=NULL;
+  gal_data_free(p->spread); p->spread=NULL;
+
   /* Re-run sigma-clipping on threads. */
   gal_threads_spin_off(multioperand_on_thread, p, p->out->size,
                        numthreads, list->minmapsize, list->quietmmap);
diff --git a/lib/options.c b/lib/options.c
index ff323d72..732ba6e3 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -3388,8 +3388,8 @@ options_as_fits_keywords_write(gal_fits_list_key_t **keys,
   void *vptr;
   int vptrfree;
   uint8_t vtype;
-  char *name, *doc;
   gal_list_str_t *tmp;
+  char *name, *doc, *strval;
 
   for(i=0; !gal_options_is_last(&options[i]); ++i)
     if( options[i].set && option_is_printable(&options[i]) )
@@ -3400,12 +3400,14 @@ options_as_fits_keywords_write(gal_fits_list_key_t 
**keys,
               tmp!=NULL; tmp=tmp->next)
             {
               /* 'name' and 'doc' have a 'const' qualifier. */
-              gal_checkset_allocate_copy(options[i].name, &name);
+              gal_checkset_allocate_copy(tmp->v, &strval);
               gal_checkset_allocate_copy(options[i].doc,  &doc);
+              gal_checkset_allocate_copy(options[i].name, &name);
               gal_fits_key_list_add(keys, GAL_TYPE_STRING, name, 1,
-                                    tmp->v, 0, doc, 1, NULL, 0);
+                                    strval, 1, doc, 1, NULL, 0);
             }
-        /* Normal types. */
+
+        /* Normal (single element) types. */
         else
           {
             /* If the option is associated with a special function for
diff --git a/lib/statistics.c b/lib/statistics.c
index 6da3470e..3593ad58 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -417,8 +417,8 @@ statistics_mad_in_sorted_no_blank(gal_data_t *sorted, 
gal_data_t *med,
     case GAL_TYPE_UINT64: type=GAL_TYPE_INT64; break;
     default:              type=GAL_TYPE_INVALID; /* Not necessary. */
     }
-  use=gal_data_copy_to_new_type(sorted, ( type!=GAL_TYPE_INVALID
-                                          ? type : sorted->type) );
+  use=gal_data_copy_to_new_type(sorted, ( type==GAL_TYPE_INVALID
+                                          ? sorted->type : type ) );
 
   /* Subtract the median from the input. */
   use=gal_arithmetic(GAL_ARITHMETIC_OP_MINUS, 1, flags, use, med);
@@ -2663,17 +2663,19 @@ statistics_clip(gal_data_t *input, float multip, float 
param,
             printf("%-5zu %-10zu %-12.5e %-12.5e\n", num+1, size, center,
                    spread);
 
-          /* If we are working by tolerance, check if we should jump out of
-             the loop. Normally, 'oldstd' should be larger than std,
-             because the possible outliers have been removed. If it is not,
-             it means that we have clipped too much and must stop anyway,
-             so we don't need an absolute value on the difference!
-
-             Note that when all the elements are identical after the clip,
-             'std' will be zero. In this case we shouldn't calculate the
-             tolerance (because it will be infinity and thus lager than the
-             requested tolerance level value). */
-          if( bytolerance && num>0 )
+          /* See if we should break out of the loop:
+             - When the spread is zero we should break out in any case (if
+               it is by tolerance or number of clips): this can happen in
+               two situtaions: when all the elements are identical after
+               the clip (resulting in both MAD and STD to be zero), or when
+               we have three numbers (for example) and two of them are the
+               same (resulting in a MAD of zero).
+             - If we are working by tolerance, normally, 'oldspread' should
+               be larger than 'spread', because the possible outliers have
+               been removed. If it is not, it means that we have clipped
+               too much and must stop anyway, so we don't need an absolute
+               value on the difference! */
+          if( spread==0 || (bytolerance && num>0) )
             if( spread==0 || ((oldspread - spread) / spread) < param )
               {
                 if(spread==0) oldspread=spread;