Changeset 29060 for branches/sc_branches/trunkTest/psphot

Timestamp:

Aug 26, 2010, 9:18:39 AM (16 years ago)

Author:

Serge CHASTEL

Message:

Merging trunk in branch

Location:

branches/sc_branches/trunkTest

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/doc/notes.20100715.txt (copied) (copied from trunk/psphot/doc/notes.20100715.txt )
• psphot/doc/psf.load.txt (copied) (copied from trunk/psphot/doc/psf.load.txt )
• psphot/doc/visual.txt (copied) (copied from trunk/psphot/doc/visual.txt )
• psphot/src/Makefile.am (modified) (1 diff)
• psphot/src/psphot.h (modified) (8 diffs)
• psphot/src/psphotAddNoise.c (modified) (3 diffs)
• psphot/src/psphotApResid.c (modified) (1 diff)
• psphot/src/psphotArguments.c (modified) (1 diff)
• psphot/src/psphotBlendFit.c (modified) (12 diffs)
• psphot/src/psphotChoosePSF.c (modified) (4 diffs)
• psphot/src/psphotEllipticalContour.c (modified) (1 diff)
• psphot/src/psphotExtendedSourceAnalysis.c (modified) (3 diffs)
• psphot/src/psphotExtendedSourceFits.c (modified) (9 diffs)
• psphot/src/psphotFindDetections.c (modified) (1 diff)
• psphot/src/psphotFitSet.c (modified) (1 diff)
• psphot/src/psphotFitSourcesLinear.c (modified) (2 diffs)
• psphot/src/psphotFitSourcesLinearStack.c (modified) (1 diff)
• psphot/src/psphotGuessModels.c (modified) (1 diff)
• psphot/src/psphotLoadSRCTEXT.c (modified) (1 diff)
• psphot/src/psphotMagnitudes.c (modified) (5 diffs)
• psphot/src/psphotModelWithPSF.c (modified) (2 diffs)
• psphot/src/psphotPSFConvModel.c (modified) (2 diffs)
• psphot/src/psphotPetrosianVisual.c (modified) (5 diffs)
• psphot/src/psphotRadiusChecks.c (modified) (9 diffs)
• psphot/src/psphotReadout.c (modified) (1 diff)
• psphot/src/psphotSetThreads.c (modified) (2 diffs)
• psphot/src/psphotSignificanceImage.c (modified) (1 diff)
• psphot/src/psphotSourceFits.c (modified) (23 diffs)
• psphot/src/psphotSourceSize.c (modified) (8 diffs)
• psphot/src/psphotSourceStats.c (modified) (4 diffs)
• psphot/src/psphotStack.c (modified) (1 diff)
• psphot/src/psphotStackMatchPSFsUtils.c (modified) (2 diffs)
• psphot/src/psphotVisual.c (modified) (46 diffs)

branches/sc_branches/trunkTest/psphot/src/Makefile.am

@@ -163 +163 @@
         psphotOutput.c                 \
         psphotFakeSources.c            \
-        psphotModelWithPSF.c           \
         psphotExtendedSourceAnalysis.c \
         psphotExtendedSourceAnalysisByObject.c \
         psphotExtendedSourceFits.c     \
         psphotKernelFromPSF.c          \
-        psphotPSFConvModel.c           \
         psphotFitSet.c                 \
         psphotSourceFreePixels.c       \

branches/sc_branches/trunkTest/psphot/src/psphot.h

@@ -12 +12 @@
 
 #define PSPHOT_RECIPE_PSF_FAKE_ALLOW "PSF.FAKE.ALLOW" // Name for recipe component permitting fake PSFs
-
-// pmPCMData : PSF Convolved Model data storage structure
-typedef struct {
-    psImage *model;
-    psArray *dmodels;
-    psImage *modelConv;
-    psArray *dmodelsConv;
-} pmPCMData;
 
 // top-level psphot functions
@@ -114 +106 @@
 bool            psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view, const char *filerule);
 bool            psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+bool            psphotExtendedSourceFits_Threaded (psThreadJob *job);
 
 bool            psphotApResid (pmConfig *config, const pmFPAview *view, const char *filerule);
@@ -188 +181 @@
 
 // functions to set the correct source pixels
-bool            psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type);
+bool            psphotInitRadiusPSF (psMetadata *recipe, pmReadout *readout);
+bool            psphotInitRadiusEXT (psMetadata *recipe, pmReadout *readout);
 
 bool            psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal);
 bool            psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal, float dR);
-bool            psphotInitRadiusEXT (psMetadata *recipe, pmModelType type);
-bool            psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal);
-float           psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal);
+bool            psphotSetRadiusFootprint (float *radius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float factor);
+bool            psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep);
 
 bool            psphotDumpMoments (psMetadata *recipe, psArray *sources);
@@ -210 +203 @@
 //  functions to support the source fitting process
 bool            psphotInitLimitsPSF (psMetadata *recipe, pmReadout *readout);
-bool            psphotInitLimitsEXT (psMetadata *recipe);
-bool            psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal);
-bool            psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal);
-bool            psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal);
-pmModel        *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal);
-psArray        *psphotFitDBL (pmReadout *readout, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
+bool            psphotInitLimitsEXT (psMetadata *recipe, pmReadout *readout);
+bool            psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
+bool            psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
+bool            psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
+pmModel        *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal);
+psArray        *psphotFitDBL (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
 
 // functions to support simultaneous multi-source fitting
@@ -252 +245 @@
 bool            psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout);
 bool            psphotVisualShowSignificance (psImage *image, float min, float max);
+bool            psphotVisualShowLogSignificance (psImage *image, float min, float max);
 bool            psphotVisualShowPeaks (pmDetections *detections);
 bool            psphotVisualShowFootprints (pmDetections *detections);
@@ -292 +286 @@
 bool psphotRadialBins (psMetadata *recipe, pmSource *source, float radiusMax, float skynoise);
 
-// structures & functions to support psf-convolved model fitting
-
-// psf-convolved model fitting
-bool psphotModelWithPSF_LMM (
-    psMinimization *min,
-    psImage *covar,
-    psVector *params,
-    psMinConstraint *constraint,
-    pmSource *source,
-    const psKernel *psf,
-    psMinimizeLMChi2Func func);
-
-psF32 psphotModelWithPSF_SetABX(
-    psImage  *alpha,
-    psVector *beta,
-    const psVector *params,
-    const psVector *paramMask,
-    pmPCMData *pcm,
-    const pmSource *source,
-    const psKernel *psf,
-    psMinimizeLMChi2Func func);
-
-pmPCMData *pmPCMDataAlloc (
-    const psVector *params,
-    const psVector *paramMask,
-    pmSource *source);
-
-psImage *pmPCMDataSaveImage (pmPCMData *pcm);
-
 int psphotKapaOpen (void);
 bool psphotKapaClose (void);
@@ -364 +329 @@
 
 bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final);
-int pmPhotObjSortBySN (const void **a, const void **b);
-int pmPhotObjSortByX (const void **a, const void **b);
 
 typedef enum {
@@ -462 +425 @@
 bool psphotStackObjectsUnifyPosition (psArray *objects);
 
+bool psphotFitSersicIndex (pmModel *model, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
+
+bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
+pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
+
 #endif

branches/sc_branches/trunkTest/psphot/src/psphotAddNoise.c

@@ -34 +34 @@
 
     bool status = false;
-    psEllipseShape oldshape;
-    psEllipseShape newshape;
-    psEllipseAxes axes;
 
     // find the currently selected readout
@@ -86 +83 @@
         if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
 
-        // select appropriate model
-        pmModel *model = pmSourceGetModel (NULL, source);
-        if (model == NULL) continue;  // model must be defined
-
-        if (add) {
-            psTrace ("psphot", 4, "adding noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
-        } else {
-            psTrace ("psphot", 4, "remove noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
-        }
-
-        psF32 *PAR = model->params->data.F32;
-
-        // save original values
-        float oldI0  = PAR[PM_PAR_I0];
-        oldshape.sx  = PAR[PM_PAR_SXX];
-        oldshape.sy  = PAR[PM_PAR_SYY];
-        oldshape.sxy = PAR[PM_PAR_SXY];
-
-        // XXX can this be done more intelligently?
-        if (oldI0 == 0.0) continue;
-        if (!isfinite(oldI0)) continue;
-
-        // increase size and height of source
-        axes = psEllipseShapeToAxes (oldshape, 20.0);
-        axes.major *= SIZE;
-        axes.minor *= SIZE;
-        newshape = psEllipseAxesToShape (axes);
-        PAR[PM_PAR_I0]  = FACTOR*oldI0;
-        PAR[PM_PAR_SXX] = newshape.sx;
-        PAR[PM_PAR_SYY] = newshape.sy;
-        PAR[PM_PAR_SXY] = newshape.sxy;
-
-        // XXX if we use pmSourceOp, the size (and possibly Io) will not be respected
-        pmSourceOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, add, maskVal, 0, 0);
-
-        // restore original values
-        PAR[PM_PAR_I0]  = oldI0;
-        PAR[PM_PAR_SXX] = oldshape.sx;
-        PAR[PM_PAR_SYY] = oldshape.sy;
-        PAR[PM_PAR_SXY] = oldshape.sxy;
+        pmSourceNoiseOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, FACTOR, SIZE, add, maskVal, 0, 0);
     }
     if (add) {
@@ -132 +90 @@
         psLogMsg ("psphot.noise", PS_LOG_INFO, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
     }
+
     return true;
 }

branches/sc_branches/trunkTest/psphot/src/psphotApResid.c

@@ -459 +459 @@
         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
 
-        bool status = pmSourceMagnitudes (source, psf, photMode, maskVal);
+        bool status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
 
         // clear the mask bit

branches/sc_branches/trunkTest/psphot/src/psphotArguments.c

@@ -176 +176 @@
     }
 
-    // visual : interactive display mode
-    if ((N = psArgumentGet (argc, argv, "-visual"))) {
-        psArgumentRemove (N, &argc, argv);
-        pmVisualSetVisual(true);
-    }
-
     // break : used from recipe throughout psphotReadout
     if ((N = psArgumentGet (argc, argv, "-break"))) {

branches/sc_branches/trunkTest/psphot/src/psphotBlendFit.c

@@ -65 +65 @@
     assert (status && fitIter > 0);
 
-    float fitTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_TOL"); // Fit tolerance
-    assert (status && isfinite(fitTol) && fitTol > 0);
+    float fitMinTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MIN_TOL"); // Fit tolerance
+    if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
+        fitMinTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance
+        if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
+            psAbort("PSF_FIT_MIN_TOL (and PSF_FIT_TOL) not defined or positive");
+        }
+    }
+
+    float fitMaxTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_TOL"); // Fit tolerance
+    if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) {
+        fitMaxTol = 1.0;
+    }
 
     bool poisson = psMetadataLookupBool(&status, recipe, "POISSON.ERRORS.PHOT.LMM"); // Poisson errors?
     assert (status);
 
+    float maxChisqDOF = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_CHISQ"); // Fit tolerance
+
     float skySig = psMetadataLookupF32(&status, recipe, "SKY_SIG");
     assert (status && isfinite(skySig) && skySig > 0);
 
     // Define source fitting parameters for extended source fits
-    pmSourceFitModelInit(fitIter, fitTol, PS_SQR(skySig), poisson);
+    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
+    fitOptions->nIter         = fitIter;
+    fitOptions->minTol        = fitMinTol;
+    fitOptions->maxTol        = fitMaxTol;
+    fitOptions->maxChisqDOF   = maxChisqDOF;
+    fitOptions->poissonErrors = poisson;
+    fitOptions->weight        = PS_SQR(skySig);
+    fitOptions->mode          = PM_SOURCE_FIT_PSF;
 
     psphotInitLimitsPSF (recipe, readout);
-    psphotInitLimitsEXT (recipe);
-    psphotInitRadiusPSF (recipe, readout->analysis, psf->type);
+    psphotInitLimitsEXT (recipe, readout);
+    psphotInitRadiusPSF (recipe, readout);
 
     // starts the timer, sets up the array of fitSets
@@ -88 +107 @@
     if (!sources->n) {
         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
+        psFree (fitOptions);
         return true;
     }
@@ -112 +132 @@
             psArrayAdd(job->args, 1, psf);
             psArrayAdd(job->args, 1, newSources); // return for new sources
+            psArrayAdd(job->args, 1, fitOptions); // default fit options
             psFree (newSources);
 
@@ -121 +142 @@
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree (fitOptions);
                 return NULL;
             }
-
-# if (0)
-            {
-                int nfit = 0;
-                int npsf = 0;
-                int next = 0;
-                int nfail = 0;
-                psArray *newSources = psArrayAllocEmpty(16);
-
-                if (!psphotBlendFit_Unthreaded (&nfit, &npsf, &next, &nfail, readout, recipe, cells->data[j], psf, newSources)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                    return NULL;
-                }
-                Nfit += nfit;
-                Npsf += npsf;
-                Next += next;
-                Nfail += nfail;
-
-                // add these back onto sources
-                for (int k = 0; k < newSources->n; k++) {
-                    psArrayAdd (sources, 16, newSources->data[k]);
-                }
-                psFree (newSources);
-            }
-# endif
         }
 
@@ -153 +150 @@
         if (!psThreadPoolWait (false)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            psFree (fitOptions);
             return NULL;
         }
@@ -163 +161 @@
             } else {
                 psScalar *scalar = NULL;
-                scalar = job->args->data[5];
+                scalar = job->args->data[6];
                 Nfit += scalar->data.S32;
-                scalar = job->args->data[6];
+                scalar = job->args->data[7];
                 Npsf += scalar->data.S32;
-                scalar = job->args->data[7];
+                scalar = job->args->data[8];
                 Next += scalar->data.S32;
-                scalar = job->args->data[8];
+                scalar = job->args->data[9];
                 Nfail += scalar->data.S32;
 
@@ -186 +184 @@
       psphotSaveImage (NULL, readout->image,  "image.v2.fits");
     }
+    psFree (fitOptions);
 
     psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit);
@@ -204 +203 @@
     psScalar *scalar = NULL;
 
-    pmReadout *readout  = job->args->data[0];
-    psMetadata *recipe  = job->args->data[1];
-    psArray *sources    = job->args->data[2];
-    pmPSF *psf          = job->args->data[3];
-    psArray *newSources = job->args->data[4];
+    pmReadout *readout             = job->args->data[0];
+    psMetadata *recipe             = job->args->data[1];
+    psArray *sources               = job->args->data[2];
+    pmPSF *psf                     = job->args->data[3];
+    psArray *newSources            = job->args->data[4];
+    pmSourceFitOptions *fitOptions = job->args->data[5];
 
     // bit-masks to test for good/bad pixels
@@ -269 +269 @@
         Nfit ++;
 
+        if (0) {
+            psF32 Mxx = source->moments->Mxx;
+            psF32 Myy = source->moments->Myy;
+            fprintf (stderr, "1: Mxx: %f, Myy: %f\n", Mxx, Myy);
+        }
+
         // try fitting PSFs or extended sources depending on source->mode
         // these functions subtract the resulting fitted source
         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-            if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
+            if (psphotFitBlob (readout, source, newSources, psf, fitOptions, maskVal, markVal)) {
                 source->type = PM_SOURCE_TYPE_EXTENDED;
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
@@ -280 +286 @@
             }
         } else {
-            if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
+            if (psphotFitBlend (readout, source, psf, fitOptions, maskVal, markVal)) {
                 source->type = PM_SOURCE_TYPE_STAR;
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
@@ -289 +295 @@
         }
 
+        if (0) {
+            psF32 Mxx = source->moments->Mxx;
+            psF32 Myy = source->moments->Myy;
+            fprintf (stderr, "2: Mxx: %f, Myy: %f\n", Mxx, Myy);
+        }
+
         psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);
         Nfail ++;
@@ -298 +310 @@
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    scalar = job->args->data[5];
+    scalar = job->args->data[6];
     scalar->data.S32 = Nfit;
 
-    scalar = job->args->data[6];
+    scalar = job->args->data[7];
     scalar->data.S32 = Npsf;
 
-    scalar = job->args->data[7];
+    scalar = job->args->data[8];
     scalar->data.S32 = Next;
 
-    scalar = job->args->data[8];
+    scalar = job->args->data[9];
     scalar->data.S32 = Nfail;
 

branches/sc_branches/trunkTest/psphot/src/psphotChoosePSF.c

@@ -74 +74 @@
 
     // structure to store user options defining the psf
-    pmPSFOptions *options = pmPSFOptionsAlloc ();
+    pmPSFOptions *options = pmPSFOptionsAlloc();
 
     // load user options from the recipe. no need to check existence -- they are
@@ -138 +138 @@
         return false;
     }
-    float fitTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance
-    if (!status || !isfinite(fitTol) || fitTol <= 0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSF_FIT_TOL is not positive");
-        return false;
-    }
-    pmSourceFitModelInit(fitIter, fitTol, PS_SQR(SKY_SIG), options->poissonErrorsPhotLMM);
-
+    float fitMinTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MIN_TOL"); // Fit tolerance
+    if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
+        fitMinTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance
+        if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSF_FIT_MIN_TOL (and PSF_FIT_TOL) not defined or positive");
+            return false;
+        }
+    }
+    float fitMaxTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MAX_TOL"); // Fit tolerance
+    if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) {
+        fitMaxTol = 1.0;
+    }
+    float maxChisqDOF = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MAX_CHISQ"); // Fit tolerance
+
+    // options which modify the behavior of the model fitting
+    options->fitOptions                = pmSourceFitOptionsAlloc();
+    options->fitOptions->nIter         = fitIter;
+    options->fitOptions->minTol        = fitMinTol;
+    options->fitOptions->maxTol        = fitMaxTol;
+    options->fitOptions->maxChisqDOF   = maxChisqDOF;
+    options->fitOptions->poissonErrors = options->poissonErrorsPhotLMM;
+    options->fitOptions->weight        = PS_SQR(SKY_SIG);
+    options->fitOptions->mode          = PM_SOURCE_FIT_PSF;
+    
     psArray *stars = psArrayAllocEmpty (sources->n);
 
@@ -227 +244 @@
 
     // try each model option listed in config
+    // pmPSFtryModel makes a local copy of the sources -- those points are not the same as those for 'sources'
     for (int i = 0; i < modelNames->n; i++) {
         char *modelName = modelNames->data[i];
@@ -304 +322 @@
 
     // unset the PSFSTAR flag for stars not used for PSF model
+    // XXX a more efficient way of achieving this would be to record a pair of arrays
+    // of the source index and the source id for the psf stars.  but that would require we do
+    // not re-sort the source list in the meanwhile
+    int nDrop = 0;
     for (int i = 0; i < try->sources->n; i++) {
         pmSource *source = try->sources->data[i];
         if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-            source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
-        }
-    }
+            // need to find this source in the original list (these are copies, not pointers)
+            for (int j = 0; j < sources->n; j++) {
+                pmSource *realSource = sources->data[j];
+                if (realSource->id != source->id) continue;
+                realSource->mode &= ~PM_SOURCE_MODE_PSFSTAR;
+                source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
+                nDrop ++;
+                break;
+            }
+        }
+    }
+    // fprintf (stderr, "drop %d stars as PSF stars\n", nDrop);
+
+    // XXX is this working?
+    // int N1 = 0;
+    // for (int i = 0; i < try->sources->n; i++) {
+    //     pmSource *source = try->sources->data[i];
+    //  fprintf (stderr, "%llx : %d\n", (long long int) source, (source->mode & PM_SOURCE_MODE_PSFSTAR));
+    //  if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
+    //      N1 ++;
+    //     }
+    // }
+    // int N2 = 0;
+    // for (int i = 0; i < sources->n; i++) {
+    //     pmSource *source = sources->data[i];
+    //  fprintf (stderr, "%llx : %d\n", (long long int) source, (source->mode & PM_SOURCE_MODE_PSFSTAR));
+    //  if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
+    //      N2 ++;
+    //     }
+    // }
+    // fprintf (stderr, "N1: %d, N2: %d\n", N1, N2);
 
     // build a PSF residual image

branches/sc_branches/trunkTest/psphot/src/psphotEllipticalContour.c

@@ -82 +82 @@
     params->data.F32[PAR_RMIN]    = Rmin;
 
-    psMinimization *myMin = psMinimizationAlloc (25, 0.001);
+    psMinimization *myMin = psMinimizationAlloc (25, 0.01, 1.00);
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
     

branches/sc_branches/trunkTest/psphot/src/psphotExtendedSourceAnalysis.c

@@ -121 +121 @@
         if (source->peak->y > AnalysisRegion.y1) continue;
 
+        // fprintf (stderr, "xsrc: %f, %f : %f\n", source->peak->xf, source->peak->yf, source->peak->SN); 
+
         // replace object in image
         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
@@ -157 +159 @@
         }
 
-# if (0)
-        // Isophotal Mags
-        if (doIsophotal) {
-            if (!psphotIsophotal (source, recipe, maskVal)) {
-                psTrace ("psphot", 5, "failed to measure isophotal mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-            } else {
-                psTrace ("psphot", 5, "measured isophotal mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-                Nisophot ++;
-                source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
-            }
-        }
-        // Kron Mags
-        if (doKron) {
-            if (!psphotKron (source, recipe, maskVal)) {
-                psTrace ("psphot", 5, "failed to measure kron mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-            } else {
-                psTrace ("psphot", 5, "measure kron mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-                Nkron ++;
-                source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
-            }
-        }
-# endif
-
         // re-subtract the object, leave local sky
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
@@ -201 +180 @@
     }
 
+    // fprintf (stderr, "xsrc : tried %ld objects\n", sources->n);
+
     return true;
 }

branches/sc_branches/trunkTest/psphot/src/psphotExtendedSourceFits.c

@@ -31 +31 @@
 // non-linear model fitting for extended sources
 bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+
+    bool status;
+    int Next = 0;
+    int Nconvolve = 0;
+    int NconvolvePass = 0;
+    int Nplain = 0;
+    int NplainPass = 0;
+
+    psTimerStart ("psphot.extended");
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->allSources;
+    psAssert (sources, "missing sources?");
+
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
+        return true;
+    }
+
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    assert (maskVal);
+
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels
+    assert (markVal);
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    // select the collection of desired models
+    psMetadata *models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS");
+    if (!status) {
+        psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n");
+        return true;
+    }
+    if (models->list->n == 0) {
+        psWarning ("extended source model fits requested but no models are specified\n");
+        return true;
+    }
+
+    // validate the model entries
+    psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
+    psMetadataItem *item = NULL;
+    while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
+
+      if (item->type != PS_DATA_METADATA) {
+        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
+      }
+
+      psMetadata *model = (psMetadata *) item->data.md;
+
+      // check on the model type
+      char *modelName = psMetadataLookupStr (&status, model, "MODEL");
+      int modelType = pmModelClassGetType (modelName);
+      if (modelType < 0) {
+        psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName);
+      }
+      psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType);
+
+      // check on the SNLIM, set a float value
+      char *SNword = psMetadataLookupStr (&status, model, "SNLIM");
+      if (!status) {
+        psAbort("SNLIM not defined for extended source model %s\n", item->name);
+      }
+      float SNlim = atof (SNword);
+      psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim);
+
+      // check on the PSF-Convolution status
+      char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED");
+      if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) {
+        psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name);
+      }
+      bool convolved = !strcasecmp (convolvedWord, "true");
+      psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved);
+    }
+    psFree (iter);
+
+    // option to limit analysis to a specific region
+    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
+    psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0);
+    *AnalysisRegion = psRegionForImage(readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined");
+
+    // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
+    int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
+    assert (status);
+
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
+    int Cx = 1, Cy = 1;
+    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
+
+    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
+
+    for (int i = 0; i < cellGroups->n; i++) {
+
+        psArray *cells = cellGroups->data[i];
+
+        for (int j = 0; j < cells->n; j++) {
+
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_EXTENDED_FIT");
+
+            psArrayAdd(job->args, 1, readout);
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, models);
+            psArrayAdd(job->args, 1, AnalysisRegion); // XXX make a pointer
+
+            PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
+            PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK);
+
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nconvolve
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NconvolvePass
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nplain
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NplainPass
+
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree(AnalysisRegion);
+                return false;
+            }
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            psFree(AnalysisRegion);
+            return false;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            } else {
+                psScalar *scalar = NULL;
+                scalar = job->args->data[7];
+                Next += scalar->data.S32;
+                scalar = job->args->data[8];
+                Nconvolve += scalar->data.S32;
+                scalar = job->args->data[9];
+                NconvolvePass += scalar->data.S32;
+                scalar = job->args->data[10];
+                Nplain += scalar->data.S32;
+                scalar = job->args->data[11];
+                NplainPass += scalar->data.S32;
+            }
+            psFree(job);
+        }
+    }
+    psFree (cellGroups);
+    psFree(AnalysisRegion);
+
+    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
+    psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
+    psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
+    return true;
+}
+
+// non-linear model fitting for extended sources
+bool psphotExtendedSourceFits_Threaded (psThreadJob *job) {
 
     bool status;
@@ -39 +221 @@
     int NplainPass = 0;
     bool savePics = false;
-
-    // find the currently selected readout
-    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
-    psAssert (file, "missing file?");
-
-    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
-    psAssert (readout, "missing readout?");
-
-    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
-    psAssert (detections, "missing detections?");
-
-    psArray *sources = detections->allSources;
-    psAssert (sources, "missing sources?");
-
-    if (!sources->n) {
-        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
-        return true;
-    }
-
-    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
-    assert (maskVal);
-
-    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels
-    assert (markVal);
-
-    // maskVal is used to test for rejected pixels, and must include markVal
-    maskVal |= markVal;
-
-    // select the collection of desired models
-    psMetadata *models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS");
-    if (!status) {
-        psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n");
-        return true;
-    }
-    if (models->list->n == 0) {
-        psWarning ("extended source model fits requested but no models are specified\n");
-        return true;
-    }
-
-    // validate the model entries
-    psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
-    psMetadataItem *item = NULL;
-    while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
-
-      if (item->type != PS_DATA_METADATA) {
-        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
-      }
-
-      psMetadata *model = (psMetadata *) item->data.md;
-
-      // check on the model type
-      char *modelName = psMetadataLookupStr (&status, model, "MODEL");
-      int modelType = pmModelClassGetType (modelName);
-      if (modelType < 0) {
-        psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName);
-      }
-      psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType);
-
-      // check on the SNLIM, set a float value
-      char *SNword = psMetadataLookupStr (&status, model, "SNLIM");
-      if (!status) {
-        psAbort("SNLIM not defined for extended source model %s\n", item->name);
-      }
-      float SNlim = atof (SNword);
-      psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim);
-
-      // check on the PSF-Convolution status
-      char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED");
-      if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) {
-        psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name);
-      }
-      bool convolved = !strcasecmp (convolvedWord, "true");
-      psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved);
-    }
-    psFree (iter);
-
-    // option to limit analysis to a specific region
-    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
-    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
-    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
-
-    // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
-    int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
-    assert (status);
-
-    // source analysis is done in S/N order (brightest first)
-    sources = psArraySort (sources, pmSourceSortBySN);
+    float radius;
+    psScalar *scalar = NULL;
+
+    // arguments: readout, sources, models, region, psfSize, maskVal, markVal
+    pmReadout *readout      = job->args->data[0];
+    psArray *sources        = job->args->data[1];
+    psMetadata *models      = job->args->data[2];
+    psRegion *region        = job->args->data[3];
+    int psfSize             = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
+
+    // Define source fitting parameters for extended source fits
+    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
+    fitOptions->mode          = PM_SOURCE_FIT_EXT;
+    // XXX for now, use the defaults for the rest:
+    // fitOptions->nIter         = fitIter;
+    // fitOptions->tol           = fitTol;
+    // fitOptions->poissonErrors = poisson;
+    // fitOptions->weight        = PS_SQR(skySig);
 
     // choose the sources of interest
@@ -140 +253 @@
 
         // XXX this should use peak?
-        if (source->peak->x < AnalysisRegion.x0) continue;
-        if (source->peak->y < AnalysisRegion.y0) continue;
-        if (source->peak->x > AnalysisRegion.x1) continue;
-        if (source->peak->y > AnalysisRegion.y1) continue;
+        if (source->peak->x < region->x0) continue;
+        if (source->peak->y < region->y0) continue;
+        if (source->peak->x > region->x1) continue;
+        if (source->peak->y > region->y1) continue;
 
         // if model is NULL, we don't have a starting guess
@@ -155 +268 @@
         Next ++;
 
+        // set the radius based on the footprint (also sets the mask pixels)
+        if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) {
+            fprintf (stderr, "skipping (1) %f, %f\n", source->peak->xf, source->peak->yf);
+            psFree (fitOptions)
+            return false;
+        }
+
+        // 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
+        if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) {
+            fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf);
+            // subtract the best fit from the object, leave local sky
+            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+            // XXX raise an error of some kind
+            continue;
+        }
+
         // save the modelFlux here in case we need to subtract it (for failure)
         psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux);
-
+        if (!modelFluxStart) {
+            pmSourceCacheModel (source, maskVal);
+            modelFluxStart = psMemIncrRefCounter (source->modelFlux);
+            if (!modelFluxStart) {
+                fprintf (stderr, "skipping (3) %f, %f\n", source->peak->xf, source->peak->yf);
+                // XXX raise an error of some kind?
+                continue;
+            }
+        }
+        
         if (savePics) {
           psphotSaveImage (NULL, readout->image, "image.xp.fits");
@@ -183 +325 @@
           assert (status);
 
+          // fprintf (stderr, "xfit: %f, %f : %f\n", source->peak->xf, source->peak->yf, source->peak->SN); 
+
           // limit selection to some SN limit
           assert (source->peak); // how can a source not have a peak?
@@ -201 +345 @@
           pmModel *modelFit = NULL;
           if (convolved) {
-              modelFit = psphotPSFConvModel (readout, source, modelType, maskVal, markVal, psfSize);
+              modelFit = psphotFitPCM (readout, source, fitOptions, modelType, maskVal, markVal, psfSize);
               if (!modelFit) {
                   psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My);
@@ -215 +359 @@
               psFree (source->modelFlux);
               source->modelFlux = NULL;
-              modelFit = psphotFitEXT (readout, source, modelType, maskVal, markVal);
+              modelFit = psphotFitEXT (readout, source, fitOptions, modelType, maskVal, markVal);
               if (!modelFit) {
                   psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
@@ -233 +377 @@
 
           // test for fit quality / result
+          modelFit->fitRadius = radius;
           psArrayAdd (source->modelFits, 4, modelFit);
 
@@ -311 +456 @@
         }
     }
-
-    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next);
-    psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
-    psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
+    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)
+    scalar = job->args->data[7];
+    scalar->data.S32 = Next;
+
+    scalar = job->args->data[8];
+    scalar->data.S32 = Nconvolve;
+
+    scalar = job->args->data[9];
+    scalar->data.S32 = NconvolvePass;
+
+    scalar = job->args->data[10];
+    scalar->data.S32 = Nplain;
+
+    scalar = job->args->data[11];
+    scalar->data.S32 = NplainPass;
+
     return true;
 }

branches/sc_branches/trunkTest/psphot/src/psphotFindDetections.c

@@ -90 +90 @@
     psphotVisualShowSignificance (significance, -1.0, PS_SQR(3.0*NSIGMA_PEAK));
 
-    // XXX getting some strange results from significance image
-    if (0) {
-        psImage *lsig = (psImage *) psUnaryOp (NULL, significance, "log");
-        psphotVisualShowSignificance (lsig, 0.0, 4.0);
-        psFree (lsig);
-    }   
+    // display the log significance image
+    psphotVisualShowLogSignificance (significance, 0.0, 4.5);
 
     // detect the peaks in the significance image

branches/sc_branches/trunkTest/psphot/src/psphotFitSet.c

@@ -24 +24 @@
     }
 
+    // Define source fitting parameters for extended source fits
+    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
+    fitOptions->mode          = PM_SOURCE_FIT_EXT;
+    // XXX for now, use the defaults for the rest:
+    // fitOptions->nIter         = fitIter;
+    // fitOptions->tol           = fitTol;
+    // fitOptions->poissonErrors = poisson;
+    // fitOptions->weight        = PS_SQR(skySig);
+
     // XXX pmSourceFitSet must cache the modelFlux?
-    pmSourceFitSet (source, modelSet, mode, maskVal);
+    pmSourceFitSet (source, modelSet, fitOptions, maskVal);
 
     // write out positive object

branches/sc_branches/trunkTest/psphot/src/psphotFitSourcesLinear.c

@@ -100 +100 @@
     if (!status) {
         SKY_FIT_LINEAR = false;
+    }
+    
+    float MIN_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MIN_VALID_FLUX");
+    if (!status) {
+        MIN_VALID_FLUX = 1e-8;
+    }
+    float MAX_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MAX_VALID_FLUX");
+    if (!status) {
+        MAX_VALID_FLUX = 1e+8;
     }
 
@@ -228 +237 @@
 
     psSparseConstraint constraint;
-    constraint.paramMin   = 0.0;
-    constraint.paramMax   = 1e8;
-    constraint.paramDelta = 1e8;
+    constraint.paramMin   = MIN_VALID_FLUX;
+    constraint.paramMax   = MAX_VALID_FLUX;
+    constraint.paramDelta = 1e7;
 
     // solve for normalization terms (need include local sky?)

branches/sc_branches/trunkTest/psphot/src/psphotFitSourcesLinearStack.c

@@ -170 +170 @@
     return true;
 }
-
-// sort by X (ascending)
-int pmPhotObjSortByX (const void **a, const void **b)
-{
-    pmPhotObj *objA = *(pmPhotObj **)a;
-    pmPhotObj *objB = *(pmPhotObj **)b;
-
-    psF32 fA = objA->x;
-    psF32 fB = objB->x;
-
-    psF32 diff = fA - fB;
-    if (diff > FLT_EPSILON) return (+1);
-    if (diff < FLT_EPSILON) return (-1);
-    return (0);
-}

branches/sc_branches/trunkTest/psphot/src/psphotGuessModels.c

@@ -80 +80 @@
 
     // setup the PSF fit radius details
-    psphotInitRadiusPSF (recipe, readout->analysis, psf->type);
+    psphotInitRadiusPSF (recipe, readout);
 
     // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)

branches/sc_branches/trunkTest/psphot/src/psphotLoadSRCTEXT.c

@@ -84 +84 @@
             source->peak->yf   = PAR[PM_PAR_YPOS]; // but we know the pixel coordinate
 
-            source->pixWeight = 1.0;
+            source->pixWeightNotBad = 1.0;
+            source->pixWeightNotPoor = 1.0;
             source->crNsigma  = 0.0;
             source->extNsigma = 0.0;

branches/sc_branches/trunkTest/psphot/src/psphotMagnitudes.c

@@ -124 +124 @@
                 return false;
             }
-
-# if (0)
-                int nap = 0;
-                if (!psphotMagnitudes_Unthreaded (&nap, cells->data[j], psf, binning, backModel, backStdev, photMode, maskVal)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                    return false;
-                }
-                Nap += nap;
-# endif
         }
 
@@ -186 +177 @@
         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
 
-        status = pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
+        status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
         if (status && isfinite(source->apMag)) Nap ++;
 
@@ -268 +259 @@
             psArrayAdd(job->args, 1, cells->data[j]); // sources
             PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
 
             if (!psThreadJobAddPending(job)) {
@@ -304 +296 @@
     psArray *sources                = job->args->data[0];
     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[1],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA);
 
     for (int i = 0; i < sources->n; i++) {
@@ -312 +305 @@
         if (model == NULL) {
           psTrace ("psphot", 3, "fail mag : no valid model");
-          source->pixWeight = NAN;
+          source->pixWeightNotBad = NAN;
+          source->pixWeightNotPoor = NAN;
           continue;
         }
 
-        status = pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
+        status = pmSourcePixelWeight (&source->pixWeightNotBad, &source->pixWeightNotPoor, model, source->maskObj, maskVal, markVal);
         if (!status) {
           psTrace ("psphot", 3, "fail to measure pixel weight");
-          source->pixWeight = NAN;
+          source->pixWeightNotBad = NAN;
+          source->pixWeightNotPoor = NAN;
           continue;
         }

branches/sc_branches/trunkTest/psphot/src/psphotModelWithPSF.c

@@ -73 +73 @@
 
     // iterate until the tolerance is reached, or give up
-    while ((min->iter < min->maxIter) && ((min->lastDelta > min->tol) || !isfinite(min->lastDelta))) {
+    while ((min->iter < min->maxIter) && ((min->lastDelta > min->minTol) || !isfinite(min->lastDelta))) {
         psTrace("psphot", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
-        psTrace("psphot", 5, "Last delta is %f.  Min->tol is %f.\n", min->lastDelta, min->tol);
+        psTrace("psphot", 5, "Last delta is %f.  Min->minTol is %f.\n", min->lastDelta, min->minTol);
 
 
@@ -166 +166 @@
     psFree(pcm);
 
-    if (min->iter == min->maxIter) {
-        psTrace("psphot", 3, "---- end (false) ----\n");
-        return(false);
-    }
-
-    psTrace("psphot", 3, "---- end (true) ----\n");
-    return(true);
+    // if the last improvement was at least as good as maxTol, accept the fit:
+    if (min->lastDelta <= min->maxTol) {
+        psTrace("psphot", 6, "---- end (true) ----\n");
+        return(true);
+    }
+    psTrace("psphot", 6, "---- end (false) ----\n");
+    return(false);
 }
 

branches/sc_branches/trunkTest/psphot/src/psphotPSFConvModel.c

@@ -4 +4 @@
 // save as static values so they may be set externally
 static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
-static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
+static psF32 PM_SOURCE_FIT_MODEL_MIN_TOL = 0.1;
+static psF32 PM_SOURCE_FIT_MODEL_MAX_TOL = 2.0;
 
 // input source has both modelPSF and modelEXT.  on successful exit, we set the
 // modelConv to contain the fitted parameters, and the modelFlux to contain the 
 // convolved model image.
-pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
+// XXX need to generalize this -- number of fitted parameters must be flexible based on the fitOptions
+
+pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
     
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -90 +94 @@
 
     // set up the minimization process
-    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
+    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_MIN_TOL, PM_SOURCE_FIT_MODEL_MAX_TOL);
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);

branches/sc_branches/trunkTest/psphot/src/psphotPetrosianVisual.c

@@ -1 +1 @@
 # include "psphotInternal.h"
-# define FORCE_VISUAL 0
 
 // this function displays representative images as the psphot analysis progresses:
@@ -54 +53 @@
     Graphdata graphdata;
 
-    // return true;
-    if (!FORCE_VISUAL && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.petro.byangle", 2)) return true;
 
     if (kapa2 == -1) {
@@ -101 +99 @@
     Graphdata graphdata;
 
-    // return true;
-    if (!FORCE_VISUAL && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.petro.radii", 2)) return true;
 
     if (kapa == -1) {
@@ -173 +170 @@
     KapaSection section;
 
-    if (!FORCE_VISUAL && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.petro.stats", 2)) return true;
 
     if (kapa2 == -1) {
@@ -311 +308 @@
     Graphdata graphdata;
 
-    if (!FORCE_VISUAL && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.petro.ellipse", 2)) return true;
 
     if (kapa == -1) {

branches/sc_branches/trunkTest/psphot/src/psphotRadiusChecks.c

@@ -8 +8 @@
                                         // and a per-object radius is calculated)
 
-bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type) {
+bool psphotInitRadiusPSF(psMetadata *recipe, pmReadout *readout) {
 
     bool status = true;
@@ -15 +15 @@
     PSF_FIT_PADDING = psMetadataLookupF32(&status, recipe, "PSF_FIT_PADDING");
 
-    PSF_FIT_RADIUS =  psMetadataLookupF32(&status, analysis, "PSF_FIT_RADIUS");
+    PSF_FIT_RADIUS =  psMetadataLookupF32(&status, readout->analysis, "PSF_FIT_RADIUS");
     if (!status) {
         PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS");
     }
 
-    PSF_APERTURE =  psMetadataLookupF32(&status, analysis, "PSF_APERTURE");
+    PSF_APERTURE =  psMetadataLookupF32(&status, readout->analysis, "PSF_APERTURE");
     if (!status) {
         PSF_APERTURE =  psMetadataLookupF32(&status, recipe, "PSF_APERTURE");
@@ -28 +28 @@
 
     if (PSF_FIT_RADIUS == 0.0) {
-        float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA");
+        float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
         if (!status) {
             gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
@@ -37 +37 @@
 
     if (PSF_APERTURE == 0.0) {
-        float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA");
+        float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
         if (!status) {
             gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
@@ -122 +122 @@
 }
 
+static float EXT_FIT_SKY_SIG;
 static float EXT_FIT_NSIGMA;
 static float EXT_FIT_PADDING;
 static float EXT_FIT_MAX_RADIUS;
 
-bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type) {
+bool psphotInitRadiusEXT (psMetadata *recipe, pmReadout *readout) {
 
     bool status;
@@ -134 +135 @@
     EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS");
 
+    float skyMean  = psMetadataLookupF32 (&status, readout->analysis, "SKY_MEAN");
+    float skyStdev = psMetadataLookupF32 (&status, readout->analysis, "SKY_STDEV");
+
+    fprintf (stderr, "sky: %f +/- %f\n", skyMean, skyStdev);
+
+    EXT_FIT_SKY_SIG = skyStdev;
+
     return true;
 }
 
 // call this function whenever you (re)-define the EXT model
-float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal) {
+bool psphotSetRadiusFootprint (float *radius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float factor) {
 
     psAssert (source, "source not defined??");
@@ -146 +154 @@
 
     // set the radius based on the footprint:
-    if (!peak->footprint) goto escape;
+    if (!peak->footprint) return false;
     pmFootprint *footprint = peak->footprint;
-    if (!footprint->spans) goto escape;
-    if (footprint->spans->n < 1) goto escape;
+    if (!footprint->spans) return false;
+    if (footprint->spans->n < 1) return false;
 
     // find the max radius
-    float radius = 0.0;
+    float rawRadius = 0.0;
     for (int j = 0; j < footprint->spans->n; j++) {
         pmSpan *span = footprint->spans->data[j];
@@ -160 +168 @@
         float dX1 = span->x1 - peak->xf;
 
-        radius = PS_MAX (radius, hypot(dY, dX0));
-        radius = PS_MAX (radius, hypot(dY, dX1));
-    }
-
-    radius += EXT_FIT_PADDING;
-    if (isnan(radius)) psAbort("error in radius");
-
-    radius = PS_MIN (radius, EXT_FIT_MAX_RADIUS);
+        rawRadius = PS_MAX (rawRadius, hypot(dY, dX0));
+        rawRadius = PS_MAX (rawRadius, hypot(dY, dX1));
+    }
+    if (isnan(rawRadius)) return false;
+    rawRadius = PS_MIN (factor*rawRadius + EXT_FIT_PADDING, EXT_FIT_MAX_RADIUS);
 
     // redefine the pixels if needed
-    pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, radius);
-
-    // set the mask to flag the excluded pixels
-    psImageKeepCircle (source->maskObj, peak->xf, peak->yf, radius, "OR", markVal);
-    return radius;
-
-escape:
-    return NAN;
-    // bool result = psphotCheckRadiusEXT (readout, source, model, markVal);
-    // return result;
+    pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, rawRadius);
+
+    // set the mask to flag the excluded pixels
+    psImageKeepCircle (source->maskObj, peak->xf, peak->yf, rawRadius, "OR", markVal);
+
+    *radius = rawRadius;
+    return true;
 }
 
 // alternative EXT radius based on model guess (for use without footprints)
-bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal) {
-
-    psAbort ("do not use this function");
+bool psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep) {
 
     psF32 *PAR = model->params->data.F32;
@@ -193 +193 @@
 
     // set the fit radius based on the object flux limit and the model
-    float rawRadius = model->modelRadius (model->params, EXT_FIT_NSIGMA*moments->dSky);
-
-    model->fitRadius = rawRadius + EXT_FIT_PADDING;
-    if (isnan(model->fitRadius)) psAbort("error in radius");
+    float flux = deep ? EXT_FIT_NSIGMA*EXT_FIT_SKY_SIG : 0.1 * model->params->data.F32[PM_PAR_I0];
+
+    float rawRadius = model->modelRadius (model->params, flux);
+    if (isnan(rawRadius)) return false;
+
+    rawRadius = PS_MIN (rawRadius + EXT_FIT_PADDING, EXT_FIT_MAX_RADIUS);
+    model->fitRadius = rawRadius;
 
     // redefine the pixels if needed
-    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);
+    pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);
 
     // set the mask to flag the excluded pixels
     psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal);
-    return status;
-}
+    return true;
+}

branches/sc_branches/trunkTest/psphot/src/psphotReadout.c

@@ -53 +53 @@
     }
 
-    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are determined and saved on
-    // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT
+    // load the psf model, if suppled.  FWHM_MAJ,FWHM_MIN,etc are determined and saved on
+    // readout->analysis. XXX Note: this function currently only works with a single
+    // PSPHOT.INPUT
     if (!psphotLoadPSF (config, view)) { // ??? need to supply 2 ?
         psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");

branches/sc_branches/trunkTest/psphot/src/psphotSetThreads.c

@@ -15 +15 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_PSF_WEIGHTS", 2);
+    task = psThreadTaskAlloc("PSPHOT_PSF_WEIGHTS", 3);
     task->function = &psphotPSFWeights_Threaded;
     psThreadTaskAdd(task);
@@ -30 +30 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_BLEND_FIT", 9);
+    task = psThreadTaskAlloc("PSPHOT_BLEND_FIT", 10);
     task->function = &psphotBlendFit_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 12);
+    task->function = &psphotExtendedSourceFits_Threaded;
     psThreadTaskAdd(task);
     psFree(task);

branches/sc_branches/trunkTest/psphot/src/psphotSignificanceImage.c

@@ -97 +97 @@
                 smooth_im->data.F32[j][i] = 0.0;
             } else {
-                smooth_im->data.F32[j][i] = factor * PS_SQR(value) / smooth_wt->data.F32[j][i];
+                float v2 = value + PS_SQR(value/1000.0);
+                smooth_im->data.F32[j][i] = factor * PS_SQR(v2) / smooth_wt->data.F32[j][i];
             }
         }

branches/sc_branches/trunkTest/psphot/src/psphotSourceFits.c

@@ -2 +2 @@
 
 // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
-// XXX this function does not call pmSourceFitModelInit : fix this?
 
 static int NfitPSF = 0;
@@ -8 +7 @@
 static int NfitDBL = 0;
 static int NfitEXT = 0;
+static int NfitPCM = 0;
 
 bool psphotFitInit (int nThreads) {
@@ -18 +18 @@
 
     psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n",
-             NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits"));
-    return true;
-}
-
-bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) {
+              NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits"));
+    return true;
+}
+
+bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
 
     float x, y, dR;
+
+    pmSourceFitOptions options = *fitOptions;
 
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -31 +33 @@
     // if this source is not a possible blend, just fit as PSF
     if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
-        bool status = psphotFitPSF (readout, source, psf, maskVal, markVal);
+        bool status = psphotFitPSF (readout, source, psf, fitOptions, maskVal, markVal);
         return status;
     }
@@ -91 +93 @@
 
     // fit PSF model
-    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal);
+    options.mode = PM_SOURCE_FIT_PSF;
+    pmSourceFitSet (source, modelSet, &options, maskVal);
 
     // clear the circular mask
@@ -154 +157 @@
 }
 
-bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) {
+bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
 
     double chiTrend;
+    pmSourceFitOptions options = *fitOptions;
 
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -171 +175 @@
 
     // fit PSF model (set/unset the pixel mask)
-    pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF, maskVal);
+    options.mode = PM_SOURCE_FIT_PSF;
+    pmSourceFitModel (source, PSF, &options, maskVal);
 
     if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
@@ -199 +204 @@
 }
 
+// save a local, static copy of the EXT model type so we don't have to lookup for each object
 static pmModelType modelTypeEXT;
 
-bool psphotInitLimitsEXT (psMetadata *recipe) {
+bool psphotInitLimitsEXT (psMetadata *recipe, pmReadout *readout) {
 
     bool status;
@@ -208 +214 @@
     char *modelNameEXT = psMetadataLookupStr (&status, recipe, "EXT_MODEL");
     modelTypeEXT = pmModelClassGetType (modelNameEXT);
-    psphotInitRadiusEXT (recipe, modelTypeEXT);
-
-    return true;
-}
-
-bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) {
-
+
+    psphotInitRadiusEXT (recipe, readout);
+
+    return true;
+}
+
+bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
+
+    float radius;
     bool okEXT, okDBL;
     float chiEXT, chiDBL;
@@ -222 +230 @@
     pmModel *EXT = NULL;
     psArray *DBL = NULL;
+    pmMoments psfMoments;
 
     // skip the source if we don't think it is extended
+    // XXX are these robust, or do we have better info from the source-size analysis??
     if (source->type == PM_SOURCE_TYPE_UNKNOWN) return false;
     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
@@ -229 +239 @@
 
     // set the radius based on the footprint (also sets the mask pixels)
-    float radius = psphotSetRadiusEXT (readout, source, markVal);
+    if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false;
 
     // XXX note that this changes the source moments that are published...
+    // XXX all published moments should use the same measurement
     // 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?
-    if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false;
+    // XXX save the psf-based moments for output
+    psfMoments = *source->moments;
+    if (!pmSourceMoments (source, radius, 0.0, 0.5, maskVal)) {
+      *source->moments = psfMoments;
+      return false;
+    }
 
     psTrace ("psphot", 5, "trying blob...\n");
@@ -244 +260 @@
     // this temporary source is used as a place-holder by the psphotEval functions below
     tmpSrc = pmSourceAlloc ();
-
-    // XXX need to handle failures better here
-    EXT = psphotFitEXT (readout, source, modelTypeEXT, maskVal, markVal);
-    if (!EXT) goto escape;
-    if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape;
-
-    okEXT = psphotEvalEXT (tmpSrc, EXT);
-    chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN;
-
-    // DBL will always be defined, but DBL->data[n] might not
-    DBL = psphotFitDBL (readout, source, maskVal, markVal);
-    if (!DBL) goto escape;
-    if (!DBL->n) goto escape;
-
-    okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
-    okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
-    // XXX should I keep / save the flags set in the eval functions?
+    {
+        // DBL will always be defined, but DBL->data[n] might not
+        DBL = psphotFitDBL (readout, source, fitOptions, maskVal, markVal);
+        if (!DBL) goto escape;
+        if (!DBL->n) goto escape;
+
+        okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
+        okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
+        // XXX should I keep / save the flags set in the eval functions?
+
+        // correct first model chisqs for flux trend
+        chiDBL = NAN;
+        ONE = DBL->data[0];
+        if (ONE) {
+            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
+            chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
+            ONE->chisqNorm = ONE->chisq / chiTrend;
+            chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison
+            ONE->fitRadius = radius;
+        }
+
+        // correct second model chisqs for flux trend
+        ONE = DBL->data[1];
+        if (ONE) {
+            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
+            chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
+            ONE->chisqNorm = ONE->chisq / chiTrend;
+            ONE->fitRadius = radius;
+        }
+    }
+
+    { 
+        // XXX need to handle failures better here
+        EXT = psphotFitEXT (readout, source, fitOptions, modelTypeEXT, maskVal, markVal);
+        if (!EXT) goto escape;
+        if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape;
+
+        okEXT = psphotEvalEXT (tmpSrc, EXT);
+        chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN;
+    }
 
     // clear the circular mask
     psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
-
-    // correct first model chisqs for flux trend
-    chiDBL = NAN;
-    ONE = DBL->data[0];
-    if (ONE) {
-        if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
-      chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
-      ONE->chisqNorm = ONE->chisq / chiTrend;
-      chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison
-    }
-
-    // correct second model chisqs for flux trend
-    ONE = DBL->data[1];
-    if (ONE) {
-        if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
-      chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
-      ONE->chisqNorm = ONE->chisq / chiTrend;
-    }
 
     psFree (tmpSrc);
@@ -298 +319 @@
 
     // both models failed; reject them both
-    // XXX -- change type flags to psf in this case and keep original moments?
+    // XXX -- change type flags to psf in this case, and make sure we subtract it?
+    // reset the psf moments
+    *source->moments = psfMoments;
+
     psFree (EXT);
     psFree (DBL);
@@ -308 +332 @@
 
     // save new model
+    // XXX save the correct radius...
     source->modelEXT = EXT;
-    source->modelEXT->fitRadius = radius;
     source->type = PM_SOURCE_TYPE_EXTENDED;
     source->mode |= PM_SOURCE_MODE_EXTMODEL;
@@ -327 +351 @@
 # endif
 
+    // reset the psf moments
+    *source->moments = psfMoments;
     return true;
 
@@ -337 +363 @@
     source->modelPSF = psMemIncrRefCounter (DBL->data[0]);
     source->mode     |= PM_SOURCE_MODE_PAIR;
-    source->modelPSF->fitRadius = radius;
 
     // copy most data from the primary source (modelEXT, blends stay NULL)
-    pmSource *newSrc = pmSourceCopy (source);
+    pmSource *newSrc = pmSourceCopyData (source);
     newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
-    newSrc->modelPSF->fitRadius = radius;
 
     // build cached models and subtract
@@ -365 +389 @@
 # endif
 
+    // reset the (original) psf moments
+    *source->moments = psfMoments;
+    *newSrc->moments = psfMoments;
+
     psArrayAdd (newSources, 100, newSrc);
     psFree (newSrc);
@@ -371 +399 @@
 
 escape:
+    // reset the psf moments
+    *source->moments = psfMoments;
     psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
     psFree (tmpSrc);
@@ -379 +409 @@
 
 // fit a double PSF source to an extended blob
-psArray *psphotFitDBL (pmReadout *readout, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) {
+psArray *psphotFitDBL (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
 
     float dx, dy;
@@ -387 +417 @@
     psEllipseMoments moments;
     psArray *modelSet;
+
+    pmSourceFitOptions options = *fitOptions;
 
     NfitDBL ++;
@@ -403 +435 @@
     axes = psEllipseMomentsToAxes (moments, 20.0);
 
+    if (isnan(axes.major)) return NULL;
+    if (isnan(axes.minor)) return NULL;
+    if (isnan(axes.theta)) return NULL;
+
     // XXX this is really arbitrary: 4 pixel separation?
     dx = 2 * cos (axes.theta);
@@ -426 +462 @@
 
     // fit PSF model (set/unset the pixel mask)
-    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal);
+    options.mode = PM_SOURCE_FIT_PSF;
+    pmSourceFitSet (source, modelSet, &options, maskVal);
     return (modelSet);
 }
 
-pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {
+pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {
+
+    if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
+        psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
+    }
+
+    pmSourceFitOptions options = *fitOptions;
 
     NfitEXT ++;
@@ -438 +481 @@
 
     // use the source moments, etc to guess basic model parameters
-    pmModel *EXT = pmSourceModelGuess (source, modelType);
-    if (!EXT) {
+    pmModel *model = pmSourceModelGuess (source, modelType);
+    if (!model) {
         psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
         return NULL;
     }
 
+    // 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);
+    }
+
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        options.mode = PM_SOURCE_FIT_NO_INDEX;
+    } else {
+        options.mode = PM_SOURCE_FIT_EXT;
+    }
+
+    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
+    pmSourceFitModel (source, model, &options, maskVal);
+    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+
+    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+    return (model);
+}
+
+pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
     if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
         psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
     }
 
+    pmSourceFitOptions options = *fitOptions;
+
+    NfitPCM ++;
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    // allocate the model
+    pmModel *model = pmModelAlloc(modelType);
+    if (!model) {
+        return NULL;
+    }
+
+    pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
+    if (!pcm) {
+        psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
+        model->flags |= PM_MODEL_STATUS_BADARGS; // XXX this is probably already set in pmPCMinit
+        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 (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)) {
+            model->flags |= PM_MODEL_STATUS_BADARGS;
+            return model;
+        }
+    }
+
+    if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
+        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
+    }
+
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        options.mode = PM_SOURCE_FIT_NO_INDEX;
+    } else {
+        options.mode = PM_SOURCE_FIT_EXT;
+    }
+    // update the pcm elements if we have changed the circumstance (options.mode or source->pixels)
+    pmPCMupdate(pcm, source, &options, model);
+
+    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
+    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+
+    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+    psFree (pcm);
+
+    return model;
+}
+
+// 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
+
+// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
+// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
+bool psphotFitSersicIndex (pmModel *model, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
+
+    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
+
+    pmSourceFitOptions options = *fitOptions;
+    
     // fit EXT (not PSF) model (set/unset the pixel mask)
-    pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT, maskVal);
-    return (EXT);
-}
+    options.mode = PM_SOURCE_FIT_NO_INDEX;
+    options.nIter = 4;
+
+    int iMin = -1;
+    float xMin = NAN;
+    float chiSquare[N_INDEX_GUESS];
+
+    for (int i = 0; i < N_INDEX_GUESS; i++) {
+        model->params->data.F32[PM_PAR_7] = indexGuess[i];
+
+        if (!model->modelGuess(model, source)) {
+            model->flags |= PM_MODEL_STATUS_BADARGS;
+            return false;
+        }
+
+        // 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);
+
+        chiSquare[i] = model->chisqNorm;
+        if (i == 0) {
+            xMin = chiSquare[i];
+            iMin = i;
+        } else {
+            if (chiSquare[i] < xMin) {
+                xMin = chiSquare[i];
+                iMin = i;
+            }
+        }
+    }
+    assert (iMin >= 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->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;
+}
+
+// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
+// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
+bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
+    pmModel *model = pcm->modelConv;
+
+    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
+
+    pmSourceFitOptions options = *fitOptions;
+    
+    // fit EXT (not PSF) model (set/unset the pixel mask)
+    options.mode = PM_SOURCE_FIT_NO_INDEX;
+    options.nIter = 4;
+
+    // update the pcm elements if we have changed the circumstance (here, options.mode)
+    pmPCMupdate(pcm, source, &options, model);
+
+    int iMin = -1;
+    float xMin = NAN;
+    float chiSquare[N_INDEX_GUESS];
+
+    for (int i = 0; i < N_INDEX_GUESS; i++) {
+        model->params->data.F32[PM_PAR_7] = indexGuess[i];
+        
+        if (!model->modelGuess(model, source)) {
+            model->flags |= PM_MODEL_STATUS_BADARGS;
+            return false;
+        }
+
+        // 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)
+        pmSourceFitModel (source, model, &options, maskVal);
+# else
+        pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
+        pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+# endif
+        // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+
+        chiSquare[i] = model->chisq;
+        if (i == 0) {
+            xMin = chiSquare[i];
+            iMin = i;
+        } else {
+            if (chiSquare[i] < xMin) {
+                xMin = chiSquare[i];
+                iMin = i;
+            }
+        }
+    }
+    assert (iMin >= 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];
+
+    pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
+
+    return true;
+}

branches/sc_branches/trunkTest/psphot/src/psphotSourceSize.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 # include <gsl/gsl_sf_gamma.h>
+
+# define KRON 1
 
 typedef struct {
@@ -168 +170 @@
     psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
 
+    psImageMaskType markVal = options->markVal;
     psImageMaskType maskVal = options->maskVal | options->markVal;
 
@@ -189 +192 @@
 
         // XXX can we test if psfMag is set and calculate only if needed?
-        pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
+        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
 
         // clear the mask bit
@@ -197 +200 @@
         pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
 
+        // XXX test: switch to kron flux
+# if (KRON)
+        float apMag = -2.5*log10(source->moments->KronFlux);
+# else
         float apMag = -2.5*log10(source->moments->Sum);
+# endif
         float dMag = source->psfMag - apMag;
 
@@ -288 +296 @@
     pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
 
+    psImageMaskType markVal = options->markVal;
     psImageMaskType maskVal = options->maskVal | options->markVal;
 
@@ -335 +344 @@
 
         // XXX can we test if psfMag is set and calculate only if needed?
-        pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
+        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
 
         // clear the mask bit
@@ -343 +352 @@
         pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
 
+# if (KRON)
+        float apMag = -2.5*log10(source->moments->KronFlux);
+# else
         float apMag = -2.5*log10(source->moments->Sum);
+# endif
         float dMag = source->psfMag - apMag;
 
         // set nSigma to include both systematic and poisson error terms
         // XXX the 'poisson error' contribution for size is probably wrong...
-        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
+        // XXX add in a hard floor on the Ap Sys Err (to be a bit generous)
+        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.025));
         float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag);
         float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag);
 
-        // partially-masked sources are more likely to be mis-measured PSFs
+        // fprintf (stderr, "%f %f : Mxx: %f, Myy: %f, dx: %f, dy: %f, psfMag: %f, apMag: %f, dMag: %f, errMag: %f, nSigmaMag: %f, nSigmaMxx: %f, nSigmaMyy: %f\n", 
+        // source->peak->xf, source->peak->yf, Mxx, Myy, source->peak->xf - source->moments->Mx, source->peak->yf - source->moments->My, 
+        // source->psfMag, apMag, dMag, source->errMag, nSigmaMAG, nSigmaMXX, nSigmaMYY);
+
+        // XXX double check on ths stuff!! partially-masked sources are more likely to be mis-measured PSFs
         float sizeBias = 1.0;
-        if (source->pixWeight < 0.9) {
+        if (source->pixWeightNotBad < 0.9) {
+            sizeBias = 3.0;
+        }
+        if (source->pixWeightNotPoor < 0.9) {
             sizeBias = 3.0;
         }
@@ -389 +410 @@
         if (isCR) {
             psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g %f\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
-                    source->peak->xf,source->peak->yf,source->pixWeight,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
+                    source->peak->xf,source->peak->yf,source->pixWeightNotBad,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
                     options->nSigmaApResid,sizeBias*options->nSigmaMoments);
             source->mode |= PM_SOURCE_MODE_DEFECT;

branches/sc_branches/trunkTest/psphot/src/psphotSourceStats.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+void pmSourceMomentsSetVerbose(bool state);
 
 // convert detections to sources and measure their basic properties (moments, local sky, sky
@@ -369 +370 @@
     maskVal |= markVal;
 
+    // XXX test : pmSourceMomentsSetVerbose(true);
+
     // threaded measurement of the sources moments
     int Nfail = 0;
@@ -408 +411 @@
         // measure basic source moments (no S/N clipping on input pixels)
         status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal);
+        // XXX moments / aperture test:
+        if (0) {
+            // clear the mask bit and set the circular mask pixels
+            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+            psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, 4.0*source->moments->Mrf, "OR", markVal);
+            
+            float apMag = NAN;
+            pmSourcePhotometryAper (&apMag, NULL, source->pixels, source->maskObj, maskVal);
+            fprintf (stderr, "apMag: %f, kronMag: %f\n", apMag, -2.5*log10(source->moments->KronFlux));
+
+            // clear the mask bit
+            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+        }           
         if (status) {
             Nmoments ++;
@@ -502 +518 @@
         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
-        psphotVisualPlotMoments (recipe, analysis, sources);
+        if (pmVisualTestLevel("psphot.moments.full", 2)) {
+            psphotVisualPlotMoments (recipe, analysis, sources);
+        }
 #endif
 

branches/sc_branches/trunkTest/psphot/src/psphotStack.c

@@ -24 +24 @@
         exit (psphotGetExitStatus());
     }
+    // Dump configuration
+    bool mdok;
+    psString dump = psMetadataLookupStr(&mdok, config->arguments, "DUMP_CONFIG"); // File for config
+    if (dump && !pmConfigDump(config, dump)) {
+        psError(psErrorCodeLast(), false, "Unable to dump configuration.");
+        exit (psphotGetExitStatus());
+    }
+
+
 
     psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete"));

branches/sc_branches/trunkTest/psphot/src/psphotStackMatchPSFsUtils.c

@@ -346 +346 @@
     if (!fake) goto escape;
 
-    dumpImage(fake, readoutSrc, index, "fake");
-    dumpImage(readoutSrc,  readoutSrc, index, "real");
+    // dumpImage(fake, readoutSrc, index, "fake");
+    // dumpImage(readoutSrc,  readoutSrc, index, "real");
 
     if (threads) pmSubtractionThreadsInit();
@@ -380 +380 @@
     }
 
-    dumpImage(readoutOut, readoutSrc, index, "conv");
-    dumpImageDiff(readoutOut, fake, readoutSrc, index, "diff");
+    // dumpImage(readoutOut, readoutSrc, index, "conv");
+    // dumpImageDiff(readoutOut, fake, readoutSrc, index, "diff");
 
     psFree(fake);

branches/sc_branches/trunkTest/psphot/src/psphotVisual.c

@@ -66 +66 @@
     int myKapa = psphotKapaChannel (channel);
     if (!(strcasecmp (overlay, "all"))) {
-      KiiEraseOverlay (myKapa, "red");
-      KiiEraseOverlay (myKapa, "green");
-      KiiEraseOverlay (myKapa, "blue");
-      KiiEraseOverlay (myKapa, "yellow");
-      return true;
+        KiiEraseOverlay (myKapa, "red");
+        KiiEraseOverlay (myKapa, "green");
+        KiiEraseOverlay (myKapa, "blue");
+        KiiEraseOverlay (myKapa, "yellow");
+        return true;
     }
     KiiEraseOverlay (myKapa, overlay);
@@ -182 +182 @@
 bool psphotVisualShowImage (pmReadout *readout) {
 
-    if (!pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.image", 1)) return true;
 
     int kapa = psphotKapaChannel (1);
@@ -199 +199 @@
     pmReadout *backgnd;
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.image.backgnd", 2)) return true;
 
     int kapa = psphotKapaChannel (1);
@@ -208 +208 @@
 
     if (file->mode == PM_FPA_MODE_INTERNAL) {
-        backgnd = file->readout;
+        backgnd = file->readout;
     } else {
-        backgnd = pmFPAviewThisReadout (view, file->fpa);
+        backgnd = pmFPAviewThisReadout (view, file->fpa);
     }
 
@@ -222 +222 @@
 bool psphotVisualShowSignificance (psImage *image, float min, float max) {
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.image.signif", 2)) return true;
 
     int kapa = psphotKapaChannel (1);
@@ -233 +233 @@
 }
 
+bool psphotVisualShowLogSignificance (psImage *image, float min, float max) {
+
+    if (!pmVisualTestLevel("psphot.image.logsignif", 3)) return true;
+
+    int kapa = psphotKapaChannel (1);
+    if (kapa == -1) return false;
+
+    psImage *lsig = (psImage *) psUnaryOp (NULL, image, "log");
+    psphotVisualRangeImage (kapa, lsig, "log-signif", 2, min, max);
+    psFree (lsig);
+
+    pmVisualAskUser(NULL);
+    return true;
+}
+
+// XXX : requires psphotVisualShowImage
 bool psphotVisualShowPeaks (pmDetections *detections) {
 
@@ -238 +254 @@
     KiiOverlay *overlay;
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.objects.peaks", 1)) return true;
 
     int kapa = psphotKapaChannel (1);
@@ -252 +268 @@
     for (int i = 0; i < peaks->n; i++) {
 
-        pmPeak *peak = peaks->data[i];
-        if (peak == NULL) continue;
-
-        overlay[Noverlay].type = KII_OVERLAY_BOX;
-        overlay[Noverlay].x = peak->xf;
-        overlay[Noverlay].y = peak->yf;
-        overlay[Noverlay].dx = 2.0;
-        overlay[Noverlay].dy = 2.0;
-        overlay[Noverlay].angle = 0.0;
-        overlay[Noverlay].text = NULL;
-        Noverlay ++;
+        pmPeak *peak = peaks->data[i];
+        if (peak == NULL) continue;
+
+        overlay[Noverlay].type = KII_OVERLAY_BOX;
+        overlay[Noverlay].x = peak->xf;
+        overlay[Noverlay].y = peak->yf;
+        overlay[Noverlay].dx = 2.0;
+        overlay[Noverlay].dy = 2.0;
+        overlay[Noverlay].angle = 0.0;
+        overlay[Noverlay].text = NULL;
+        Noverlay ++;
     }
 
@@ -272 +288 @@
 }
 
+// XXX : requires psphotVisualShowImage
 bool psphotVisualShowFootprints (pmDetections *detections) {
 
@@ -277 +294 @@
     KiiOverlay *overlay;
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.objects.footprints", 3)) return true;
 
     int kapa = psphotKapaChannel (1);
@@ -292 +309 @@
     for (int i = 0; i < footprints->n; i++) {
 
-        pmSpan *span = NULL;
-
-        pmFootprint *footprint = footprints->data[i];
-        if (footprint == NULL) continue;
-        if (footprint->spans == NULL) continue;
-        if (footprint->spans->n < 1) continue;
-
-        // draw the top
-        // XXX need to allow top (and bottom) to have more than one span
-        span = footprint->spans->data[0];
-        overlay[Noverlay].type = KII_OVERLAY_LINE;
-        overlay[Noverlay].x = span->x0;
-        overlay[Noverlay].y = span->y;
-        overlay[Noverlay].dx = span->x1 - span->x0;
-        overlay[Noverlay].dy = 0;
-        overlay[Noverlay].angle = 0.0;
-        overlay[Noverlay].text = NULL;
-        Noverlay ++;
-        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
-
-        int ys = span->y;
-        int x0s = span->x0;
-        int x1s = span->x1;
-
-        // draw the outer span edges
-        for (int j = 1; j < footprint->spans->n; j++) {
-            pmSpan *span1 = footprint->spans->data[j];
-
-            int ye = span1->y;
-            int x0e = span1->x0;
-            int x1e = span1->x1;
-
-            // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right)
-            // find all of the spans in this row and generate x0e, x01:
-            for (int k = j + 1; k < footprint->spans->n; k++) {
-                pmSpan *span2 = footprint->spans->data[k];
-                if (span2->y > span1->y) break;
-                x0e = PS_MIN (x0e, span2->x0);
-                x1e = PS_MAX (x1e, span2->x1);
-                j++;
-            }
-
-            overlay[Noverlay].type = KII_OVERLAY_LINE;
-            overlay[Noverlay].x = x0s;
-            overlay[Noverlay].y = ys;
-            overlay[Noverlay].dx = x0e - x0s;
-            overlay[Noverlay].dy = ye - ys;
-            overlay[Noverlay].angle = 0.0;
-            overlay[Noverlay].text = NULL;
-            Noverlay ++;
-            CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
-
-            overlay[Noverlay].type = KII_OVERLAY_LINE;
-            overlay[Noverlay].x = x1s;
-            overlay[Noverlay].y = ys;
-            overlay[Noverlay].dx = x1e - x1s;
-            overlay[Noverlay].dy = ye - ys;
-            overlay[Noverlay].angle = 0.0;
-            overlay[Noverlay].text = NULL;
-            Noverlay ++;
-            CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
-
-            ys = ye;
-            x0s = x0e;
-            x1s = x1e;
-        }
-
-        // draw the bottom
-        span = footprint->spans->data[footprint->spans->n - 1];
-        overlay[Noverlay].type = KII_OVERLAY_LINE;
-        overlay[Noverlay].x = span->x0;
-        overlay[Noverlay].y = span->y;
-        overlay[Noverlay].dx = span->x1 - span->x0;
-        overlay[Noverlay].dy = 0;
-        overlay[Noverlay].angle = 0.0;
-        overlay[Noverlay].text = NULL;
-        Noverlay ++;
-        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
+        pmSpan *span = NULL;
+
+        pmFootprint *footprint = footprints->data[i];
+        if (footprint == NULL) continue;
+        if (footprint->spans == NULL) continue;
+        if (footprint->spans->n < 1) continue;
+
+        // draw the top
+        // XXX need to allow top (and bottom) to have more than one span
+        span = footprint->spans->data[0];
+        overlay[Noverlay].type = KII_OVERLAY_LINE;
+        overlay[Noverlay].x = span->x0;
+        overlay[Noverlay].y = span->y;
+        overlay[Noverlay].dx = span->x1 - span->x0;
+        overlay[Noverlay].dy = 0;
+        overlay[Noverlay].angle = 0.0;
+        overlay[Noverlay].text = NULL;
+        Noverlay ++;
+        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
+
+        int ys = span->y;
+        int x0s = span->x0;
+        int x1s = span->x1;
+
+        // draw the outer span edges
+        for (int j = 1; j < footprint->spans->n; j++) {
+            pmSpan *span1 = footprint->spans->data[j];
+
+            int ye = span1->y;
+            int x0e = span1->x0;
+            int x1e = span1->x1;
+
+            // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right)
+            // find all of the spans in this row and generate x0e, x01:
+            for (int k = j + 1; k < footprint->spans->n; k++) {
+                pmSpan *span2 = footprint->spans->data[k];
+                if (span2->y > span1->y) break;
+                x0e = PS_MIN (x0e, span2->x0);
+                x1e = PS_MAX (x1e, span2->x1);
+                j++;
+            }
+
+            overlay[Noverlay].type = KII_OVERLAY_LINE;
+            overlay[Noverlay].x = x0s;
+            overlay[Noverlay].y = ys;
+            overlay[Noverlay].dx = x0e - x0s;
+            overlay[Noverlay].dy = ye - ys;
+            overlay[Noverlay].angle = 0.0;
+            overlay[Noverlay].text = NULL;
+            Noverlay ++;
+            CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
+
+            overlay[Noverlay].type = KII_OVERLAY_LINE;
+            overlay[Noverlay].x = x1s;
+            overlay[Noverlay].y = ys;
+            overlay[Noverlay].dx = x1e - x1s;
+            overlay[Noverlay].dy = ye - ys;
+            overlay[Noverlay].angle = 0.0;
+            overlay[Noverlay].text = NULL;
+            Noverlay ++;
+            CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
+
+            ys = ye;
+            x0s = x0e;
+            x1s = x1e;
+        }
+
+        // draw the bottom
+        span = footprint->spans->data[footprint->spans->n - 1];
+        overlay[Noverlay].type = KII_OVERLAY_LINE;
+        overlay[Noverlay].x = span->x0;
+        overlay[Noverlay].y = span->y;
+        overlay[Noverlay].dx = span->x1 - span->x0;
+        overlay[Noverlay].dy = 0;
+        overlay[Noverlay].angle = 0.0;
+        overlay[Noverlay].text = NULL;
+        Noverlay ++;
+        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     }
 
@@ -379 +396 @@
 }
 
+// XXX : requires psphotVisualShowImage
 bool psphotVisualShowMoments (psArray *sources) {
 
@@ -387 +405 @@
     psEllipseAxes axes;
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.objects.moments", 2)) return true;
 
     int kapa = psphotKapaChannel (1);
@@ -401 +419 @@
     for (int i = 0; i < sources->n; i++) {
 
-        pmSource *source = sources->data[i];
-        if (source == NULL) continue;
-
-        pmMoments *moments = source->moments;
-        if (moments == NULL) continue;
-
-        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
-        overlay[Noverlay].x = moments->Mx;
-        overlay[Noverlay].y = moments->My;
-
-        emoments.x2 = moments->Mxx;
-        emoments.xy = moments->Mxy;
-        emoments.y2 = moments->Myy;
-
-        axes = psEllipseMomentsToAxes (emoments, 20.0);
-
-        overlay[Noverlay].dx = 2.0*axes.major;
-        overlay[Noverlay].dy = 2.0*axes.minor;
-
-        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
-
-        overlay[Noverlay].text = NULL;
-        Noverlay ++;
+        pmSource *source = sources->data[i];
+        if (source == NULL) continue;
+
+        pmMoments *moments = source->moments;
+        if (moments == NULL) continue;
+
+        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
+        overlay[Noverlay].x = moments->Mx;
+        overlay[Noverlay].y = moments->My;
+
+        emoments.x2 = moments->Mxx;
+        emoments.xy = moments->Mxy;
+        emoments.y2 = moments->Myy;
+
+        axes = psEllipseMomentsToAxes (emoments, 20.0);
+
+        overlay[Noverlay].dx = 2.0*axes.major;
+        overlay[Noverlay].dy = 2.0*axes.minor;
+
+        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
+
+        overlay[Noverlay].text = NULL;
+        Noverlay ++;
     }
 
@@ -439 +457 @@
     KapaSection section;
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.moments", 1)) return true;
 
     int myKapa = psphotKapaChannel (2);
@@ -456 +474 @@
     float Ymin = 1000.0, Ymax = 0.0;
     {
-        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
-        for (int n = 0; n < nRegions; n++) {
-
-            char regionName[64];
-            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
-            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
-
-            float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
-            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
-            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
-            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
-
-            float X0 = psfX - 4.0*psfdX;
-            float X1 = psfX + 4.0*psfdX;
-            float Y0 = psfY - 4.0*psfdY;
-            float Y1 = psfY + 4.0*psfdY;
-
-            if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
-            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
-            if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
-            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
-        }
+        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
+        for (int n = 0; n < nRegions; n++) {
+
+            char regionName[64];
+            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
+            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
+
+            float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
+            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
+            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
+            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
+
+            float X0 = psfX - 4.0*psfdX;
+            float X1 = psfX + 4.0*psfdX;
+            float Y0 = psfY - 4.0*psfdY;
+            float Y1 = psfY + 4.0*psfdY;
+
+            if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
+            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
+            if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
+            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
+        }
     }
     Xmin = PS_MAX(Xmin, -0.1);
@@ -498 +516 @@
     int nF = 0;
     for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-        if (source->moments == NULL)
-            continue;
-
-        xFaint->data.F32[nF] = source->moments->Mxx;
-        yFaint->data.F32[nF] = source->moments->Myy;
-        mFaint->data.F32[nF] = -2.5*log10(source->moments->Sum);
-        nF++;
-
-        // XXX make this a user-defined cutoff
-        if (source->moments->SN < SN_LIM)
-            continue;
-
-        xBright->data.F32[nB] = source->moments->Mxx;
-        yBright->data.F32[nB] = source->moments->Myy;
-        mBright->data.F32[nB] = -2.5*log10(source->moments->Sum);
-        nB++;
+        pmSource *source = sources->data[i];
+        if (source->moments == NULL)
+            continue;
+
+        xFaint->data.F32[nF] = source->moments->Mxx;
+        yFaint->data.F32[nF] = source->moments->Myy;
+        mFaint->data.F32[nF] = -2.5*log10(source->moments->Sum);
+        nF++;
+
+        // XXX make this a user-defined cutoff
+        if (source->moments->SN < SN_LIM)
+            continue;
+
+        xBright->data.F32[nB] = source->moments->Mxx;
+        yBright->data.F32[nB] = source->moments->Myy;
+        mBright->data.F32[nB] = -2.5*log10(source->moments->Sum);
+        nB++;
     }
     xFaint->n = nF;
@@ -652 +670 @@
     // draw N circles to outline the clumps
     {
-        KapaSelectSection (myKapa, "MxxMyy");
-
-        // draw a circle centered on psfX,Y with size of the psf limit
-        psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
-        psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
-
-        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
-        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
-
-        graphdata.color = KapaColorByName ("blue");
-        graphdata.style = 0;
-
-        graphdata.xmin = Xmin;
-        graphdata.ymin = Ymin;
-        graphdata.xmax = Xmax;
-        graphdata.ymax = Ymax;
-        KapaSetLimits (myKapa, &graphdata);
-
-        for (int n = 0; n < nRegions; n++) {
-
-            char regionName[64];
-            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
-            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
-
-            float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
-            float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
-            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
-            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
-            float Rx = psfdX * PSF_CLUMP_NSIGMA;
-            float Ry = psfdY * PSF_CLUMP_NSIGMA;
-
-            for (int i = 0; i < xLimit->n; i++) {
-                xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX;
-                yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
-            }
-            KapaPrepPlot (myKapa, xLimit->n, &graphdata);
-            KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x");
-            KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y");
-        }
-        psFree (xLimit);
-        psFree (yLimit);
+        KapaSelectSection (myKapa, "MxxMyy");
+
+        // draw a circle centered on psfX,Y with size of the psf limit
+        psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
+        psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
+
+        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
+        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
+
+        graphdata.color = KapaColorByName ("blue");
+        graphdata.style = 0;
+
+        graphdata.xmin = Xmin;
+        graphdata.ymin = Ymin;
+        graphdata.xmax = Xmax;
+        graphdata.ymax = Ymax;
+        KapaSetLimits (myKapa, &graphdata);
+
+        for (int n = 0; n < nRegions; n++) {
+
+            char regionName[64];
+            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
+            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
+
+            float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
+            float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
+            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
+            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
+            float Rx = psfdX * PSF_CLUMP_NSIGMA;
+            float Ry = psfdY * PSF_CLUMP_NSIGMA;
+
+            for (int i = 0; i < xLimit->n; i++) {
+                xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX;
+                yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
+            }
+            KapaPrepPlot (myKapa, xLimit->n, &graphdata);
+            KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x");
+            KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y");
+        }
+        psFree (xLimit);
+        psFree (yLimit);
     }
 
@@ -721 +739 @@
     for (int i = 0; i < sources->n; i++) {
 
-        pmSource *source = sources->data[i];
-        if (source == NULL) continue;
-
-        if (source->type != type) continue;
-        if (mode && !(source->mode & mode)) continue;
-
-        pmMoments *moments = source->moments;
-        if (moments == NULL) continue;
-
-        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
-        overlay[Noverlay].x = moments->Mx;
-        overlay[Noverlay].y = moments->My;
-
-        emoments.x2 = moments->Mxx;
-        emoments.y2 = moments->Myy;
-        emoments.xy = moments->Mxy;
-
-        axes = psEllipseMomentsToAxes (emoments, 20.0);
-
-        overlay[Noverlay].dx = 2.0*axes.major;
-        overlay[Noverlay].dy = 2.0*axes.minor;
-        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
-        overlay[Noverlay].text = NULL;
-        Noverlay ++;
+        pmSource *source = sources->data[i];
+        if (source == NULL) continue;
+
+        if (source->type != type) continue;
+        if (mode && !(source->mode & mode)) continue;
+
+        pmMoments *moments = source->moments;
+        if (moments == NULL) continue;
+
+        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
+        overlay[Noverlay].x = moments->Mx;
+        overlay[Noverlay].y = moments->My;
+
+        emoments.x2 = moments->Mxx;
+        emoments.y2 = moments->Myy;
+        emoments.xy = moments->Mxy;
+
+        axes = psEllipseMomentsToAxes (emoments, 20.0);
+
+        overlay[Noverlay].dx = 2.0*axes.major;
+        overlay[Noverlay].dy = 2.0*axes.minor;
+        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
+        overlay[Noverlay].text = NULL;
+        Noverlay ++;
     }
 
@@ -753 +771 @@
 }
 
+// XXX : requires psphotVisualShowImage
 bool psphotVisualShowRoughClass (psArray *sources) {
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.objects.size", 3)) return true;
 
     int myKapa = psphotKapaChannel (1);
@@ -776 +795 @@
 bool psphotVisualShowPSFModel (pmReadout *readout, pmPSF *psf) {
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.psf.model", 1)) return true;
 
     int myKapa = psphotKapaChannel (3);
@@ -797 +816 @@
     // generate a fake model at each of the 3x3 image grid positions
     for (int x = -2; x <= +2; x ++) {
-        for (int y = -2; y <= +2; y ++) {
-            // use the center of the center pixel of the image
-            float xc = (int)((0.5 + 0.225*x)*readout->image->numCols) + readout->image->col0 + 0.5;
-            float yc = (int)((0.5 + 0.225*y)*readout->image->numRows) + readout->image->row0 + 0.5;
-
-            // assign the x and y coords to the image center
-            // create an object with center intensity of 1000
-            modelRef->params->data.F32[PM_PAR_SKY] = 0;
-            modelRef->params->data.F32[PM_PAR_I0] = 1000;
-            modelRef->params->data.F32[PM_PAR_XPOS] = xc;
-            modelRef->params->data.F32[PM_PAR_YPOS] = yc;
-
-            // create modelPSF from this model
-            pmModel *model = pmModelFromPSF (modelRef, psf);
-            if (!model) continue;
-
-            // place the reference object in the image center
-            // no need to mask the source here
-            // XXX should we measure this for the analytical model only or the full model?
-            pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
-            pmModelAddWithOffset (funMosaic, NULL, model, PM_MODEL_OP_FUNC | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
-            pmModelAddWithOffset (resMosaic, NULL, model, PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1 | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
-            psFree (model);
-        }
+        for (int y = -2; y <= +2; y ++) {
+            // use the center of the center pixel of the image
+            float xc = (int)((0.5 + 0.225*x)*readout->image->numCols) + readout->image->col0 + 0.5;
+            float yc = (int)((0.5 + 0.225*y)*readout->image->numRows) + readout->image->row0 + 0.5;
+
+            // assign the x and y coords to the image center
+            // create an object with center intensity of 1000
+            modelRef->params->data.F32[PM_PAR_SKY] = 0;
+            modelRef->params->data.F32[PM_PAR_I0] = 1000;
+            modelRef->params->data.F32[PM_PAR_XPOS] = xc;
+            modelRef->params->data.F32[PM_PAR_YPOS] = yc;
+
+            // create modelPSF from this model
+            pmModel *model = pmModelFromPSF (modelRef, psf);
+            if (!model) continue;
+
+            // place the reference object in the image center
+            // no need to mask the source here
+            // XXX should we measure this for the analytical model only or the full model?
+            pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
+            pmModelAddWithOffset (funMosaic, NULL, model, PM_MODEL_OP_FUNC | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
+            pmModelAddWithOffset (resMosaic, NULL, model, PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1 | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
+            psFree (model);
+        }
     }
 
@@ -842 +861 @@
     bool status;
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.psf.stars", 2)) return true;
 
     int myKapa = psphotKapaChannel (3);
@@ -871 +890 @@
     for (int i = 0; i < sources->n; i++) {
 
-        pmSource *source = sources->data[i];
-
-        bool keep = false;
-        keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
-        if (!keep) continue;
-
-        // how does this subimage get placed into the output image?
-        // DX = source->pixels->numCols
-        // DY = source->pixels->numRows
-
-        if (dX + DX > NX) {
-            // too wide for the rest of this row
-            if (dX == 0) {
-                // alone on this row
-                NY += DY;
-                dX = 0;
-                dY = 0;
-            } else {
-                // start the next row
-                NY += dY;
-                dX = DX;
-                dY = DY;
-            }
-        } else {
-            // extend this row
-            dX += DX;
-            dY = PS_MAX (dY, DY);
-        }
+        pmSource *source = sources->data[i];
+
+        bool keep = false;
+        keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
+        if (!keep) continue;
+
+        // how does this subimage get placed into the output image?
+        // DX = source->pixels->numCols
+        // DY = source->pixels->numRows
+
+        if (dX + DX > NX) {
+            // too wide for the rest of this row
+            if (dX == 0) {
+                // alone on this row
+                NY += DY;
+                dX = 0;
+                dY = 0;
+            } else {
+                // start the next row
+                NY += dY;
+                dX = DX;
+                dY = DY;
+            }
+        } else {
+            // extend this row
+            dX += DX;
+            dY = PS_MAX (dY, DY);
+        }
     }
     NY += DY;
@@ -918 +937 @@
     for (int i = 0; i < sources->n; i++) {
 
-        pmSource *source = sources->data[i];
-
-        bool keep = false;
-        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
-            nPSF ++;
-            keep = true;
-        }
-        if (!keep) continue;
-
-        if (Xo + DX > NX) {
-            // too wide for the rest of this row
-            if (Xo == 0) {
-                // place source alone on this row
-                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
-                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
-
-                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
-
-                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-
-                Yo += DY;
-                Xo = 0;
-                dY = 0;
-            } else {
-                // start the next row
-                Yo += dY;
-                Xo = 0;
-
-                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
-                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
-
-                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
-
-                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-
-                Xo = DX;
-                dY = DY;
-            }
-        } else {
-            // extend this row
-            bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
-            if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-            psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
-
-            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-            psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
-            if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-
-            Xo += DX;
-            dY = PS_MAX (dY, DY);
-        }
+        pmSource *source = sources->data[i];
+
+        bool keep = false;
+        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
+            nPSF ++;
+            keep = true;
+        }
+        if (!keep) continue;
+
+        if (Xo + DX > NX) {
+            // too wide for the rest of this row
+            if (Xo == 0) {
+                // place source alone on this row
+                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
+
+                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
+
+                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+
+                Yo += DY;
+                Xo = 0;
+                dY = 0;
+            } else {
+                // start the next row
+                Yo += dY;
+                Xo = 0;
+
+                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
+
+                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
+
+                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+
+                Xo = DX;
+                dY = DY;
+            }
+        } else {
+            // extend this row
+            bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+            if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+            psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
+
+            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+            psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
+            if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+
+            Xo += DX;
+            dY = PS_MAX (dY, DY);
+        }
     }
 
@@ -994 +1013 @@
     bool status;
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.psf.sat", 3)) return true;
 
     int myKapa = psphotKapaChannel (3);
@@ -1023 +1042 @@
     for (int i = 0; i < sources->n; i++) {
 
-        pmSource *source = sources->data[i];
-
-        // only show "real" saturated stars (not defects)
-        if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;;
-        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;;
-
-        // how does this subimage get placed into the output image?
-        // DX = source->pixels->numCols
-        // DY = source->pixels->numRows
-
-        if (dX + DX > NX) {
-            // too wide for the rest of this row
-            if (dX == 0) {
-                // alone on this row
-                NY += DY;
-                dX = 0;
-                dY = 0;
-            } else {
-                // start the next row
-                NY += dY;
-                dX = DX;
-                dY = DY;
-            }
-        } else {
-            // extend this row
-            dX += DX;
-            dY = PS_MAX (dY, DY);
-        }
+        pmSource *source = sources->data[i];
+
+        // only show "real" saturated stars (not defects)
+        if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;;
+        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;;
+
+        // how does this subimage get placed into the output image?
+        // DX = source->pixels->numCols
+        // DY = source->pixels->numRows
+
+        if (dX + DX > NX) {
+            // too wide for the rest of this row
+            if (dX == 0) {
+                // alone on this row
+                NY += DY;
+                dX = 0;
+                dY = 0;
+            } else {
+                // start the next row
+                NY += dY;
+                dX = DX;
+                dY = DY;
+            }
+        } else {
+            // extend this row
+            dX += DX;
+            dY = PS_MAX (dY, DY);
+        }
     }
     NY += DY;
@@ -1068 +1087 @@
     for (int i = 0; i < sources->n; i++) {
 
-        pmSource *source = sources->data[i];
-
-        // only show "real" saturated stars (not defects)
-        if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;;
-        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;;
-        nSAT ++;
-
-        if (Xo + DX > NX) {
-            // too wide for the rest of this row
-            if (Xo == 0) {
-                // place source alone on this row
-                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
-                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
-                if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-
-                Yo += DY;
-                Xo = 0;
-                dY = 0;
-            } else {
-                // start the next row
-                Yo += dY;
-                Xo = 0;
-
-                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
-                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
-                if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-
-                Xo = DX;
-                dY = DY;
-            }
-        } else {
-            // extend this row
-            bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
-            if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-            psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
-            if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-
-            Xo += DX;
-            dY = PS_MAX (dY, DY);
-        }
+        pmSource *source = sources->data[i];
+
+        // only show "real" saturated stars (not defects)
+        if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;;
+        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;;
+        nSAT ++;
+
+        if (Xo + DX > NX) {
+            // too wide for the rest of this row
+            if (Xo == 0) {
+                // place source alone on this row
+                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
+                if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+
+                Yo += DY;
+                Xo = 0;
+                dY = 0;
+            } else {
+                // start the next row
+                Yo += dY;
+                Xo = 0;
+
+                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
+                if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+
+                Xo = DX;
+                dY = DY;
+            }
+        } else {
+            // extend this row
+            bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+            if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+            psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
+            if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+
+            Xo += DX;
+            dY = PS_MAX (dY, DY);
+        }
     }
 
@@ -1151 +1170 @@
     float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
     for (int iy = 0; iy < source->pixels->numRows; iy++) {
-        for (int ix = 0; ix < source->pixels->numCols; ix++) {
-            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
-                rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
-                // rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
-                Rb->data.F32[nb] = log10(rb->data.F32[nb]);
-                fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]);
-                nb++;
-            } else {
-                rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
-                // rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ;
-                Rg->data.F32[ng] = log10(rg->data.F32[ng]);
-                fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]);
-                ng++;
-            }
-        }
+        for (int ix = 0; ix < source->pixels->numCols; ix++) {
+            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
+                rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
+                // rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
+                Rb->data.F32[nb] = log10(rb->data.F32[nb]);
+                fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]);
+                nb++;
+            } else {
+                rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
+                // rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ;
+                Rg->data.F32[ng] = log10(rg->data.F32[ng]);
+                fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]);
+                ng++;
+            }
+        }
     }
 
@@ -1358 +1377 @@
     KapaSection section;  // put the positive profile in one and the residuals in another?
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.profiles", 3)) return true;
 
     int myKapa = psphotKapaChannel (2);
@@ -1394 +1413 @@
     for (int i = 0; i < sources->n; i++) {
 
-        pmSource *source = sources->data[i];
-        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
-
-        psphotVisualPlotRadialProfile (myKapa, source, maskVal);
-
-        // pause and wait for user input:
-        // continue, save (provide name), ??
-        char key[10];
-        fprintf (stdout, "[e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : ");
-        if (!fgets(key, 8, stdin)) {
-            psWarning("Unable to read option");
-        }
-        if (key[0] == 'e') {
-            KapaClearPlots (myKapa);
-        }
-        if (key[0] == 's') {
-            break;
-        }
+        pmSource *source = sources->data[i];
+        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+
+        psphotVisualPlotRadialProfile (myKapa, source, maskVal);
+
+        // pause and wait for user input:
+        // continue, save (provide name), ??
+        char key[10];
+        fprintf (stdout, "[e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : ");
+        if (!fgets(key, 8, stdin)) {
+            psWarning("Unable to read option");
+        }
+        if (key[0] == 'e') {
+            KapaClearPlots (myKapa);
+        }
+        if (key[0] == 's') {
+            break;
+        }
     }
 
@@ -1429 +1448 @@
     return true;
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.objects.flags", 3)) return true;
 
     int myKapa = psphotKapaChannel (1);
@@ -1445 +1464 @@
     for (int i = 0; i < sources->n; i++) {
 
-        float Xo, Yo, Rmaj, Rmin, cs, sn;
-
-        pmSource *source = sources->data[i];
-        if (source == NULL) continue;
-
-        pmMoments *moments = source->moments;
-        if (0) {
-            emoments.x2 = moments->Mxx;
-            emoments.y2 = moments->Myy;
-            emoments.xy = moments->Mxy;
-            Xo = moments->Mx;
-            Yo = moments->My;
-
-            axes = psEllipseMomentsToAxes (emoments, 20.0);
-            Rmaj = 2.0*axes.major;
-            Rmin = 2.0*axes.minor;
-            cs = cos(axes.theta);
-            sn = sin(axes.theta);
-        } else {
-            Rmaj = Rmin = 5.0;
-            cs = 1.0;
-            sn = 0.0;
-            Xo = source->peak->xf;
-            Yo = source->peak->yf;
-        }
-
-        unsigned short int flagMask = 0x01;
-        for (int j = 0; j < 8; j++) {
-            if (source->mode & flagMask) {
-                overlayE[NoverlayE].type = KII_OVERLAY_LINE;
-                overlayE[NoverlayE].x = Xo;
-                overlayE[NoverlayE].y = Yo;
-
-                float phi = j*M_PI/4.0;
-                overlayE[NoverlayE].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
-                overlayE[NoverlayE].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
-                overlayE[NoverlayE].angle = 0;
-                overlayE[NoverlayE].text = NULL;
-                NoverlayE ++;
-                CHECK_REALLOCATE (overlayE, KiiOverlay, NOVERLAYE, NoverlayE, 100);
-            }
-            flagMask <<= 1;
-
-            if (source->mode & flagMask) {
-                overlayO[NoverlayO].type = KII_OVERLAY_LINE;
-                overlayO[NoverlayO].x = Xo + 1;
-                overlayO[NoverlayO].y = Yo;
-
-                float phi = j*M_PI/4.0;
-                overlayO[NoverlayO].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
-                overlayO[NoverlayO].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
-                overlayO[NoverlayO].angle = 0;
-                overlayO[NoverlayO].text = NULL;
-                NoverlayO ++;
-                CHECK_REALLOCATE (overlayO, KiiOverlay, NOVERLAYO, NoverlayO, 100);
-            }
-            flagMask <<= 1;
-        }
+        float Xo, Yo, Rmaj, Rmin, cs, sn;
+
+        pmSource *source = sources->data[i];
+        if (source == NULL) continue;
+
+        pmMoments *moments = source->moments;
+        if (0) {
+            emoments.x2 = moments->Mxx;
+            emoments.y2 = moments->Myy;
+            emoments.xy = moments->Mxy;
+            Xo = moments->Mx;
+            Yo = moments->My;
+
+            axes = psEllipseMomentsToAxes (emoments, 20.0);
+            Rmaj = 2.0*axes.major;
+            Rmin = 2.0*axes.minor;
+            cs = cos(axes.theta);
+            sn = sin(axes.theta);
+        } else {
+            Rmaj = Rmin = 5.0;
+            cs = 1.0;
+            sn = 0.0;
+            Xo = source->peak->xf;
+            Yo = source->peak->yf;
+        }
+
+        unsigned short int flagMask = 0x01;
+        for (int j = 0; j < 8; j++) {
+            if (source->mode & flagMask) {
+                overlayE[NoverlayE].type = KII_OVERLAY_LINE;
+                overlayE[NoverlayE].x = Xo;
+                overlayE[NoverlayE].y = Yo;
+
+                float phi = j*M_PI/4.0;
+                overlayE[NoverlayE].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
+                overlayE[NoverlayE].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
+                overlayE[NoverlayE].angle = 0;
+                overlayE[NoverlayE].text = NULL;
+                NoverlayE ++;
+                CHECK_REALLOCATE (overlayE, KiiOverlay, NOVERLAYE, NoverlayE, 100);
+            }
+            flagMask <<= 1;
+
+            if (source->mode & flagMask) {
+                overlayO[NoverlayO].type = KII_OVERLAY_LINE;
+                overlayO[NoverlayO].x = Xo + 1;
+                overlayO[NoverlayO].y = Yo;
+
+                float phi = j*M_PI/4.0;
+                overlayO[NoverlayO].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
+                overlayO[NoverlayO].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
+                overlayO[NoverlayO].angle = 0;
+                overlayO[NoverlayO].text = NULL;
+                NoverlayO ++;
+                CHECK_REALLOCATE (overlayO, KiiOverlay, NOVERLAYO, NoverlayO, 100);
+            }
+            flagMask <<= 1;
+        }
     }
 
@@ -1531 +1550 @@
     for (int i = 0; i < sources->n; i++) {
 
-        pmSource *source = sources->data[i];
-        if (source == NULL) continue;
-
-        if (mode) {
-            if (keep) {
-                if (!(source->mode & mode)) continue;
-            } else {
-                if (source->mode & mode) continue;
-            }
-        }
-
-        pmMoments *moments = source->moments;
-        if (moments == NULL) continue;
-
-        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
-        overlay[Noverlay].x = moments->Mx;
-        overlay[Noverlay].y = moments->My;
-
-        emoments.x2 = moments->Mxx;
-        emoments.y2 = moments->Myy;
-        emoments.xy = moments->Mxy;
-
-        axes = psEllipseMomentsToAxes (emoments, 20.0);
-
-        overlay[Noverlay].dx = scale*2.0*axes.major;
-        overlay[Noverlay].dy = scale*2.0*axes.minor;
-        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
-        overlay[Noverlay].text = NULL;
-        Noverlay ++;
+        pmSource *source = sources->data[i];
+        if (source == NULL) continue;
+
+        if (mode) {
+            if (keep) {
+                if (!(source->mode & mode)) continue;
+            } else {
+                if (source->mode & mode) continue;
+            }
+        }
+
+        pmMoments *moments = source->moments;
+        if (moments == NULL) continue;
+
+        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
+        overlay[Noverlay].x = moments->Mx;
+        overlay[Noverlay].y = moments->My;
+
+        emoments.x2 = moments->Mxx;
+        emoments.y2 = moments->Myy;
+        emoments.xy = moments->Mxy;
+
+        axes = psEllipseMomentsToAxes (emoments, 20.0);
+
+        overlay[Noverlay].dx = scale*2.0*axes.major;
+        overlay[Noverlay].dy = scale*2.0*axes.minor;
+        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
+        overlay[Noverlay].text = NULL;
+        Noverlay ++;
     }
 
@@ -1570 +1589 @@
 bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) {
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.objects.size", 2)) return true;
 
     int myKapa = psphotKapaChannel (1);
@@ -1597 +1616 @@
     KapaSection section;
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.size", 2)) return true;
 
     int myKapa = psphotKapaChannel (2);
@@ -1612 +1631 @@
     float Ymin = 1000.0, Ymax = 0.0;
     {
-        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
-        for (int n = 0; n < nRegions; n++) {
-
-            char regionName[64];
-            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
-            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
-
-            float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
-            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
-            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
-            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
-
-            float X0 = psfX - 10.0*psfdX;
-            float X1 = psfX + 10.0*psfdX;
-            float Y0 = psfY - 10.0*psfdY;
-            float Y1 = psfY + 10.0*psfdY;
-
-            if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
-            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
-            if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
-            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
-        }
+        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
+        for (int n = 0; n < nRegions; n++) {
+
+            char regionName[64];
+            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
+            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
+
+            float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
+            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
+            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
+            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
+
+            float X0 = psfX - 10.0*psfdX;
+            float X1 = psfX + 10.0*psfdX;
+            float Y0 = psfY - 10.0*psfdY;
+            float Y1 = psfY + 10.0*psfdY;
+
+            if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
+            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
+            if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
+            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
+        }
     }
     Xmin = PS_MAX(Xmin, -0.1);
@@ -1677 +1696 @@
     int nCR  = 0;
     for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-        if (source->moments == NULL) continue;
+        pmSource *source = sources->data[i];
+        if (source->moments == NULL) continue;
 
         // only plot the measured sources...
-        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) continue;
-
-        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) {
-            xCR->data.F32[nCR] = source->moments->Mxx;
-            yCR->data.F32[nCR] = source->moments->Myy;
-            mCR->data.F32[nCR] = -2.5*log10(source->moments->Sum);
-            sCR->data.F32[nCR] = source->extNsigma;
-            nCR++;
-        }
-        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
-            xSAT->data.F32[nSAT] = source->moments->Mxx;
-            ySAT->data.F32[nSAT] = source->moments->Myy;
-            mSAT->data.F32[nSAT] = -2.5*log10(source->moments->Sum);
-            sSAT->data.F32[nSAT] = source->extNsigma;
-            nSAT++;
-        }
-        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-            xEXT->data.F32[nEXT] = source->moments->Mxx;
-            yEXT->data.F32[nEXT] = source->moments->Myy;
-            mEXT->data.F32[nEXT] = -2.5*log10(source->moments->Sum);
-            sEXT->data.F32[nEXT] = source->extNsigma;
-            nEXT++;
-            continue;
-        }
-        if (source->mode & PM_SOURCE_MODE_DEFECT) {
-            xDEF->data.F32[nDEF] = source->moments->Mxx;
-            yDEF->data.F32[nDEF] = source->moments->Myy;
-            mDEF->data.F32[nDEF] = -2.5*log10(source->moments->Sum);
-            sDEF->data.F32[nDEF] = source->extNsigma;
-            nDEF++;
-            continue;
-        }
-        if (source->errMag > 0.1) {
-            xLOW->data.F32[nLOW] = source->moments->Mxx;
-            yLOW->data.F32[nLOW] = source->moments->Myy;
-            mLOW->data.F32[nLOW] = -2.5*log10(source->moments->Sum);
-            sLOW->data.F32[nLOW] = source->extNsigma;
-            nLOW++;
-            continue;
-        }
-        xPSF->data.F32[nPSF] = source->moments->Mxx;
-        yPSF->data.F32[nPSF] = source->moments->Myy;
-        mPSF->data.F32[nPSF] = -2.5*log10(source->moments->Sum);
-        sPSF->data.F32[nPSF] = source->extNsigma;
-        nPSF++;
+        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) continue;
+
+        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) {
+            xCR->data.F32[nCR] = source->moments->Mxx;
+            yCR->data.F32[nCR] = source->moments->Myy;
+            mCR->data.F32[nCR] = -2.5*log10(source->moments->Sum);
+            sCR->data.F32[nCR] = source->extNsigma;
+            nCR++;
+        }
+        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
+            xSAT->data.F32[nSAT] = source->moments->Mxx;
+            ySAT->data.F32[nSAT] = source->moments->Myy;
+            mSAT->data.F32[nSAT] = -2.5*log10(source->moments->Sum);
+            sSAT->data.F32[nSAT] = source->extNsigma;
+            nSAT++;
+        }
+        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+            xEXT->data.F32[nEXT] = source->moments->Mxx;
+            yEXT->data.F32[nEXT] = source->moments->Myy;
+            mEXT->data.F32[nEXT] = -2.5*log10(source->moments->Sum);
+            sEXT->data.F32[nEXT] = source->extNsigma;
+            nEXT++;
+            continue;
+        }
+        if (source->mode & PM_SOURCE_MODE_DEFECT) {
+            xDEF->data.F32[nDEF] = source->moments->Mxx;
+            yDEF->data.F32[nDEF] = source->moments->Myy;
+            mDEF->data.F32[nDEF] = -2.5*log10(source->moments->Sum);
+            sDEF->data.F32[nDEF] = source->extNsigma;
+            nDEF++;
+            continue;
+        }
+        if (source->errMag > 0.1) {
+            xLOW->data.F32[nLOW] = source->moments->Mxx;
+            yLOW->data.F32[nLOW] = source->moments->Myy;
+            mLOW->data.F32[nLOW] = -2.5*log10(source->moments->Sum);
+            sLOW->data.F32[nLOW] = source->extNsigma;
+            nLOW++;
+            continue;
+        }
+        xPSF->data.F32[nPSF] = source->moments->Mxx;
+        yPSF->data.F32[nPSF] = source->moments->Myy;
+        mPSF->data.F32[nPSF] = -2.5*log10(source->moments->Sum);
+        sPSF->data.F32[nPSF] = source->extNsigma;
+        nPSF++;
     }
 
@@ -2029 +2048 @@
     // draw N circles to outline the clumps
     {
-        KapaSelectSection (myKapa, "MxxMyy");
-
-        // draw a circle centered on psfX,Y with size of the psf limit
-        psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
-        psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
-
-        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
-        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
-
-        graphdata.color = KapaColorByName ("blue");
-        graphdata.style = 0;
-
-        graphdata.xmin = Xmin;
-        graphdata.ymin = Ymin;
-        graphdata.xmax = Xmax;
-        graphdata.ymax = Ymax;
-        KapaSetLimits (myKapa, &graphdata);
-
-        for (int n = 0; n < nRegions; n++) {
-
-            char regionName[64];
-            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
-            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
-
-            float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
-            float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
-            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
-            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
-            float Rx = psfdX * PSF_CLUMP_NSIGMA;
-            float Ry = psfdY * PSF_CLUMP_NSIGMA;
-
-            for (int i = 0; i < xLimit->n; i++) {
-                xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX;
-                yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
-            }
-            KapaPrepPlot (myKapa, xLimit->n, &graphdata);
-            KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x");
-            KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y");
-        }
-        psFree (xLimit);
-        psFree (yLimit);
+        KapaSelectSection (myKapa, "MxxMyy");
+
+        // draw a circle centered on psfX,Y with size of the psf limit
+        psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
+        psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
+
+        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
+        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
+
+        graphdata.color = KapaColorByName ("blue");
+        graphdata.style = 0;
+
+        graphdata.xmin = Xmin;
+        graphdata.ymin = Ymin;
+        graphdata.xmax = Xmax;
+        graphdata.ymax = Ymax;
+        KapaSetLimits (myKapa, &graphdata);
+
+        for (int n = 0; n < nRegions; n++) {
+
+            char regionName[64];
+            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
+            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
+
+            float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
+            float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
+            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
+            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
+            float Rx = psfdX * PSF_CLUMP_NSIGMA;
+            float Ry = psfdY * PSF_CLUMP_NSIGMA;
+
+            for (int i = 0; i < xLimit->n; i++) {
+                xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX;
+                yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
+            }
+            KapaPrepPlot (myKapa, xLimit->n, &graphdata);
+            KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x");
+            KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y");
+        }
+        psFree (xLimit);
+        psFree (yLimit);
     }
 
@@ -2108 +2127 @@
 bool psphotVisualShowResidualImage (pmReadout *readout) {
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.image.resid", 2)) return true;
 
     int myKapa = psphotKapaChannel (1);
@@ -2124 +2143 @@
     float lineX[2], lineY[2];
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.apresid", 1)) return true;
 
     int myKapa = psphotKapaChannel (2);
@@ -2144 +2163 @@
     int n = 0;
     for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-        if (!source) continue;
-        if (source->type != PM_SOURCE_TYPE_STAR) continue;
-        if (!isfinite (source->apMag)) continue;
-        if (!isfinite (source->psfMag)) continue;
-
-        x->data.F32[n] = source->psfMag;
-        y->data.F32[n] = source->apMag - source->psfMag;
-        dy->data.F32[n] = source->errMag;
-        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
-        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
-        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]);
-        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]);
-
-        n++;
+        pmSource *source = sources->data[i];
+        if (!source) continue;
+        if (source->type != PM_SOURCE_TYPE_STAR) continue;
+        if (!isfinite (source->apMag)) continue;
+        if (!isfinite (source->psfMag)) continue;
+
+        x->data.F32[n] = source->psfMag;
+        y->data.F32[n] = source->apMag - source->psfMag;
+        dy->data.F32[n] = source->errMag;
+        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
+        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
+        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]);
+        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]);
+
+        n++;
     }
     x->n = y->n = dy->n = n;
@@ -2242 +2261 @@
     Graphdata graphdata;
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.chisq", 1)) return true;
 
     int myKapa = psphotKapaChannel (2);
@@ -2263 +2282 @@
     int n = 0;
     for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-        if (!source) continue;
-        if (source->type != PM_SOURCE_TYPE_STAR) continue;
-        if (!source->moments) continue;
-        if (!isfinite(source->moments->Sum)) continue;
-        if (!source->modelPSF) continue;
-        if (!isfinite(source->modelPSF->chisq)) continue;
-
-        x->data.F32[n] = -2.5*log10(source->moments->Sum);
-        y->data.F32[n] = source->modelPSF->chisq / source->modelPSF->nDOF;
-        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
-        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
-        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]);
-        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]);
-
-        fprintf (f, "%d %d %f %f\n", i, n, x->data.F32[n], y->data.F32[n]);
-
-        n++;
+        pmSource *source = sources->data[i];
+        if (!source) continue;
+        if (source->type != PM_SOURCE_TYPE_STAR) continue;
+        if (!source->moments) continue;
+        if (!isfinite(source->moments->Sum)) continue;
+        if (!source->modelPSF) continue;
+        if (!isfinite(source->modelPSF->chisq)) continue;
+
+        x->data.F32[n] = -2.5*log10(source->moments->Sum);
+        y->data.F32[n] = source->modelPSF->chisq / source->modelPSF->nDOF;
+        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
+        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
+        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]);
+        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]);
+
+        fprintf (f, "%d %d %f %f\n", i, n, x->data.F32[n], y->data.F32[n]);
+
+        n++;
     }
     x->n = y->n = n;
@@ -2326 +2345 @@
     KiiOverlay *overlay;
 
-    if (!DEBUG && !pmVisualIsVisual()) return true;
+    if (!pmVisualTestLevel("psphot.objects.petro", 2)) return true;
 
     int kapa = psphotKapaChannel (1);
@@ -2336 +2355 @@
 
     for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-
-        if (!source) continue;
-        if (!source->extpars) continue;
-        if (!source->extpars->petProfile) continue;
-
-        float petrosianRadius = source->extpars->petrosianRadius;
+        pmSource *source = sources->data[i];
+
+        if (!source) continue;
+        if (!source->extpars) continue;
+        if (!source->extpars->petProfile) continue;
+
+        float petrosianRadius = source->extpars->petrosianRadius;
         psEllipseAxes *axes = &source->extpars->axes;
 
-        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
-        overlay[Noverlay].x = source->peak->xf;
-        overlay[Noverlay].y = source->peak->yf;
-        overlay[Noverlay].dx = 1.0*petrosianRadius;
-        overlay[Noverlay].dy = 1.0*petrosianRadius*axes->minor/axes->major;
-        overlay[Noverlay].angle = axes->theta * PS_DEG_RAD;
-        overlay[Noverlay].text = NULL;
-        Noverlay ++;
-        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
-
-        // overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
-        // overlay[Noverlay].x = source->peak->xf;
-        // overlay[Noverlay].y = source->peak->yf;
-        // overlay[Noverlay].dx = 2.0*petrosianRadius;
-        // overlay[Noverlay].dy = 2.0*petrosianRadius*axes->minor/axes->major;
-        // overlay[Noverlay].angle = axes->theta * PS_DEG_RAD;
-        // overlay[Noverlay].text = NULL;
-        // Noverlay ++;
-        // CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
+        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
+        overlay[Noverlay].x = source->peak->xf;
+        overlay[Noverlay].y = source->peak->yf;
+        overlay[Noverlay].dx = 1.0*petrosianRadius;
+        overlay[Noverlay].dy = 1.0*petrosianRadius*axes->minor/axes->major;
+        overlay[Noverlay].angle = axes->theta * PS_DEG_RAD;
+        overlay[Noverlay].text = NULL;
+        Noverlay ++;
+        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
+
+        // overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
+        // overlay[Noverlay].x = source->peak->xf;
+        // overlay[Noverlay].y = source->peak->yf;
+        // overlay[Noverlay].dx = 2.0*petrosianRadius;
+        // overlay[Noverlay].dy = 2.0*petrosianRadius*axes->minor/axes->major;
+        // overlay[Noverlay].angle = axes->theta * PS_DEG_RAD;
+        // overlay[Noverlay].text = NULL;
+        // Noverlay ++;
+        // CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     }