Changeset 30975 for branches/eam_branches/ipp-20110213/psphot/src/psphotFindPeaks.c

Timestamp:

Mar 18, 2011, 2:10:35 PM (15 years ago)

Author:

eugene

Message:

we had 3 items conceptually equivalent to the brightness of a peak: "value", "flux", and "SN"; Ive modified pmPeak to carry explicitly named elements "rawFlux", "smoothFlux", and "detValue". Ive also added the corresponding variance values for rawFlux and smoothFlux. This lets us choose independently of the peak detection process whether flux comparisons are done relative to the smoothed or unsmoothed image. It also simplifies tests of the peak flux compared to something else

File:

• branches/eam_branches/ipp-20110213/psphot/src/psphotFindPeaks.c (modified) (4 diffs)

branches/eam_branches/ipp-20110213/psphot/src/psphotFindPeaks.c

@@ -4 +4 @@
 // image must be constructed to represent (S/N)^2.  If nMax is non-zero, only return a maximum
 // of nMax peaks
-psArray *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax) {
+psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax) {
 
     bool status = false;
@@ -11 +11 @@
 
     // find the peaks in the smoothed image
-    psArray *peaks = pmPeaksInImage (significance, threshold);
+    // NOTE : significance->variance actually carries the detection S/N image
+    psArray *peaks = pmPeaksInImage (significance->variance, threshold);
     if (peaks == NULL) {
         // we only get a NULL peaks array due to a programming or config error. 
@@ -34 +35 @@
     for (int i = 0; i < peaks->n; i++) {
         pmPeak *peak = peaks->data[i];
-        peak->SN = sqrt(peak->value);
-        peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0];
-        // if (peak->flux / peak->value > 5.0/12.0) {
-        //     psWarning ("odd peak levels (1)");
-        // }
-        // if (peak->value / peak->flux > 5*12.0) {
-        //     psWarning ("odd peak levels (2)");
-        // }
+        peak->rawFlux = readout->image->data.F32[peak->y-row0][peak->x-col0];
+        peak->rawFluxStdev = sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]);
+        peak->smoothFlux = significance->image->data.F32[peak->y-row0][peak->x-col0];
+        peak->smoothFluxStdev = peak->smoothFlux / sqrt(significance->variance->data.F32[peak->y-row0][peak->x-col0]);
+        // NOTE smoothFluxStdev is actually (sqrt(variance) / covar_factor)
+
+        // do we need this or not?
+        // peak->SN = sqrt(peak->detValue);
+
         if (readout->variance && isfinite (peak->dx)) {
             peak->dx *= sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]);
@@ -51 +53 @@
 
     // limit the total number of returned peaks as specified
-    psArraySort (peaks, pmPeakSortBySN);
+    psArraySort (peaks, pmPeaksSortByRawFluxDescend);
     if (nMax && (peaks->n > nMax)) {
         psArray *tmpPeaks = psArrayAllocEmpty (nMax);