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

Timestamp:

Aug 9, 2010, 8:27:57 AM (16 years ago)

Author:

eugene

Message:

separate penalties for kernels for image 1 and 2, convolve second moments with fwhm to get penalties; separate normalization windows for images 1 and 2; update visual to have facilities like psphot visual; curve of growth analysis for norm window

File:

• branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionEquation.c (modified) (18 diffs)

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

@@ -38 +38 @@
                                   const psImage *polyValues, // Spatial polynomial values
                                   int footprint, // (Half-)Size of stamp
-                                  int normWindow, // Window (half-)size for normalisation measurement
+                                  int normWindow1, // Window (half-)size for normalisation measurement
+                                  int normWindow2, // Window (half-)size for normalisation measurement
                                   const pmSubtractionEquationCalculationMode mode
                                   )
@@ -184 +185 @@
             double one = 1.0;
 
-            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
+            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) {
                 normI1 += ref;
+            }
+            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) {
                 normI2 += in;
             }
@@ -214 +217 @@
 
     *norm = normI2 / normI1;
+
+    fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm);
 
     if (mode & PM_SUBTRACTION_EQUATION_NORM) {
@@ -262 +267 @@
                                       const psImage *polyValues, // Spatial polynomial values
                                       int footprint, // (Half-)Size of stamp
-                                      int normWindow, // Window (half-)size for normalisation measurement
+                                      int normWindow1, // Window (half-)size for normalisation measurement
+                                      int normWindow2, // Window (half-)size for normalisation measurement
                                       const pmSubtractionEquationCalculationMode mode
                                       )
@@ -492 +498 @@
             double i1i2 = i1 * i2;
 
-            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
+            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) {
                 normI1 += i1;
+            }
+            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) {
                 normI2 += i2;
             }
@@ -522 +530 @@
 
     *norm = normI2 / normI1;
+    fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm);
 
     if (mode & PM_SUBTRACTION_EQUATION_NORM) {
@@ -559 +568 @@
 // Add in penalty term to least-squares vector
 bool calculatePenalty(psImage *matrix,                     // Matrix to which to add in penalty term
-                             psVector *vector,                    // Vector to which to add in penalty term
-                             const pmSubtractionKernels *kernels, // Kernel parameters
-                             float norm                           // Normalisation
-    )
+                      psVector *vector,                    // Vector to which to add in penalty term
+                      const pmSubtractionKernels *kernels, // Kernel parameters
+                      float norm                           // Normalisation
+  )
 {
     if (kernels->penalty == 0.0) {
@@ -568 +577 @@
     }
 
-    psVector *penalties = kernels->penalties; // Penalties for each kernel component
+    psVector *penalties1 = kernels->penalties1; // Penalties for each kernel component (input)
+    psVector *penalties2 = kernels->penalties2; // Penalties for each kernel component (ref)
+
     int spatialOrder = kernels->spatialOrder; // Order of spatial variations
     int numKernels = kernels->num; // Number of kernel components
@@ -588 +599 @@
             for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
                 // Contribution to chi^2: a_i^2 P_i
-                psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
-                matrix->data.F64[index][index] += norm * penalties->data.F32[i];
+                psAssert(isfinite(penalties1->data.F32[i]), "Invalid penalty");
+                fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], norm * penalties1->data.F32[i], norm, penalties1->data.F32[i]);
+                matrix->data.F64[index][index] += norm * penalties1->data.F32[i];
                 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-                    matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i];
+                    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]);
+                    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)
                     // penalties scale with second moments
@@ -682 +695 @@
 
     pmSubtractionStampList *stamps = job->args->data[0]; // List of stamps
-    const pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
+    pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
     int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
     pmSubtractionEquationCalculationMode mode  = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model
@@ -689 +702 @@
 }
 
-bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
+bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels,
                                          int index, const pmSubtractionEquationCalculationMode mode)
 {
@@ -778 +791 @@
         status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1,
                                        weight, window, stamp->convolutions1, kernels,
-                                       polyValues, footprint, stamps->normWindow, mode);
+                                       polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode);
         break;
       case PM_SUBTRACTION_MODE_2:
         status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2,
                                        weight, window, stamp->convolutions2, kernels,
-                                       polyValues, footprint, stamps->normWindow, mode);
+                                       polyValues, footprint, stamps->normWindow2, stamps->normWindow1, mode);
         break;
       case PM_SUBTRACTION_MODE_DUAL:
@@ -789 +802 @@
                                            stamp->image1, stamp->image2,
                                            weight, window, stamp->convolutions1, stamp->convolutions2,
-                                           kernels, polyValues, footprint, stamps->normWindow, mode);
+                                           kernels, polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode);
         break;
       default:
@@ -830 +843 @@
 }
 
-bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
+bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels,
                                     const pmSubtractionEquationCalculationMode mode)
 {
@@ -996 +1009 @@
             }
 
+            // double normValue = 1.0;
             double normValue = stats->robustMedian;
             // double bgValue = 0.0;
@@ -1023 +1037 @@
         }
 # endif
+
+#if (1)
+        for (int i = 0; i < solution->n; i++) {
+            fprintf(stderr, "Single solution %d: %lf\n", i, solution->data.F64[i]);
+        }
+#endif
 
         if (!kernels->solution1) {
@@ -1096 +1116 @@
 
         int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 1000.0);
+        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 10000.0);
 #endif
 
@@ -1177 +1197 @@
 
 
-#ifdef TESTING
+#if (1)
         for (int i = 0; i < solution->n; i++) {
             fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]);