Changeset 26609

Timestamp:

Jan 14, 2010, 11:15:29 AM (16 years ago)

Author:

Paul Price

Message:

Had wrong sign for penalty function.

File:

• branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (modified) (7 diffs)

branches/eam_branches/20091201/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_WEIGHT                      // Include weight (1/variance) in equation?
@@ -487 +487 @@
                 // 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];
+                matrix->data.F64[index][index] += norm * penalties->data.F32[i];
             }
         }
@@ -1337 +1337 @@
     float peak = 0.0;
     for (int y = - footprint; y <= footprint; y++) {
-        for (int x = - footprint; x <= footprint; x++) {
-            sum += 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
-            peak = PS_MAX(peak, 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm));
-        }
+        for (int x = - footprint; x <= footprint; x++) {
+            sum += 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
+            peak = PS_MAX(peak, 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm));
+        }
     }
 
@@ -1353 +1353 @@
 
     for (int y = - footprint; y <= footprint; y++) {
-        for (int x = - footprint; x <= footprint; x++) {
-            float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
-            if (dflux < 0.02*sum) continue;
-            dflux1 += residual->kernel[y][x];
-            dflux2 += PS_SQR(residual->kernel[y][x]);
-            dmax = PS_MAX(residual->kernel[y][x], dmax);
-            dmin = PS_MIN(residual->kernel[y][x], dmin);
-            npix ++;
-        }
+        for (int x = - footprint; x <= footprint; x++) {
+            float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
+            if (dflux < 0.02*sum) continue;
+            dflux1 += residual->kernel[y][x];
+            dflux2 += PS_SQR(residual->kernel[y][x]);
+            dmax = PS_MAX(residual->kernel[y][x], dmax);
+            dmin = PS_MIN(residual->kernel[y][x], dmin);
+            npix ++;
+        }
     }
     float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
@@ -1398 +1398 @@
     psVector *keepStamps  = psVectorAlloc(stamps->num, PS_TYPE_S32);
     psVectorInit (keepStamps, 0);
-    { 
-        psVector *flux  = psVectorAlloc(stamps->num, PS_TYPE_F32);
-        psVectorInit (flux, 0.0);
-        
-        for (int i = 0; i < stamps->num; i++) {
-            pmSubtractionStamp *stamp = stamps->stamps->data[i];
-            if (!isfinite(stamp->flux)) continue;
-            flux->data.F32[i] = stamp->flux;
-        }
-
-        psVector *index = psVectorSortIndex(NULL, flux);
-        for (int i = 0; (i < stamps->num) && (i < 9); i++) {
-            int n = stamps->num - i - 1;
-            keepStamps->data.S32[index->data.S32[n]] = 1;
-            fprintf (stderr, "keeping %d (%d of %d)\n", index->data.S32[n], n, 9);
-        }
-        psFree (flux);
-        psFree (index);
-
-        // this function is called multiple times in the iteration, but 
-        // we only know after the interation is done if we will try again.
-        // therefore we must save the sample each time, and blow away the old one
-        // if it exists.
-        psFree (kernels->sampleStamps);
-        kernels->sampleStamps = psArrayAllocEmpty(9);
+    {
+        psVector *flux  = psVectorAlloc(stamps->num, PS_TYPE_F32);
+        psVectorInit (flux, 0.0);
+
+        for (int i = 0; i < stamps->num; i++) {
+            pmSubtractionStamp *stamp = stamps->stamps->data[i];
+            if (!isfinite(stamp->flux)) continue;
+            flux->data.F32[i] = stamp->flux;
+        }
+
+        psVector *index = psVectorSortIndex(NULL, flux);
+        for (int i = 0; (i < stamps->num) && (i < 9); i++) {
+            int n = stamps->num - i - 1;
+            keepStamps->data.S32[index->data.S32[n]] = 1;
+            fprintf (stderr, "keeping %d (%d of %d)\n", index->data.S32[n], n, 9);
+        }
+        psFree (flux);
+        psFree (index);
+
+        // this function is called multiple times in the iteration, but
+        // we only know after the interation is done if we will try again.
+        // therefore we must save the sample each time, and blow away the old one
+        // if it exists.
+        psFree (kernels->sampleStamps);
+        kernels->sampleStamps = psArrayAllocEmpty(9);
     }
 
@@ -1498 +1498 @@
             }
 
-            if (keepStamps->data.S32[i]) {
-                psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32);
-                psArrayAdd (kernels->sampleStamps, 9, sample);
-                psFree (sample);
-            }
-
-            pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint);
+            if (keepStamps->data.S32[i]) {
+                psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32);
+                psArrayAdd (kernels->sampleStamps, 9, sample);
+                psFree (sample);
+            }
+
+            pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint);
 
         } else {
@@ -1534 +1534 @@
                 }
             }
-            if (keepStamps->data.S32[i]) {
-                psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32);
-                psArrayAdd (kernels->sampleStamps, 9, sample);
-                psFree (sample);
-            }
-
-            pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint);
+            if (keepStamps->data.S32[i]) {
+                psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32);
+                psArrayAdd (kernels->sampleStamps, 9, sample);
+                psFree (sample);
+            }
+
+            pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint);
         }