Changeset 31990

Timestamp:

Aug 3, 2011, 10:30:02 AM (15 years ago)

Author:

eugene

Message:

add new window function using first radial moment; use new sersic guess function

Location:

branches/eam_branches/ipp-20110710/psphot/src

Files:

• psphot.h (modified) (1 diff)
• psphotExtendedSourceFits.c (modified) (15 diffs)
• psphotRadiusChecks.c (modified) (1 diff)
• psphotSersicModelClass.c (modified) (5 diffs)
• psphotSourceFits.c (modified) (19 diffs)

branches/eam_branches/ipp-20110710/psphot/src/psphot.h

@@ -469 +469 @@
 psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc);
 
+bool psphotSersicModelClassGuessPCM (pmPCMdata *pcm, pmSource *source);
+void psphotSersicModelClassInit ();
+void psphotSersicModelClassCleanup ();
+
+bool psphotSetRadiusMoments (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal);
+
 #endif

branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c

@@ -45 +45 @@
     int Nfail = 0;
 
+    psphotSersicModelClassInit();
+
     psTimerStart ("psphot.extended");
 
@@ -103 +105 @@
 
       if (item->type != PS_DATA_METADATA) {
+        // XXX we could cull the bad entries or build a validated model folder
         psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name);
-        // XXX we could cull the bad entries or build a validated model folder
       }
 
@@ -247 +249 @@
     psFree(AnalysisRegion);
 
+    psphotSersicModelClassCleanup();
+
     psLogMsg ("psphot", PS_LOG_INFO, "extended source model fits: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
@@ -265 +269 @@
     int Nplain = 0;
     int NplainPass = 0;
-    bool savePics = false;
-    float radius;
+    float fitRadius, windowRadius;
     psScalar *scalar = NULL;
     pmMoments psfMoments;
@@ -279 +282 @@
     psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
 
-    // pthread_t tid = pthread_self();     // Thread identifier
-
     // Define source fitting parameters for extended source fits
     pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
@@ -309 +310 @@
         if (source->peak->y > region->y1) continue;
 
-        // if model is NULL, we don't have a starting guess
-        // XXX use the parameters guessed from moments
-        // if (source->modelEXT == NULL) continue;
-
-        // fprintf (stderr, "fit %d,%d in thread %d\n", source->peak->x, source->peak->y, (int) tid);
-
         // replace object in image
         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
@@ -321 +316 @@
         Next ++;
 
-        // set the radius based on the footprint (also sets the mask pixels)
-        if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) {
+        // set the radius based on the first radial moment (also sets the mask pixels)
+        if (!psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal)) {
             fprintf (stderr, "skipping (1) %f, %f\n", source->peak->xf, source->peak->yf);
             pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-            // XXX raise an error of some kind
+            // XXX raise an error of some kind or set a flag bit
             continue;
         }
 
-        // XXX note that this changes the source moments that are published...
-        // recalculate the source moments using the larger extended-source moments radius
-        // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
-        // this uses the footprint to judge both radius and aperture?
-        // XXX save the psf-based moments for output
+        // Recalculate the source moments using the larger extended-source moments radius.  At
+        // this stage, skip Gaussian windowing, and do not clip pixels by S/N.  Save the
+        // psf-based moments for output
         psfMoments = *source->moments;
-        if (!pmSourceMoments (source, radius, 0.0, 0.0, 0.0, maskVal)) {
+        if (!pmSourceMoments (source, windowRadius, 0.0, 0.0, 0.0, maskVal & PS_NOT_IMAGE_MASK(markVal))) {
             // subtract the best fit from the object, leave local sky
             fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf);
@@ -344 +337 @@
         }
 
-        // save the modelFlux here in case we need to subtract it (for failure)
+        // save the PSF-based modelFlux here in case we need to subtract it (for failure).  If
+        // it does not exist, attempt to generate it.  Give up if we cannot even do that.
         psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux);
         if (!modelFluxStart) {
@@ -358 +352 @@
         }
         
-        if (savePics) {
-          psphotSaveImage (NULL, readout->image, "image.xp.fits");
-        }
-
         // array to store the pointers to the model flux images while the models are being fitted
         psArray *modelFluxes = psArrayAllocEmpty (models->list->n);
@@ -368 +358 @@
         if (source->modelFits == NULL) {
             source->modelFits = psArrayAllocEmpty (models->list->n);
+        }
+        // extFitPars are the non-parametric data measured for this model fit (moments, kron, etc)
+        if (source->extFitPars == NULL) {
+            source->extFitPars = psArrayAllocEmpty (models->list->n);
         }
 
@@ -397 +391 @@
           bool convolved = psMetadataLookupBool (&status, model, "PSF_CONVOLVED_VALUE");
           assert (status);
-
-          // XXX psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6);
-          // XXX psTraceSetLevel ("psphot.psphotModelWithPSF_LMM", 6);
-
-          // XXX this does not make sense in a threaded context
-          psTimerStart ("psphot.extended.fit");
 
           // fit the model as convolved or not
@@ -439 +427 @@
           }
 
-          psLogMsg ("psphot", PS_LOG_INFO, "model fit: %7.5f sec, %5d npix, %2d iter, %s type\n", psTimerMark ("psphot.extended.fit"), modelFit->nPix, modelFit->nIter, pmModelClassGetName (modelFit->type));
-
-          // save each of the model flux images and store the best
+          // save each of the model flux images for now, and (below) store the best
           psArrayAdd (modelFluxes, 4, source->modelFlux);
 
+          pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
+          extFitPars->Mxx = source->moments->Mxx;
+          extFitPars->Mxy = source->moments->Mxy;
+          extFitPars->Myy = source->moments->Myy;
+          extFitPars->Mrf = source->moments->Mrf;
+          extFitPars->Mrh = source->moments->Mrh;
+          extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
+          
+          // save kron mag, but assign apMag & psfMag on output (not yet calculated)
+          extFitPars->krMag = -2.5*log10(source->moments->KronFlux);
+
+          psArrayAdd (source->extFitPars, 4, extFitPars);
+          psFree (extFitPars);
+
           // test for fit quality / result
-          modelFit->fitRadius = radius;
+          modelFit->fitRadius = fitRadius;
           psArrayAdd (source->modelFits, 4, modelFit);
 
@@ -470 +470 @@
         }
 
+        // clear the circular mask
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
+
         if (minModel == -1) {
-          // re-subtract the object, leave local sky
+          // no valid extended fit; re-subtract the object, leave local sky
           psTrace ("psphot", 5, "failed to fit extended source model to object at %f, %f", source->moments->Mx, source->moments->My);
 
@@ -483 +486 @@
           psFree (modelFluxes);
 
-          if (savePics) {
-              psphotSaveImage (NULL, readout->image, "image.xp.fits");
-              char key[10];
-              fprintf (stdout, "continue? ");
-              if (!fgets (key, 8, stdin)) {
-                  psWarning("Couldn't read anything.");
-              } else if (key[0] == 'n') {
-                  savePics = false;
-              }
-          }
           continue;
         }
@@ -516 +509 @@
         psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq);
         psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-
-        if (savePics) {
-          psphotSaveImage (NULL, readout->image, "image.xm.fits");
-          char key[10];
-          fprintf (stdout, "continue? ");
-          if (!fgets (key, 8, stdin)) {
-              psWarning("Couldn't read anything.");
-          } else if (key[0] == 'n') {
-              savePics = false;
-          }
-        }
     }
     psFree (fitOptions);
-
-    // fprintf (stderr, "xfit : tried %ld objects\n", sources->n);
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)

branches/eam_branches/ipp-20110710/psphot/src/psphotRadiusChecks.c

@@ -142 +142 @@
 }
 
+# define MIN_WINDOW 5.0
+# define SCALE1 5.0
+# define SCALE2 12.5
+
+// call this function whenever you (re)-define the EXT model
+bool psphotSetRadiusMoments (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal) {
+
+    psAssert (source, "source not defined??");
+    psAssert (source->moments, "moments not defined??");
+
+    *fitRadius = SCALE1 * source->moments->Mrf;
+    *fitRadius = PS_MIN (PS_MAX(*fitRadius, MIN_WINDOW), EXT_FIT_MAX_RADIUS);
+
+    *windowRadius = SCALE2 * source->moments->Mrf;
+    *windowRadius = PS_MIN (PS_MAX(*windowRadius, 2.5*MIN_WINDOW), 2.5*EXT_FIT_MAX_RADIUS);
+
+    // redefine the pixels if needed
+    pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, *windowRadius);
+
+    // set the mask to flag the excluded pixels
+    psImageKeepCircle (source->maskObj, source->peak->xf, source->peak->yf, *fitRadius, "OR", markVal);
+
+    return true;
+}
+
 // call this function whenever you (re)-define the EXT model
 bool psphotSetRadiusFootprint (float *radius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float factor) {

branches/eam_branches/ipp-20110710/psphot/src/psphotSersicModelClass.c

@@ -505 +505 @@
 }
 
-void psphotWriteGuess(pmModel *model, float dKronMag, float PSFratio, float ASPratio);
+# define TIMING false
 
 bool psphotSersicModelClassGuessPCM (pmPCMdata *pcm, pmSource *source) {
+
+    float t1, t2, t4, t5;
+    if (TIMING) { psTimerStart ("SersicGuess"); }
 
     psF32 *PAR = pcm->modelConv->params->data.F32;
@@ -549 +552 @@
 
     // find the containing model class:
+
+    if (TIMING) { t1 = psTimerMark ("SersicGuess"); }
 
     psphotSersicModelClass *class = NULL;    
@@ -561 +566 @@
     }
         
+    if (TIMING) { t2 = psTimerMark ("SersicGuess"); }
+
     psAssert (class, "PSFratio and ASPratio must be in range");
 
@@ -610 +617 @@
     PAR[PM_PAR_I0]  = 1.0;
 
+    if (TIMING) { t4 = psTimerMark ("SersicGuess"); }
+
     // set the normalization by linear fit between model and data
     psphotSersicModelNorm (pcm, source);
+
+    if (TIMING) { t5 = psTimerMark ("SersicGuess"); }
+
+    if (TIMING) {
+        fprintf (stderr, "guess, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", t1, t2, t4, t5);
+    }
 
     // float flux = pcm->modelConv->modelFlux(pcm->modelConv->params);
@@ -617 +632 @@
     // float Io = pow(10.0, -0.4*(totalMag - normMag));
     // PAR[PM_PAR_I0] = Io;
-
-    psphotWriteGuess(pcm->modelConv, dKronMag, PSFratio, ASPratio);
 
     return true;

branches/eam_branches/ipp-20110710/psphot/src/psphotSourceFits.c

@@ -211 +211 @@
 bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
 
-    float radius;
+    float fitRadius, windowRadius;
     bool okEXT, okDBL;
     pmModel *ONE = NULL;
@@ -225 +225 @@
     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
 
+# define TEST_X -540.0
+# define TEST_Y 540.0
+   
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        fprintf (stderr, "test galaxy\n");
+    }
+
+# undef TEST_X
+# undef TEST_Y
+
     // set the radius based on the footprint (also sets the mask pixels)
-    if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false;
+    if (!psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal)) return false;
+    // fprintf (stderr, "rad: %6.1f %6.1f  | %5.2f %5.2f %5.2f  ", source->peak->xf, source->peak->yf, source->moments->Mrf, fitRadius, windowRadius);
 
     // XXX note that this changes the source moments that are published...
@@ -234 +245 @@
     // this uses the footprint to judge both radius and aperture?
     // XXX save the psf-based moments for output
+    // XXX do not limit to pixels in the marked circle
     psfMoments = *source->moments;
-    if (!pmSourceMoments (source, radius, 0.0, 0.5, 0.0, maskVal)) {
+    if (!pmSourceMoments (source, windowRadius, 0.0, 0.5, 0.0, maskVal & PS_NOT_IMAGE_MASK(markVal))) {
       *source->moments = psfMoments;
+      fprintf (stderr, "skip (bad moments)\n");
       return false;
     }
+    // fprintf (stderr, "r1: %f\n", source->moments->Mrf);
+
+    // XXX note that we do not re-try to measure the moments if the resulting Mrf is large compared 
+    // to the initial value
 
     psTrace ("psphot", 5, "trying blob...\n");
@@ -263 +280 @@
         ONE = DBL->data[0];
         if (ONE) {
-            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
+            psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
             chiDBL = ONE->chisqNorm; // save chisq for double-star/galaxy comparison
-            ONE->fitRadius = radius;
+            ONE->fitRadius = fitRadius;
         }
 
@@ -271 +288 @@
         ONE = DBL->data[1];
         if (ONE) {
-            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
-            ONE->fitRadius = radius;
+            psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
+            ONE->fitRadius = fitRadius;
         }
     }
@@ -284 +301 @@
         okEXT = psphotEvalEXT (tmpSrc, EXT);
         chiEXT = EXT ? EXT->chisqNorm : NAN;
+        EXT->fitRadius = fitRadius;
     }
 
@@ -476 +494 @@
     // for sersic models, use a grid search to choose an index, then float the params there
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-        // for the test fits, use a somewhat smaller radius
-        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
         psphotFitSersicIndex (model, readout, source, fitOptions, maskVal, markVal);
-    }
-
-    if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
-        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
     }
 
@@ -499 +511 @@
 }
 
+# define TIMING 0
+
 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
 
@@ -517 +531 @@
         return NULL;
     }
+
+# define TEST_X -540.0
+# define TEST_Y 540.0
+   
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
+    }
+
+    float t1, t2, t4, t5;
+    if (TIMING) { psTimerStart ("psphotFitPCM"); }
 
     pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
@@ -524 +548 @@
         return model;
     }
-
-    // use the source moments, etc to guess basic model parameters
-    if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
-        psFree (pcm);
-        model->flags |= PM_MODEL_STATUS_BADARGS;
-        return model;
-    }
-
-    // for sersic models, use a grid search to choose an index, then float the params there
+    if (TIMING) { t1 = psTimerMark ("psphotFitPCM"); }
+
+    // get the guess for sersic models 
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-        // for the test fits, use a somewhat smaller radius
-        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
-
-        if (!psphotFitSersicIndexPCM (pcm, readout, source, fitOptions, maskVal, markVal, psfSize)) {
+        // use the source moments, etc to guess basic model parameters
+        if (!psphotSersicModelClassGuessPCM (pcm, source)) {
             psFree (pcm);
             model->flags |= PM_MODEL_STATUS_BADARGS;
             return model;
         }
-    }
-
-    if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
-        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
-    }
+    } else {
+        // use the source moments, etc to guess basic model parameters
+        if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
+            psFree (pcm);
+            model->flags |= PM_MODEL_STATUS_BADARGS;
+            return model;
+        }
+    }
+
+    if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); }
 
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
@@ -555 +576 @@
     // update the pcm elements if we have changed the circumstance (options.mode or source->pixels)
     pmPCMupdate(pcm, source, &options, model);
+    if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
-    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+    if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }
+
+    if (TIMING) {
+        fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, t1, t2, t4, t5);
+    }
+
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
+    }
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
@@ -566 +596 @@
 }
 
+# undef TEST_X
+# undef TEST_Y
+
 // note that these should be 1/2n of the standard sersic index
-float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
-# define N_INDEX_GUESS 6
+// float indexGuess[] = {0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0};
+// float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
+float indexGuess[] = {1.0, 2.0, 3.0, 4.0};
+# define N_INDEX_GUESS 4
 
 // A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
@@ -586 +621 @@
     float chiSquare[N_INDEX_GUESS];
 
+# define TEST_X -540.0
+# define TEST_Y 540.0
+   
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
+    }
+
     for (int i = 0; i < N_INDEX_GUESS; i++) {
-        model->params->data.F32[PM_PAR_7] = indexGuess[i];
+        model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[i];
 
         if (!model->modelGuess(model, source)) {
@@ -594 +636 @@
         }
 
-        // each time we change the model guess, we need to adjust the radius
-        // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
-        if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
-            psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
-        }
-        
         pmSourceFitModel (source, model, &options, maskVal);
-        // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+        fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
 
         chiSquare[i] = model->chisqNorm;
@@ -616 +652 @@
     assert (iMin >= 0);
 
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
+    }
+
     model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
-    model->params->data.F32[PM_PAR_7] = indexGuess[iMin];
+    model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[iMin];
     model->modelGuess(model, source);
-
-    // each time we change the model guess, we need to adjust the radius
-    // if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
-    //  psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal);
-    // }
 
     return true;
@@ -648 +683 @@
     float xMin = NAN;
     float chiSquare[N_INDEX_GUESS];
+
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
+    }
 
     for (int i = 0; i < N_INDEX_GUESS; i++) {
@@ -657 +696 @@
         }
 
-        // each time we change the model guess, we need to adjust the radius
-        // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
-        if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
-            psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
-        }
-        
 # if (0)
+        // this block is to test the relative speed of straight and PCM fits
         pmSourceFitModel (source, model, &options, maskVal);
 # else
@@ -669 +703 @@
         pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
 # endif
+        fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
         // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
 
@@ -684 +719 @@
     assert (iMin >= 0);
     
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
+    }
+
     model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
     model->params->data.F32[PM_PAR_7] = indexGuess[iMin];