Changeset 26330

Timestamp:

Dec 3, 2009, 2:30:10 AM (16 years ago)

Author:

eugene

Message:

initial mods to calculate normalization and kernel coeffs independently (how about background?)

Location:

branches/eam_branches/20091201/psModules/src/imcombine

Files:

• pmSubtractionEquation.c (modified) (10 diffs)
• pmSubtractionEquation.h (modified) (1 diff)
• pmSubtractionMatch.c (modified) (2 diffs)
• pmSubtractionThreads.c (modified) (1 diff)

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c

@@ -35 +35 @@
                                   const pmSubtractionKernels *kernels, // Kernels
                                   const psImage *polyValues, // Spatial polynomial values
-                                  int footprint // (Half-)Size of stamp
+                                  int footprint, // (Half-)Size of stamp
+                                  const bool normOnly
                                   )
 {
@@ -53 +54 @@
 
     // Evaluate polynomial-polynomial terms
+    // XXX we can skip this for normOnly
     for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) {
         for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) {
@@ -67 +69 @@
 
 
-    for (int i = 0; i < numKernels; i++) {
-        psKernel *iConv = convolutions->data[i]; // Convolution for index i
-        for (int j = i; j < numKernels; j++) {
-            psKernel *jConv = convolutions->data[j]; // Convolution for index j
-
-            double sumCC = 0.0;         // Sum of convolution products
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = - footprint; x <= footprint; x++) {
-                    double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
+    // if we only want to calculate the normalization (and background model?)
+    // then we should assign these values in the matrix and vector to 1 (diagonal) or 0 (off-diagonal)
+    // XXX need to disable normalization calculation if !normOnly
+    if (normOnly) {
+        psImageInit(matrix, 0.0);
+        psVectorInit(vector, 1.0);
+        for (int i = 0; i < matrix->numCols; i++) {
+            matrix->data.F64[i][i] = 1.0;
+        }
+    } else {
+        for (int i = 0; i < numKernels; i++) {
+            psKernel *iConv = convolutions->data[i]; // Convolution for index i
+            for (int j = i; j < numKernels; j++) {
+                psKernel *jConv = convolutions->data[j]; // Convolution for index j
+
+                double sumCC = 0.0;         // Sum of convolution products
+                for (int y = - footprint; y <= footprint; y++) {
+                    for (int x = - footprint; x <= footprint; x++) {
+                        double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
 #ifdef USE_WEIGHT
-                    cc *= weight->kernel[y][x];
+                        cc *= weight->kernel[y][x];
+#endif
+#ifdef USE_WINDOW
+                        if (window) {
+                            cc *= window->kernel[y][x];
+                        }
+#endif
+                        sumCC += cc;
+                    }
+                }
+
+                // Spatial variation
+                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
+                        double value = sumCC * poly2[iTerm][jTerm];
+                        matrix->data.F64[iIndex][jIndex] = value;
+                        matrix->data.F64[jIndex][iIndex] = value;
+                    }
+                }
+            }
+
+            double sumRC = 0.0;             // Sum of the reference-convolution products
+            double sumIC = 0.0;             // Sum of the input-convolution products
+            double sumC = 0.0;              // Sum of the convolution
+            for (int y = - footprint; y <= footprint; y++) {
+                for (int x = - footprint; x <= footprint; x++) {
+                    float conv = iConv->kernel[y][x];
+                    float in = input->kernel[y][x];
+                    float ref = reference->kernel[y][x];
+                    double ic = in * conv;
+                    double rc = ref * conv;
+                    double c = conv;
+#ifdef USE_WEIGHT
+                    float wtVal = weight->kernel[y][x];
+                    ic *= wtVal;
+                    rc *= wtVal;
+                    c *= wtVal;
 #endif
 #ifdef USE_WINDOW
                     if (window) {
-                        cc *= window->kernel[y][x];
+                        float winVal = window->kernel[y][x];
+                        ic *= winVal;
+                        rc *= winVal;
+                        c  *= winVal;
                     }
 #endif
-                    sumCC += cc;
-                }
-            }
-
-            // Spatial variation
-            for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
-                for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
-                    double value = sumCC * poly2[iTerm][jTerm];
-                    matrix->data.F64[iIndex][jIndex] = value;
-                    matrix->data.F64[jIndex][iIndex] = value;
-                }
-            }
-        }
-
-        double sumRC = 0.0;             // Sum of the reference-convolution products
-        double sumIC = 0.0;             // Sum of the input-convolution products
-        double sumC = 0.0;              // Sum of the convolution
-        for (int y = - footprint; y <= footprint; y++) {
-            for (int x = - footprint; x <= footprint; x++) {
-                float conv = iConv->kernel[y][x];
-                float in = input->kernel[y][x];
-                float ref = reference->kernel[y][x];
-                double ic = in * conv;
-                double rc = ref * conv;
-                double c = conv;
-#ifdef USE_WEIGHT
-                float wtVal = weight->kernel[y][x];
-                ic *= wtVal;
-                rc *= wtVal;
-                c *= wtVal;
-#endif
-#ifdef USE_WINDOW
-                if (window) {
-                    float winVal = window->kernel[y][x];
-                    ic *= winVal;
-                    rc *= winVal;
-                    c  *= winVal;
+                    sumIC += ic;
+                    sumRC += rc;
+                    sumC += c;
                 }
-#endif
-                sumIC += ic;
-                sumRC += rc;
-                sumC += c;
-            }
-        }
-        // Spatial variation
-        for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
-            double normTerm = sumRC * poly[iTerm];
-            double bgTerm = sumC * poly[iTerm];
-            matrix->data.F64[iIndex][normIndex] = normTerm;
-            matrix->data.F64[normIndex][iIndex] = normTerm;
-            matrix->data.F64[iIndex][bgIndex] = bgTerm;
-            matrix->data.F64[bgIndex][iIndex] = bgTerm;
-            vector->data.F64[iIndex] = sumIC * poly[iTerm];
-        }
-    }
+            }
+            // Spatial variation
+            for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+                double normTerm = sumRC * poly[iTerm];
+                double bgTerm = sumC * poly[iTerm];
+                matrix->data.F64[iIndex][normIndex] = normTerm;
+                matrix->data.F64[normIndex][iIndex] = normTerm;
+                matrix->data.F64[iIndex][bgIndex] = bgTerm;
+                matrix->data.F64[bgIndex][iIndex] = bgTerm;
+                vector->data.F64[iIndex] = sumIC * poly[iTerm];
+            }
+        }
+    }
+
+    // XXX need to disable normalization calculation if !normOnly
 
     double sumRR = 0.0;                 // Sum of the reference product
@@ -574 +589 @@
     int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
 
-    return pmSubtractionCalculateEquationStamp(stamps, kernels, index);
+    return pmSubtractionCalculateEquationStamp(stamps, kernels, index, normOnly);
 }
 
 bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
-                                         int index)
+                                         int index, bool normOnly)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
@@ -591 +606 @@
     int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
 
+    // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial).  
+    // = \sum_i^N_Gaussians [(order + 1) * (order + 2) / 2], eg for 1 Gauss and 1st order, numKernels = 3
+
     // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the
     // number of coefficients for the spatial polynomial, normalisation and a constant background offset.
@@ -598 +616 @@
     psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE, "We only operate on stamps with this state.");
 
-    // Generate convolutions
+    // Generate convolutions: these are generated once and saved
     if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", index);
@@ -648 +666 @@
         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1,
                                        stamp->weight, stamps->window, stamp->convolutions1, kernels,
-                                       polyValues, footprint);
+                                       polyValues, footprint, normOnly);
         break;
       case PM_SUBTRACTION_MODE_2:
         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2,
                                        stamp->weight, stamps->window, stamp->convolutions2, kernels,
-                                       polyValues, footprint);
+                                       polyValues, footprint, normOnly);
         break;
       case PM_SUBTRACTION_MODE_DUAL:
+        psAbort ("dual is disabled: need to add window and normOnly");
         if (new) {
             stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64);
@@ -739 +758 @@
 }
 
-bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels)
+bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, const bool normOnly)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
@@ -759 +778 @@
             psArrayAdd(job->args, 1, (pmSubtractionKernels*)kernels); // Casting away const to put on array
             PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, normOnly, PS_TYPE_BOOL);
             if (!psThreadJobAddPending(job)) {
                 psFree(job);
@@ -765 +785 @@
             psFree(job);
         } else {
-            pmSubtractionCalculateEquationStamp(stamps, kernels, i);
+            pmSubtractionCalculateEquationStamp(stamps, kernels, i, normOnly);
         }
     }

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.h

@@ -12 +12 @@
 bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, ///< Stamps
                                          const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                         int index ///< Index of stamp
+                                         int index, ///< Index of stamp
+                                         const bool normOnly
                                     );
 
 /// Calculate the least-squares equation to match the image quality
 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps
-                                    const pmSubtractionKernels *kernels ///< Kernel parameters
+                                    const pmSubtractionKernels *kernels, ///< Kernel parameters
+                                    const bool normOnly
                                     );
 

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c

@@ -560 +560 @@
                 }
 
+                // XXX step 1: calculate normalization
+                psTrace("psModules.imcombine", 3, "Calculating equation...\n");
+                if (!pmSubtractionCalculateEquation(stamps, kernels, true)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+
+                // XXX step 2: calculate only other parameters
+# if (0)
                 psTrace("psModules.imcombine", 3, "Calculating equation...\n");
                 if (!pmSubtractionCalculateEquation(stamps, kernels)) {
@@ -565 +574 @@
                     goto MATCH_ERROR;
                 }
+# endif
 
                 memCheck("  calculate equation");

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionThreads.c

@@ -69 +69 @@
 
     {
-        psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION", 3);
+        psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION", 4);
         task->function = &pmSubtractionCalculateEquationThread;
         psThreadTaskAdd(task);