Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 42087)
+++ /trunk/psphot/src/psphot.h	(revision 42088)
@@ -216,5 +216,5 @@
 // used by psphotFindDetections
 pmReadout      *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal);
-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            psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal);
 psErrorCode     psphotCullPeaks(const pmReadout *readout, const pmReadout *signifRO, const psMetadata *recipe, psArray *footprints);
Index: /trunk/psphot/src/psphotFindDetections.c
===================================================================
--- /trunk/psphot/src/psphotFindDetections.c	(revision 42087)
+++ /trunk/psphot/src/psphotFindDetections.c	(revision 42088)
@@ -112,5 +112,5 @@
     // detect the peaks in the significance image
     int totalPeaks = 0;
-    detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX, &totalPeaks);
+    detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX, &totalPeaks, firstPass);
     psMetadataAddF32  (readout->analysis, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold);
     if (!detections->peaks) {
Index: /trunk/psphot/src/psphotFindPeaks.c
===================================================================
--- /trunk/psphot/src/psphotFindPeaks.c	(revision 42087)
+++ /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) {
Index: /trunk/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- /trunk/psphot/src/psphotFitSourcesLinear.c	(revision 42087)
+++ /trunk/psphot/src/psphotFitSourcesLinear.c	(revision 42088)
@@ -217,4 +217,12 @@
     int Nsat = 0;
 
+    // track number kept at several stages of the analysis
+    int Nstep_0 = 0;
+    int Nstep_1 = 0;
+    int Nstep_2 = 0;
+    int Nstep_3 = 0;
+    int Nstep_4 = 0;
+    int Nstep_5 = 0;
+
     // select the sources which will be used for the fitting analysis
     for (int i = 0; i < sources->n; i++) {
@@ -246,4 +254,6 @@
         if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
 
+	Nstep_0++;
+
 	// XXX count saturated stars
         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
@@ -255,4 +265,5 @@
             if (!pmSourceCacheModel (source, maskVal)) continue;
         }
+	Nstep_1++;
 
         // save the original coords
@@ -265,4 +276,5 @@
         if (x > AnalysisRegion.x1) continue;
         if (y > AnalysisRegion.y1) continue;
+	Nstep_2++;
 
 	// check the integral of the model : is it large enough?
@@ -283,4 +295,6 @@
             continue;
         }
+	Nstep_3++;
+
 	bool isPSF = false;
         pmModel *model = pmSourceGetModel (&isPSF, source);
@@ -306,5 +320,8 @@
         } 
         if (modelSum  < cutModelSum * normFlux) continue;
+	Nstep_4++;
+
         if (maskedSum < cutMaskedSum * normFlux) continue;
+	Nstep_5++;
 
 	// clear the 'mark' pixels and remask on the fit aperture 
@@ -322,4 +339,5 @@
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), sources->n);
+    psLogMsg ("psphot.ensemble", PS_LOG_INFO, "Number of sources kept at each stage: %d %d %d %d %d %d\n", Nstep_0, Nstep_1, Nstep_2, Nstep_3, Nstep_4, Nstep_5);
 
     if (fitSources->n == 0) {
Index: /trunk/psphot/src/psphotSignificanceImage.c
===================================================================
--- /trunk/psphot/src/psphotSignificanceImage.c	(revision 42087)
+++ /trunk/psphot/src/psphotSignificanceImage.c	(revision 42088)
@@ -105,4 +105,5 @@
 	      //	      float v2 = value + PS_SQR(value/100.0);
 	      // CZW 2013-06-20: I don't think this hack was helping.
+	      // EAM 2021-12-21: see note below about divots in the significance image
 	      float v2 = value;
 	      smooth_wt->data.F32[j][i] = factor * PS_SQR(v2) / smooth_wt->data.F32[j][i];
@@ -122,4 +123,13 @@
     }
     psImageConvolveSetThreads(oldThreads);
+
+    // We now have the significance image and the signal image.  In some cases (e.g.,
+    // stacks), the variance on pixels in the cores of stars is elevated compared to pure
+    // poisson statistics.  In this case, especially at high signal levels, the ratio of
+    // signal / noise in the core of the star can be lower than the surrounding ring of
+    // pixels.  This results in a divot in the center of the star in the significance
+    // image.  the apparent peak of the significance is then not centered on the star and
+    // chaos ensues.  A possible fix is to use the signal image for signficance for
+    // the high S/N detection pass.
 
     pmReadout *significanceRO = pmReadoutAlloc(NULL);
Index: /trunk/psphot/src/psphotSourceFits.c
===================================================================
--- /trunk/psphot/src/psphotSourceFits.c	(revision 42087)
+++ /trunk/psphot/src/psphotSourceFits.c	(revision 42088)
@@ -499,5 +499,7 @@
 
 # if (PS_TRACE_ON) 
-    if (psTraceGetLevel ("psphot") >= 6) {
+    if ( (source->peak->xf > 5500) &&
+	 (source->peak->yf < 500) &&
+	 psTraceGetLevel ("psphot") >= 6) {
 
       // Moments-based shapes parameters
@@ -546,14 +548,14 @@
 
 # if (PS_TRACE_ON) 
-    if (psTraceGetLevel ("psphot") >= 5) {
-      if (psTraceGetLevel ("psphot") >= 6) {
-	fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
-      }
+    if ((source->peak->xf > 5500) &&
+	(source->peak->yf < 500) &&
+	psTraceGetLevel ("psphot") >= 6) {
+      fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
       fprintf (stderr, "--- fitted values ---\n");
       for (int i = 0; i < model->params->n; i++) {
 	fprintf (stderr, "par %d: %f\n", i, model->params->data.F32[i]);
       }
-      psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
-    }
+    }
+    psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
 # endif