Changeset 37804 for branches/eam_branches/ipp-20140904/psLib/src/math/psMinimizePolyFit.c

Timestamp:

Jan 11, 2015, 2:05:35 PM (12 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ipp-20140904

Files:

• . (modified) (1 prop)
• psLib/src/math/psMinimizePolyFit.c (modified) (5 diffs, 1 prop)

branches/eam_branches/ipp-20140904/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;
+    bool median_zero = 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))||
@@ -721 +729 @@
           sigma = 1.0;
         }
-        if ((median == 0.0)&&(sigma == 1.0)) {
-          scale = false;
+        if (fabs(median) < 1e-10) {
+          // CZW 2014-10-21: This median is small and close to zero.  This can cause issues with the
+          // scaling,
+          median_zero = true;
+          if ((sigma == 1.0)||(fabs(sigma) <= 1e-10)) {
+            // Don't bother scaling if sigma is unity (it's already scaled) or if the sigma calculation has gone wrong.
+            scale = false;
+          }
         }
+
         // 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));
@@ -733 +748 @@
 #endif
 
-        
         result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, z64);
         psFree(z64); // Done with this.
@@ -758 +772 @@
 #endif
         }
+
         for (psS32 i = 0; i <= poly->nX; i++) {
           poly->coeff[i] = 0.0;
           poly->coeffErr[i] = 0.0;
-          
-          for (psS32 j = 0; j <= poly->nX; j++) {
+          if (median_zero) {  // If the median is zero, the obtained solution just needs to be scaled by the sigma values.
+            poly->coeff[i]    = Zcoeff[i] * pow(1.0 / sigma,i);
+            poly->coeffErr[i] = ZcoeffErr[i] * pow(1.0 / sigma,i);
+          }
+          else { // Otherwise, do the correct transformations by expanding the (x-m)/s terms.
+            for (psS32 j = 0; j <= poly->nX; j++) {
 #if (CZW)
-            printf("        %d %d %f %f %f %f => %f\n",
-                   i,j,Zcoeff[j],
-                   pow(1.0 / sigma,j) * pow(-1,j - i),
-                   pow(median,j - i),
-                   1.0 * psBinomialCoeff(j,i),
-                   Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j  -i) * pow(median,j - i) * 1.0 * psBinomialCoeff(j,i)
-                   );
+              printf("        %d %d %f %f %f %f => %f\n",
+                     i,j,Zcoeff[j],
+                     pow(1.0 / sigma,j) * pow(-1,j - i),
+                     pow(median,j - i),
+                     1.0 * psBinomialCoeff(j,i),
+                     Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j  -i) * pow(median,j - i) * 1.0 * psBinomialCoeff(j,i)
+                     );
 #endif
-            poly->coeff[i] += Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - i) * psBinomialCoeff(j,i);
-            poly->coeffErr[i] += pow(ZcoeffErr[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - 1) * psBinomialCoeff(j,i),2);
+              poly->coeff[i] += Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - i) * psBinomialCoeff(j,i);
+              poly->coeffErr[i] += pow(ZcoeffErr[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - 1) * psBinomialCoeff(j,i),2);
+            }
+            poly->coeffErr[i] = sqrt(poly->coeffErr[i]);
           }
-          poly->coeffErr[i] = sqrt(poly->coeffErr[i]);
 #if (CZW)
           printf("poly1d: unscaled parameters: %d %f %f\n",
@@ -784 +804 @@
         psFree(ZcoeffErr);
 
-      }
+      } // End scaling block.
         
         break;