Changeset 29004 for trunk/psphot/src

Timestamp:

Aug 20, 2010, 1:14:11 PM (16 years ago)

Author:

eugene

Message:

Many changes:

• psphot-related issues

1) added a new feature to psLib/src/fft/psImageFFT to allow for a pre-FFTed kernel generated for a specific image dimension which can then be applied to an arbitrary number of images for convolution -- this avoids the kernel FFT for every image convolution.
2) updated the psMinimizeLMM to include 2 levels of tolerance: minTol and maxTol: if the minimization reaches minTol, it stops. if it only reaches maxTol within maxIter, it stops and is considered successful, if it fails to reach maxTol, then it fails. This allows us to accept fits that are actually acceptable (within error) even if they are not as close as ideally possible.
3) add new stat: psfWeightNotPoor (vs psfWeightNotBad) : the first gives the fraction of psf-weighted unmasked pixels considering any mask bits (except the internal 'mark' bit), while the former considers only 'bad' mask bits -- these are written to QF_PSF and QF_PSF_PERFECT in the CMF files.
4) define user-set parameters for max and min valid flux in the linear fit analysis: Note this was the cause of the non-negative fluxes in forced photometry -- the min limit was hard-wired to 0.0.
5) significance image is now constructed as (image + (image/1000)^{2) / sqrt(variance) so that bright sources have well-defined peaks (other wise, they become somewhat flat-topped as image = sqrt(variance), leading to ill-defined peaks).
6) modification of the visualization functions to accept facility and level values akin to psTrace
7) modification of the source fitting APIs to allow thread- and source- independent fit options (iterations, tolerances, etc)
8) use Kron magnitude as test for source size (CR vs EXT) instead of PSF-based aperture
9) set a min systematic error in the aperture mags when used for size classification significance
10) threaded the psf-convolved model (PCM) fitting process (extended source fits)
11) for extended sources, adjust the radius based on the footprint and re-calculate moments for guess; save the psf-based moments for output in psf table
12) for Sersic fitting, choose the best index with a grid search using only a few iterations; iterate fully on the selected value.
13) same for Sersic fitting using the PCM fitting process
14) fixed a bug in which the PSF candidate stars which failed in the psf fitting were not correctly unmarked as being PSF stars
15) define galaxy fit radius based on sky stdev (global, not local)
16) move the PCM code to psModules
17) save and report the raw aperture mag in addition to the curve-of-growth corrected one (PS1_V3)
18) save and report the Kron parameters and higher-order moments
19) some psModule header reorganization for code clarity
20) fixed a bug in which the diff stats were counting 'marked' pixels (outside area of interest) as 'masked'.
21) save the radial profile aperture sizes in the headers
22) better guess for Sersic parameters based on moments (requires setting the functional form so that the scale length is right)}

• ppSub-related:

1) ensure masked pixels are NANed in output diff image
2) add code to flag detections if they have bright positive neighbors (+ output of these in PS1_DV2)
3) define separate penalties for each image (based on their fwhm values) (this requires the penalties to be calculated later in the code).
4) define separate apertures for each image for flux normalization
5) choose aperture based on curve-of-growth (was based on fixed fraction of full aperture flux, and thus noisy)
6) some fine tuning of the penalty factor (this still seems arbitrary, and results are somewhat sensitive to the right value)

Location:

trunk/psphot

Files:

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

trunk/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       \

trunk/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

trunk/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;
 }

trunk/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

trunk/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"))) {

trunk/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
+    assert (status && isfinite(fitMinTol) && fitMinTol > 0);
+
+    float fitMaxTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_TOL"); // Fit tolerance
+    assert (status && isfinite(fitMaxTol) && fitMaxTol > 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 +100 @@
     if (!sources->n) {
         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
+        psFree (fitOptions);
         return true;
     }
@@ -112 +125 @@
             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 +135 @@
             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 +143 @@
         if (!psThreadPoolWait (false)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            psFree (fitOptions);
             return NULL;
         }
@@ -163 +154 @@
             } 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 +177 @@
       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 +196 @@
     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 +262 @@
         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 +279 @@
             }
         } 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 +288 @@
         }
 
+        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 +303 @@
 
     // 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;
 

trunk/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

trunk/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);
     

trunk/psphot/src/psphotExtendedSourceAnalysis.c

@@ -157 +157 @@
         }
 
-# 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);

trunk/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
@@ -154 +267 @@
         }
         Next ++;
+
+        // set the radius based on the footprint (also sets the mask pixels)
+        if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) {
+            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)) {
+            // 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)
@@ -201 +332 @@
           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 +346 @@
               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 +364 @@
 
           // test for fit quality / result
+          modelFit->fitRadius = radius;
           psArrayAdd (source->modelFits, 4, modelFit);
 
@@ -311 +443 @@
         }
     }
-
-    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);
+
+    // 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;
 }

trunk/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

trunk/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

trunk/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?)

trunk/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);
-}

trunk/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)

trunk/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;

trunk/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;
         }

trunk/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);
 }
 

trunk/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);

trunk/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) {

trunk/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;
+}

trunk/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");

trunk/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);

trunk/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];
             }
         }

trunk/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);
+        }
+        
+        pmSourceFitModel (source, model, &options, maskVal);
+        fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+
+        // pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
+        // pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+
+        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;
+}

trunk/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;

trunk/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
 

trunk/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);
     }