Changeset 28894

Timestamp:

Aug 11, 2010, 9:18:20 AM (16 years ago)

Author:

eugene

Message:

some test prints are on; kernel second moments added in quadrature to image second moment

Location:

branches/eam_branches/ipp-20100621/psModules/src/imcombine

Files:

• pmSubtractionEquation.c (modified) (2 diffs)
• pmSubtractionKernels.c (modified) (2 diffs)
• pmSubtractionStamps.c (modified) (4 diffs)
• pmSubtractionStamps.h (modified) (1 diff)
• pmSubtractionVisual.c (modified) (1 diff)

branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionEquation.c

@@ -603 +603 @@
                 matrix->data.F64[index][index] += norm * penalties1->data.F32[i];
                 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-                    fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index + numParams + 2][index + numParams + 2], norm * penalties2->data.F32[i], norm, penalties2->data.F32[i]);
+                    fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], 
+                             matrix->data.F64[index + numParams + 2][index + numParams + 2], norm * penalties2->data.F32[i], norm, penalties2->data.F32[i]);
                     matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties2->data.F32[i];                        
                     // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
@@ -1116 +1117 @@
 
         int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 10000.0);
+        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 100.0);
 #endif
 

branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionKernels.c

@@ -218 +218 @@
         if (uOrder % 2 == 0 && vOrder % 2 == 0) {
             // Even functions: normalise to unit sum and subtract null pixel so that sum is zero
-            scale2D = 1.0 / fabs(sum);
+            // Re-normalize 
+            // scale2D  = 1.0 / fabs(sum);
+            scale2D  = 1.0 / sqrt(sum2);
             zeroNull = true;
         } else {
@@ -252 +254 @@
 
     if (zeroNull) {
-        preCalc->kernel->kernel[0][0] -= 1.0;
-    }
-
-#if 0
+        // preCalc->kernel->kernel[0][0] -= 1.0;
+        preCalc->kernel->kernel[0][0] -= sum / sqrt (sum2);
+    }
+
+#if 1
     {
-        double sum = 0.0;   // Sum of kernel component
+        double Sum = 0.0;   // Sum of kernel component
+        double Sum2 = 0.0;   // Sum of kernel component
         float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
         for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
             for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
-                sum += preCalc->kernel->kernel[v][u];
+                double value = preCalc->kernel->kernel[v][u];
+                Sum += value;
+                Sum2 += PS_SQR(value);
                 min = PS_MIN(preCalc->kernel->kernel[v][u], min);
                 max = PS_MAX(preCalc->kernel->kernel[v][u], max);
             }
         }
-        fprintf(stderr, "%d mod: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, sum, preCalc->kernel->kernel[0][0], min, max, scale2D);
+        fprintf(stderr, "%d sum: %lf, sum2: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, Sum, Sum2, preCalc->kernel->kernel[0][0], min, max, scale2D);
     }
 #endif

branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.c

@@ -891 +891 @@
     pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
 
+    bool zeroNull = false;
+    int uOrder = kernels->u->data.S32[index];
+    int vOrder = kernels->v->data.S32[index];
+    if (uOrder % 2 == 0 && vOrder % 2 == 0) zeroNull = true;
+
     if (EMPIRICAL) {
         psKernel *convolution1 = stamp->convolutions1->data[index];
-        penalty1 = pmSubtractionKernelPenaltySingle(convolution1);
+        penalty1 = pmSubtractionKernelPenaltySingle(convolution1, zeroNull);
 
         psKernel *convolution2 = stamp->convolutions2->data[index];
-        penalty2 = pmSubtractionKernelPenaltySingle(convolution2);
+        penalty2 = pmSubtractionKernelPenaltySingle(convolution2, zeroNull);
     } else {
         pmSubtractionKernelPreCalc *kernel = kernels->preCalc->data[index];
-        float M2 = pmSubtractionKernelPenaltySingle(kernel->kernel);
-
-        penalty1 = M2 * PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment 
-        penalty2 = M2 * PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+        float M2 = pmSubtractionKernelPenaltySingle(kernel->kernel, zeroNull);
+
+        penalty1 = M2 + PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+        penalty2 = M2 + PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment 
     }
     kernels->penalties1->data.F32[index] = kernels->penalty * penalty1;
@@ -915 +920 @@
 }
 
-float pmSubtractionKernelPenaltySingle(psKernel *kernel)
+float pmSubtractionKernelPenaltySingle(psKernel *kernel, bool zeroNull)
 {
     // Calculate moments
@@ -924 +929 @@
         for (int u = kernel->xMin; u <= kernel->xMax; u++) {
             double value = kernel->kernel[v][u];
+            if (false && zeroNull && (u == 0) && (v == 0)) {
+                value += 1.0;
+            }
             double value2 = PS_SQR(value);
             sum += value;
@@ -934 +942 @@
     penalty *= 1.0 / sum2;
 
-    if (0) {
-        fprintf(stderr, "min: %lf, max: %lf, moment: %lf\n", min, max, penalty);
+    if (1) {
+        fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2);
         // psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty);
     }

branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.h

@@ -200 +200 @@
 bool pmSubtractionKernelPenaltiesStamp(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels);
 bool pmSubtractionKernelPenalties(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels, int index);
-float pmSubtractionKernelPenaltySingle(psKernel *kernel);
+float pmSubtractionKernelPenaltySingle(psKernel *kernel, bool zeroNull);
 
 #endif

branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionVisual.c

@@ -234 +234 @@
 bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels) {
 
-    if (!pmVisualTestLevel("ppsub.kern.final", 1)) return true;
+    if (!pmVisualTestLevel("ppsub.kernels.final", 1)) return true;
 
     if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {