Changeset 25833

Timestamp:

Oct 13, 2009, 7:54:28 PM (17 years ago)

Author:

Paul Price

Message:

Attempting to fix dual convolution. Starting off by consolidating and making code plain to do forwards convolution. Adding some small optimisations: it's slightly quicker to generate the solvedKernel using the pre-computed full image kernel, but for convolving stamps, faster to use the separate x and y kernels.

Location:

branches/pap/psModules/src/imcombine

Files:

• pmSubtraction.c (modified) (4 diffs)
• pmSubtractionEquation.c (modified) (4 diffs)
• pmSubtractionKernels.c (modified) (3 diffs)

branches/pap/psModules/src/imcombine/pmSubtraction.c

@@ -135 +135 @@
           case PM_SUBTRACTION_KERNEL_ISIS: {
               psArray *preCalc = kernels->preCalc->data[i]; // Precalculated values
+#if 0
               psVector *xKernel = preCalc->data[0]; // Kernel in x
               psVector *yKernel = preCalc->data[1]; // Kernel in y
               // Iterating over the kernel
               for (int y = 0, v = -size; v <= size; y++, v++) {
+                  float yValue = value * yKernel->data.F32[y];
                   for (int x = 0, u = -size; u <= size; x++, u++) {
-                      kernel->kernel[v][u] +=  value * xKernel->data.F32[x] * yKernel->data.F32[y];
+                      kernel->kernel[v][u] +=  yValue * xKernel->data.F32[x];
                   }
               }
@@ -147 +149 @@
                   kernel->kernel[0][0] -= value;
               }
+#else
+              psKernel *k = preCalc->data[2]; // Kernel image
+              for (int v = -size; v <= size; v++) {
+                  for (int u = -size; u <= size; u++) {
+                      kernel->kernel[v][u] +=  value * k->kernel[v][u];
+                  }
+              }
+#endif
               break;
           }
@@ -597 +607 @@
       case PM_SUBTRACTION_KERNEL_ISIS: {
           psArray *preCalc = kernels->preCalc->data[index]; // Precalculated values
+#if 1
           psVector *xKernel = preCalc->data[0]; // Kernel in x
           psVector *yKernel = preCalc->data[1]; // Kernel in y
@@ -640 +651 @@
           }
           return convolved;
+#else
+          return p_pmSubtractionConvolveStampPrecalc(image, preCalc->data[2]);
+#endif
+
       }
       case PM_SUBTRACTION_KERNEL_RINGS: {

branches/pap/psModules/src/imcombine/pmSubtractionEquation.c

@@ -15 +15 @@
 #include "pmSubtractionVisual.h"
 
-// #define TESTING                         // TESTING output for debugging; may not work with threads!
+//#define TESTING                         // TESTING output for debugging; may not work with threads!
 
 #define USE_VARIANCE                    // Include variance in equation?
@@ -22 +22 @@
 // Private (file-static) functions
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+#if 1
+// Calculate the least-squares matrix and vector
+static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated
+                                  psVector *vector, // Least-squares vector, updated
+                                  const psKernel *input, // Input image (target)
+                                  const psKernel *reference, // Reference image (convolution source)
+                                  const psKernel *variance,  // Variance image
+                                  const psArray *convolutions,         // Convolutions for each kernel
+                                  const pmSubtractionKernels *kernels, // Kernels
+                                  const psImage *polyValues, // Spatial polynomial values
+                                  int footprint // (Half-)Size of stamp
+                                  )
+{
+    // (I - R * sum_i a_i k_i - g) (R * k_j) = 0
+    // I C_j = sum_i C_i C_j
+
+    // Background: C_i = 1.0
+    // Normalisation: C_i = R
+
+    int numKernels = kernels->num;                      // Number of kernels
+    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+    int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
+    int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
+    int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
+    double poly[numPoly];                                 // Polynomial terms
+    double poly2[numPoly][numPoly];                       // Polynomial-polynomial values
+
+    // Evaluate polynomial-polynomial terms
+    for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) {
+        for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) {
+            double iPoly = polyValues->data.F64[iyOrder][ixOrder]; // Value of polynomial
+            poly[iIndex] = iPoly;
+            for (int jyOrder = 0, jIndex = 0; jyOrder <= spatialOrder; jyOrder++) {
+                for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; jxOrder++, jIndex++) {
+                    double jPoly = polyValues->data.F64[jyOrder][jxOrder];
+                    poly2[iIndex][jIndex] = iPoly * jPoly;
+                }
+            }
+        }
+    }
+
+
+    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_VARIANCE
+                    cc /= variance->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_VARIANCE
+                float varVal = 1.0 / variance->kernel[y][x];
+                ic *= varVal;
+                rc *= varVal;
+                c *= varVal;
+#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];
+        }
+    }
+
+    double sumRR = 0.0;                 // Sum of the reference product
+    double sumIR = 0.0;                 // Sum of the input-reference product
+    double sum1 = 0.0;                  // Sum of the background
+    double sumR = 0.0;                  // Sum of the reference
+    double sumI = 0.0;                  // Sum of the input
+    for (int y = - footprint; y <= footprint; y++) {
+        for (int x = - footprint; x <= footprint; x++) {
+            double in = input->kernel[y][x];
+            double ref = reference->kernel[y][x];
+            double ir = in * ref;
+            double rr = PS_SQR(ref);
+#ifdef USE_VARIANCE
+            float varVal = 1.0 / variance->kernel[y][x];
+            rr *= varVal;
+            ir *= varVal;
+            in *= varVal;
+            ref *= varVal;
+#endif
+            sumRR += rr;
+            sumIR += ir;
+            sumR += ref;
+            sumI += in;
+            sum1 += 1.0;
+        }
+    }
+    matrix->data.F64[normIndex][normIndex] = sumRR;
+    matrix->data.F64[bgIndex][bgIndex] = sum1;
+    matrix->data.F64[normIndex][bgIndex] = matrix->data.F64[bgIndex][normIndex] = sumR;
+    vector->data.F64[normIndex] = sumIR;
+    vector->data.F64[bgIndex] = sumI;
+
+    return true;
+}
+
+#if 0
+// Calculate the least-squares matrix and vector for dual convolution
+static bool calculateDualMatrixVector(psImage *matrix1, // Least-squares matrix, updated
+                                      psVector *vector1, // Least-squares vector, updated
+                                      psImage *matrix2,  // Least-squares matrix, updated
+                                      psVector *vector2, // Least-squares vector, updated
+                                      psImage *matrixX,  // Cross-matrix
+                                      const psKernel *input, // Input image (target)
+                                      const psKernel *reference, // Reference image (convolution source)
+                                      const psKernel *variance,  // Variance image
+                                      const psArray *convolutions,         // Convolutions for each kernel
+                                      const pmSubtractionKernels *kernels, // Kernels
+                                      const psImage *polyValues, // Spatial polynomial values
+                                      int footprint // (Half-)Size of stamp
+                                      )
+{
+    // (I - R * sum_i a_i k_i - g) (R * k_j) = 0
+    // I C_j = sum_i C_i C_j
+
+    // Background: C_i = 1.0
+    // Normalisation: C_i = R
+
+    int numKernels = kernels->num;                      // Number of kernels
+    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+    int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
+    int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
+    int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
+    double poly[numPoly];                                 // Polynomial terms
+    double poly2[numPoly][numPoly];                       // Polynomial-polynomial values
+
+    // Evaluate polynomial-polynomial terms
+    for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) {
+        for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) {
+            double iPoly = polyValues->data.F64[iyOrder][ixOrder]; // Value of polynomial
+            poly[iIndex] = iPoly;
+            for (int jyOrder = 0, jIndex = 0; jyOrder <= spatialOrder; jyOrder++) {
+                for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; jxOrder++, jIndex++) {
+                    double jPoly = polyValues->data.F64[jyOrder][jxOrder];
+                    poly2[iIndex][jIndex] = iPoly * jPoly;
+                }
+            }
+        }
+    }
+
+
+    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_VARIANCE
+                    cc /= variance->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_VARIANCE
+                float varVal = 1.0 / variance->kernel[y][x];
+                ic *= varVal;
+                rc *= varVal;
+                c *= varVal;
+#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];
+        }
+    }
+
+    double sumRR = 0.0;                 // Sum of the reference product
+    double sumIR = 0.0;                 // Sum of the input-reference product
+    double sum1 = 0.0;                  // Sum of the background
+    double sumR = 0.0;                  // Sum of the reference
+    double sumI = 0.0;                  // Sum of the input
+    for (int y = - footprint; y <= footprint; y++) {
+        for (int x = - footprint; x <= footprint; x++) {
+            double in = input->kernel[y][x];
+            double ref = reference->kernel[y][x];
+            double ir = in * ref;
+            double rr = PS_SQR(ref);
+#ifdef USE_VARIANCE
+            float varVal = 1.0 / variance->kernel[y][x];
+            rr *= varVal;
+            ir *= varVal;
+            in *= varVal;
+            ref *= varVal;
+#endif
+            sumRR += rr;
+            sumIR += ir;
+            sumR += ref;
+            sumI += in;
+            sum1 += 1.0;
+        }
+    }
+    matrix->data.F64[normIndex][normIndex] = sumRR;
+    matrix->data.F64[bgIndex][bgIndex] = sum1;
+    matrix->data.F64[normIndex][bgIndex] = matrix->data.F64[bgIndex][normIndex] = sumR;
+    vector->data.F64[normIndex] = sumIR;
+    vector->data.F64[bgIndex] = sumI;
+
+    return true;
+}
+#endif
+#endif
+
 
 // Calculate the sum over a stamp product
@@ -582 +860 @@
     switch (kernels->mode) {
       case PM_SUBTRACTION_MODE_1:
+#if 0
         status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
                                  stamp->variance, polyValues, footprint, true);
         status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1,
                                   stamp->image2, stamp->variance, polyValues, footprint, true);
+#else
+        status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1,
+                                       stamp->variance, stamp->convolutions1, kernels, polyValues,
+                                       footprint);
+#endif
         break;
       case PM_SUBTRACTION_MODE_2:
+#if 0
         status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2,
                                  stamp->variance, polyValues, footprint, true);
         status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2,
                                   stamp->image1, stamp->variance, polyValues, footprint, true);
+#else
+        status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2,
+                                       stamp->variance, stamp->convolutions2, kernels, polyValues,
+                                       footprint);
+#endif
         break;
       case PM_SUBTRACTION_MODE_DUAL:
@@ -629 +919 @@
 
 #ifdef TESTING
-    if (psTraceGetLevel("psModules.imcombine.equation") >= 10) {
+    {
         psString matrixName = NULL;
         psStringAppend(&matrixName, "matrix1_%d.fits", index);

branches/pap/psModules/src/imcombine/pmSubtractionKernels.c

@@ -140 +140 @@
         for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
             for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
-                psArray *preCalc = psArrayAlloc(2); // Array to hold precalculated values
+                psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values
                 psVector *xKernel = preCalc->data[0] = subtractionKernelISIS(sigma, uOrder, size); // x Kernel
                 psVector *yKernel = preCalc->data[1] = subtractionKernelISIS(sigma, vOrder, size); // y Kernel
+                psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size);      // Kernel
 
                 // Calculate moments
@@ -149 +150 @@
                     for (int u = -size, x = 0; u <= size; u++, x++) {
                         double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel
+                        kernel->kernel[v][u] = value;
                         moment += value * (PS_SQR(u) + PS_SQR(v));
                     }
@@ -164 +166 @@
                     psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
                     psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
+                    psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32));
+                    kernel->kernel[0][0] -= 1.0;
                     moment *= PS_SQR(sum);
                 }