Changeset 18239 for branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionEquation.c

Timestamp:

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

Author:

Paul Price

Message:

Adding penalties for large kernel components.

File:

• branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionEquation.c (modified) (5 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