Changeset 37404 for trunk/psLib/src/math/psMinimizePolyFit.c

Timestamp:

Sep 19, 2014, 4:33:54 PM (12 years ago)

Author:

eugene

Message:

do not include masked elements when finding rescaling sigma

File:

• trunk/psLib/src/math/psMinimizePolyFit.c (modified) (2 diffs)

trunk/psLib/src/math/psMinimizePolyFit.c

@@ -700 +700 @@
 
     // Define values that may be used by PS_POLYNOMIAL_ORD form.
+    // these scaling values put the dynamic range of the data into some vaguely sensible location
     bool scale = false;
     double median = 0.0;
     double sigma = 1.0;
-    psVector *sorted = NULL;
-    psVector *z64 = NULL;
-    
+
     switch (poly->type) {
     case PS_POLYNOMIAL_ORD:
       if ((f64->n < 10000)&&(poly->nX > 1)) {
         scale = true;
-        sorted = psVectorSort(NULL,x64);
+
+        // generate a subset of the unmasked values:
+        psVector *tmp = psVectorAllocEmpty (x64->n, PS_TYPE_F64);
+        for (int itmp = 0; itmp < x64->n; itmp++) {
+          if (mask && (mask->data.PS_TYPE_VECTOR_MASK_DATA[itmp] && maskValue)) continue;
+          psVectorAppend (tmp, x64->data.F64[itmp]);
+        }
+        psVector *sorted = psVectorSort(NULL,tmp);
         median = sorted->data.F64[sorted->n / 2];
         // CZW: I'm not bothering to scale this because it doesn't really matter.
         sigma  = (sorted->data.F64[3 * sorted->n / 4] - sorted->data.F64[sorted->n / 4]); 
         psFree(sorted);
+        psFree(tmp);
 
         if ((!isfinite(median))||
@@ -725 +732 @@
         }
         // I can't see a way to not clobber x if it's already F64, so make a copy.x
-        z64 = psVectorCopy(NULL,x64,PS_TYPE_F64);
+        psVector *z64 = psVectorCopy(NULL,x64,PS_TYPE_F64);
         psBinaryOp(z64,z64,"-",psScalarAlloc(median,PS_TYPE_F64));
         psBinaryOp(z64,z64,"/",psScalarAlloc(sigma,PS_TYPE_F64));