Changeset 18239

Timestamp:

Jun 20, 2008, 10:00:01 AM (18 years ago)

Author:

Paul Price

Message:

Adding penalties for large kernel components.

Location:

branches/pap_branch_080617/psModules/src/imcombine

Files:

• pmSubtractionEquation.c (modified) (5 diffs)
• pmSubtractionKernels.c (modified) (13 diffs)
• pmSubtractionKernels.h (modified) (1 diff)
• pmSubtractionMatch.c (modified) (2 diffs)

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

@@ -382 +382 @@
 }
 
+
+// Add in penalty term to least-squares vector
+static bool calculatePenalty(int numPixels, // Number of pixels; for normalisation
+                             psVector *vector, // Vector to which to add in penalty term
+                             const pmSubtractionKernels *kernels // Kernel parameters
+    )
+{
+    if (kernels->penalty == 0.0) {
+        return true;
+    }
+
+    float penalty = numPixels * kernels->penalty; // Penalty value
+    psVector *penalties = kernels->penalties; // Penalties for each kernel component
+    int spatialOrder = kernels->spatialOrder; // Order of spatial variations
+    int numKernels = kernels->num; // Number of kernel components
+    for (int i = 0; i < numKernels; i++) {
+        for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
+            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
+                vector->data.F64[index] -= penalty * penalties->data.F32[i];
+            }
+        }
+    }
+
+    return true;
+}
+
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 // Semi-public functions
-// XXX EAM these cannot be inline :: talk to Josh
+// XXX We might like to define these functions as "extern inline" but gcc currently doesn't handle this in c99
+// mode.  See http://gcc.gnu.org/ml/gcc/2006-11/msg00006.html
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
@@ -675 +702 @@
         psVectorInit(sumVector, 0.0);
         psImageInit(sumMatrix, 0.0);
+        int numStamps = 0;              // Number of good stamps
         for (int i = 0; i < stamps->num; i++) {
             pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
@@ -680 +708 @@
                 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1);
                 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1);
-            }
-        }
+                numStamps++;
+            }
+        }
+        calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector, kernels);
 
         psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
@@ -716 +746 @@
         psVectorInit(sumVector2, 0.0);
 
+        int numStamps = 0;              // Number of good stamps
         for (int i = 0; i < stamps->num; i++) {
             pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
@@ -724 +755 @@
                 (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1);
                 (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2);
-            }
-        }
+                numStamps++;
+            }
+        }
+        calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector1, kernels);
+        calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector2, kernels);
 
 #if 0

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

@@ -24 +24 @@
     psFree(kernels->vStop);
     psFree(kernels->preCalc);
+    psFree(kernels->penalties);
     psFree(kernels->solution1);
     psFree(kernels->solution2);
@@ -60 +61 @@
     kernels->v = psVectorRealloc(kernels->v, start + numNew);
     kernels->preCalc = psArrayRealloc(kernels->preCalc, start + numNew);
+    kernels->penalties = psVectorRealloc(kernels->penalties, start + numNew);
     kernels->inner = start;
 
@@ -74 +76 @@
             kernels->v->data.S32[index] = v;
             kernels->preCalc->data[index] = NULL;
+            kernels->penalties->data.F32[index] = PS_SQR(u) + PS_SQR(v);
 
             psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
@@ -122 +125 @@
                 psKernel *preCalc = psKernelAlloc(-size, size, -size, size);
                 double sum = 0.0;       // Normalisation
+                double moment = 0.0;    // Moment, for penalty
                 for (int v = -size; v <= size; v++) {
                     for (int u = -size; u <= size; u++) {
                         sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) *
                             expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma));
+                        moment += preCalc->kernel[v][u] * sqrtf(PS_SQR(u) + PS_SQR(v));
                     }
                 }
@@ -146 +151 @@
                 }
                 kernels->preCalc->data[index] = preCalc;
+                kernels->penalties->data.F32[index] = moment;
 
                 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index,
@@ -173 +179 @@
     kernels->widths = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     kernels->preCalc = psArrayAlloc(numBasisFunctions);
+    kernels->penalty = 0.0;
+    kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     kernels->uStop = NULL;
     kernels->vStop = NULL;
@@ -324 +332 @@
     psFree(widths);
 
+    psWarning("Kernel penalty for dual-convolution is not configured for SPAM kernels.");
+    psVectorInit(kernels->penalties, 0.0);
+
     return kernels;
 }
@@ -413 +424 @@
     psFree(start);
     psFree(stop);
+
+    psWarning("Kernel penalty for dual-convolution is not configured for FRIES kernels.");
+    psVectorInit(kernels->penalties, 0.0);
 
     return kernels;
@@ -520 +534 @@
                 psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords
                 psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial
+                double moment = 0.0;    // Moment, for penalty
 
                 if (i == 0) {
@@ -527 +542 @@
                     uCoords->n = vCoords->n = poly->n = 1;
                     radiusLast = 0;
+                    moment = 0.0;
                 } else {
                     int j = 0;          // Index for data
@@ -546 +562 @@
                                     poly->data.F32[j] = polyVal;
                                     norm += polyVal;
+                                    moment += polyVal * sqrtf(PS_SQR(u) + PS_SQR(v));
 
                                     psVectorExtend(uCoords, RINGS_BUFFER, 1);
@@ -571 +588 @@
                         psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
                     }
+                    moment /= norm;
                 }
 
@@ -578 +596 @@
                 kernels->u->data.S32[index] = uOrder;
                 kernels->v->data.S32[index] = vOrder;
+                kernels->penalties->data.F32[index] = moment;
 
                 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index,

branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionKernels.h

@@ -33 +33 @@
     psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
     psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS)
+    float penalty;                      ///< Penalty value
+    psVector *penalties;                ///< Penalty for each kernel component
     int size;                           ///< The half-size of the kernel
     int inner;                          ///< The size of an inner region

branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionMatch.c

@@ -430 +430 @@
                 // Generate image with convolution kernels
                 int fullSize = 2 * size + 1 + 1;    // Full size of kernel
-                psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32);
+                int imageSize = (2 * KERNEL_MOSAIC + 1) * fullSize;
+                psImage *convKernels = psImageAlloc((mode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) *
+                                                    imageSize - 1 +
+                                                    (mode == PM_SUBTRACTION_MODE_DUAL ? 4 : 0),
+                                                    imageSize - 1, PS_TYPE_F32);
                 psImageInit(convKernels, NAN);
                 for (int j = -KERNEL_MOSAIC; j <= KERNEL_MOSAIC; j++) {
@@ -451 +455 @@
                         }
                         psFree(kernel);
+
+                        if (mode == PM_SUBTRACTION_MODE_DUAL) {
+                            kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC,
+                                                              (float)j / (float)KERNEL_MOSAIC,
+                                                              true); // Image of the kernel
+                            if (!kernel) {
+                                psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
+                                psFree(convKernels);
+                                goto MATCH_ERROR;
+                            }
+
+                            if (psImageOverlaySection(convKernels, kernel,
+                                                      (2 * KERNEL_MOSAIC + 1 + i + KERNEL_MOSAIC) * fullSize +
+                                                      4,
+                                                      (j + KERNEL_MOSAIC) * fullSize, "=") == 0) {
+                                psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
+                                psFree(kernel);
+                                psFree(convKernels);
+                                goto MATCH_ERROR;
+                            }
+                            psFree(kernel);
+                        }
+
                     }
                 }