Changeset 30383 for branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c

Timestamp:

Jan 26, 2011, 5:21:42 PM (15 years ago)

Author:

eugene

Message:

report the diff chisq as well as the residual chisq; truncate SVD at dynamic range of 1e10; fix the scaling process

File:

• branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c (modified) (5 diffs)

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c

@@ -476 +476 @@
                         const psArray *sources, const char *stampsName,
                         pmSubtractionKernelsType type, int size, int spatialOrder,
-                        const psVector *isisWidths, const psVector *isisOrders,
+                        psVector *isisWidths, const psVector *isisOrders,
                         int inner, int ringsOrder, int binning, float penalty,
                         bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
@@ -695 +695 @@
             // check on the kernel scaling -- if the kron-based radial moments are very different, adjust to match them
             { 
-                float fwhm1;
-                float fwhm2;
-                pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
-                psAssert(isfinite(fwhm1), "fwhm 1 not set");
-                psAssert(isfinite(fwhm2), "fwhm 2 not set");
+                // float fwhm1;
+                // float fwhm2;
+                // pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
+                // psAssert(isfinite(fwhm1), "fwhm 1 not set");
+                // psAssert(isfinite(fwhm2), "fwhm 2 not set");
 
                 // XXX this is BAD: depends on the relationship below:
@@ -706 +706 @@
                 float radMoment1 = stamps->normWindow1 / 2.75;
                 float radMoment2 = stamps->normWindow2 / 2.75;
-
-                float maxFWHM = PS_MAX(fwhm1, fwhm2);
-                float maxRadial = PS_MAX(radMoment1, radMoment2);
+                pmSubtractionParamsScale(NULL, NULL, isisWidths, radMoment1, radMoment2);
+
+                // float maxFWHM = PS_MAX(fwhm1, fwhm2);
+                // float maxRadial = PS_MAX(radMoment1, radMoment2);
                 
                 // if (fabs(2.0*(maxFWHM - maxRadial)/(maxFWHM + maxRadial)) > 0.25) {
-                if (1) {
-
-                    float scale = maxRadial / maxFWHM;
-                    psLogMsg ("psModules.imcombine", PS_LOG_INFO, "Kron and FWHM scales are quite different, re-scale by %f to use Kron", scale);
-                    
-                    for (int i = 0; i < isisWidths->n; i++) {
-                        isisWidths->data.F32[i] *= scale;
-                    }
-                }
+                // if (1) {
+                // 
+                //     float scale = maxRadial / maxFWHM;
+                //     psLogMsg ("psModules.imcombine", PS_LOG_INFO, "Kron and FWHM scales are quite different, re-scale by %f to use Kron", scale);
+                //     
+                //     for (int i = 0; i < isisWidths->n; i++) {
+                //      isisWidths->data.F32[i] *= scale;
+                //     }
+                // }
             }
 
@@ -1252 +1253 @@
 }
 
-
-bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths,
-                              float scaleRef, float scaleMin, float scaleMax)
+static float scaleRefOption = NAN;
+static float scaleMinOption = NAN;
+static float scaleMaxOption = NAN;
+static bool  scaleOption = false;
+
+bool pmSubtractionParamScaleOptions(bool scale, float scaleRef, float scaleMin, float scaleMax) { 
+
+    if (scale) {
+        PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false);
+        PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false);
+        PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false);
+        PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false);
+    }
+
+    scaleRefOption = scaleRef;
+    scaleMinOption = scaleMin;
+    scaleMaxOption = scaleMax;
+    scaleOption = scale;
+    
+    return true;
+}
+
+bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, float fwhm1, float fwhm2)
 {
-    PS_ASSERT_PTR_NON_NULL(kernelSize, false);
-    PS_ASSERT_PTR_NON_NULL(stampSize, false);
+    // PS_ASSERT_PTR_NON_NULL(kernelSize, false);
+    // PS_ASSERT_PTR_NON_NULL(stampSize, false);
     PS_ASSERT_VECTOR_NON_NULL(widths, false);
     PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false);
-    PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false);
-    PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false);
-    PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false);
-    PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false);
-
-    float fwhm1;
-    float fwhm2;
-
-    pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
-    psAssert(isfinite(fwhm1), "fwhm 1 not set");
-    psAssert(isfinite(fwhm2), "fwhm 2 not set");
+
+    if (!scaleOption) return true;
+
+    // pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
+    // psAssert(isfinite(fwhm1), "fwhm 1 not set");
+    // psAssert(isfinite(fwhm2), "fwhm 2 not set");
     
     // float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference
-    float scale = PS_MAX(fwhm1, fwhm2) / scaleRef;      // Scaling factor
-
-    if (isfinite(scaleMin) && scale < scaleMin) {
-        scale = scaleMin;
-    }
-    if (isfinite(scaleMax) && scale > scaleMax) {
-        scale = scaleMax;
+    float scale = PS_MAX(fwhm1, fwhm2) / scaleRefOption;      // Scaling factor
+
+    if (isfinite(scaleMinOption) && scale < scaleMinOption) {
+        scale = scaleMinOption;
+    }
+    if (isfinite(scaleMaxOption) && scale > scaleMaxOption) {
+        scale = scaleMaxOption;
     }
 
@@ -1285 +1301 @@
         widths->data.F32[i] *= scale;
     }
-    *kernelSize = *kernelSize * scale + 0.5;
-    *stampSize = *stampSize * scale + 0.5;
-
-    psLogMsg("psModules.imcombine", PS_LOG_INFO,
-             "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize);
+    if (kernelSize) {
+        *kernelSize = *kernelSize * scale + 0.5;
+    }
+    if (stampSize) {
+        *stampSize = *stampSize * scale + 0.5;
+    }
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Scaling kernel parameters by %f", scale);
+    if (kernelSize) psLogMsg("psModules.imcombine", PS_LOG_INFO, " modified kernel size %d", *kernelSize);
+    if (stampSize) psLogMsg("psModules.imcombine", PS_LOG_INFO, " modified stamp size %d", *stampSize);
 
     return true;