Changeset 26593

Timestamp:

Jan 14, 2010, 10:10:51 AM (16 years ago)

Author:

eugene

Message:

clean up residual stats calculation; add sample residual image support; change normalization to use rms not max for odd terms

Location:

branches/eam_branches/20091201/psModules/src/imcombine

Files:

• pmSubtractionAnalysis.c (modified) (1 diff)
• pmSubtractionEquation.c (modified) (4 diffs)
• pmSubtractionKernels.c (modified) (3 diffs)
• pmSubtractionKernels.h (modified) (1 diff)

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

@@ -117 +117 @@
     }
 
+    // sample difference images
+    { 
+        psMetadataAddArray(analysis, PS_LIST_TAIL, "SUBTRACTION.SAMPLE.STAMP.SET", PS_META_DUPLICATE_OK, "Sample Difference Stamps", kernels->sampleStamps);
+    }
 
 #ifdef TESTING

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

@@ -1331 +1331 @@
 }
 
+bool pmSubtractionResidualStats(psVector *fSigRes, psVector *fMaxRes, psVector *fMinRes, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) {
+
+    // XXX measure some useful stats on the residuals
+    float sum = 0.0;
+    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));
+        }
+    }
+
+    // only count pixels with more than X% of the source flux
+    // calculate stdev(dflux)
+    float dflux1 = 0.0;
+    float dflux2 = 0.0;
+    int npix = 0;
+
+    float dmax = 0.0;
+    float dmin = 0.0;
+
+    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 ++;
+        }
+    }
+    float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
+    // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak);
+    psVectorAppend(fSigRes, sigma/sum);
+    psVectorAppend(fMaxRes, dmax/peak);
+    psVectorAppend(fMinRes, dmin/peak);
+    return true;
+}
+
 psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps,
                                            pmSubtractionKernels *kernels)
@@ -1353 +1393 @@
     psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+
+    // we want to save the residual images for the 9 brightest stamps.
+    // identify the 9 brightest stamps
+    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);
+    }
 
     for (int i = 0; i < stamps->num; i++) {
@@ -1427 +1498 @@
             }
 
-            // XXX measure some useful stats on the residuals
-            float sum = 0.0;
-            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));
-                }
-            }
-
-            // only count pixels with more than X% of the source flux
-            // calculate stdev(dflux)
-            float dflux1 = 0.0;
-            float dflux2 = 0.0;
-            int npix = 0;
-
-            float dmax = 0.0;
-            float dmin = 0.0;
-
-            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 ++;
-                }
-            }
-            float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
-            // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak);
-            psVectorAppend(fSigRes, sigma/sum);
-            psVectorAppend(fMaxRes, dmax/peak);
-            psVectorAppend(fMinRes, dmin/peak);
+            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 {
             // Dual convolution
@@ -1490 +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);
         }
 

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

@@ -29 +29 @@
     psFree(kernels->solution1);
     psFree(kernels->solution2);
+    psFree(kernels->sampleStamps);
 }
 
@@ -226 +227 @@
     if (!forceZeroNull && ((uOrder % 2) || (vOrder % 2))) {
         // Odd function
-        scale2D = 1.0 / sqrt(sum2);
+        scale2D = 1.0 / sqrt(sum2);
     }
 
@@ -573 +574 @@
     kernels->rms = NAN;
     kernels->numStamps = 0;
+    kernels->sampleStamps = NULL;
 
     kernels->fSigResMean  = NAN;

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

@@ -55 +55 @@
     float fMinResMean;                  ///< mean fractional negative swing in residuals
     float fMinResStdev;                 ///< stdev of fractional negative swing in residuals
+    psArray *sampleStamps;              ///< array of brightest set of stamps for output visualizations
 } pmSubtractionKernels;