Changeset 30624

Timestamp:

Feb 13, 2011, 12:33:05 PM (15 years ago)

Author:

eugene

Message:

adding some of the metadata needed by PSPS to output headers; skip models with few valid pixels (eg, only the edge is showing); only do the linear fit on pixels within the fit radius; modify psf-match auto-scaling process; enable multiple target psfs for matched-psf aperture photometry; speed up analysis of the radial apertures

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• doc/notes.20101222.txt (copied) (copied from branches/eam_branches/ipp-20101205/psphot/doc/notes.20101222.txt )
• doc/stack.txt (modified) (1 diff)
• src/Makefile.am (modified) (1 diff)
• src/psphot.h (modified) (6 diffs)
• src/psphotApResid.c (modified) (1 prop)
• src/psphotBlendFit.c (modified) (1 diff)
• src/psphotCheckStarDistribution.c (modified) (3 diffs)
• src/psphotChoosePSF.c (modified) (16 diffs)
• src/psphotExtendedSourceAnalysis.c (modified) (1 diff)
• src/psphotExtendedSourceAnalysisByObject.c (modified) (1 diff)
• src/psphotExtendedSourceFits.c (modified) (11 diffs)
• src/psphotFindDetections.c (modified) (2 diffs)
• src/psphotFindFootprints.c (modified) (1 diff)
• src/psphotFitSourcesLinear.c (modified) (4 diffs)
• src/psphotFitSourcesLinearStack.c (modified) (1 diff)
• src/psphotImageQuality.c (modified) (2 diffs)
• src/psphotMakeFluxScale.c (modified) (1 prop)
• src/psphotMakePSFReadout.c (modified) (1 diff)
• src/psphotMakeResiduals.c (modified) (2 diffs)
• src/psphotMergeSources.c (modified) (1 diff)
• src/psphotOutput.c (modified) (2 diffs)
• src/psphotPetrosianAnalysis.c (modified) (1 diff)
• src/psphotPetrosianProfile.c (modified) (1 diff)
• src/psphotRadialApertures.c (modified) (8 diffs)
• src/psphotRadialAperturesByObject.c (modified) (6 diffs)
• src/psphotReadout.c (modified) (1 diff)
• src/psphotReadoutFindPSF.c (modified) (1 diff)
• src/psphotReadoutKnownSources.c (modified) (1 diff)
• src/psphotReplaceUnfit.c (modified) (2 diffs)
• src/psphotSetThreads.c (modified) (1 diff)
• src/psphotSourceFits.c (modified) (2 diffs)
• src/psphotSourceMatch.c (modified) (5 diffs)
• src/psphotSourceSize.c (modified) (1 diff)
• src/psphotStackArguments.c (modified) (1 diff)
• src/psphotStackChisqImage.c (modified) (4 diffs)
• src/psphotStackImageLoop.c (modified) (5 diffs)
• src/psphotStackMatchPSFs.c (modified) (3 diffs)
• src/psphotStackMatchPSFsNext.c (copied) (copied from branches/eam_branches/ipp-20101205/psphot/src/psphotStackMatchPSFsNext.c )
• src/psphotStackMatchPSFsUtils.c (modified) (2 diffs)
• src/psphotStackPSF.c (modified) (2 diffs)
• src/psphotStackParseCamera.c (modified) (6 diffs)
• src/psphotStackReadout.c (modified) (7 diffs)
• src/psphotVisual.c (modified) (4 diffs)
• test/tap_psphot_deteff.pro (copied) (copied from branches/eam_branches/ipp-20101205/psphot/test/tap_psphot_deteff.pro )
• test/tap_psphot_diff.pro (copied) (copied from branches/eam_branches/ipp-20101205/psphot/test/tap_psphot_diff.pro )
• test/tap_psphot_diffsuite.pro (copied) (copied from branches/eam_branches/ipp-20101205/psphot/test/tap_psphot_diffsuite.pro )
• test/tap_psphot_stack.pro (copied) (copied from branches/eam_branches/ipp-20101205/psphot/test/tap_psphot_stack.pro )

trunk/psphot/doc/stack.txt

@@ +1 @@
+
+20101221 
+
+  psphotStackReadout is now correctly subtracting the PSF models from
+  the images as it measures the radial aperture fluxes. 
+
+  Some issues:
+
+  * the source properties get buggered up by the radial aperture flux
+    analysis (we need to re-determine the psf, re-create the source
+    pixels, and re-fit the models (linearly) to subtract them
+    correctly).
+
+  * the standard analysis sequence is not doing a second pass
+
+  * the choice of the best model is ambiguous
+
+  * the radial aperture (and source addition / subtraction) is only
+    using the PSF model
+
+  * make sure psphotEfficiency actually subtracts the inserted fake
+    sources
+
+  * XXX how does the PSF-convolved model work with changing psf models?
+
+20101207
+
+  header fields for PSPS:
+
+  * stack needs to count # of inputs
+  * propagate the input header to the output image (DONE)
+  * save PSF parameters in header 
+  * stack type is only known to the launcher
+
+  * what are we doing for the stack calibrations??
+
+skycellID    : code
+surveyID     : code
+filterID     : code
+stackMetaID  : stack_id 
+photoCalID   : photcode -> number
+magSat       : FSATUR       saturation magnitude level ** not correctly set
+completMag   : FLIMIT       95% completion level in mag ** not correctly set 
+stackTypeID  : STK_TYPE     stack type identifier ** deep stack, nightly stack, best IQ stack?
+refImageID   :              identifier of image used as reference for analysis
+subtrImageID : N/A (stack)  identifier of image subtracted to generate difference image
+analVer      :              analysis version index  ** index for tess_id + skycell_id + filter?
+nP2Images    : NINPUTS      number of P2 images contributing to this cell ** missing from input stack
+astroScat    :              astrometric scatter for chip ** measure scatter on stack creation?
+photoScat    :              photometric scatter for chip ** internal scatter?  
+nAstroRef    :              number of astrometric reference sources ** (connected to above)
+nPhoRef      :              number of photometric reference sources ** same
+psfFwhm      :              PSF full width at half maximum ** 0.5(FWHM_MAJ + FWHM_MIN)
+psfmodelID   : * PSFMODEL   PSF model identifier ** save the PSF model name in the header
+psfWidMajor  : * FWHM_MAJ     PSF parameters
+psfWidMinor  : * FWHM_MIN     PSF parameters
+psfTheta     : * ANGLE        PSF parameters
+psfExtra1    : * PSF_EXT1   PSF parameters ** (at field center?)
+psfExtra2    : * PSF_EXT2   PSF parameters ** (not set for all models)
+photoZero    :              local derived photometric zero point 
+photoColor   :              local derived photometric color term
+ctype1       : * CTYPE1             name of astrometric projection in RA ** propagate from input stacks
+ctype2       : * CTYPE2             name of astrometric projection in DEC
+crval1       : * CRVAL1             RA corresponding to reference pixel
+crval2       : * CRVAL2             DEC corresponding to reference pixel
+crpix1       : * CRPIX1             reference pixel value for RA
+crpix2       : * CRPIX2             reference pixel value for DEC
+cdelt1       : * CDELT1             scale factor for RA
+cdelt2       : * CDELT2             scale factor for DEC
+pc001001     : * PC001001           elements of rotation/Dcale matrix
+pc001002     : * PC001002           elements of rotation/Dcale matrix
+pc002001     : * PC002001           elements of rotation/Dcale matrix
+pc002002     : * PC002002           elements of rotation/Dcale matrix
+polyOrder    : * NPLYTERM           polynomial order of astrometry fit ** default to 1
+pca1x3y0     : * PCA1X3Y0     polynomial coefficients for the astrometric fit
+pca1x2y1     : * PCA1X2Y1     polynomial coefficients for the astrometric fit
+pca1x1y2     : * PCA1X1Y2     polynomial coefficients for the astrometric fit
+pca1x0y3     : * PCA1X0Y3     polynomial coefficients for the astrometric fit
+pca1x2y0     : * PCA1X2Y0     polynomial coefficients for the astrometric fit
+pca1x1y1     : * PCA1X1Y1     polynomial coefficients for the astrometric fit
+pca1x0y2     : * PCA1X0Y2     polynomial coefficients for the astrometric fit
+pca2x3y0     : * PCA2X3Y0     polynomial coefficients for the astrometric fit
+pca2x2y1     : * PCA2X2Y1     polynomial coefficients for the astrometric fit
+pca2x1y2     : * PCA2X1Y2     polynomial coefficients for the astrometric fit
+pca2x0y3     : * PCA2X0Y3     polynomial coefficients for the astrometric fit
+pca2x2y0     : * PCA2X2Y0     polynomial coefficients for the astrometric fit
+pca2x1y1     : * PCA2X1Y1     polynomial coefficients for the astrometric fit
+pca2x0y2     : * PCA2X0Y2     polynomial coefficients for the astrometric fit
+calibModNum  :              calibration modification number
+dataRelease  :              Data release
 
 20100506:

trunk/psphot/src/Makefile.am

@@ -98 +98 @@
         psphotStackMatchPSFsUtils.c   \
         psphotStackMatchPSFsPrepare.c \
+        psphotStackMatchPSFsNext.c    \
         psphotStackOptions.c          \
         psphotStackObjects.c          \

trunk/psphot/src/psphot.h

@@ -75 +75 @@
 bool            psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
 
-bool            psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule);
-bool            psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+bool            psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule, bool newSources);
+bool            psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool newSources);
 
 bool            psphotGuessModels (pmConfig *config, const pmFPAview *view, const char *filerule);
@@ -260 +260 @@
 bool            psphotVisualShowSourceSize (pmReadout *readout, psArray *sources);
 bool            psphotVisualPlotSourceSizeAlt (psMetadata *recipe, psMetadata *analysis, psArray *sources);
-bool            psphotVisualShowResidualImage (pmReadout *readout);
+bool            psphotVisualShowResidualImage (pmReadout *readout, bool reshow);
+bool            psphotVisualShowObjectRegions (pmReadout *readout, psMetadata *recipe, psArray *sources);
 bool            psphotVisualPlotApResid (psArray *sources, float mean, float error, bool useApMag);
 bool            psphotVisualPlotChisq (psArray *sources);
@@ -312 +313 @@
 
 int psphotFileruleCount(const pmConfig *config, const char *filerule);
+bool psphotFileruleCountSet(const pmConfig *config, const char *filerule, int num);
 
 bool psphotAddKnownSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *inSources);
@@ -403 +405 @@
 bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
 bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
+bool psphotStackMatchPSFsNext (bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize);
+bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, psVector *fwhmValues, int lastSize);
 
 // psphotStackMatchPSFsUtils
@@ -426 +430 @@
 bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule);
 bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax);
+bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry);
 
 bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);
-bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);
+bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule, int nMatchedPSF);
 
 bool psphotStackObjectsUnifyPosition (psArray *objects);
@@ -440 +444 @@
 bool psphotCleanInputs (pmConfig *config, const pmFPAview *view, const char *filerule);
 
+bool psphotResetModels (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotResetModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+
+bool psphotRedefinePixels (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotRedefinePixelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+
+bool psphotSourceChildren (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc);
+bool psphotSourceChildrenReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
+psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc);
+
 #endif

trunk/psphot/src/psphotBlendFit.c

@@ -187 +187 @@
     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);
 
-    psphotVisualShowResidualImage (readout);
+    psphotVisualShowResidualImage (readout, false);
+    psphotVisualShowObjectRegions (readout, recipe, sources);
     psphotVisualShowFlags (sources);
 

trunk/psphot/src/psphotCheckStarDistribution.c

@@ -34 +34 @@
         pmSource *source = stars->data[i];
         if (source->peak == NULL) continue;
-        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+        if (!(source->tmpFlags & PM_SOURCE_TMPF_CANDIDATE_PSFSTAR)) continue;
 
         int binX = psImageBinningGetRuffX (binning, source->peak->xf);
@@ -85 +85 @@
         if (source->peak == NULL) continue;
         if (source->moments == NULL) continue;
-        if (source->mode & PM_SOURCE_MODE_PSFSTAR) continue;
+        if (source->tmpFlags & PM_SOURCE_TMPF_CANDIDATE_PSFSTAR) continue;
         if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
         if (source->type != PM_SOURCE_TYPE_STAR) continue;
@@ -97 +97 @@
         if (y > Ye) continue;
 
-        source->mode |= PM_SOURCE_MODE_PSFSTAR;
+        source->tmpFlags |= PM_SOURCE_TMPF_CANDIDATE_PSFSTAR;
         psArrayAdd (stars, 200, source);
 

trunk/psphot/src/psphotChoosePSF.c

@@ -2 +2 @@
 
 // generate a PSF model for inputs without PSF models already loaded
-bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule)
+bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule, bool newSources)
 {
     bool status = true;
@@ -19 +19 @@
     for (int i = 0; i < num; i++) {
         if (i == chisqNum) continue; // skip chisq image
-        if (!psphotChoosePSFReadout (config, view, filerule, i, recipe)) {
+        if (!psphotChoosePSFReadout (config, view, filerule, i, recipe, newSources)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for %s entry %d", filerule, i);
             return false;
@@ -28 +28 @@
 
 // try PSF models and select best option
-bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool newSources) {
 
     bool status;
@@ -50 +50 @@
     psAssert (detections, "missing detections?");
 
-    psArray *sources = detections->newSources;
+    psArray *sources = newSources ? detections->newSources : detections->allSources;
     psAssert (sources, "missing sources?");
 
@@ -166 +166 @@
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
-        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
-            // keep NSTARS PSF stars, unmark the rest
+        if (source->tmpFlags & PM_SOURCE_TMPF_CANDIDATE_PSFSTAR) {
+            // keep NSTARS PSF stars
             if (stars->n < NSTARS) {
                 psArrayAdd (stars, 200, source);
-            } else {
-                source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
             }
         }
@@ -214 +212 @@
         psFree(options);
 
-        // unset the PSFSTAR flags (none are used):
-        for (int i = 0; i < sources->n; i++) {
-            pmSource *source = sources->data[i];
-            source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
-        }
+        // no sources are used as PSF stars
 
         // XXX set sxx, etc from FWHM in recipe
@@ -292 +286 @@
         psLogMsg ("psphot.pspsf", PS_LOG_INFO, "Using guess PSF model");
 
-        // unset the PSFSTAR flags (none are used):
-        for (int i = 0; i < sources->n; i++) {
-            pmSource *source = sources->data[i];
-            source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
-        }
+        // no sources are used as PSF stars
 
         // XXX set sxx, etc from FWHM in recipe
@@ -322 +312 @@
     pmPSFtry *try = models->data[bestN];
 
-    // 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;
+    // set the PSFSTAR flag for stars used for the PSF model
+    int nKeep = 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]) {
-            // 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);
+        if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
+        source->mode |= PM_SOURCE_MODE_PSFSTAR;
+
+        // this source was used: find the parent and set its PSFSTAR flag
+        pmSource *realSource = source->parent;
+        psAssert (realSource, "pmPSFtryAlloc should have set the parent pointers");
+        realSource->mode |= PM_SOURCE_MODE_PSFSTAR;
+        nKeep ++;
+    }
+    psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "used %d of %ld candidate PSF objects\n", nKeep, try->sources->n);
 
     // build a PSF residual image
+    // we need to use the 'try' set because they have the associated models defined
     if (!psphotMakeResiduals (try->sources, recipe, try->psf, maskVal)) {
         psError(PSPHOT_ERR_PSF, false, "Unable to construct residual table for PSF");
@@ -386 +352 @@
     }
 
-    // XXX test dump of psf star data and psf-subtracted image
+    // test dump of psf star data and psf-subtracted image
     if (psTraceGetLevel("psphot.psfstars") > 5) {
         psphotDumpPSFStars (readout, try, options->fitRadius, maskVal, markVal);
@@ -413 +379 @@
         return false;
     }
-    psFree (psf); // XXX double-check
-
-    // move into psphotChoosePSF
+    psFree (psf);
+
     psphotVisualShowPSFModel (readout, psf);
 
@@ -436 +401 @@
     psVector *fwhmMajor = psVectorAllocEmpty (100, PS_DATA_F32);
     psVector *fwhmMinor = psVectorAllocEmpty (100, PS_DATA_F32);
+    psVector *psfExtra1 = psVectorAllocEmpty (100, PS_DATA_F32);
+    psVector *psfExtra2 = psVectorAllocEmpty (100, PS_DATA_F32);
 
     for (float ix = -0.4; ix <= +0.4; ix += 0.1) {
@@ -459 +426 @@
             shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
             axes = psEllipseShapeToAxes (shape, 20.0);
-            psFree (modelPSF);
 
             float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
             if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) {
                 fprintf (stderr, "!");
+                psFree (modelPSF);
                 continue;
             }
             psVectorAppend (fwhmMajor, FWHM_MAJOR);
             psVectorAppend (fwhmMinor, FWHM_MINOR);
+
+            if (modelPSF->params->n >= 7) {
+              psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]);
+            }
+            if (modelPSF->params->n >= 8) {
+              psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]);
+            }
+            psFree (modelPSF);
         }
     }
@@ -474 +449 @@
     if (!psVectorStats (stats, fwhmMajor, NULL, NULL, 0)) {
         psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MAJOR");
-        return false;
+        goto escape;
     }
 
@@ -486 +461 @@
     if (!psVectorStats (stats, fwhmMinor, NULL, NULL, 0)) {
         psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MINOR");
-        return false;
+        goto escape;
     }
     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MIN",   PS_META_REPLACE, "PSF FWHM Minor axis (mean)", stats->sampleMean);
@@ -494 +469 @@
 
     float fwhmMin = stats->sampleMean;  // FWHM on minor axis
+
     if (readout->parent) {
+
+        // we now have 2 possible measurements of the seeing : the PSF model based version and
+        // the Moments based value we need to define a definitive "CHIP.SEEING" value.  This
+        // value is used by the PSF-matching programs (ppSub and ppStack).  The moments-based
+        // value is probably more representative of the value needed by those tools (it is also
+        // more stable?)
+        bool status = false;
+        float fwhmMajorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1");
+        float fwhmMinorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1");
+
         pmChip *chip = readout->parent->parent; // Parent chip
         psAssert(chip, "Cell should be attached to a chip.");
         psMetadataItem *item = psMetadataLookup(chip->concepts, "CHIP.SEEING"); // Item with chip
-        item->data.F32 = 0.5 * (fwhmMaj + fwhmMin);
+        item->data.F32 = 0.5 * (fwhmMajorMoments + fwhmMinorMoments);
+
+        psLogMsg ("psphot", PS_LOG_DETAIL, "fwhm (psf): %f,%f (moments): %f,%f", fwhmMaj, fwhmMin, fwhmMajorMoments, fwhmMinorMoments);
+    }
+
+    if (psfExtra1->n) {
+      if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1");
+        goto escape;
+      }
+      psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean);
+    }
+
+    if (psfExtra2->n) {
+      if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2");
+        goto escape;
+      }
+      psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean);
     }
 
     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
     psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars);
-    psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", true);
+
+    char *psfModelName = pmModelClassGetName(psf->type);
+    psMetadataAddStr(readout->analysis,  PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "PSF Model Name", psfModelName);
+    psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSF_OK",   PS_META_REPLACE, "Valid PSF Model?", true);
+
+    int nParams = pmModelClassParameterCount(psf->type);
+    psMetadataAddS32(readout->analysis, PS_LIST_TAIL, "PSF_NPAR",   PS_META_REPLACE, "Number of PSF parameters", nParams);
 
     psFree (fwhmMajor);
     psFree (fwhmMinor);
+    psFree (psfExtra1);
+    psFree (psfExtra2);
     psFree (stats);
     return true;
+
+escape:
+    psFree (fwhmMajor);
+    psFree (fwhmMinor);
+    psFree (psfExtra1);
+    psFree (psfExtra2);
+    psFree (stats);
+    return false;
 }
 
@@ -568 +588 @@
     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_LQ",   PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", 0);
     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_UQ",   PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", 0);
-    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           FWHM_T);
-    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", 0);
-    psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", false);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE",      PS_META_REPLACE, "PSF angle",           FWHM_T);
+    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR",   PS_META_REPLACE, "Number of stars used to make PSF", 0);
+    psMetadataAddStr(readout->analysis,  PS_LIST_TAIL, "PSFMODEL",   PS_META_REPLACE, "PSF Model Name", "NONE");
+    psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSF_OK",     PS_META_REPLACE, "Valid PSF Model?", false);
 
     return true;

trunk/psphot/src/psphotExtendedSourceAnalysis.c

@@ -173 +173 @@
     psLogMsg ("psphot", PS_LOG_INFO, "  %d kron\n", Nkron);
 
-    psphotVisualShowResidualImage (readout);
+    psphotVisualShowResidualImage (readout, false);
 
     if (doPetrosian) {

trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c

@@ -53 +53 @@
         pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
         psAssert (readout, "missing readout?");
+
+        psLogMsg("psphot", PS_LOG_INFO, "petrosians for image %d", i);
+        psphotVisualShowImage(readout);
 
         readouts->data[i] = psMemIncrRefCounter(readout);

trunk/psphot/src/psphotExtendedSourceFits.c

@@ -18 +18 @@
     int num = psphotFileruleCount(config, filerule);
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
         if (!psphotExtendedSourceFitsReadout (config, view, filerule, i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on to fit extended sources for %s entry %d", filerule, i);
@@ -37 +42 @@
     int Nplain = 0;
     int NplainPass = 0;
+    int Nfaint = 0;
 
     psTimerStart ("psphot.extended");
@@ -46 +52 @@
     pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     psAssert (readout, "missing readout?");
+
+    psLogMsg("psphot", PS_LOG_INFO, "extended source fits for image %d", index);
+    psphotVisualShowImage(readout);
 
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
@@ -167 +176 @@
             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
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfain
 
             if (false && !psThreadJobAddPending(job)) {
@@ -189 +199 @@
                 scalar = job->args->data[11];
                 NplainPass += scalar->data.S32;
+                scalar = job->args->data[12];
+                Nfaint += scalar->data.S32;
                 psFree(job);
             }
@@ -217 +229 @@
                 scalar = job->args->data[11];
                 NplainPass += scalar->data.S32;
+                scalar = job->args->data[12];
+                Nfaint += scalar->data.S32;
             }
             psFree(job);
@@ -227 +241 @@
     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);
+    psLogMsg ("psphot", PS_LOG_INFO, "  %d too faint to fit\n", Nfaint);
     return true;
 }
@@ -238 +253 @@
     int NconvolvePass = 0;
     int Nplain = 0;
+    int Nfaint = 0;
     int NplainPass = 0;
     bool savePics = false;
@@ -326 +342 @@
 
         // array to store the pointers to the model flux images while the models are being fitted
-        psArray *modelFluxes = psArrayAllocEmpty (4);
+        psArray *modelFluxes = psArrayAllocEmpty (models->list->n);
 
         // allocate the array to store the model fits
         if (source->modelFits == NULL) {
-            source->modelFits = psArrayAllocEmpty (4);
+            source->modelFits = psArrayAllocEmpty (models->list->n);
         }
 
@@ -350 +366 @@
           // limit selection to some SN limit
           assert (source->peak); // how can a source not have a peak?
-          if (source->peak->SN < SNlim) continue;
+          if (source->peak->SN < SNlim) {
+              Nfaint ++;
+              continue;
+          }
 
           // check on the model type
@@ -497 +516 @@
     scalar->data.S32 = NplainPass;
 
+    scalar = job->args->data[12];
+    scalar->data.S32 = Nfaint;
+
     return true;
 }

trunk/psphot/src/psphotFindDetections.c

@@ -86 +86 @@
     psImage *significance = psphotSignificanceImage (readout, recipe, maskVal);
 
-    // display the significance image
-    psphotVisualShowSignificance (significance, -1.0, PS_SQR(3.0*NSIGMA_PEAK));
-
     // display the log significance image
     psphotVisualShowLogSignificance (significance, 0.0, 4.5);
+
+    // display the significance image
+    psphotVisualShowSignificance (significance, 0.98*threshold, 1.02*threshold);
 
     // detect the peaks in the significance image
@@ -108 +108 @@
     }
 
+    // XXX do a second (or third?) pass with rebinning (to detected more extended sources)
+
     psFree (significance);
 

trunk/psphot/src/psphotFindFootprints.c

@@ -21 +21 @@
     psArray *footprints = pmFootprintsFind (significance, threshold, npixMin);
 
-    pmFootprintsAssignPeaks(footprints, detections->peaks);
-    // XXX handle the error conditions here
+    if (pmFootprintsAssignPeaks(footprints, detections->peaks) != PS_ERR_NONE) {
+        psAbort ("inconsistent peaks and footprints");
+    }
 
     // footprints now owns the peaks; after culling (below), we will rebuild the peaks array

trunk/psphot/src/psphotFitSourcesLinear.c

@@ -46 +46 @@
             return false;
         }
+
+        psphotVisualShowResidualImage (readout, (num > 0)); 
+        psphotVisualShowPeaks (detections);
+        psphotVisualShowObjectRegions (readout, recipe, sources);
     }
     return true;
@@ -151 +155 @@
         if (x > AnalysisRegion.x1) continue;
         if (y > AnalysisRegion.y1) continue;
+
+        // check the integral of the model : is it large enough?
+        float modelSum = 0.0;
+        for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
+            for (int ix = 0; ix < source->modelFlux->numCols; ix++) {
+                modelSum += source->modelFlux->data.F32[iy][ix];
+            }
+        }
+        if (modelSum < 0.5) continue; // skip sources with no model constraint (somewhat arbitrary limit)
+        // if (modelSum < 0.01) continue; // skip sources with no model constraint (somewhat arbitrary limit)
+        if (modelSum < 0.8) {
+            fprintf (stderr, "low-sig model @ %f, %f (%f sum, %f peak)\n",
+                     source->peak->xf, source->peak->yf, modelSum, source->peak->flux);
+        }
+
+        pmModel *model = pmSourceGetModel (NULL, source);
+
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, model->fitRadius, "OR", markVal);
 
         source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
@@ -270 +293 @@
         model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
         model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
-        // XXX is the value of 'errors' modified by the sky fit?
 
         // subtract object
@@ -297 +319 @@
     psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
 
-    psphotVisualShowResidualImage (readout);
     psphotVisualPlotChisq (sources);
     // psphotVisualShowFlags (sources);

trunk/psphot/src/psphotFitSourcesLinearStack.c

@@ -176 +176 @@
 
     psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
-
     return true;
 }

trunk/psphot/src/psphotImageQuality.c

@@ -79 +79 @@
         // select by PSFSTAR mode?
         // ??
-        if (source->type != PM_SOURCE_TYPE_STAR || !(source->mode & PM_SOURCE_MODE_PSFSTAR)) {
+        if (source->type != PM_SOURCE_TYPE_STAR || !(source->tmpFlags & PM_SOURCE_TMPF_CANDIDATE_PSFSTAR)) {
             psTrace("psphot", 10, "Ignoring source for image quality because not a good star");
             continue;
@@ -133 +133 @@
     if (num == 0) {
         psLogMsg ("psphot", PS_LOG_INFO, "no valid sources for image quality, skipping");
+        psFree(FWHM_MAJOR);
+        psFree(FWHM_MINOR);
+        psFree(M2);
+        psFree(M2c);
+        psFree(M2s);
+        psFree(M3);
+        psFree(M4);
+
         return true;
     }

trunk/psphot/src/psphotMakePSFReadout.c

@@ -60 +60 @@
     // Use bright stellar objects to measure PSF. If we do not have enough stars to generate
     // the PSF, build one from the SEEING guess and model class
-    if (!psphotChoosePSF (config, view, filerule)) {
+    if (!psphotChoosePSF (config, view, filerule, true)) {
         psLogMsg ("psphot", 3, "failure to construct a psf model");
         return psphotReadoutCleanup (config, view, filerule);

trunk/psphot/src/psphotMakeResiduals.c

@@ -136 +136 @@
         psFree (interp);
     }
+    xSize = PS_MIN(xSize, 2*radiusMax+3);
+    ySize = PS_MIN(ySize, 2*radiusMax+3);
+
     pmResiduals *resid = pmResidualsAlloc (xSize, ySize, xBin, yBin);
     psImageInit (resid->mask, 0);
@@ -316 +319 @@
     psFree (B);
 
-    psLogMsg ("psphot.pspsf", PS_LOG_MINUTIA, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("psphot.residuals"));
+    psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("psphot.residuals"));
 
     psFree (xC);

trunk/psphot/src/psphotMergeSources.c

@@ -513 +513 @@
     }
 
-    return true;
-}
-
+    // loop over the sources, redefine their pixels to point at the new filerule image,
+    // copy the source data, and add a reference back to the original source
+    
+
+    return true;
+}
+
+// create source children from ruleSrc for ruleOut
+bool psphotSourceChildren (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc)
+{
+    bool status = true;
+
+    int num = psphotFileruleCount(config, ruleSrc);
+
+    // skip the chisq image because it is a duplicate of the detection version
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+        if (!psphotSourceChildrenReadout (config, view, ruleOut, ruleSrc, i)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to copy sources from %s to %s entry %d", ruleSrc, ruleOut, i);
+            return false;
+        }
+    }
+    return true;
+}
+
+// create source children from ruleSrc for ruleOut for this entry.  XXX currently, this is only
+// used by psphotStackReadout (sources go on allSources so that psphotChoosePSF can be called
+// repeatedly)
+bool psphotSourceChildrenReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) {
+
+    bool status;
+
+    // find the currently selected readout
+    pmFPAfile *fileSrc = pmFPAfileSelectSingle(config->files, ruleSrc, index); // File of interest
+    psAssert (fileSrc, "missing file?");
+
+    pmReadout *readoutSrc = pmFPAviewThisReadout(view, fileSrc->fpa);
+    psAssert (readoutSrc, "missing readout?");
+
+    pmDetections *detectionsSrc = psMetadataLookupPtr (&status, readoutSrc->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detectionsSrc, "missing detections?");
+
+    psArray *sourcesSrc = detectionsSrc->allSources;
+    psAssert (sourcesSrc, "missing sources?");
+
+    // find the currently selected readout
+    pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, ruleOut, index); // File of interest
+    psAssert (fileOut, "missing file?");
+
+    pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa);
+    psAssert (readoutOut, "missing readout?");
+
+    // generate a new detection structure for the output filerule
+    pmDetections *detectionsOut = psMetadataLookupPtr (&status, readoutOut->analysis, "PSPHOT.DETECTIONS");
+    if (!detectionsOut) {
+        detectionsOut = pmDetectionsAlloc();
+        detectionsOut->allSources = psArrayAllocEmpty (100);
+        // save detections on the readout->analysis
+        if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detectionsOut)) {
+            psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout");
+            return false;
+        }
+        psFree(detectionsOut); // a copy remains on the analysis metadata
+    }
+
+    // loop over the sources, redefine their pixels to point at the new filerule image,
+    // copy the source data, and add a reference back to the original source
+    
+    // copy the sources from sourceSrcs to the new detection structure
+    for (int i = 0; i < sourcesSrc->n; i++) {
+      pmSource *sourceSrc = sourcesSrc->data[i];
+
+      pmSource *sourceOut = pmSourceCopy(sourceSrc);
+      sourceOut->parent = sourceSrc;
+      
+      // keep the original source flags
+      sourceOut->type     = sourceSrc->type;
+      sourceOut->mode     = sourceSrc->mode;
+      sourceOut->mode2    = sourceSrc->mode2;
+      sourceOut->tmpFlags = sourceSrc->tmpFlags;
+
+      // does this copy all model data? (NO)
+      sourceOut->modelPSF = pmModelCopy(sourceSrc->modelPSF);
+      sourceOut->modelEXT = pmModelCopy(sourceSrc->modelEXT);
+
+      if (sourceSrc->modelFits) {
+          sourceOut->modelFits = psArrayAlloc(sourceSrc->modelFits->n);
+          for (int j = 0; j < sourceSrc->modelFits->n; j++) {
+              sourceOut->modelFits->data[j] = pmModelCopy(sourceSrc->modelFits->data[j]);
+          }
+      }
+
+      // drop the references to the original image pixels:
+      pmSourceFreePixels (sourceOut);
+
+      // allocate image, weight, mask for the new image for each peak
+      pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius);
+
+      // child sources have not been subtracted in this image, but this flag may be raised if
+      // they were subtracted in the parent's image
+      sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+
+      psArrayAdd (detectionsOut->allSources, 100, sourceOut);
+      psFree (sourceOut);
+    }
+    psLogMsg ("psphot", 3, "%ld known sources supplied", detectionsOut->allSources->n);
+
+
+    return true;
+}
+
+// create source children associated with 'filerule' from the objectsSrc.  XXX currently, this
+// is only used by psphotStackReadout (sources go on allSources so that psphotChoosePSF can be
+// called repeatedly)
+psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc) {
+
+    bool status;
+
+    int nImages = psphotFileruleCount(config, filerule);
+
+    // generate look-up arrays for detections and readouts
+    psArray *detArrays = psArrayAlloc(nImages);
+    psArray *readouts = psArrayAlloc(nImages);
+
+    for (int i = 0; i < nImages; i++) {
+
+        // find the currently selected readout
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // 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");
+        if (!detections) {
+            detections = pmDetectionsAlloc();
+            detections->allSources = psArrayAllocEmpty (100);
+            detections->peaks = psArrayAllocEmpty (100);
+            // save detections on the readout->analysis
+            if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detections)) {
+                psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout");
+                return NULL;
+            }
+            psFree(detections); // a copy remains on the analysis metadata
+            psAssert (detections, "missing detections?");
+        }
+        detArrays->data[i] = psMemIncrRefCounter(detections);
+        readouts->data[i] = psMemIncrRefCounter(readout);
+    }
+
+    psArray *objectsOut = psArrayAlloc(objectsSrc->n);
+
+    // copy all sources for each object
+    for (int k = 0; k < objectsSrc->n; k++) {
+
+        pmPhotObj *objectSrc = objectsSrc->data[k];
+        if (!objectSrc) continue;
+        if (!objectSrc->sources) continue;
+
+        pmPhotObj *objectOut = pmPhotObjAlloc();
+        objectsOut->data[k] = objectOut;
+
+        objectOut->SN = objectSrc->SN;
+        objectOut->x  = objectSrc->x;
+        objectOut->y  = objectSrc->y;
+        
+        objectOut->sources = psArrayAlloc(objectSrc->sources->n);
+
+        // copy the sources from sourceSrcs to the new detection structure
+        // loop over the sources, redefine their pixels to point at the new filerule image,
+        // copy the source data, and add a reference back to the original source
+        for (int i = 0; i < objectSrc->sources->n; i++) {
+
+            pmSource *sourceSrc = objectSrc->sources->data[i];
+
+            pmSource *sourceOut = pmSourceCopy(sourceSrc);
+            objectOut->sources->data[i] = sourceOut;
+
+            sourceOut->parent = sourceSrc;
+
+            // keep the original source flags
+            sourceOut->type     = sourceSrc->type;
+            sourceOut->mode     = sourceSrc->mode;
+            sourceOut->mode2    = sourceSrc->mode2;
+            sourceOut->tmpFlags = sourceSrc->tmpFlags;
+
+            // does this copy all model data? (NO)
+            sourceOut->modelPSF = pmModelCopy(sourceSrc->modelPSF);
+            sourceOut->modelEXT = pmModelCopy(sourceSrc->modelEXT);
+
+            if (sourceSrc->modelFits) {
+                sourceOut->modelFits = psArrayAlloc(sourceSrc->modelFits->n);
+                for (int j = 0; j < sourceSrc->modelFits->n; j++) {
+                    sourceOut->modelFits->data[j] = pmModelCopy(sourceSrc->modelFits->data[j]);
+                }
+            }
+
+            // drop the references to the original image pixels:
+            pmSourceFreePixels (sourceOut);
+
+            // set the output readotu
+            int index = sourceOut->imageID;
+            if (index >= readouts->n) continue; // skip the sources generated by the chisq image
+            pmReadout *readout = readouts->data[index];
+
+            // allocate image, weight, mask for the new image for each peak
+            pmSourceRedefinePixels (sourceOut, readout, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius);
+
+            // child sources have not been subtracted in this image, but this flag may be raised if
+            // they were subtracted in the parent's image
+            sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+
+            // set the output detections:
+            pmDetections *detectionsOut = detArrays->data[index];
+            psArrayAdd (detectionsOut->allSources, 100, sourceOut);
+            psArrayAdd (detectionsOut->peaks, 100, sourceOut->peak);
+        }
+    }
+
+    for (int i = 0; i < nImages; i++) {
+        pmDetections *detections = detArrays->data[i];
+        psLogMsg ("psphot", 3, "%ld source children for image %d", detections->allSources->n, i);
+    }
+
+    psFree (detArrays);
+    psFree (readouts);
+
+    return objectsOut;
+}
+

trunk/psphot/src/psphotOutput.c

@@ -15 +15 @@
     psFree (name);
     return num;
+}
+
+// convert filerule to filerule.NUM and look up in the config->arguments metadata
+bool psphotFileruleCountSet(const pmConfig *config, const char *filerule, int num) {
+
+    psString name = NULL;
+    psStringAppend(&name, "%s.NUM", filerule);
+
+    bool status = psMetadataAddS32(config->arguments, PS_LIST_TAIL, name, PS_META_REPLACE, "", num);
+    psFree (name);
+
+    return status;
 }
 
@@ -261 +273 @@
     psMetadataItemSupplement (&status, header, analysis, "ANGLE");
 
+    psMetadataItemSupplement (&status, header, analysis, "PSFMODEL");
+    psMetadataItemSupplement (&status, header, analysis, "PSF_OK");
+
     // Image Quality measurements
     psMetadataItemSupplement (&status, header, analysis, "IQ_NSTAR");

trunk/psphot/src/psphotPetrosianAnalysis.c

@@ -64 +64 @@
     }
 
-    psphotVisualShowResidualImage (readout);
+    psphotVisualShowResidualImage (readout, false);
     return true;
 }

trunk/psphot/src/psphotPetrosianProfile.c

@@ -68 +68 @@
     // XXX this will only work in the psphot context, not the psphotPetrosianStudy...
     // XXX add the petrosian to the pmSource structure...
-    // psphotVisualShowResidualImage (readout);
     psphotVisualShowPetrosian (source, petrosian);
 

trunk/psphot/src/psphotRadialApertures.c

@@ -102 +102 @@
         if (source->peak->x > AnalysisRegion.x1) continue;
         if (source->peak->y > AnalysisRegion.y1) continue;
+
+        // allocate pmSourceExtendedParameters, if not already defined
+        if (!source->radialAper) {
+            source->radialAper = psArrayAlloc(1);
+        }
 
         // replace object in image
@@ -116 +121 @@
         pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
 
-        if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax)) {
+        if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, 0)) {
             psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
         } else {
@@ -130 +135 @@
 }
 
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax) {
-
-    // allocate pmSourceExtendedParameters, if not already defined
-    if (!source->radial) {
-        source->radial = pmSourceRadialAperturesAlloc ();
+bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
+
+    // if we are a child source, save the results to the parent source radial aperture array
+    psArray *radialAperSet = source->radialAper;
+    if (source->parent) {
+        radialAperSet = source->parent->radialAper;
+    }
+    psAssert(radialAperSet, "this should be defined before calling");
+    psAssert(radialAperSet->data[entry] == NULL, "why is this already defined?");
+
+    pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc ();
+    radialAperSet->data[entry] = radialAper;
+
+    // storage for the derived pixel values
+    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
+
+    // outer-most radius for initial truncation
+    float Rmax  = aperRadii->data.F32[aperRadii->n - 1];
+    float Rmax2 = PS_SQR(Rmax);
+
+    // store the R^2 values for the apertures
+    psVector *aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32);
+    for (int i = 0; i < aperRadii->n; i++) {
+        aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]);
+    }
+
+    // center of the apertures
+    float xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
+    float yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
+
+    // one pass through the pixels to select the valid pixels and calculate R^2
+    for (int iy = 0; iy < source->pixels->numRows; iy++) {
+
+        float yDiff = iy - yCM;
+        if (fabs(yDiff) > Rmax) continue;
+
+        float *vPix = source->pixels->data.F32[iy];
+        float *vWgt = source->variance->data.F32[iy];
+        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy];
+
+        for (int ix = 0; ix < source->pixels->numCols; ix++, vPix++, vWgt++) {
+
+            if (vMsk) {
+                if (*vMsk & maskVal) {
+                    vMsk++;
+                    continue;
+                }
+                vMsk++;
+            }
+            if (isnan(*vPix)) continue;
+
+            float xDiff = ix - xCM;
+            if (fabs(xDiff) > Rmax) continue;
+
+            // radius is just a function of (xDiff, yDiff)
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            if (r2 > Rmax2) continue;
+
+            psVectorAppend(pixRadius2, r2);
+            psVectorAppend(pixFlux, *vPix);
+            psVectorAppend(pixVar, *vWgt);
+        }
+    }
+
+    psVector *flux    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
+    psVector *fluxErr = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
+    psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
+
+    psVectorInit (flux,    0.0);
+    psVectorInit (fluxErr, 0.0);
+    psVectorInit (fill,    0.0);
+
+    float *rPix2 = pixRadius2->data.F32;
+    for (int i = 0; i < pixRadius2->n; i++, rPix2++) {
+
+        float *aRad2 = aperRadii2->data.F32;
+        for (int j = 0; (*rPix2 < *aRad2) && (j < aperRadii2->n); j++, aRad2++) {
+            flux->data.F32[j]    += pixFlux->data.F32[i];
+            fluxErr->data.F32[j] += pixVar->data.F32[i];
+            fill->data.F32[j]    += 1.0;
+        }
+    }
+
+    for (int i = 0; i < flux->n; i++) {
+        // calculate the total flux for bin 'nOut'
+        float Area = M_PI*aperRadii2->data.F32[i];
+        fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
+        fill->data.F32[i] /= Area;
+        psTrace ("psphot", 5, "radial bins: %3d  %5.1f : %8.1f +/- %7.2f : %4.2f %6.1f\n", 
+                 i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], fill->data.F32[i], Area);
     }
     
-    psVector *radius  = psVectorAllocEmpty(100, PS_TYPE_F32);
+    radialAper->flux = flux;
+    radialAper->fluxErr = fluxErr;
+    radialAper->fill = fill;
+
+    psFree (aperRadii2);
+    psFree (pixRadius2);
+    psFree (pixFlux);
+    psFree (pixVar);
+
+    return true;
+}
+
+static int nCalls = 0;
+static int nPass = 0;
+static int nPix = 0;
+
+bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry) {
+
+    psAssert(source->radialAper->data[entry] == NULL, "why is this already defined?");
+
+    pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc ();
+    source->radialAper->data[entry] = radialAper;
+
+    psVector *pixRadius  = psVectorAllocEmpty(100, PS_TYPE_F32);
     psVector *pixFlux = psVectorAllocEmpty(100, PS_TYPE_F32);
     psVector *pixVar  = psVectorAllocEmpty(100, PS_TYPE_F32);
@@ -144 +259 @@
         for (int ix = 0; ix < source->pixels->numCols; ix++) {
 
-            // 0.5 PIX: get radius as a function of pixel coord
+            // 0.5 PIX: get pixRadius as a function of pixel coord
             float x = ix + 0.5 - source->peak->xf + source->pixels->col0;
             float y = iy + 0.5 - source->peak->yf + source->pixels->row0;
@@ -150 +265 @@
             float r = hypot(x, y);
 
-            psVectorAppend(radius, r);
+            psVectorAppend(pixRadius, r);
             psVectorAppend(pixFlux, source->pixels->data.F32[iy][ix]);
             psVectorAppend(pixVar, source->variance->data.F32[iy][ix]);
-        }
-    }
-    psphotRadialAperturesSortFlux(radius, pixFlux, pixVar);
+            nPix ++;
+            // if (nPix % 10000 == 0) {fprintf (stderr, "?");}
+        }
+    }
+    psphotRadialAperturesSortFlux(pixRadius, pixFlux, pixVar);
 
     psVector *flux    = psVectorAllocEmpty(radMax->n, PS_TYPE_F32); // surface brightness of radial bin
@@ -174 +291 @@
 
     // XXX assume (or enforce) that the bins are contiguous and non-overlapping (Rmax[i] = Rmin[i+1])
-    for (int i = 0; !done && (i < radius->n); i++) {
-        if (radius->data.F32[i] > Rmax) {
+    for (int i = 0; !done && (i < pixRadius->n); i++) {
+        if (pixRadius->data.F32[i] > Rmax) {
             // calculate the total flux for bin 'nOut'
             float Area = M_PI*PS_SQR(Rmax);
@@ -185 +302 @@
                      nOut, radMax->data.F32[nOut], flux->data.F32[nOut], fluxErr->data.F32[nOut], fill->data.F32[nOut], Area);
 
+            nPass ++;
+            // if (nPass % 1000 == 0) {fprintf (stderr, "!");}
+
             nOut ++;
             if (nOut >= radMax->n) break;
@@ -195 +315 @@
     flux->n = fluxErr->n = fill->n = nOut;
     
-    psFree(source->radial->flux);
-    psFree(source->radial->fluxErr);
-    psFree(source->radial->fill);
-
-    source->radial->flux = flux;
-    source->radial->fluxErr = fluxErr;
-    source->radial->fill = fill;
-
-    psFree (radius);
+    radialAper->flux = flux;
+    radialAper->fluxErr = fluxErr;
+    radialAper->fill = fill;
+
+    psFree (pixRadius);
     psFree (pixFlux);
     psFree (pixVar);
 
+    nCalls ++;
+    // if (nCalls % 100 == 0) {fprintf (stderr, "*");}
     return true;
 }

trunk/psphot/src/psphotRadialAperturesByObject.c

@@ -3 +3 @@
 // aperture-like measurements for extended sources
 // flux in simple, circular apertures
-bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule) {
+
+// **** it looks like this function will re-point the source pixels at the specified FILERULE
+// **** I need to distinguish PSF-matched images from raw
+// **** save (somewhere) the PSF-matched PSF values
+
+// this function measures the radial aperture fluxes for the set of readouts.  this function
+// may be called multiple times (presumably with different matched PSF sizes).  we must have
+// already added an entry to the readout->analysis identifying the FWHM of this version.
+
+bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule, int nMatchedPSF) {
 
     bool status;
@@ -28 +37 @@
     psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
     psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
+    float outerRadius = radMax->data.F32[radMax->n - 1];
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -39 +49 @@
     float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
 
+    // how many target PSFs do we want?
+    int nPSFsizes = 0;
+    { 
+        psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
+        if (status) {
+            nPSFsizes = 1;
+        } else {
+            psVector *fwhmValues = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
+            psAssert (status, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
+            nPSFsizes = fwhmValues->n;
+        }
+    }
+    
     // source analysis is done in S/N order (brightest first)
     objects = psArraySort (objects, pmPhotObjSortBySN);
@@ -52 +75 @@
         pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
         psAssert (readout, "missing readout?");
+
+        psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
+        if (!fwhmValues) {
+            psError (PSPHOT_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout");
+            return false;
+        }
+        if (fwhmValues->n != nMatchedPSF + 1) {
+            psError (PSPHOT_ERR_CONFIG, true, "convolved or measured FWHM sequence is inconsistent this readout");
+            return false;
+        }
+        psLogMsg ("psphot", PS_LOG_DETAIL, "PSF FWHM of %s : %f pixels\n", file->name, fwhmValues->data.F32[nMatchedPSF]);
 
         readouts->data[i] = psMemIncrRefCounter(readout);
@@ -77 +111 @@
             if (source->peak->SN < SN_LIM) continue;
 
+            int index = source->imageID;
+            if (index >= readouts->n) continue; // skip the sources generated by the chisq image
+            pmReadout *readout = readouts->data[index];
+
+            // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
+            // psphotVisualShowImage(readout);
+
+            // allocate pmSourceExtendedParameters, if not already defined
+            if (source->parent) {
+                if (!source->parent->radialAper) {
+                    source->parent->radialAper = psArrayAlloc(nPSFsizes);
+                }
+            } else {
+                if (!source->radialAper) {
+                    source->radialAper = psArrayAlloc(nPSFsizes);
+                }
+            }
+
             // replace object in image
             if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
                 pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
             }
+
+            // we need to change the view for the radial aperture analysis, but we want to recover exactly
+            // the original view; the following elements get destroyed by pmSourceRedefinePixels so save them:
+            psImage *oldMaskObj   = psMemIncrRefCounter(source->maskObj);
+            psImage *oldModelFlux = psMemIncrRefCounter(source->modelFlux);
+            psImage *oldPSFimage  = psMemIncrRefCounter(source->psfImage);
+            psRegion oldRegion    = source->region;
+
             Nradial ++;
 
-            int index = source->imageID;
-            pmReadout *readout = readouts->data[index];
+            // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
+            // psphotVisualShowImage(readout);
 
             // force source image to be a bit larger...
-            float radius = source->peak->xf - source->pixels->col0;
-            radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
-            radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
-            radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
-            pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
+            // float radius = source->peak->xf - source->pixels->col0;
+            // radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
+            // radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
+            // radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
+            pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
 
-            if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax)) {
+            if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, nMatchedPSF)) {
                 psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
             } else {
@@ -99 +159 @@
             }
 
+            pmSourceRedefinePixelsByRegion (source, readout, oldRegion);
+            psFree(source->maskObj);   source->maskObj   = oldMaskObj;
+            psFree(source->modelFlux); source->modelFlux = oldModelFlux;
+            psFree(source->psfImage);  source->psfImage  = oldPSFimage;
+
             // re-subtract the object, leave local sky
             pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+
+            // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
+            // psphotVisualShowImage(readout);
         }
     }

trunk/psphot/src/psphotReadout.c

@@ -112 +112 @@
     // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
     // this step is skipped
-    if (!psphotChoosePSF (config, view, filerule)) { // pass 1
+    if (!psphotChoosePSF (config, view, filerule, true)) { // pass 1
         psLogMsg ("psphot", 3, "failure to construct a psf model");
         return psphotReadoutCleanup (config, view, filerule);

trunk/psphot/src/psphotReadoutFindPSF.c

@@ -49 +49 @@
     }
 
-    if (!psphotChoosePSF(config, view, filerule)) {
+    if (!psphotChoosePSF(config, view, filerule, true)) {
         psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model");
         return psphotReadoutCleanup (config, view, filerule);

trunk/psphot/src/psphotReadoutKnownSources.c

@@ -43 +43 @@
     }
 
-    if (!psphotChoosePSF (config, view, filerule)) {
+    if (!psphotChoosePSF (config, view, filerule, true)) {
         psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model");
         return psphotReadoutCleanup (config, view, filerule);

trunk/psphot/src/psphotReplaceUnfit.c

@@ -75 +75 @@
       pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     }
+
+    psphotVisualShowImage(readout);
     psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
     return true;
@@ -101 +103 @@
     return true;
 }
+
+// modify the sources to point at the corresponding pixels for the given filerule
+bool psphotRedefinePixels (pmConfig *config, const pmFPAview *view, const char *filerule)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (!psphotRedefinePixelsReadout (config, view, filerule, i, recipe)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i);
+            return false;
+        }
+    }
+    return true;
+}
+
+bool psphotRedefinePixelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+
+    bool status;
+    pmSource *source;
+
+    psTimerStart ("psphot.replace");
+
+    // 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?");
+
+    // XXX the sources have already been copied (merge into here?)
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->allSources;
+    psAssert (sources, "missing sources?");
+
+    // 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
+    psAssert (maskVal, "missing mask value?");
+
+    for (int i = 0; i < sources->n; i++) {
+      source = sources->data[i];
+
+      // sources have not yet been subtracted in this image (but this flag may be raised)
+      source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+
+      float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+      float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+      float radius = source->modelPSF->fitRadius;
+
+      // force a redefine to this image
+      pmSourceFreePixels(source);
+      pmSourceRedefinePixels (source, readout, Xo, Yo, radius);
+    }
+    return true;
+}
+
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotResetModels (pmConfig *config, const pmFPAview *view, const char *filerule)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (!psphotResetModelsReadout (config, view, filerule, i, recipe)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i);
+            return false;
+        }
+    }
+    return true;
+}
+
+bool psphotResetModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+
+    bool status;
+    pmSource *source;
+
+    psTimerStart ("psphot.replace");
+
+    // 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?");
+
+    // XXX the sources have already been copied (merge into here?)
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->allSources;
+    psAssert (sources, "missing sources?");
+
+    // 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
+    psAssert (maskVal, "missing mask value?");
+
+    // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
+    int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
+    assert (status);
+
+    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+    psAssert (psf, "missing psf?");
+
+    for (int i = 0; i < sources->n; i++) {
+      source = sources->data[i];
+
+      // *** we need to cache the 'best' model, and we have 3 cases:
+      // 1) model is the psf model --> generate from the new psf
+      // 2) model is an unconvolved extended model --> just cache the copy (not perfect)
+      // 3) model is a convolved extended model --> re-generate
+
+      // use the 'best' model to cache the model (PSF or EXT : EXT may point at one of modelFits
+      bool isPSF = false;
+      pmModel *model = pmSourceGetModel(&isPSF, source);
+      float radius = model->fitRadius; // save for future use below
+
+      // regenerate the PSF if the model is a PSF, or if we need the PSF for a PCM
+      if (isPSF || model->isPCM) {
+          // the guess central intensity comes from the peak:
+          float Io = source->peak->flux;
+          float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+          float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+
+          // generate a model for this object with Io = 1.0
+          pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io);
+          if (modelPSF == NULL) {
+              psWarning ("Failed to determine PSF model for source at (%f,%f), skipping", Xo, Yo);
+              continue;
+          }
+
+          // set the source PSF model
+          psFree (source->modelPSF);
+          source->modelPSF = modelPSF;
+          source->modelPSF->fitRadius = radius;
+      }
+
+      if (model->isPCM) {
+          psAssert(false, "this section is not complete");
+
+          pmSourceCachePSF (source, maskVal);
+
+          psKernel *psfKernel = pmPCMkernelFromPSF(source, psfSize);
+          if (!psfKernel) { 
+              psWarning ("no psf kernel");
+          }
+
+          // generate an image of the right size
+          psImage *rawModelFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
+          psImageInit (rawModelFlux, 0.0);
+          
+          // insert the model image normalized to 1.0
+          pmModelAdd (rawModelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
+
+          psImageConvolveFFT (source->modelFlux, rawModelFlux, NULL, 0, psfKernel);
+          
+          psFree (psfKernel);
+          psFree (rawModelFlux);
+      }
+
+      pmSourceCacheModel (source, maskVal);  // ALLOC x14 (!)
+    }
+
+    psLogMsg ("psphot.replace", PS_LOG_INFO, "subtracted models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
+    return true;
+}
+

trunk/psphot/src/psphotSetThreads.c

@@ -35 +35 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 12);
+    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 13);
     task->function = &psphotExtendedSourceFits_Threaded;
     psThreadTaskAdd(task);

trunk/psphot/src/psphotSourceFits.c

@@ -365 +365 @@
 
     // copy most data from the primary source (modelEXT, blends stay NULL)
-    pmSource *newSrc = pmSourceCopyData (source);
+    pmSource *newSrc = pmSourceCopy (source);
     newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
 
@@ -385 +385 @@
             psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelPSF->params->data.F32[i], source->modelPSF->dparams->data.F32[i]);
         }
-        psphotVisualShowResidualImage (readout);
+        psphotVisualShowResidualImage (readout, false);
     }
 # endif

trunk/psphot/src/psphotSourceMatch.c

@@ -187 +187 @@
         pmPhotObj *obj = objects->data[i]; 
 
+        // we will find the input source with the max number of spans and reproduce that footprint
+        int nSpansMax = 0;
+        int iSpansMax = -1;
+
         // mark the images for which sources have been found
         psVectorInit (found, 0);
@@ -196 +200 @@
             psAssert (index < found->n, "invalid index");
 
+            if (src->peak && src->peak->footprint && src->peak->footprint->nspans > nSpansMax) {
+                nSpansMax = src->peak->footprint->nspans;
+                iSpansMax = j;
+            }
+
             found->data.U8[index] = 1;
+        }
+
+        // we make a copy of the largest footprint; this will be used for all new sources associated with this object
+        pmFootprint *footprint = NULL;
+        if (iSpansMax != -1) { // copy the footprint info
+            pmSource *src = obj->sources->data[iSpansMax]; 
+            psAssert(src->peak, "source does not exist?");
+            psAssert(src->peak->footprint, "footprint does not exist");
+            psAssert(src->peak->footprint->nspans == nSpansMax, "wrong footprint?");
+            
+            // we only care about the spans, do not worry about the image of this footprint
+            footprint = pmFootprintCopyData(src->peak->footprint, NULL);
         }
 
@@ -219 +240 @@
             peak->dy = NAN;
             
-            // XXX assign to a footprint?
-
+            // assign to a footprint on this readout->image
+            peak->footprint = pmFootprintCopyData(footprint, readout->image);
+
+            // the peak does not claim ownership of the footprint (it does not free it). save a copy of this 
+            // footprint on detections->footprints so we can free it later
+            psArrayAdd(detections->footprints, 100, peak->footprint); 
+            psFree (peak->footprint);
+            
             // create a new source
             pmSource *source = pmSourceAlloc();
@@ -227 +254 @@
 
             // add the peak
-            source->peak = psMemIncrRefCounter(peak);
+            source->peak = peak;
 
             // allocate space for moments
@@ -239 +266 @@
             psArrayAdd (detections->newSources, 100, source);
             psFree (source);
-            psFree (peak);
-        }
+        }
+        psFree (footprint);
     }
 

trunk/psphot/src/psphotSourceSize.c

@@ -499 +499 @@
 
         // clear the mask bit and set the circular mask pixels
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
-
+        // psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
+        // psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
         pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
 
         // clear the mask bit
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
+        // psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
 
         // re-subtract the object, leave local sky

trunk/psphot/src/psphotStackArguments.c

@@ -41 +41 @@
     if ((N = psArgumentGet (argc, argv, "-break"))) {
         if (argc <= N+1) {
-          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
+          psErrorStackPrint(stderr, "Expected to see an argument for -break");
           exit(PS_EXIT_CONFIG_ERROR);
         }

trunk/psphot/src/psphotStackChisqImage.c

@@ -6 +6 @@
 
 // XXX supply filename or keep PSPHOT.INPUT fixed?
-bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view, const char *ruleDet, const char *ruleCnv)
+bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view, const char *ruleDet, const char *ruleSrc)
 {
     psTimerStart ("psphot.chisq.image");
@@ -27 +27 @@
 
     psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.CHISQ.NUM", PS_META_REPLACE, "", num);
+
+    // we need to increment the counter for ruleDet and ruleSrc:
     num++;
-    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", num);
 
     // save the resulting image in the 'detection' set
@@ -35 +36 @@
         return false;
     }
+    psphotFileruleCountSet(config, ruleDet, num);
 
-    // save the resulting image in the 'convolved' set
-    if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, ruleCnv, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) {
-        psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files");
-        return false;
+    // also save the resulting image in the 'source' set (analysis set)
+    if (strcmp(ruleDet, ruleSrc)) {
+        if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, ruleSrc, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) {
+            psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files");
+            return false;
+        }
+        psphotFileruleCountSet(config, ruleSrc, num);
     }
 
@@ -119 +124 @@
     psAssert (status, "programming error: must define PSPHOT.CHISQ.NUM");
 
-    int inputNum = psphotFileruleCount(config, "PSPHOT.INPUT");
+    int inputNum = psphotFileruleCount(config, filerule);
 
     pmFPAfileRemoveSingle (config->files, filerule, chisqNum);
 
     inputNum --;
-    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", inputNum);
+    psphotFileruleCountSet(config, filerule, inputNum);
 
     return true;

trunk/psphot/src/psphotStackImageLoop.c

@@ -101 +101 @@
 */
 
-# define UPDATE_HEADER 0
+# define UPDATE_HEADER 1
 
 bool GetAstrometryFPA (pmConfig *config, pmFPAview *view) {
@@ -174 +174 @@
         pmChip *outChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output
 
-# if (UPDATE_HEADER)
         pmHDU *outHDU = pmFPAviewThisHDU (view, output->fpa);
         if (!outHDU) {
-            pmFPAAddSourceFromView(output->fpa, "name", view, output->format);
+            pmFPAAddSourceFromView(output->fpa, view, output->format);
             outHDU = pmFPAviewThisHDU (view, output->fpa);
             psAssert (outHDU, "failed to make HDU");
         }
-# endif
+        if (!outHDU->header) {
+          outHDU->header = psMetadataAlloc();
+        }
 
         if (bilevelAstrometry) {
@@ -188 +189 @@
                 continue;
             }
-# if (UPDATE_HEADER)
             if (!pmAstromWriteBilevelChip(outHDU->header, outChip, WCS_NONLIN_TOL)) {
                 psWarning("Unable to generate WCS header.");
                 continue;
             }
-# endif
         } else {
             // we use a default FPA pixel scale of 1.0
@@ -200 +199 @@
                 continue;
             }
-# if (UPDATE_HEADER)
-            if (UPDATE_HEADER && !pmAstromWriteWCS(outHDU->header, output->fpa, outChip, WCS_NONLIN_TOL)) {
+            if (!pmAstromWriteWCS(outHDU->header, output->fpa, outChip, WCS_NONLIN_TOL)) {
                 psWarning("Unable to generate WCS header.");
                 continue;
             }
-# endif
         }
     }
@@ -225 +222 @@
         psAssert (output, "missing file?");
 
-# if (UPDATE_HEADER)
         pmHDU *PHU = pmFPAviewThisPHU(view, output->fpa);
         if (!PHU) {
-            pmFPAAddSourceFromView(output->fpa, "name", view, output->format);
+            pmFPAAddSourceFromView(output->fpa, view, output->format);
             PHU = pmFPAviewThisPHU (view, output->fpa);
             psAssert (PHU, "failed to make PHU");
         }
+        if (!PHU->header) {
+          PHU->header = psMetadataAlloc();
+        }
 
         if (!pmAstromWriteBilevelMosaic(PHU->header, output->fpa, WCS_NONLIN_TOL)) {
             psWarning("Unable to generate WCS header.");
         }
-# endif
     }
 

trunk/psphot/src/psphotStackMatchPSFs.c

@@ -65 +65 @@
 bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) {
 
+    psImageMaskType maskValue;
+    psImageMaskType markValue;
+
+    // get the PSPHOT.MASK value from the config
+    if (!pmConfigMaskSetBits (&maskValue, &markValue, config)) {
+        psError (PS_ERR_UNKNOWN, true, "Unable to define the mask bit values");
+        return false;
+    }
+
     pmFPAfile *fileSrc = psphotStackGetConvolveSource(config, options, index);
     if (!fileSrc) {
@@ -84 +93 @@
 
     // set NAN pixels to 'SAT'
-    // XXX replace this is pmReadoutMaskInvalid?
-    psImageMaskType maskVal = pmConfigMaskGet("SAT", config);
-    if (!pmReadoutMaskNonfinite(readoutSrc, maskVal)) {
+    psImageMaskType maskSat = pmConfigMaskGet("SAT", config);
+    if (!pmReadoutMaskInvalid(readoutSrc, maskValue, maskSat)) {
         psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in readout.");
         return false;
@@ -95 +103 @@
         matchKernel(config, readoutOut, readoutSrc, options, index);
         saveMatchData(readoutOut, options, index);
-        // renormKernel(readoutCnv, options, index);
-    } else {
-        // only match the flux (NO! not for multi-filter, at least!)
-        // XXX do not generate readoutCnv in this case?
-        // float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation
-        // psBinaryOp(readoutRaw->image, readoutRaw->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
-        // psBinaryOp(readoutRaw->variance, readoutRaw->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32));
     }
-
     rescaleData(readoutOut, config, options, index);
 
-    // dumpImage(readoutOut, readoutSrc, index, "convolved");
+    // save the output fwhm values in the readout->analysis.  we may have / will have multiple output PSF sizes,
+    // so we save this in a vector.  if the vector is not yet defined, create it
+    bool mdok = false;
+    psVector *fwhmValues = psMetadataLookupVector(&mdok, readoutOut->analysis, "STACK.PSF.FWHM.VALUES");
+    if (!fwhmValues) {
+        fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
+        psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
+        psFree(fwhmValues); // drops the extra copy
+    }
+    psVectorAppend(fwhmValues, options->targetSeeing);
 
     return true;
 }
-
-
-# if (0)
-// Read previously produced kernel
-if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
-    loadKernel(config, readoutCnv, options, index);
-} else {
-    matchKernel(config, readoutCnv, readoutRaw, options, index);
-}
-# endif

trunk/psphot/src/psphotStackMatchPSFsUtils.c

@@ -365 +365 @@
         pmSubtractionSetFWHMs(options->inputSeeing->data.F32[index], options->targetSeeing);
 
-        if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, scaleRef, scaleMin, scaleMax)) {
-            psError(psErrorCodeLast(), false, "Unable to scale kernel parameters");
-            goto escape;
-        }
+        pmSubtractionParamScaleOptions(scale, scaleRef, scaleMin, scaleMax);
+
+        // if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, scaleRef, scaleMin, scaleMax)) {
+        //     psError(psErrorCodeLast(), false, "Unable to scale kernel parameters");
+        //     goto escape;
+        // }
 
         if (!pmSubtractionMatch(NULL, readoutOut, fake, readoutSrc, footprint, stride, regionSize, spacing, threshold, stampSources, stampsName, type, size, order, widthsCopy, orders, inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) {
@@ -383 +385 @@
         goto escape;
     }
+
+    // save the PSF on the new readout->analysis:
+    // if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot psf model", options->psf)) {
+    //     psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout");
+    //     return false;
+    // }
 
     // dumpImage(readoutOut, readoutSrc, index, "conv");

trunk/psphot/src/psphotStackPSF.c

@@ -12 +12 @@
     bool autoPSF = psMetadataLookupBool (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.AUTO");
 
+    // Get the recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe
+    psAssert(recipe, "We've thrown an error on this before.");
+        
+    char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF
+
     if (autoPSF) {
-        // Get the recipe values
-        psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe
-        psAssert(recipe, "We've thrown an error on this before.");
-
         int psfInstances = psMetadataLookupS32(NULL, recipe, "PSF.INSTANCES"); // Number of instances for PSF
         float psfRadius = psMetadataLookupF32(NULL, recipe, "PSF.RADIUS"); // Radius for PSF
-        const char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF
         int psfOrder = psMetadataLookupS32(NULL, recipe, "PSF.ORDER"); // Spatial order for PSF
 
@@ -49 +50 @@
 
         float targetFWHM = psMetadataLookupF32 (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
-        psAssert (isfinite(targetFWHM), "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
+        if (!mdok) {
+            psVector *fwhmValues = psMetadataLookupVector(&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
+            psAssert (mdok, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
+            targetFWHM = fwhmValues->data.F32[0];
+        }
 
         float Sxx = sqrt(2.0)*targetFWHM / 2.35;
 
         // XXX probably should make the model type (and par 7) optional from recipe
-        psf = pmPSFBuildSimple("PS_MODEL_PS1_V1", Sxx, Sxx, 0.0, 1.0);
+        psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0);
         if (!psf) {
             psError(PSPHOT_ERR_PSF, false, "Unable to build dummy PSF.");

trunk/psphot/src/psphotStackParseCamera.c

@@ -15 +15 @@
     }
 
+    int nRaw = 0;
+    int nCnv = 0;
     int nInputs = inputs->list->n;
     for (int i = 0; i < nInputs; i++) {
@@ -57 +59 @@
                 }
             }
+            nRaw ++;
         }
 
@@ -88 +91 @@
                 }
             }
+            nCnv ++;
         }
 
@@ -95 +99 @@
         }
 
+        // XXX what if they do not match in length
+        if (nCnv && nRaw) {
+            if (nCnv != nRaw) {
+                psError (PSPHOT_ERR_CONFIG, true, "if both RAW and CNV images are supplied, the number must match");
+                return false;
+            }
+        }
+
         psString sources = psMetadataLookupStr(&status, input, "SOURCES"); // Name of mask
-        // pmFPAfile *srcInputFile = rawInputFile ? rawInputFile : cnvInputFile;
         if (sources && strlen(sources) > 0) {
             if (!defineFile(config, NULL, "PSPHOT.STACK.SOURCES", sources, PM_FPA_FILE_CMF)) {
@@ -107 +118 @@
         // XXX output of these files should be optional
         {
+            // pmFPAfile *srcInputFile = rawInputFile ? rawInputFile : cnvInputFile;
             pmFPAfile *outputImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.STACK.OUTPUT.IMAGE");
             if (!outputImage) {
@@ -150 +162 @@
     }
     psMetadataRemoveKey(config->arguments, "FILENAMES");
+    psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.STACK.INPUT.RAW.NUM", PS_META_REPLACE, "number of inputs", nRaw);
+    psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.STACK.INPUT.CNV.NUM", PS_META_REPLACE, "number of inputs", nCnv);
+    psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.STACK.OUTPUT.IMAGE.NUM", PS_META_REPLACE, "number of inputs", nInputs);
     psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", nInputs);
 

trunk/psphot/src/psphotStackReadout.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
+// we have 3 possible real filesets:
 # define STACK_RAW "PSPHOT.STACK.INPUT.RAW"
-# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"
+# define STACK_CNV "PSPHOT.STACK.INPUT.CNV"
+# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"  /* the psf-matched image */
+
+// we have 3 files on which we operate:
+// DET (detection image)       : nominally RAW (optionally CNV?)
+// SRC (source analysis image) : nominally CNV (optionally RAW)
+// OUT (psf-matched images)    : always OUT
+
+bool psphotStackVisualFilerule(pmConfig *config, const pmFPAview *view, const char *filerule) {
+
+    bool status = false;
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+
+        // find the currently selected readout
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // 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?");
+
+        psphotVisualShowResidualImage (readout, true);
+        psphotVisualShowObjectRegions (readout, recipe, sources);
+    }
+    return true;
+}
 
 bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
@@ -20 +58 @@
     PS_ASSERT_PTR_NON_NULL (breakPt, false);
 
+    // XXX or do I set OUT to be a pmFPAfile pointing at the input of interest?
+    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
+    char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV;
+    char *STACK_DET = STACK_RAW; // XXX optionally allow this to be CNV?
+
     // we have 3 relevant files: RAW, CNV, OUT 
 
     // set the photcode for each image
-    if (!psphotAddPhotcode (config, view, STACK_OUT)) {
+    if (!psphotAddPhotcode (config, view, STACK_SRC)) {
         psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
         return false;
@@ -30 +73 @@
     // Generate the mask and weight images
     // XXX this should be done before we perform the convolutions
-    if (!psphotSetMaskAndVariance (config, view, STACK_RAW)) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+    if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "NOTHING")) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
@@ -40 +83 @@
     // XXX I think this is not defined correctly for an array of images.
     // XXX probably need to subtract the model (same model?) for both RAW and OUT
-    if (!psphotModelBackground (config, view, STACK_RAW)) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
-    }
-    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+    if (!psphotModelBackground (config, view, STACK_DET)) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+    if (!psphotSubtractBackground (config, view, STACK_DET)) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "BACKMDL")) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
-    }
-
-    // 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
-    if (!psphotLoadPSF (config, view, STACK_RAW)) {
-        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
-    }
-
-    if (!psphotStackChisqImage(config, view, STACK_RAW, STACK_OUT)) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "CHISQ")) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     // find the detections (by peak and/or footprint) in the image.
     // This finds the detections on Chisq image as well as the individuals
-    if (!psphotFindDetections (config, view, STACK_RAW, true)) { // pass 1
+    if (!psphotFindDetections (config, view, STACK_DET, true)) { // pass 1
         // this only happens if we had an error in psphotFindDetections
         psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
-    }
-
-    // copy the detections from RAW to OUT
-    if (!psphotCopySources (config, view, STACK_OUT, STACK_RAW)) {
-        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    // copy the detections from DET to SRC 
+    if (strcmp(STACK_SRC, STACK_DET)) {
+        if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
+            psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
+            return psphotReadoutCleanup (config, view, STACK_SRC);
+        }
     }
 
     // construct sources and measure basic stats (saved on detections->newSources)
-    // only run this on detections from the input images, not chisq image
-    if (!psphotSourceStats (config, view, STACK_OUT, true)) { // pass 1
+    if (!psphotSourceStats (config, view, STACK_SRC, true)) { // pass 1
         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    if (!strcasecmp (breakPt, "TEST1")) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     // generate the objects (object unify the sources from the different images)
-    psArray *objects = psphotMatchSources (config, view, STACK_OUT);
+    // XXX this could just match the detections for the chisq image, and not bother measuring the
+    // source stats in that case...
+    psArray *objects = psphotMatchSources (config, view, STACK_SRC);
+
+    if (!strcasecmp (breakPt, "TEST2")) {
+        psFree(objects);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
 
     // construct sources for the newly-generated sources (from other images)
-    if (!psphotSourceStats (config, view, STACK_OUT, false)) { // pass 1
+    if (!psphotSourceStats (config, view, STACK_SRC, false)) { // pass 1
+        psFree(objects);
         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
@@ -98 +147 @@
     // if (!psphotDeblendSatstars (config, view)) {
     //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
-    //     return psphotReadoutCleanup (config, view, STACK_OUT);
+    //     return psphotReadoutCleanup (config, view, STACK_SRC);
     // }
 
@@ -104 +153 @@
     // if (!psphotBasicDeblend (config, view)) {
     //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
-    //     return psphotReadoutCleanup (config, view, STACK_OUT);
+    //     return psphotReadoutCleanup (config, view, STACK_SRC);
     // }
 
     // classify sources based on moments, brightness
     // only run this on detections from the input images, not chisq image
-    if (!psphotRoughClass (config, view, STACK_OUT)) {
+    if (!psphotRoughClass (config, view, STACK_SRC)) {
+        psFree(objects);
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     // only run this on detections from the input images, not chisq image
-    if (!psphotImageQuality (config, view, STACK_OUT)) { // pass 1
+    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
+        psFree(objects);
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "MOMENTS")) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        psFree(objects);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
     // this step is skipped
-    if (!psphotChoosePSF (config, view, STACK_OUT)) { // pass 1
+    if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
+        psFree(objects);
         psLogMsg ("psphot", 3, "failure to construct a psf model");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "PSFMODEL")) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        psFree(objects);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
-    psphotGuessModels (config, view, STACK_OUT);
+    psphotGuessModels (config, view, STACK_SRC);
 
     // merge the newly selected sources into the existing list
     // NOTE: merge OLD and NEW
-    psphotMergeSources (config, view, STACK_OUT);
+    psphotMergeSources (config, view, STACK_SRC);
 
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
     psphotFitSourcesLinearStack (config, objects, FALSE);
+    psphotStackVisualFilerule(config, view, STACK_SRC);
 
     // identify CRs and extended sources
-    psphotSourceSize (config, view, STACK_OUT, TRUE);
+    psphotSourceSize (config, view, STACK_SRC, TRUE);
+
+    // XXX do we want to do a preliminary (unconvolved) model fit here, and then 
+    // do a second detection pass? (like standard psphot)
 
     // measure aperture photometry corrections
-    if (!psphotApResid (config, view, STACK_OUT)) {
+    if (!psphotApResid (config, view, STACK_SRC)) {
         psFree (objects);
         psLogMsg ("psphot", 3, "failed on psphotApResid");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     psphotStackObjectsUnifyPosition (objects);
 
-    // measure circular, radial apertures (objects sorted by S/N)
-    psphotRadialAperturesByObject (config, objects, view, STACK_OUT); 
-
     // measure elliptical apertures, petrosians (objects sorted by S/N)
-    psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_OUT); // pass 1 (detections->allSources)
+    psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
 
     // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
-    psphotExtendedSourceFits (config, view, STACK_OUT); // pass 1 (detections->allSources)
+    psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
 
     // calculate source magnitudes
-    psphotMagnitudes(config, view, STACK_OUT);
+    psphotMagnitudes(config, view, STACK_SRC);
+
+    // create source children for the OUT filerule (for radial aperture photometry) 
+    psArray *objectsRadial = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
+    if (!objectsRadial) {
+        psFree(objects);
+        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    bool smoothAgain = true;
+    for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) {
+
+        // re-measure the PSF for the smoothed image (using entries in 'allSources')
+        psphotChoosePSF (config, view, STACK_OUT, false);
+
+        // this is necessary to update the models based on the new PSF
+        psphotResetModels (config, view, STACK_OUT);
+
+        // this is necessary to get the right normalization for the new models
+        psphotFitSourcesLinear (config, view, STACK_OUT, false);
+
+        // measure circular, radial apertures (objects sorted by S/N)
+        psphotRadialAperturesByObject (config, objectsRadial, view, STACK_OUT, nMatchedPSF); 
+
+        // replace the flux in the image so it is returned to its original state
+        psphotReplaceAllSources (config, view, STACK_OUT);
+
+        // smooth to the next FWHM, or set 'smoothAgain' to false if no more 
+        psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF);
+    }
 
     if (0 && !psphotEfficiency(config, view, STACK_OUT)) {
@@ -176 +261 @@
 
     // replace background in residual image
-    psphotSkyReplace (config, view, STACK_RAW);
+    psphotSkyReplace (config, view, STACK_DET);
 
     // drop the references to the image pixels held by each source
-    psphotSourceFreePixels (config, view, STACK_OUT);
-
-    // remove chisq image from config->file:PSPHOT.INPUT (why?)
-    psphotStackRemoveChisqFromInputs(config, STACK_RAW);
+    // psphotSourceFreePixels (config, view, STACK_OUT);
+    psphotSourceFreePixels (config, view, STACK_SRC);
+
+    // remove chisq image from config->file:PSPHOT.INPUT
+    psphotStackRemoveChisqFromInputs(config, STACK_DET);
+    if (strcmp(STACK_SRC, STACK_DET)) {
+        psphotStackRemoveChisqFromInputs(config, STACK_SRC);
+    }
 
     psFree (objects);
+    psFree (objectsRadial);
 
     // create the exported-metadata and free local data
-    return psphotReadoutCleanup (config, view, STACK_OUT);
+    return psphotReadoutCleanup (config, view, STACK_SRC);
 }
 
+/* here is the process:
+
+ * we have three(*) images:
+ * RAW : unconvolved image stack 
+ * CNV : input convolved image
+
+ * OUT : psf-matched output image (there may be more than one of
+ * these.  we will generate the first matched image by selecting the
+ * target PSF and doing a full psf-maching process (as used by ppStack
+ * and ppSub).  But, additional target output files should use a
+ * simple gaussian convolution kernel determind from therms of the
+ * current and the target).
+
+ * the output should be / could be one of the matched images, but not
+ * all.  should we ensure the first gets written out, and ot save the
+ * others (or only optionally).
+
+ * by default, we probably only sve the cmf ffile outputs.
+
+ * load the RAW image (unconvolved stacks)
+ * add photcode to the output headers / readout->analysis
+ * generate mask and variance image (this is probably never needed in
+   practice: we always load an input mask & var.
+ * generate & subtract a model background for ?? (RAW? CNV? OUT? all?)
+ * load a PSF (probably not yet working)
+
+ * generate the CHISQ image from the RAW input images (why save on OUT?)
+
+ * find detections on RAW
+
+ * copy detections to OUT
+
+ * generate source stats (moments) for OUT
+
+ * match sources across inputs (on OUT?)
+
+ * generate source stats for the new constructions
+
+ * rough class (star, galaxy, cosmic, etc)
+
+ * Image quality
+
+ * generate PSF
+
+ * guess models
+
+ * merge sources (new -> old)
+
+ * linear fit to the psf
+
+ * find ApResid
+
+ * assign common positions
+
+ * radial apertures (** this should be on the PSF-matched images
+
+ * extended analysis (elliptical profile & petrosian)
+
+ * extended fits (sersic, etc)
+
+ * psphot magnitudes
+
+
+ ******
+
+ the above is all wrong:  first, we should be doing the full
+ morphology analysis (ExtendedAnalysis & ExtendedFits) on the CNV or
+ RAW image (as desired optionally), etc.
+
+ In the discussion below, 'BST' (best) means optionally RAW or CNV
+
+ * detection : RAW & CHISQ (of RAW)
+ * moments : used by psf analysis & classification (BST)
+ * rough class : uses moments, not pixels
+ * image quality : uses moments as well
+ * generate PSF : (BST)
+ * guess models (BST)
+ * linear fit (BST)
+ * find ApResid (BST) -- uses sources not pixels
+ * extended analysis (BST)
+ * extended fits (BST)
+ * detection efficiency (BST)
+ 
+ * somehow need to copy the sources so they point at the pixels on the
+ * OUT image 
+
+ * foreach target PSF
+   * radial aperture
+   * convolve to next target PSF
+
+   * somehow need to organize the output file to have the values from
+   * the different PSFs in separate tables (with header info to
+   * specify the size of that PSF)
+
+   */

trunk/psphot/src/psphotVisual.c

@@ -93 +93 @@
     strcpy (coords.ctype, "RA---TAN");
 
-    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
-    if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) {
-        fprintf (stderr, "failed to get background values\n");
-        return false;
-    }
-
     image.Nx = inImage->numCols;
     image.Ny = inImage->numRows;
@@ -125 +118 @@
     free (image.data2d);
 
+    return true;
+}
+
+bool psphotVisualShowObjectRegions (pmReadout *readout, psMetadata *recipe, psArray *sources) {
+
+    KiiImage image;
+    KapaImageData data;
+    Coords coords;
+
+    bool status = false;
+
+    // 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 (maskVal);
+
+    maskVal |= markVal;
+
+    if (!pmVisualTestLevel("psphot.image.objects", 2)) return true;
+
+    int kapa = psphotKapaChannel (1);
+    if (kapa == -1) return false;
+
+    strcpy (coords.ctype, "RA---TAN");
+
+    psImage *inImage = readout->image;
+    psImage *inMask = readout->mask;
+    image.Nx = inImage->numCols;
+    image.Ny = inImage->numRows;
+    
+    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
+    if (!psImageBackground(stats, NULL, inImage, inMask, 0xffff, rng)) {
+        fprintf (stderr, "failed to get background values\n");
+        return false;
+    }
+
+    ALLOCATE (image.data2d, float *, image.Ny);
+    for (int iy = 0; iy < image.Ny; iy++) {
+        ALLOCATE (image.data2d[iy], float, image.Nx);
+        for (int ix = 0; ix < image.Nx; ix++) {
+            image.data2d[iy][ix] = 0;
+        }
+    }
+
+    // loop over sources and set unmasked pixels to 0
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+        if (source == NULL) continue;
+
+        psImage *mask = source->maskObj;
+        if (mask == NULL) continue;
+
+        for (int iy = 0; iy < mask->numRows; iy++) {
+            int jy = iy + mask->row0;
+            if (jy < 0) continue;
+            if (jy >= inImage->numRows) continue;
+            for (int ix = 0; ix < mask->numCols; ix++) {
+                int jx = ix + mask->col0;
+                if (jx < 0) continue;
+                if (jx >= inImage->numCols) continue;
+
+                if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue;
+                image.data2d[jy][jx] = 1;
+            }
+        }
+    }
+
+    for (int iy = 0; iy < image.Ny; iy++) {
+        for (int ix = 0; ix < image.Nx; ix++) {
+            image.data2d[iy][ix] = (image.data2d[iy][ix] == 0.0) ? NAN : inImage->data.F32[iy][ix];
+        }
+    }
+
+    strcpy (data.name, "maskObj");
+    strcpy (data.file, "maskObj");
+    // data.zero = 0.0;
+    // data.range = 1.0;
+    data.zero = stats->robustMedian - stats->robustStdev;
+    data.range = 5*stats->robustStdev;
+    data.logflux = 0;
+
+    KiiSetChannel (kapa, 2);
+    KiiNewPicture2D (kapa, &image, &data, &coords);
+
+    for (int iy = 0; iy < image.Ny; iy++) {
+        free (image.data2d[iy]);
+    }
+    free (image.data2d);
+
     psFree (stats);
     psFree (rng);
-
+    
+    pmVisualAskUser(NULL);
     return true;
 }
@@ -2401 +2488 @@
 }
 
-bool psphotVisualShowResidualImage (pmReadout *readout) {
+// option to redo variance since in some cases we may have displayed a different image in the meanwhile
+bool psphotVisualShowResidualImage (pmReadout *readout, bool reshow) {
 
     if (!pmVisualTestLevel("psphot.image.resid", 2)) return true;
@@ -2408 +2496 @@
     if (myKapa == -1) return false;
 
-    psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", 1);
+    if (reshow) {
+        psphotVisualShowMask (myKapa, readout->mask, "mask", 2);
+        psphotVisualScaleImage (myKapa, readout->variance, readout->mask, "variance", 1);
+    }
+    psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", 0);
 
     pmVisualAskUser(NULL);