Index: trunk/psphot/src/psphotFindPeaks.c
===================================================================
--- trunk/psphot/src/psphotFindPeaks.c	(revision 33761)
+++ trunk/psphot/src/psphotFindPeaks.c	(revision 42088)
@@ -4,5 +4,5 @@
 // image must be constructed to represent (S/N)^2.  If nMax is non-zero, only return a maximum
 // of nMax peaks
-psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks) {
+psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks, bool firstPass) {
 
     bool status = false;
@@ -10,13 +10,74 @@
     psTimerStart ("psphot.peaks");
 
-    // find the peaks in the smoothed image
-    // 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. 
-	// this will result in a failure.
-        psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
-        return NULL;
+    psArray *peaks = NULL;
+
+    bool useSignalImage = psMetadataLookupF32 (&status, recipe, "PEAKS_USE_SIGNAL_IMAGE"); PS_ASSERT (status, NULL);
+
+    if (firstPass && useSignalImage) {
+	// find the approximate smoothed signal image that corresponds to the desired threshold
+
+	// find all pixels within 25% of the threshold:
+	float minThresh = threshold * 0.80;
+	float maxThresh = threshold * 1.25;
+
+	psImage *smooth_sn = significance->variance;
+	psImage *smooth_im = significance->image;
+
+	psVector *SNvalues = psVectorAllocEmpty (1000, PS_TYPE_F32);
+	for (int iy = 0; iy < smooth_im->numRows; iy++) {
+	    for (int ix = 0; ix < smooth_im->numCols; ix++) {
+		// select all valid pixels within the S/N range
+		if (!isfinite(smooth_sn->data.F32[iy][ix])) continue;
+		if (!isfinite(smooth_im->data.F32[iy][ix])) continue;
+		if (smooth_sn->data.F32[iy][ix] < minThresh) continue;
+		if (smooth_sn->data.F32[iy][ix] > maxThresh) continue;
+		psVectorAppend (SNvalues, smooth_im->data.F32[iy][ix]);
+            }
+        }
+
+	// what is the median of the selected values?
+	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+	if (!psVectorStats(stats, SNvalues, NULL, NULL, 0)) {
+	    // psVectorStats will only fail due to a programming error (e.g., invalid vector type)
+	    psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
+	    psFree (SNvalues);
+	    return NULL;
+	}
+	psFree (SNvalues);
+
+	if (!isfinite(stats->sampleMedian)) {
+	    // could not determine relationship between SN threshold and image values
+	    // XXX fall-back could be the standard analysis above.
+	    psLogMsg ("psphot", PS_LOG_INFO, "failure to map SN to image values");
+	    psFree (stats);
+	    goto use_significance;
+	}
+	
+	float alt_threshold = stats->sampleMedian;
+	psFree (stats);
+
+	peaks = pmPeaksInImage (significance->image, alt_threshold);
+	if (peaks == NULL) {
+	    // we only get a NULL peaks array due to a programming or config error. 
+	    // this will result in a failure.
+	    psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
+	    return NULL;
+	}
     }
+
+use_significance:
+    if (!peaks) {
+	// peaks is NULL at this point if we did not use signal image (by choice or by failure).
+	// find the peaks in the smoothed image
+	// NOTE : significance->variance actually carries the detection S/N image
+	peaks = pmPeaksInImage (significance->variance, threshold);
+	if (peaks == NULL) {
+	    // we only get a NULL peaks array due to a programming or config error. 
+	    // this will result in a failure.
+	    psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
+	    return NULL;
+	}
+    }
+
     // return the total number of peaks found before the nMax limit is applied
     if (totalPeaks) {