Changeset 33415 for branches/meh_branches/ppstack_test/psphot

Timestamp:

Mar 5, 2012, 5:19:48 PM (14 years ago)

Author:

mhuber

Message:

merging latest r33407 trunk changes to meh_branches/ppstack_test

Location:

branches/meh_branches/ppstack_test

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/doc/sourcepixels.txt (copied) (copied from trunk/psphot/doc/sourcepixels.txt )
• psphot/src (modified) (2 props)
• psphot/src/Makefile.am (modified) (9 diffs)
• psphot/src/psphot.c (modified) (1 diff)
• psphot/src/psphot.h (modified) (11 diffs)
• psphot/src/psphotAddNoise.c (modified) (1 diff)
• psphot/src/psphotApResid.c (modified) (4 diffs)
• psphot/src/psphotBlendFit.c (modified) (7 diffs)
• psphot/src/psphotChoosePSF.c (modified) (4 diffs)
• psphot/src/psphotEfficiency.c (modified) (11 diffs)
• psphot/src/psphotEllipticalContour.c (modified) (3 diffs)
• psphot/src/psphotExtendedSourceAnalysis.c (modified) (8 diffs)
• psphot/src/psphotExtendedSourceFits.c (modified) (18 diffs)
• psphot/src/psphotFake.c (modified) (1 diff)
• psphot/src/psphotFindDetections.c (modified) (1 diff)
• psphot/src/psphotFindFootprints.c (modified) (2 diffs)
• psphot/src/psphotFindPeaks.c (modified) (1 diff)
• psphot/src/psphotFitSourcesLinear.c (modified) (8 diffs)
• psphot/src/psphotFitSourcesLinearStack.c (modified) (4 diffs)
• psphot/src/psphotFootprintSaddles.c (copied) (copied from trunk/psphot/src/psphotFootprintSaddles.c )
• psphot/src/psphotForcedReadout.c (modified) (1 diff)
• psphot/src/psphotGuessModels.c (modified) (5 diffs)
• psphot/src/psphotImageLoop.c (modified) (1 diff)
• psphot/src/psphotImageQuality.c (modified) (1 diff)
• psphot/src/psphotKronIterate.c (copied) (copied from trunk/psphot/src/psphotKronIterate.c )
• psphot/src/psphotKronMasked.c (modified) (13 diffs)
• psphot/src/psphotLoadSRCTEXT.c (modified) (1 diff)
• psphot/src/psphotMagnitudes.c (modified) (4 diffs)
• psphot/src/psphotMergeSources.c (modified) (10 diffs)
• psphot/src/psphotModelBackground.c (modified) (6 diffs)
• psphot/src/psphotModelTest.c (modified) (1 diff)
• psphot/src/psphotModelTestArguments.c (copied) (copied from trunk/psphot/src/psphotModelTestArguments.c )
• psphot/src/psphotModelTestReadout.c (copied) (copied from trunk/psphot/src/psphotModelTestReadout.c )
• psphot/src/psphotOutput.c (modified) (3 diffs)
• psphot/src/psphotPSFConvModel.c (modified) (1 diff)
• psphot/src/psphotPetrosianProfile.c (modified) (1 diff)
• psphot/src/psphotPetrosianStats.c (modified) (4 diffs)
• psphot/src/psphotRadialApertures.c (modified) (16 diffs)
• psphot/src/psphotRadialAperturesByObject.c (modified) (1 diff)
• psphot/src/psphotRadialProfile.c (modified) (1 diff)
• psphot/src/psphotRadialProfileWings.c (copied) (copied from trunk/psphot/src/psphotRadialProfileWings.c )
• psphot/src/psphotRadiusChecks.c (modified) (2 diffs)
• psphot/src/psphotReadout.c (modified) (14 diffs)
• psphot/src/psphotReadoutFindPSF.c (modified) (2 diffs)
• psphot/src/psphotReadoutForcedKnownSources.c (modified) (1 diff)
• psphot/src/psphotReadoutKnownSources.c (modified) (1 diff)
• psphot/src/psphotReadoutMinimal.c (modified) (1 diff)
• psphot/src/psphotReplaceUnfit.c (modified) (9 diffs)
• psphot/src/psphotRoughClass.c (modified) (2 diffs)
• psphot/src/psphotSersicGuess.00.00.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.00.00.h )
• psphot/src/psphotSersicGuess.00.01.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.00.01.h )
• psphot/src/psphotSersicGuess.00.02.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.00.02.h )
• psphot/src/psphotSersicGuess.00.03.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.00.03.h )
• psphot/src/psphotSersicGuess.01.00.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.01.00.h )
• psphot/src/psphotSersicGuess.01.01.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.01.01.h )
• psphot/src/psphotSersicGuess.01.02.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.01.02.h )
• psphot/src/psphotSersicGuess.01.03.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.01.03.h )
• psphot/src/psphotSersicGuess.02.00.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.02.00.h )
• psphot/src/psphotSersicGuess.02.01.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.02.01.h )
• psphot/src/psphotSersicGuess.02.02.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.02.02.h )
• psphot/src/psphotSersicGuess.02.03.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.02.03.h )
• psphot/src/psphotSersicGuess.03.00.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.03.00.h )
• psphot/src/psphotSersicGuess.03.01.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.03.01.h )
• psphot/src/psphotSersicGuess.03.02.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.03.02.h )
• psphot/src/psphotSersicGuess.03.03.h (copied) (copied from trunk/psphot/src/psphotSersicGuess.03.03.h )
• psphot/src/psphotSersicModelClass.c (copied) (copied from trunk/psphot/src/psphotSersicModelClass.c )
• psphot/src/psphotSetThreads.c (modified) (2 diffs)
• psphot/src/psphotSourceFits.c (modified) (27 diffs)
• psphot/src/psphotSourceMatch.c (modified) (11 diffs)
• psphot/src/psphotSourceSize.c (modified) (5 diffs)
• psphot/src/psphotSourceStats.c (modified) (7 diffs)
• psphot/src/psphotStackChisqImage.c (modified) (1 diff)
• psphot/src/psphotStackImageLoop.c (modified) (5 diffs, 1 prop)
• psphot/src/psphotStackMatchPSFs.c (modified) (4 diffs)
• psphot/src/psphotStackMatchPSFsNext.c (modified) (7 diffs)
• psphot/src/psphotStackMatchPSFsPrepare.c (modified) (1 diff)
• psphot/src/psphotStackMatchPSFsUtils.c (modified) (5 diffs)
• psphot/src/psphotStackObjects.c (modified) (1 diff)
• psphot/src/psphotStackOptions.c (modified) (4 diffs)
• psphot/src/psphotStackPSF.c (modified) (3 diffs)
• psphot/src/psphotStackReadout.c (modified) (10 diffs)
• psphot/src/psphotSubtractBackground.c (modified) (1 diff)
• psphot/src/psphotTest.c (modified) (1 diff)
• psphot/test/tap_dense.pro (copied) (copied from trunk/psphot/test/tap_dense.pro )

branches/meh_branches/ppstack_test/psphot/src

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

branches/meh_branches/ppstack_test/psphot/src/Makefile.am

@@ -25 +25 @@
 libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 
-bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotStack
+bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotStack psphotModelTest
 # bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy 
 
@@ -43 +43 @@
 psphotStack_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 psphotStack_LDADD = libpsphot.la
+
+psphotModelTest_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+psphotModelTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
+psphotModelTest_LDADD = libpsphot.la
 
 # psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
@@ -102 +106 @@
         psphotCleanup.c
 
+# a psphot-variant that simply generates the PSF model
+psphotModelTest_SOURCES = \
+        psphotModelTest.c            \
+        psphotModelTestArguments.c   \
+        psphotParseCamera.c        \
+        psphotImageLoop.c   \
+        psphotMosaicChip.c         \
+        psphotCleanup.c
+
 # psphotTest_SOURCES = \
 #         psphotTest.c
@@ -117 +130 @@
         psphotVisual.c                 \
         psphotCullPeaks.c              \
+        psphotFootprintSaddles.c       \
         psphotVersion.c                \
         psphotModelGroupInit.c         \
@@ -129 +143 @@
         psphotMakePSFReadout.c         \
         psphotModelBackground.c        \
+        psphotModelTestReadout.c       \
         psphotMaskBackground.c         \
         psphotSubtractBackground.c     \
@@ -157 +172 @@
         psphotExtendedSourceAnalysisByObject.c \
         psphotExtendedSourceFits.c     \
+        psphotSersicModelClass.c       \
         psphotKernelFromPSF.c          \
         psphotFitSet.c                 \
@@ -168 +184 @@
         psphotRadialPlot.c             \
         psphotKronMasked.c             \
+        psphotKronIterate.c            \
+        psphotRadialProfileWings.c     \
         psphotDeblendSatstars.c        \
         psphotMosaicSubimage.c         \
@@ -184 +202 @@
         psphotRadialBins.c             \
         psphotRadialApertures.c        \
-        psphotRadialAperturesByObject.c \
         psphotPetrosian.c              \
         psphotPetrosianRadialBins.c    \
@@ -190 +207 @@
         psphotPetrosianVisual.c        \
         psphotEfficiency.c
-
-# XXX need to fix this for the new apis
-#       psphotModelTest.c              
 
 # re-instate these

branches/meh_branches/ppstack_test/psphot/src/psphot.c

@@ -27 +27 @@
     }
 
-    psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete"));
+    psLogMsg ("psphot", PS_LOG_WARN, "complete psphot run: %f sec\n", psTimerMark ("complete"));
 
     psErrorCode exit_status = psphotGetExitStatus();

branches/meh_branches/ppstack_test/psphot/src/psphot.h

@@ -19 +19 @@
     PSPHOT_FORCED,
     PSPHOT_MAKE_PSF,
+    PSPHOT_MODEL_TEST,
 } psphotImageLoopMode;
 
@@ -107 +108 @@
 bool            psphotBlendFit_Threaded (psThreadJob *job);
 
-bool            psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule);
-bool            psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+bool            psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState);
+bool            psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool ignoreState);
 
 bool            psphotAddNoise (pmConfig *config, const pmFPAview *view, const char *filerule);
@@ -117 +118 @@
 bool            psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule);
 bool            psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+bool            psphotExtendedSourceAnalysis_Threaded (psThreadJob *job);
 
 bool            psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view, const char *filerule);
@@ -174 +176 @@
 
 // in psphotReplaceUnfit.c:
-bool            psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe);
+bool            psphotRemoveAllSourcesByArray (const psArray *sources, const psMetadata *recipe);
+bool            psphotRemoveAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState);
+bool            psphotRemoveAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState);
 bool            psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal);
 
@@ -327 +331 @@
 bool psphotMakePSFReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
 
+pmConfig *psphotModelTestArguments(int argc, char **argv);
+bool psphotModelTestReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
+
 int psphotFileruleCount(const pmConfig *config, const char *filerule);
 bool psphotFileruleCountSet(const pmConfig *config, const char *filerule, int num);
@@ -372 +379 @@
     bool convolve;                      // Convolve images?
     psphotStackConvolveSource convolveSource;
-    // psArray *convImages, *convMasks, *convVariances; // Filenames for the temporary convolved images
-
-    // bool matchZPs;                      // Adjust relative fluxes based on transparency analysis?
-    // bool photometry;                    // Perform photometry?
-    // psMetadata *stats;                  // Statistics for output
-    // FILE *statsFile;                    // File to which to write statistics
-    // psArray *origImages, *origMasks, *origVariances; // Filenames of the original images
-    // psArray *origCovars;                // Original covariances matrices
-    // int quality;                        // Bad data quality flag
 
     // Prepare
@@ -387 +385 @@
     psVector *inputMask;                // Mask for inputs
 
-    float targetSeeing;                 // Target seeing FWHM
+    psVector *targetSeeing;             // Target seeing FWHMs
     psArray *sourceLists;               // Individual lists of sources for matching
     psVector *norm;                     // Normalisation for each image
     psArray *psfs;
-
-    // psVector *exposures;                // Exposure times
-    // float sumExposure;                  // Sum of exposure times
-    // float zp;                           // Zero point for output
-    // psVector *inputMask;                // Mask for inputs
-    // psArray *sources;                   // Matched sources
 
     // Convolve
@@ -402 +394 @@
     psArray *regions;                   // PSF-matching regions --- required in the stacking
     psVector *matchChi2;                // chi^2 for stamps from matching
-    psVector *weightings;               // Combination weightings for images (1/noise^2)
-    // psArray *cells;                     // Cells for convolved images --- a handle for reading again
-    // int numCols, numRows;               // Size of image
-    // psArray *convCovars;                // Convolved covariance matrices
-
-    // Combine initial
-    // pmReadout *outRO;                   // Output readout
-    // pmReadout *expRO;                   // Exposure readout
-    // psArray *inspect;                   // Array of arrays of pixels to inspect
-
-    // Rejection
-    // psArray *rejected;                  // Rejected pixels
 } psphotStackOptions;
 
@@ -420 +400 @@
 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);
+
+bool psphotStackMatchPSFsNext (pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize);
+bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, int lastSize);
+int psphotStackMatchPSFsEntries (pmConfig *config, const pmFPAview *view, const char *filerule);
 
 // psphotStackMatchPSFsUtils
-psVector *SetOptWidths (bool *optimum, psMetadata *recipe);
-pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize);
-bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index);
-bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index);
+// psVector *SetOptWidths (bool *optimum, psMetadata *recipe);
+// pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize);
+// bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index);
+// bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index);
 bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index);
 bool matchKernel(pmConfig *config, pmReadout *cnv, pmReadout *raw, psphotStackOptions *options, int index);
-bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname);
-bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname);
-bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index);
-
-pmPSF *psphotStackPSF(const pmConfig *config, int numCols, int numRows, const psArray *psfs, const psVector *inputMask);
+// bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname);
+// bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname);
+// bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index);
+
+bool psphotStackRenormaliseVariance(const pmConfig *config, pmReadout *readout);
+
+bool psphotStackPSF(const pmConfig *config, psphotStackOptions *options);
 
 psphotStackOptions *psphotStackOptionsAlloc (int num);
@@ -443 +427 @@
 bool psphotCopySourcesReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
 
-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, psImageMaskType maskVal, const psVector *radMax, int entry);
+bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule, int entry);
+bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, int entry);
+bool psphotRadialApertures_Threaded (psThreadJob *job);
+bool psphotRadialApertureSource (pmSource *source, pmReadout *readout, int entry, psVector *pixRadius2, psVector *pixFlux, psVector *pixVer);
+// bool psphotRadialApertureSource (pmSource *source, int entry);
 
 bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);
@@ -469 +455 @@
 psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc);
 
+bool psphotSourceParents (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc);
+bool psphotSourceParentsReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
+
+bool psphotCopyPeaks (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc);
+bool psphotCopyPeaksReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
+
+bool psphotSersicModelClassGuessPCM (pmPCMdata *pcm, pmSource *source);
+void psphotSersicModelClassInit ();
+void psphotSersicModelClassCleanup ();
+
+bool psphotSetRadiusMoments (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal);
+bool psphotSetRadiusMomentsExact (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal);
+
+bool psphotFootprintSaddles(pmReadout *readout, psArray *footprints);
+bool psphotMaskFootprint (pmReadout *readout, pmSource *source, psImageMaskType markVal);
+
+bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf);
+bool psphotKronIterate_Threaded (psThreadJob *job);
+
+bool psphotRadialProfileWings (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotRadialProfileWingsReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources);
+bool psphotRadialProfileWings_Threaded (psThreadJob *job);
+
+bool psphotStackObjectsSelectForAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
+
 #endif

branches/meh_branches/ppstack_test/psphot/src/psphotAddNoise.c

@@ -85 +85 @@
     }
     if (add) {
-        psLogMsg ("psphot.noise", PS_LOG_INFO, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
+        psLogMsg ("psphot.noise", PS_LOG_WARN, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
     } else {
-        psLogMsg ("psphot.noise", PS_LOG_INFO, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
+        psLogMsg ("psphot.noise", PS_LOG_WARN, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
     }
 

branches/meh_branches/ppstack_test/psphot/src/psphotApResid.c

@@ -1 +1 @@
 # include "psphotInternal.h"
-# define DEBUG
+// # define DEBUG
 
 # define SKIPSTAR(MSG) { psTrace ("psphot", 3, "invalid : %s", MSG); continue; }
@@ -9 +9 @@
 {
     bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Aperture Residuals ---");
 
     // select the appropriate recipe information
@@ -151 +154 @@
 
         // wait for the threads to finish and manage results
-        if (!psThreadPoolWait (false)) {
+        if (!psThreadPoolWait (false, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            psFree(cellGroups);
             return false;
         }
@@ -347 +351 @@
 
     psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
-    psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid"));
+    psLogMsg ("psphot.apresid", PS_LOG_WARN, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid"));
 
     psFree (xPos);

branches/meh_branches/ppstack_test/psphot/src/psphotBlendFit.c

@@ -5 +5 @@
 {
     bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Fit Sources (Non-Linear) ---");
 
     // select the appropriate recipe information
@@ -148 +151 @@
 
         // wait for the threads to finish and manage results
-        if (!psThreadPoolWait (false)) {
+        if (!psThreadPoolWait (false, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
             psFree (fitOptions);
@@ -186 +189 @@
     psFree (fitOptions);
 
-    psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit);
+    psLogMsg ("psphot.psphotBlendFit", PS_LOG_WARN, "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, false);
@@ -233 +236 @@
         pmSource *source = sources->data[i];
 
+# if (0)
+# define TEST_X 34
+# define TEST_Y 28
+   
+        if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+            fprintf (stderr, "test object\n");
+        }
+
+# undef TEST_X
+# undef TEST_Y
+# endif
+
         // skip non-astronomical objects (very likely defects)
         if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
@@ -246 +261 @@
 
         // limit selection to some SN limit
-        if (sqrt(source->peak->detValue) < FIT_SN_LIM) continue;
-
+        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+            if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) continue;
+        } else {
+            if (sqrt(source->peak->detValue) < FIT_SN_LIM) continue;
+        }
         // exclude sources outside optional analysis region
         if (source->peak->xf < AnalysisRegion.x0) continue;
@@ -267 +285 @@
         }
 
-        // replace object in image
+        // replace object in image & remove excess noise
         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
             pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
@@ -305 +323 @@
         Nfail ++;
 
-        // re-subtract the object, leave local sky
+        // re-subtract the object, leave local sky, re-bump noise
         pmSourceCacheModel (source, maskVal);
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);

branches/meh_branches/ppstack_test/psphot/src/psphotChoosePSF.c

@@ -5 +5 @@
 {
     bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Choose PSF ---");
 
     // select the appropriate recipe information
@@ -376 +379 @@
 
     char *modelName = pmModelClassGetName (psf->type);
-    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "select psf model: %f sec\n", psTimerMark ("psphot.choose.psf"));
+    psLogMsg ("psphot.pspsf", PS_LOG_WARN, "select psf model: %f sec\n", psTimerMark ("psphot.choose.psf"));
     psLogMsg ("psphot.pspsf", PS_LOG_INFO, "psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
 
@@ -443 +446 @@
             psVectorAppend (fwhmMinor, FWHM_MINOR);
 
-            if (modelPSF->params->n >= 7) {
+            if (modelPSF->params->n > 7) {
               psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]);
             }
-            if (modelPSF->params->n >= 8) {
+            if (modelPSF->params->n > 8) {
               psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]);
             }
@@ -517 +520 @@
     }
 
-    // psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
+    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);
 

branches/meh_branches/ppstack_test/psphot/src/psphotEfficiency.c

@@ -99 +99 @@
 
 /// Generate a fake image and add it in to the existing readout
-static bool effGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources
+static pmReadout *effGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources
                         const pmReadout *ro,            // Readout of interest
                         const pmPSF *psf,               // Point-spread function
@@ -152 +152 @@
         psFree(xAll);
         psFree(yAll);
-        return false;
+        return NULL;
     }
     psFree(magAll);
@@ -161 +161 @@
 
     psBinaryOp(ro->image, ro->image, "+", fakeRO->image);
-    psFree(fakeRO);
-
-    return true;
+
+    // return the readout so we can subtract it later
+    return fakeRO;
 }
 
@@ -169 +169 @@
 {
     bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Efficiency ---");
 
     // select the appropriate recipe information
@@ -252 +255 @@
 
     // remove all sources, adding noise for subtracted sources
-    psphotRemoveAllSources(realSources, recipe);
+    psphotRemoveAllSourcesByArray(realSources, recipe);
 
 #if TESTING
@@ -290 +293 @@
 
     psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row
-    if (!effGenerate(&xFake, &yFake, readout, psf, magOffsets,
-                     numSources, magLim, radius, minFlux)) {
+    pmReadout *fakeRO = effGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, radius, minFlux);
+    if (!fakeRO) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources");
         psFree(xFake);
@@ -416 +419 @@
     psFree(significance);
 
+    // psphotFitSourcesLinearReadout subtracts the model fits
     if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
@@ -427 +431 @@
     psf->ApTrend = NULL;
 
+    // measure the magnitudes and fluxes for the sources
     if (!psphotMagnitudesReadout(config, recipe, view, readout, fakeSourcesAll, psf)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources.");
@@ -433 +438 @@
         psf->ApTrend = apTrend; // Casting away const!
         return false;
+    }
+
+    // replace the subtracted model fits
+    for (int i = 0; i < fakeSourcesAll->n; i++) {
+        pmSource *source = fakeSourcesAll->data[i];
+
+        // replace other sources?
+        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
+        
+        pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     }
     psFree(fakeSourcesAll);
@@ -515 +530 @@
     psFree(fakeSources);
 
+    // subtract the faked sources from the original image
+    psBinaryOp(readout->image, readout->image, "-", fakeRO->image);
+    psFree(fakeRO);
+
     pmDetEff *de = pmDetEffAlloc(magLim, numSources, numBins); // Detection efficiency
     de->magOffsets = psVectorCopy(NULL, magOffsets, PS_TYPE_F32);
@@ -529 +548 @@
     psFree(de);
 
-    psLogMsg("psphot", PS_LOG_INFO, "Detection efficiency: %lf sec\n", psTimerClear("psphot.fake"));
+    psLogMsg("psphot", PS_LOG_WARN, "Detection efficiency: %lf sec\n", psTimerClear("psphot.fake"));
 
     return true;

branches/meh_branches/ppstack_test/psphot/src/psphotEllipticalContour.c

@@ -127 +127 @@
 psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord) {
 
-    static int pass = 0;
-
     psF32 *par = params->data.F32;
 
@@ -145 +143 @@
 
     // value is X
-    // if (coord->data.F32[1] == 0) {
-    if (pass == 0) {
-        pass = 1;
-
+    if (coord->data.F32[1] < 0.5) {
         float value = par[PAR_RMIN]*cs_alpha*r;
 
@@ -161 +156 @@
 
     // value is Y
-    // if (coord->data.F32[1] == 1) {
-    if (pass == 1) {
-        pass = 0;
-
+    if (coord->data.F32[1] > 0.5) {
         float value = par[PAR_RMIN]*sn_alpha*r;
 

branches/meh_branches/ppstack_test/psphot/src/psphotExtendedSourceAnalysis.c

@@ -2 +2 @@
 
 // measure the elliptical radial profile and use this to measure the petrosian parameters for the sources
-// XXX this function needs to be threaded
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -8 +7 @@
 {
     bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Extended Source Analysis (Petrosians) ---");
 
     // select the appropriate recipe information
@@ -59 +61 @@
     }
 
-    // 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);
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
 
     // get the sky noise from the background analysis; if this is missing, get the user-supplied value
@@ -70 +74 @@
     }
 
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, pmSourceSortByFlux);
+
+    // option to limit analysis to a specific region
+    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
+    psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0);
+    *AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined");
+
+    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
+    int Cx = 1, Cy = 1;
+    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
+
+    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
+
+    for (int i = 0; i < cellGroups->n; i++) {
+
+        psArray *cells = cellGroups->data[i];
+
+        for (int j = 0; j < cells->n; j++) {
+
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_EXTENDED_ANALYSIS");
+
+            psArrayAdd(job->args, 1, readout);
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, AnalysisRegion);
+            psArrayAdd(job->args, 1, recipe);
+
+            PS_ARRAY_ADD_SCALAR(job->args, skynoise, PS_TYPE_F32);
+
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli
+
+// set this to 0 to run without threading
+# if (1)            
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree(AnalysisRegion);
+                return false;
+            } 
+# else
+            if (!psphotExtendedSourceAnalysis_Threaded(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree(AnalysisRegion);
+                return false;
+            }
+            psScalar *scalar = NULL;
+            scalar = job->args->data[5];
+            Next += scalar->data.S32;
+            scalar = job->args->data[6];
+            Npetro += scalar->data.S32;
+            scalar = job->args->data[7];
+            Nannuli += scalar->data.S32;
+            psFree(job);
+# endif
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false, true)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            psFree(AnalysisRegion);
+            return false;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            } else {
+                psScalar *scalar = NULL;
+                scalar = job->args->data[5];
+                Next += scalar->data.S32;
+                scalar = job->args->data[6];
+                Npetro += scalar->data.S32;
+                scalar = job->args->data[7];
+                Nannuli += scalar->data.S32;
+            }
+            psFree(job);
+        }
+    }
+    psFree (cellGroups);
+    psFree(AnalysisRegion);
+
+    psLogMsg ("psphot", PS_LOG_WARN, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
+    psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
+    psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
+
+    psphotVisualShowResidualImage (readout, false);
+
+    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
+    if (doPetrosian) {
+        psphotVisualShowPetrosians (sources);
+    }
+
+    return true;
+}
+
+bool psphotExtendedSourceAnalysis_Threaded (psThreadJob *job) {
+
+    bool status;
+
+    int Next = 0;
+    int Npetro = 0;
+    int Nannuli = 0;
+
+    // arguments: readout, sources, models, region, psfSize, maskVal, markVal
+    pmReadout *readout      = job->args->data[0];
+    psArray *sources        = job->args->data[1];
+    psRegion *region        = job->args->data[2];
+    psMetadata *recipe      = job->args->data[3];
+
+    float skynoise          = PS_SCALAR_VALUE(job->args->data[4],F32);
+
     // S/N limit to perform full non-linear fits
     float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
@@ -78 +198 @@
     bool doPetroStars   = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS");
 
-    // source analysis is done in S/N order (brightest first)
-    sources = psArraySort (sources, pmSourceSortByFlux);
-
-    // option to limit analysis to a specific region
-    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
-    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
-    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
+    // 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);
 
     // choose the sources of interest
@@ -90 +206 @@
 
         pmSource *source = sources->data[i];
+
+        // if we have checked the source validity on the basis of the object set, then 
+        // we either skip these tests below or we skip the source completely
+        if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) continue;
+        if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_KEEP) goto keepSource;
 
         // skip PSF-like and non-astronomical objects
@@ -104 +225 @@
 
         // limit selection to some SN limit
-        assert (source->peak); // how can a source not have a peak?
-        if (sqrt(source->peak->detValue) < SN_LIM) continue;
-
-        // limit selection by analysis region
-        if (source->peak->x < AnalysisRegion.x0) continue;
-        if (source->peak->y < AnalysisRegion.y0) continue;
-        if (source->peak->x > AnalysisRegion.x1) continue;
-        if (source->peak->y > AnalysisRegion.y1) continue;
+        // assert (source->peak); // how can a source not have a peak?
+        // limit selection to some SN limit
+        bool skipSource = false;
+        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+            skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr);
+        } else {
+            skipSource = (sqrt(source->peak->detValue) < SN_LIM);
+        }
+        if (skipSource) continue;
+
+        // limit selection by analysis region (this automatically apply
+        if (source->peak->x < region->x0) continue;
+        if (source->peak->y < region->y0) continue;
+        if (source->peak->x > region->x1) continue;
+        if (source->peak->y > region->y1) continue;
+
+    keepSource:
 
         // replace object in image
@@ -159 +289 @@
     }
 
-    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
-    psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
-    psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
-
-    psphotVisualShowResidualImage (readout, false);
-
-    if (doPetrosian) {
-        psphotVisualShowPetrosians (sources);
-    }
-
-    // fprintf (stderr, "xsrc : tried %ld objects\n", sources->n);
+    psScalar *scalar = NULL;
+
+    // change the value of a scalar on the array (wrap this and put it in psArray.h)
+    scalar = job->args->data[5];
+    scalar->data.S32 = Next;
+
+    scalar = job->args->data[6];
+    scalar->data.S32 = Npetro;
+
+    scalar = job->args->data[7];
+    scalar->data.S32 = Nannuli;
 
     return true;

branches/meh_branches/ppstack_test/psphot/src/psphotExtendedSourceFits.c

@@ -6 +6 @@
     bool status = true;
 
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Extended Source Fits ---");
+
     // select the appropriate recipe information
     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
@@ -12 +15 @@
     // perform full extended source non-linear fits?
     if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS")) {
-        psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");
+        psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source fits\n");
         return true;
     }
@@ -45 +48 @@
     int Nfail = 0;
 
+    psphotSersicModelClassInit();
+
     psTimerStart ("psphot.extended");
 
@@ -103 +108 @@
 
       if (item->type != PS_DATA_METADATA) {
+        // XXX we could cull the bad entries or build a validated model folder
         psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name);
-        // XXX we could cull the bad entries or build a validated model folder
       }
 
@@ -213 +218 @@
 
         // wait for the threads to finish and manage results
-        if (!psThreadPoolWait (false)) {
+        if (!psThreadPoolWait (false, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
             psFree(AnalysisRegion);
@@ -247 +252 @@
     psFree(AnalysisRegion);
 
-    psLogMsg ("psphot", PS_LOG_INFO, "extended source model fits: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
+    psphotSersicModelClassCleanup();
+
+    psphotVisualShowResidualImage (readout, false);
+
+    psLogMsg ("psphot", PS_LOG_WARN, "extended source model fits: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
@@ -265 +274 @@
     int Nplain = 0;
     int NplainPass = 0;
-    bool savePics = false;
-    float radius;
+    float fitRadius, windowRadius;
     psScalar *scalar = NULL;
-    pmMoments psfMoments;
 
     // arguments: readout, sources, models, region, psfSize, maskVal, markVal
@@ -275 +282 @@
     psMetadata *models      = job->args->data[2];
     psRegion *region        = job->args->data[3];
-    int psfSize             = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
+    int psfSize             = PS_SCALAR_VALUE(job->args->data[4],S32);
     psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
     psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
-
-    // pthread_t tid = pthread_self();     // Thread identifier
 
     // Define source fitting parameters for extended source fits
@@ -299 +304 @@
 
         // skip PSF-like and non-astronomical objects
-        if (source->type == PM_SOURCE_TYPE_STAR) continue;
+        if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
@@ -309 +314 @@
         if (source->peak->y > region->y1) continue;
 
-        // if model is NULL, we don't have a starting guess
-        // XXX use the parameters guessed from moments
-        // if (source->modelEXT == NULL) continue;
-
-        // fprintf (stderr, "fit %d,%d in thread %d\n", source->peak->x, source->peak->y, (int) tid);
-
         // replace object in image
         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
@@ -321 +320 @@
         Next ++;
 
-        // set the radius based on the footprint (also sets the mask pixels)
-        if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) {
-            fprintf (stderr, "skipping (1) %f, %f\n", source->peak->xf, source->peak->yf);
-            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-            // XXX raise an error of some kind
-            continue;
-        }
-
-        // XXX note that this changes the source moments that are published...
-        // recalculate the source moments using the larger extended-source moments radius
-        // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
-        // this uses the footprint to judge both radius and aperture?
-        // XXX save the psf-based moments for output
-        psfMoments = *source->moments;
-        if (!pmSourceMoments (source, radius, 0.0, 0.0, 0.0, maskVal)) {
-            // subtract the best fit from the object, leave local sky
-            fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf);
-            *source->moments = psfMoments;
-            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-            // XXX raise an error flag of some kind
-            continue;
-        }
-
-        // save the modelFlux here in case we need to subtract it (for failure)
-        psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux);
-        if (!modelFluxStart) {
-            pmSourceCacheModel (source, maskVal);
-            modelFluxStart = psMemIncrRefCounter (source->modelFlux);
-            if (!modelFluxStart) {
-                fprintf (stderr, "skipping (3) %f, %f\n", source->peak->xf, source->peak->yf);
-                *source->moments = psfMoments;
-                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-                // XXX raise an error of some kind?
-                continue;
-            }
-        }
-        
-        if (savePics) {
-          psphotSaveImage (NULL, readout->image, "image.xp.fits");
-        }
-
-        // array to store the pointers to the model flux images while the models are being fitted
-        psArray *modelFluxes = psArrayAllocEmpty (models->list->n);
+        // set the fit radius based on the first radial moment (also sets the mask pixels)
+        psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // NOTE : 6 allocs
+
+        // UPDATE : we have changed the moments calculation.  There is now an iteration within 
+        // psphotKronMasked to determine moments appropriate for a larger object.  The values
+        // Mrf, KronFlux, and KronFluxErr are calculated for the iterated radius.  The other
+        // values are left at the psf-based values.
 
         // allocate the array to store the model fits
@@ -380 +343 @@
 
           // check the SNlim and skip model if source is too faint
-          float SNlim = psMetadataLookupF32 (&status, model, "SNLIM_VALUE");
+          float FIT_SN_LIM = psMetadataLookupF32 (&status, model, "SNLIM_VALUE");
           assert (status);
 
           // limit selection to some SN limit
           // assert (source->peak); // how can a source not have a peak?
-          if (sqrt(source->peak->detValue) < SNlim) {
+          // limit selection to some SN limit
+          bool skipSource = false;
+          if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+              skipSource = (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr);
+          } else {
+              skipSource = (sqrt(source->peak->detValue) < FIT_SN_LIM);
+          }
+          if (skipSource) {
               Nfaint ++;
               continue;
@@ -398 +368 @@
           assert (status);
 
-          // XXX psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6);
-          // XXX psTraceSetLevel ("psphot.psphotModelWithPSF_LMM", 6);
-
           // fit the model as convolved or not
           pmModel *modelFit = NULL;
           if (convolved) {
-              modelFit = psphotFitPCM (readout, source, fitOptions, modelType, maskVal, markVal, psfSize);
+              // NOTE : 4 more allocs to here
+              modelFit = psphotFitPCM (readout, source, fitOptions, modelType, maskVal, markVal, psfSize); // NOTE : 2313 allocs in here
               if (!modelFit) {
                   psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My);
@@ -410 +378 @@
                   continue;
               }
-              psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq);
+              psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", 
+                       source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
               Nconvolve ++;
               if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
@@ -425 +394 @@
                   continue;
               }
-              pmSourceCacheModel (source, maskVal);
-              psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq);
+              psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", 
+                       source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
               Nplain ++;
               if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
@@ -434 +403 @@
           }
 
-          // save each of the model flux images and store the best
-          psArrayAdd (modelFluxes, 4, source->modelFlux);
+          // XXX deprecate?
+          // XXX pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
+          // XXX extFitPars->Mxx = source->moments->Mxx;
+          // XXX extFitPars->Mxy = source->moments->Mxy;
+          // XXX extFitPars->Myy = source->moments->Myy;
+          // XXX extFitPars->Mrf = source->moments->Mrf;
+          // XXX extFitPars->Mrh = source->moments->Mrh;
+          // XXX extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
+          
+          // save kron mag, but assign apMag & psfMag on output (not yet calculated)
+          // XXX extFitPars->krMag = -2.5*log10(source->moments->KronFlux);
+
+          // psArrayAdd (source->extFitPars, 4, extFitPars);
+          // psFree (extFitPars);
 
           // test for fit quality / result
-          modelFit->fitRadius = radius;
+          modelFit->fitRadius = fitRadius;
           psArrayAdd (source->modelFits, 4, modelFit);
 
@@ -463 +444 @@
         }
 
+        // clear the circular mask
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
+
         if (minModel == -1) {
-          // re-subtract the object, leave local sky
+          // no valid extended fit; re-subtract the object, leave local sky
           psTrace ("psphot", 5, "failed to fit extended source model to object at %f, %f", source->moments->Mx, source->moments->My);
 
-          // replace original model, subtract it
-          psFree (source->modelFlux);
-          source->modelFlux = modelFluxStart;
-
-          *source->moments = psfMoments;
+          // ensure the modelEXT is cached
+          pmSourceCacheModel (source, maskVal);
+
           pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
-          psFree (modelFluxes);
-
-          if (savePics) {
-              psphotSaveImage (NULL, readout->image, "image.xp.fits");
-              char key[10];
-              fprintf (stdout, "continue? ");
-              if (!fgets (key, 8, stdin)) {
-                  psWarning("Couldn't read anything.");
-              } else if (key[0] == 'n') {
-                  savePics = false;
-              }
-          }
           continue;
         }
@@ -493 +463 @@
         source->modelEXT = psMemIncrRefCounter (source->modelFits->data[minModel]);
 
-        // save the modelFlux for the best model
-        psFree (source->modelFlux);
-        source->modelFlux = psMemIncrRefCounter (modelFluxes->data[minModel]);
-
-        // replace the original moments
-        *source->moments = psfMoments;
+        // adjust the window so the subtraction covers the faint wings
+        psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal);
+
+        // cache the model flux
+        if (source->modelEXT->isPCM) {
+            pmPCMCacheModel (source, maskVal, psfSize);
+        } else {
+            pmSourceCacheModel (source, maskVal);
+        }
 
         // subtract the best fit from the object, leave local sky
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
-        // the initial model flux is no longer needed
-        psFree (modelFluxStart);
-        psFree (modelFluxes);
-
         psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq);
         psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-
-        if (savePics) {
-          psphotSaveImage (NULL, readout->image, "image.xm.fits");
-          char key[10];
-          fprintf (stdout, "continue? ");
-          if (!fgets (key, 8, stdin)) {
-              psWarning("Couldn't read anything.");
-          } else if (key[0] == 'n') {
-              savePics = false;
-          }
-        }
     }
     psFree (fitOptions);
-
-    // fprintf (stderr, "xfit : tried %ld objects\n", sources->n);
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)

branches/meh_branches/ppstack_test/psphot/src/psphotFake.c

@@ -166 +166 @@
 
     // remove all sources, adding noise for subtracted sources
-    psphotRemoveAllSources(realSources, recipe);
+    psphotRemoveAllSourcesByArray(realSources, recipe);
     psphotAddNoise(readout, realSources, recipe);
 

branches/meh_branches/ppstack_test/psphot/src/psphotFindDetections.c

@@ -7 +7 @@
 {
     bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Find Detections ---");
 
     // select the appropriate recipe information

branches/meh_branches/ppstack_test/psphot/src/psphotFindFootprints.c

@@ -35 +35 @@
         psArray *tmp = pmFootprintArrayGrow(footprints, growRadius);
         psImageConvolveSetThreads(oldThreads);
-        psLogMsg ("psphot", PS_LOG_MINUTIA, "grow footprint coverage by %d pixels, %ld footprints become %ld footprints: %f sec\n", growRadius, footprints->n, tmp->n, psTimerMark ("psphot.footprints"));
-        psFree(footprints);
-        footprints = tmp;
+        psLogMsg ("psphot", PS_LOG_MINUTIA, "grow footprint coverage by %d pixels, %ld footprints become %ld footprints: %f sec\n", growRadius, footprints->n, tmp ? tmp->n : 0, psTimerMark ("psphot.footprints"));
+        if (tmp) {
+            psFree(footprints);
+            footprints = tmp;
+        } else {
+            psLogMsg ("psphot", PS_LOG_WARN, "pmFootprintArray grow returned NULL\n");
+        }
+            
     }
 
@@ -58 +63 @@
     psphotCullPeaks(readout, significance, recipe, detections->footprints);
     detections->peaks = pmFootprintArrayToPeaks(detections->footprints);
-    psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot.footprints"));
+
+    // psphotFootprintSaddles (readout, detections->footprints);
+    // 
+    // int nSaddle = 0;
+    // for (int i = 0; i < detections->peaks->n; i++) {
+    //  pmPeak *peak = detections->peaks->data[i];
+    //  
+    //  if (peak->saddlePoints) nSaddle += peak->saddlePoints->n;
+    // }
+    // fprintf (stderr, "%d saddle points for %ld peaks\n", nSaddle, detections->peaks->n);
+
+    psLogMsg ("psphot", PS_LOG_WARN, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot.footprints"));
 
     return detections;

branches/meh_branches/ppstack_test/psphot/src/psphotFindPeaks.c

@@ -68 +68 @@
         pmPeaksWriteText (peaks, output);
     }
-    psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot.peaks"));
+    psLogMsg ("psphot", PS_LOG_WARN, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot.peaks"));
 
     return peaks;

branches/meh_branches/ppstack_test/psphot/src/psphotFitSourcesLinear.c

@@ -17 +17 @@
     bool status = true;
 
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Fit Source (Linear) ---");
+
     // select the appropriate recipe information
     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
@@ -23 +26 @@
     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
 
         // find the currently selected readout
@@ -126 +134 @@
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
+
+        // psAssert (source->peak, "source without peak??");
+        // psAssert (source->peak->footprint, "peak without footprint??");
 
         // turn this bit off and turn it on again if we pass this test
@@ -190 +201 @@
         pmModel *model = pmSourceGetModel (&isPSF, source);
 
+        // clear the 'mark' pixels and remask on the fit aperture 
         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, model->fitRadius, "OR", markVal);
@@ -203 +215 @@
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), sources->n);
-
-    // fprintf (stderr, "****** Nsat : %d ********\n", Nsat);
 
     if (fitSources->n == 0) {
@@ -330 +340 @@
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
 
-    // XXXX **** philosophical question:
-    // we measure bright objects in three passes: 1) linear fit; 2) non-linear fit; 3) linear fit:
-    // should retain the chisq and errors from the intermediate non-linear fit?
-    // the non-linear fit provides better values for the position errors, and for
-    // extended sources, the shape errors
+    // Philosophical question: we measure bright objects in three passes: 1) linear fit; 2)
+    // non-linear fit; 3) linear fit: should we retain the chisq and errors from the
+    // intermediate non-linear fit?  the non-linear fit provides better values for the position
+    // errors, and for extended sources, the shape errors
 
     // adjust I0 for fitSources and subtract
@@ -349 +358 @@
         model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
 
-        // subtract object
+        // clear the 'mark' pixels so the subtraction covers the full window
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+
+        // subtract object & add noise
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     }
@@ -374 +386 @@
     psFree (border);
 
-    psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
+    psLogMsg ("psphot.ensemble", PS_LOG_WARN, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
 
     psphotVisualPlotChisq (sources);

branches/meh_branches/ppstack_test/psphot/src/psphotFitSourcesLinearStack.c

@@ -64 +64 @@
           // turn this bit off and turn it on again if we keep this source
           source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;
+
+          // skip non-astronomical objects (very likely defects)
+          if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+          if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+
+          // do not include CRs in the full ensemble fit
+          if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
+
+          // do not include MOMENTS_FAILURES in the fit
+          if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+
+          if (final) {
+              if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
+          } else {
+              // if (source->mode & PM_SOURCE_MODE_BLEND) continue;
+          }
 
           // generate model for sources without, or skip if we can't
@@ -83 +99 @@
           }
 
-          source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
+          bool isPSF = false;
+          pmModel *model = pmSourceGetModel (&isPSF, source);
+
+          // clear the 'mark' pixels and remask on the fit aperture 
+          psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+          psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, model->fitRadius, "OR", markVal);
+
+          // we call this function multiple times. for the first time, we have only PSF models for all objects
+          // the second time has extended sources.  If we ever fit the PSF model, we should raise this bit
+          source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
+          if (isPSF) {
+              source->mode |= PM_SOURCE_MODE_PSFMODEL;
+          }         
           psArrayAdd (fitSources, 100, source);
         }
@@ -166 +194 @@
         model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
         model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
+
+        // clear the 'mark' pixels so the subtraction covers the full window
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
 
         // subtract object
@@ -190 +221 @@
     psFree (errors);
 
-    psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
+    psLogMsg ("psphot.ensemble", PS_LOG_WARN, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
     return true;
 }

branches/meh_branches/ppstack_test/psphot/src/psphotForcedReadout.c

@@ -55 +55 @@
     psphotLoadExtSources (config, view, filerule);
 
-    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
-    psphotGuessModels (config, view, filerule);
-
     // merge the newly selected sources into the existing list
     // NOTE: merge OLD and NEW
     psphotMergeSources (config, view, filerule);
+
+    // Construct an initial model for each object, set the radius to fitRadius, set circular
+    // fit mask.  NOTE: only applied to sources without guess models
+    psphotGuessModels (config, view, filerule);
 
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)

branches/meh_branches/ppstack_test/psphot/src/psphotGuessModels.c

@@ -46 +46 @@
     psAssert (detections, "missing detections?");
 
-    psArray *sources = detections->newSources;
+    psArray *sources = detections->allSources;
     psAssert (sources, "missing sources?");
 
@@ -112 +112 @@
         // wait here for the threaded jobs to finish
         // fprintf (stderr, "wait for threads (%d, %d)\n", jx, jy);
-        if (!psThreadPoolWait (false)) {
+        if (!psThreadPoolWait (false, true)) {
+            // harvest our jobs
+            psFree(cellGroups);
             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
             return false;
@@ -141 +143 @@
     psFree (cellGroups);
 
-    psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.models"));
+    psLogMsg ("psphot.models", PS_LOG_WARN, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.models"));
     return true;
 }
@@ -159 +161 @@
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
+
+        // do not redo sources already guessed
+        if (source->tmpFlags & PM_SOURCE_TMPF_MODEL_GUESS) continue;
 
         // this is used to mark sources for which the model is measured. We check later that
@@ -169 +174 @@
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
         if (!source->peak) continue;
+
+        // psAssert (source->peak->footprint, "peak without footprint??");
 
         nSrc ++;

branches/meh_branches/ppstack_test/psphot/src/psphotImageLoop.c

@@ -129 +129 @@
                     }
                     break;
+                  case PSPHOT_MODEL_TEST:
+                    if (!psphotModelTestReadout (config, view, "PSPHOT.INPUT")) {
+                        psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
+                        psFree (view);
+                        return false;
+                    }
+                    break;
                 }
             }

branches/meh_branches/ppstack_test/psphot/src/psphotImageQuality.c

@@ -289 +289 @@
 #endif
 
-    psLogMsg ("psphot", PS_LOG_INFO, "Image Quality Stats from %ld psf stars : FWHM (major, minor) [pixels]: %f, %f\n",
+    psLogMsg ("psphot", PS_LOG_WARN, "Image Quality Stats from %ld psf stars : FWHM (major, minor) [pixels]: %f, %f\n",
               M2->n, fwhm_major, fwhm_minor);
 

branches/meh_branches/ppstack_test/psphot/src/psphotKronMasked.c

@@ -1 +1 @@
 # include "psphotInternal.h"
-
+# ifndef ROUND
+# define ROUND(X) ((int) ((X) + 0.5*SIGN(X)))
+# endif
+
+bool psphotKronMask (pmSource *source, psImage *kronMask, psImageMaskType markVal);
 bool psphotKronMag (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal);
 
@@ -45 +49 @@
 }
 
+int psphotKapaChannel (int channel);
+bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel);
+bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
+
 bool psphotKronMaskedReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) {
 
@@ -72 +80 @@
     }
 
+    float EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS");
+
     // bit-masks to test for good/bad pixels
     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
@@ -83 +93 @@
     maskVal |= markVal;
 
+    // psphotSaveImage (NULL, readout->mask, "kron.unmasked.fits");
+
+    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY;
+
+    // XXX tmp visualization
+    // int kapa = psphotKapaChannel (1);
+
     // generate the mask image: increment counter for every source overlapping the pixel
-    psImage *kronMask = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_U8);
+    psImage *kronMask = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
     psImageInit (kronMask, 0);
-    int Nx = kronMask->numCols;
-    int Ny = kronMask->numRows;
     for (int i = 0; i < sources->n; i++) {
 
@@ -96 +111 @@
         if (!source->moments) continue;
 
-        // replace object in image
-        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
-            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-        }
-
-        psphotKronMag (source, RADIUS, MIN_KRON_RADIUS, maskVal);
-
-        // re-subtract the object, leave local sky
-        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-
-        continue;
-
-        // XXX skip this code
-        // generate the pixel masks
-        // int Xo = source->moments->Mx;
-        // int Yo = source->moments->My;
-        float dX = source->moments->Mx - source->peak->xf;
-        float dY = source->moments->My - source->peak->yf;
-        float dR = hypot(dX, dY);
-        
-        float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
-        float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
-        
-        int Kr = 2.5*source->moments->Mrf;
-        int Kr2 = Kr*Kr;
-        
-        for (int iy = Yo - Kr; iy < Yo + Kr + 1; iy++) {
-            if (iy < 0) continue;
-            if (iy >= Ny) continue;
-            for (int ix = Xo - Kr; ix < Xo + Kr + 1; ix++) {
-                if (ix < 0) continue;
-                if (ix >= Nx) continue;
-
-                if (PS_SQR(ix - Xo) + PS_SQR(iy - Yo) > Kr2) continue;
-                if (kronMask->data.U8[iy][ix] < 0xff) {
-                    kronMask->data.U8[iy][ix] ++;
-                }
-            }
-        }
-    }           
+        // replace object in image
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+        }
+
+        // save the PSF-based values
+        // source->moments->KronFluxPSF    = source->moments->KronFlux;
+        // source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
+        // source->moments->KronRadiusPSF  = source->moments->Mrf;
+
+        // iterate to the window radius
+        float windowRadius = RADIUS;
+        for (int j = 0; j < 4; j++) {
+            // XXX use some S/N criterion to limit?
+            // if (source->moments->KronFlux / source->moments->KronFluxErr > 10.0) {
+            // windowRadius = PS_MAX (RADIUS, 10*source->moments->Mrf);
+            // }
+
+            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
+            pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
+
+            // mask the pixels not contained by the footprint
+            // if (psphotMaskFootprint (readout, source, markVal)) {
+            //  // psphotVisualShowMask (kapa, source->maskObj, "mask", 2);
+            //  // psphotVisualRangeImage (kapa, source->pixels, "image", 1, -50, 300);
+            //  // fprintf (stderr, "masked\n");
+            // }            
+
+            // mask the pixels associated with near neighbors
+            if (psphotKronMask (source, kronMask, markVal)) {
+                // psphotVisualShowMask (kapa, source->maskObj, "mask", 2);
+                // psphotVisualRangeImage (kapa, source->pixels, "image", 1, -50, 300);
+                // fprintf (stderr, "masked\n");
+            }       
+            // this function populates moments->Mrf,KronFlux,KronFluxErr
+            // XXX what about KronInner, KronOuter, etc?
+            psphotKronMag (source, windowRadius, MIN_KRON_RADIUS, maskVal);
+            windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
+            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+        }
+        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
+        float kmag = -2.5*log10(source->moments->KronFlux);
+        if (source->psfMag - kmag > 0.25) {
+            // psphotVisualShowMask (kapa, source->maskObj, "mask", 1);
+            // psphotVisualRangeImage (kapa, source->pixels, "image", 0, -50, 300);
+            // fprintf (stderr, "psf: %f, kron: %f, dmag: %f\n", source->psfMag, kmag, source->psfMag - kmag);
+            // fprintf (stderr, "continue\n");
+        }
+        // re-subtract the object, leave local sky
+        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+    }
+
+    // psphotSaveImage (NULL, readout->mask, "kron.masked.fits");
     // psphotSaveImage (NULL, kronMask, "kronmask.fits");
+
     psLogMsg ("psphot.kron", PS_LOG_DETAIL, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
 
@@ -142 +171 @@
     return true;
 }
+
+# define WEIGHTED 0
 
 bool psphotKronMag (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal) {
@@ -150 +181 @@
     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
 
-    psF32 Sum = 0.0;
-    psF32 Var = 0.0;
     psF32 R2 = PS_SQR(radius);
+
+# if (WEIGHTED)
+    float rsigma2 = 0.5 / PS_SQR(radius / 4.0); // use a Gaussian window with sigma = R_window / 2
+# endif
 
     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
@@ -168 +201 @@
     psF32 RS = 0.0;
 
-# if (1)
+    // the peak position is less accurate but less subject to extreme deviations
     float dX = source->moments->Mx - source->peak->xf;
     float dY = source->moments->My - source->peak->yf;
@@ -174 +207 @@
     float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
     float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
-# else
-    float Xo = source->moments->Mx;
-    float Yo = source->moments->My;
-# endif
 
     // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
@@ -210 +239 @@
             if (r2 > R2) continue;
 
+# if (WEIGHTED)
+            float z = r2 * rsigma2;
+            assert (z >= 0.0);
+            float weight  = exp(-z);
+            psF32 pDiff = *vPix * weight;
+# else
             psF32 pDiff = *vPix;
-
-            Sum += pDiff;
-
-            // Kron Flux uses the 1st radial moment (NOT Gaussian windowed)
+# endif
+
+            // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
             psF32 rf = pDiff * sqrt(r2);
             psF32 rs = pDiff;
@@ -234 +268 @@
 
     int nKronPix = 0;
-    Sum = Var = 0.0;
+    float Sum = 0.0;
+    float Var = 0.0;
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
@@ -263 +298 @@
             if (r2 > radKron2) continue;
 
+# if (WEIGHTED)
+            float z = r2 * rsigma2;
+            assert (z >= 0.0);
+            float weight  = exp(-z);
+            psF32 pDiff = *vPix * weight;
+            psF32 wDiff = *vWgt * weight;
+# else
             psF32 pDiff = *vPix;
             psF32 wDiff = *vWgt;
+# endif
 
             Sum += pDiff;
@@ -272 +315 @@
     }
 
-    source->moments->Mrh = Mrf;
-
-    // XXX for a test, save the old values here:
-    source->moments->KronFouter    = source->moments->KronFlux;
-    source->moments->KronCoreErr = source->moments->KronFluxErr;
-
-    source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
-    source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
+    source->moments->Mrf = Mrf;
+    source->moments->KronFlux    = Sum;
+    source->moments->KronFluxErr = sqrt(Var);
+    // source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
+    // source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
 
     return true;
 }
+
+enum {Y_P, Y_M, X_P, X_M};
+
+bool psphotKronMask (pmSource *source, psImage *kronMask, psImageMaskType markVal) {
+
+    int X, Y, e, dXs, dYs;
+
+    if (!source->peak) return false;
+    if (!source->peak->saddlePoints) return false;
+    if (!source->peak->saddlePoints->n) return false;
+
+    // return true;
+
+    int xMax = source->maskObj->numCols;
+    int yMax = source->maskObj->numRows;
+    // int xOff = source->maskObj->col0;
+    // int yOff = source->maskObj->row0;
+
+    int Xo = ROUND(source->peak->xf) - source->maskObj->col0;
+    int Yo = ROUND(source->peak->yf) - source->maskObj->row0;
+
+    for (int i = 0; i < source->peak->saddlePoints->n; i++) {
+
+        psVector *saddlePoint = source->peak->saddlePoints->data[i];
+        int Xs = ROUND(saddlePoint->data.S32[0]) - source->maskObj->col0;
+        int Ys = ROUND(saddlePoint->data.S32[1]) - source->maskObj->row0;
+        
+        // fprintf (stderr, "mask %d,%d @ %d,%d (%d,%d)\n", Xo, Yo, Xs, Ys, saddlePoint->data.S32[0], saddlePoint->data.S32[1]);
+
+        // We want to mask the pixels between the edge of the image and the line perpendicular
+        // to the connecting line.  Call dX,dY = (Xs-Xo,Ys-Yo).  There are 4 cases: 
+        // +y : dY > 0, fabs(dY) > fabs(dX)
+        // -y : dY < 0, fabs(dY) > fabs(dX)
+        // +x : dX > 0, fabs(dX) > fabs(dY)
+        // -x : dX < 0, fabs(dX) > fabs(dY)
+
+        int dX = Xs - Xo;
+        int dY = Ys - Yo;
+
+        // +y : dY > 0, fabs(dY) > fabs(dX)
+        // -y : dY < 0, fabs(dY) > fabs(dX)
+        if (fabs(dY) > fabs(dX)) {
+            // points to right of saddle
+            Y = Ys;
+            e = 0;
+            dXs = (dY > 0) ? +dY : -dY;  // this forces dXs > 0
+            dYs = (dY > 0) ? -dX : +dX;
+            for (X = Xs; X < xMax; X++) {
+                if (dY > 0) {
+                    for (int Ym = Y; Ym < yMax; Ym++) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+                    }
+                } else {
+                    for (int Ym = Y; Ym >= 0; Ym--) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+                    }
+                }
+                e += dYs;
+                int e2 = 2 * e;
+                if (e2 > dXs) {
+                    Y ++;
+                    e -= dXs;
+                }
+                if (e2 < -dXs) {
+                    Y --;
+                    e += dXs;
+                }
+                if (Y >= yMax) break;
+                if (Y < 0) break;
+            }
+            // points to left of saddle
+            Y = Ys;
+            e = 0;
+            for (X = Xs; X >= 0; X--) {
+                if (dY > 0) {
+                    for (int Ym = Y; Ym < yMax; Ym++) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+                    }
+                } else {
+                    for (int Ym = Y; Ym >= 0; Ym--) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+                    }
+                }
+                e -= dYs;
+                int e2 = 2 * e;
+                if (e2 > dXs) {
+                    Y ++;
+                    e -= dXs;
+                }
+                if (e2 < -dXs) {
+                    Y --;
+                    e += dXs;
+                }
+                if (Y >= yMax) break;
+                if (Y < 0) break;
+            }
+        } else {
+            // points to right of saddle
+            X = Xs;
+            e = 0;
+            dXs = (dX > 0) ? -dY : +dY;
+            dYs = (dX > 0) ? +dX : -dX;  // this forces dYs > 0
+            for (Y = Ys; Y < yMax; Y++) {
+                if (dX > 0) {
+                    for (int Xm = X; Xm < xMax; Xm++) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+                    }
+                } else {
+                    for (int Xm = X; Xm >= 0; Xm--) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+                    }
+                }
+                e += dXs;
+                int e2 = 2 * e;
+                if (e2 > dYs) {
+                    X ++;
+                    e -= dYs;
+                }
+                if (e2 < -dYs) {
+                    X --;
+                    e += dYs;
+                }
+                if (X >= yMax) break;
+                if (X < 0) break;
+            }
+            // points to left of saddle
+            X = Xs;
+            e = 0;
+            for (Y = Ys; Y >= 0; Y--) {
+                if (dX > 0) {
+                    for (int Xm = X; Xm < xMax; Xm++) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+                    }
+                } else {
+                    for (int Xm = X; Xm >= 0; Xm--) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+                    }
+                }
+                e -= dXs;
+                int e2 = 2 * e;
+                if (e2 > dYs) {
+                    X ++;
+                    e -= dYs;
+                }
+                if (e2 < -dYs) {
+                    X --;
+                    e += dYs;
+                }
+                if (X >= yMax) break;
+                if (X < 0) break;
+            }
+        }
+    }
+    return true;
+}

branches/meh_branches/ppstack_test/psphot/src/psphotLoadSRCTEXT.c

@@ -73 +73 @@
             dPAR[PM_PAR_I0]   = 0.0;
 
-            pmPSF_AxesToModel (PAR, axes);
+            pmPSF_AxesToModel (PAR, axes, modelType);
 
             float peakFlux    = 1.0;

branches/meh_branches/ppstack_test/psphot/src/psphotMagnitudes.c

@@ -4 +4 @@
 {
     bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Magnitudes ---");
 
     // select the appropriate recipe information
@@ -126 +129 @@
 
         // wait for the threads to finish and manage results
-        if (!psThreadPoolWait (false)) {
+        if (!psThreadPoolWait (false, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
             return false;
@@ -146 +149 @@
     psFree (cellGroups);
 
-    psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "measure magnitudes : %f sec for %ld objects (%d with apertures)\n", psTimerMark ("psphot.mags"), sources->n, Nap);
+    psLogMsg ("psphot.magnitudes", PS_LOG_WARN, "measure magnitudes : %f sec for %ld objects (%d with apertures)\n", psTimerMark ("psphot.mags"), sources->n, Nap);
     return true;
 }
@@ -306 +309 @@
 
         // wait for the threads to finish and manage results
-        if (!psThreadPoolWait (false)) {
+        if (!psThreadPoolWait (false, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
             return false;

branches/meh_branches/ppstack_test/psphot/src/psphotMergeSources.c

@@ -520 +520 @@
 }
 
+// copy the newPeaks from the detections of one pmFPAfile to another
+bool psphotCopyPeaks (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 (!psphotCopyPeaksReadout (config, view, ruleOut, ruleSrc, i)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to copy peaks from %s to %s entry %d", ruleSrc, ruleOut, i);
+            return false;
+        }
+    }
+    return true;
+}
+
+// add newly detected peaks to the existing list of sources
+bool psphotCopyPeaksReadout (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 *detections = psMetadataLookupPtr (&status, readoutSrc->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    // 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");
+    psAssert (detectionsOut, "missing PSPHOT.DETECTIONS?");
+
+    psAssert (detectionsOut->peaks, "programming error");
+    psAssert (!detectionsOut->oldPeaks, "programming error");
+    psAssert (detections->peaks, "programming error");
+
+    // save the OUT existing peaks on oldPeaks
+    detectionsOut->oldPeaks = detectionsOut->peaks;
+    detectionsOut->peaks = psArrayAllocEmpty(detections->peaks->n);
+
+    for (int i = 0; i < detections->peaks->n; i++) {
+        psAssert (detections->peaks->data[i], "programming error");
+        pmPeak *peak = pmPeakAlloc (0, 0, 0.0, PM_PEAK_LONE);
+        pmPeakCopy(peak, detections->peaks->data[i]);
+        psArrayAdd (detectionsOut->peaks, 100, peak);
+        psFree (peak);
+    }
+    return true;
+}
+
+// create source parents children from ruleSrc for ruleOut for orphans
+bool psphotSourceParents (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 (!psphotSourceParentsReadout (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 parents from ruleSrc for ruleOut for orphaned children for this readout.  
+bool psphotSourceParentsReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) {
+
+    bool status;
+    int nParents = 0;
+    int nNonOrphans = 0;
+
+    // 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");
+    psAssert (detectionsOut, "missing PSPHOT.DETECTIONS?");
+
+    // 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];
+      if (sourceSrc->parent) {
+          nNonOrphans ++;
+          continue; // Not an orphan
+      }
+
+      pmSource *sourceOut = pmSourceCopy(sourceSrc);
+      sourceOut->parent = sourceSrc;
+      
+      // keep the original source flags
+      sourceOut->seq      = sourceSrc->seq;
+      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
+      if (sourceOut->modelPSF) {
+        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;
+
+      nParents ++;
+      psArrayAdd (detectionsOut->allSources, 100, sourceOut);
+      psFree (sourceOut);
+    }
+    psLogMsg ("psphot", 3, "%d parents created, %d unorphaned children, %ld input vs %ld output", nParents, nNonOrphans, sourcesSrc->n, detectionsOut->allSources->n);
+
+    return true;
+}
+
 // create source children from ruleSrc for ruleOut
 bool psphotSourceChildren (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc)
@@ -542 +714 @@
 }
 
-// create source children from ruleSrc for ruleOut for this entry.  XXX currently, this is only
+// Create source children from ruleSrc for ruleOut for this entry.  Currently, this is only
 // used by psphotStackReadout (sources go on allSources so that psphotChoosePSF can be called
-// repeatedly)
+// repeatedly).
 bool psphotSourceChildrenReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) {
 
@@ -569 +741 @@
     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
+    // replace any existing DETECTION container on readoutOut->analysis with the new one
+    pmDetections *detectionsOut = pmDetectionsAlloc();
+    detectionsOut->allSources = psArrayAllocEmpty (100);
+    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;
     }
 
@@ -593 +760 @@
       
       // keep the original source flags
+      sourceOut->seq      = sourceSrc->seq;
       sourceOut->type     = sourceSrc->type;
       sourceOut->mode     = sourceSrc->mode;
@@ -612 +780 @@
       pmSourceFreePixels (sourceOut);
 
+      // XXX do we need to skip the Chisq image sources?
+
       // allocate image, weight, mask for the new image for each peak
-      if (sourceOut->modelPSF) {
-        pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius);
-      }
+      pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->windowRadius);
 
       // child sources have not been subtracted in this image, but this flag may be raised if
@@ -624 +792 @@
       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)
+    psLogMsg ("psphot", 3, "created %ld children", detectionsOut->allSources->n);
+
+    psFree(detectionsOut); // a copy remains on the analysis metadata
+
+    return true;
+}
+
+// create source children associated with 'filerule' from the objectsSrc.  returns a new object
+// array containing the child sources.  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) {
 
@@ -652 +821 @@
         psAssert (readout, "missing readout?");
 
+        // create DETECTIONS containers for each image, in case one lacks it
         pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
         if (!detections) {
@@ -665 +835 @@
             psAssert (detections, "missing detections?");
         }
+
+        // we need to save the new sources on the detection arrays of the appropriate image
         detArrays->data[i] = psMemIncrRefCounter(detections);
         readouts->data[i] = psMemIncrRefCounter(readout);
@@ -695 +867 @@
 
             pmSource *sourceOut = pmSourceCopy(sourceSrc);
+            sourceOut->parent = sourceSrc;
+
+            // save on the output object array at the same location
             objectOut->sources->data[i] = sourceOut;
-
-            sourceOut->parent = sourceSrc;
 
             // keep the original source flags and sequence ID (if set)
@@ -720 +893 @@
             pmSourceFreePixels (sourceOut);
 
-            // set the output readotu
+            // set the output readout
             int index = sourceOut->imageID;
             if (index >= readouts->n) continue; // skip the sources generated by the chisq image

branches/meh_branches/ppstack_test/psphot/src/psphotModelBackground.c

@@ -145 +145 @@
     psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
 
+    psVector *dQ = psVectorAllocEmpty (100, PS_TYPE_F32);
     psF32 **modelData = model->data.F32;
     psF32 **modelStdevData = modelStdev->data.F32;
@@ -216 +217 @@
                 }
                 modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth);
+                // fprintf (stderr, "dQ : %f - %f - %f = %f\n", stats->robustLQ, stats->robustMedian, stats->robustUQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
+                psVectorAppend (dQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
 
                 // supply sample to plotting routing
@@ -322 +325 @@
     }
 
-    psLogMsg ("psphot", PS_LOG_INFO, "built background image: %f sec\n", psTimerMark ("psphot.background"));
-
-    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value);
+    psLogMsg ("psphot", PS_LOG_WARN, "built background image: %f sec\n", psTimerMark ("psphot.background"));
+
+    psStats *statsDQ = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+    psVectorStats (statsDQ, dQ, NULL, NULL, 0);
+
+    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN",  PS_META_REPLACE, "sky mean", Value);
     psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev);
-    psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f stdev %f", Value, ValueStdev);
+    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_DQ",    PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian);
+    psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian);
 
     // measure image and background stats and save for later output
@@ -334 +341 @@
                                       PS_STAT_MAX);
     psImageStats (statsBck, model, NULL, 0);
-    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN", PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN",  PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DEV", PS_META_REPLACE, "sky stdev",               ValueStdev);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DQ",  PS_META_REPLACE, "sky quartile slope",      statsDQ->sampleMedian);
     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max);
     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min);
@@ -343 +352 @@
               statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev);
 
+    psFree(statsDQ);
+    psFree(dQ);
+
     psFree(stats);
     psFree(statsBck);
@@ -401 +413 @@
     int num = psphotFileruleCount(config, filerule);
 
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Model Background ---");
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {

branches/meh_branches/ppstack_test/psphot/src/psphotModelTest.c

@@ -1 @@
-# include "psphotInternal.h"
-# define PM_SOURCE_FIT_PSF_X_EXT PM_SOURCE_FIT_PSF_AND_SKY
+# include "psphotStandAlone.h"
 
-// XXX add more test information?
-bool psphotModelTest (pmConfig *config, const pmFPAview *view, const char *filerule, psMetadata *recipe) {
+int main (int argc, char **argv) {
 
-    bool status;
-    int modelType = -1;
-    float obsMag, fitMag, value;
-    char name[64];
-    pmPSF *psf = NULL;
-    pmSourceFitMode fitMode;
+    psMemInit();
+    psTimerStart ("complete");
+    pmErrorRegister();                  // register psModule's error codes/messages
+    psphotInit();
 
-    // bit-masks to test for good/bad pixels
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
-    assert (maskVal);
+    // load command-line arguments, options, and system config data
+    pmConfig *config = psphotModelTestArguments (argc, argv);
+    assert(config);
 
-    // bit-mask to mark pixels not used in analysis
-    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
-    assert (markVal);
+    psphotVersionPrint();
 
-    // maskVal is used to test for rejected pixels, and must include markVal
-    maskVal |= markVal;
-
-    // run model fitting tests on a single source?
-    if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false;
-
-    psTimerStart ("modelTest");
-
-    // find the currently selected readout
-    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
-    PS_ASSERT_PTR_NON_NULL (readout, false);
-
-    // use poissonian errors or local-sky errors
-    bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, filerule);
-    if (!status) POISSON_ERRORS = true;
-    pmSourceFitModelInit (15, 0.1, 1.0, POISSON_ERRORS);
-
-    // find the various fitting parameters (try test values first)
-    float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS");
-    if (!status || !isfinite(INNER)) {
-        INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
-    }
-    float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS");
-    if (!status || !isfinite(OUTER)) {
-        OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
-    }
-    float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS");
-    if (!status || !isfinite(RADIUS)) {
-        RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
-    }
-    float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS");
-    if (!status || !isfinite(mRADIUS)) {
-        mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
+    // load input data (config and images (signal, noise, mask)
+    if (!psphotParseCamera (config)) {
+        psErrorStackPrint(stderr, "Error setting up the camera\n");
+        exit (psphotGetExitStatus());
     }
 
-    // define the source of interest
-    float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
-    float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
-    if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined");
-
-    // what fitting mode to use?
-    fitMode = PM_SOURCE_FIT_EXT;
-    char *fitModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE");
-    if (fitModeWord && !strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF;
-    if (fitModeWord && !strcasecmp (fitModeWord, "CONV")) fitMode = PM_SOURCE_FIT_PSF_X_EXT;
-    if (fitModeWord && !strcasecmp (fitModeWord, "DEFAULT")) fitMode = PM_SOURCE_FIT_EXT;
-
-    // construct the source structures
-    pmSource *source = pmSourceAlloc();
-    source->peak = pmPeakAlloc (xObj, yObj, 0, 0);
-    pmSourceDefinePixels (source, readout, xObj, yObj, OUTER);
-
-    // in fitMode, psf sets the model type
-    if (fitMode == PM_SOURCE_FIT_PSF) {
-        psf = psphotLoadPSF (config, view, recipe);
-        if (!psf) psAbort("PSF_INPUT_FILE not supplied");
-        modelType = psf->type;
-        source->type = PM_SOURCE_TYPE_STAR;
-    }
-    if (fitMode == PM_SOURCE_FIT_EXT) {
-        // find the model: supplied by user or first in the PSF_MODEL list
-        char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
-        if (!status || !strcasecmp (modelName, "DEFAULT")) {
-            // get the list pointers for the PSF_MODEL entries
-
-            psList *list = NULL;
-            psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
-            if (mdi == NULL) psAbort("missing PSF_MODEL selection");
-            if (mdi->type == PS_DATA_STRING) {
-                list = psListAlloc(NULL);
-                psListAdd (list, PS_LIST_HEAD, mdi);
-            } else {
-                if (mdi->type != PS_DATA_METADATA_MULTI) psAbort("missing PSF_MODEL selection");
-                list = psMemIncrRefCounter(mdi->data.list);
-            }
-
-            // take the first list element
-            psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
-            modelName = item->data.V;
-        }
-        modelType = pmModelClassGetType (modelName);
-        if (modelType < 0) psAbort("unknown model %s", modelName);
-        source->type = PM_SOURCE_TYPE_EXTENDED;
-    }
-    if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
-        // we need to load BOTH a psf and an ext model
-        psf = psphotLoadPSF (config, view, recipe);
-        if (!psf) psAbort("PSF_INPUT_FILE not supplied");
-
-        // find the model: supplied by user or first in the PSF_MODEL list
-        char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
-        if (!status || !strcasecmp (modelName, "DEFAULT")) {
-            // get the list pointers for the PSF_MODEL entries
-
-            psList *list = NULL;
-            psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
-            if (mdi == NULL) psAbort("missing PSF_MODEL selection");
-            if (mdi->type == PS_DATA_STRING) {
-                list = psListAlloc(NULL);
-                psListAdd (list, PS_LIST_HEAD, mdi);
-            } else {
-                if (mdi->type != PS_DATA_METADATA_MULTI) psAbort("missing PSF_MODEL selection");
-                list = psMemIncrRefCounter(mdi->data.list);
-            }
-
-            // take the first list element
-            psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
-            modelName = item->data.V;
-        }
-        modelType = pmModelClassGetType (modelName);
-        if (modelType < 0) psAbort("unknown model %s", modelName);
-        source->type = PM_SOURCE_TYPE_EXTENDED;
+    // call psphot for each readout
+    if (!psphotImageLoop (config, PSPHOT_MODEL_TEST)) {
+        psErrorStackPrint(stderr, "Error in the psphot image loop\n");
+        exit (psphotGetExitStatus());
     }
 
-    // find the local sky
-    status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
-    if (!status) psAbort("pmSourceLocalSky error");
+    psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete"));
 
-    // get the source moments
-    status = pmSourceMoments (source, mRADIUS, 0.0, 1.0, maskVal);
-    if (!status) psAbort("psSourceMoments error");
-    source->peak->value = source->moments->Peak;
+    psErrorCode exit_status = psphotGetExitStatus();
+    psphotCleanup (config);
+    exit (exit_status);
+}
 
-    fprintf (stderr, "sum: %f @ (%f, %f)\n", source->moments->Sum, source->moments->Mx, source->moments->My);
-    fprintf (stderr, "moments: %f, %f - %f\n", source->moments->Mxx, source->moments->Myy, source->moments->Mxy);
-
-    psEllipseMoments moments;
-    moments.x2 = source->moments->Mxx;
-    moments.y2 = source->moments->Myy;
-    moments.xy = source->moments->Mxy;
-    psEllipseAxes axes = psEllipseMomentsToAxes (moments, 20.0);
-
-    fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
-
-    // get the initial model parameter guess
-    pmModel *model = pmSourceModelGuess (source, modelType);
-    source->modelEXT = model;
-
-    // if any parameters are defined by the user, take those values
-    int nParams = pmModelClassParameterCount (modelType);
-    psF32 *params = model->params->data.F32;
-    params[PM_PAR_XPOS] = xObj; // XXX use the user-supplied value,
-    params[PM_PAR_YPOS] = yObj; // XXX or use the centroid
-    for (int i = 0; i < nParams; i++) {
-        if (i == PM_PAR_XPOS) continue;
-        if (i == PM_PAR_YPOS) continue;
-
-        sprintf (name, "TEST_FIT_PAR%d", i);
-        value = psMetadataLookupF32 (&status, recipe, name);
-        if (status && isfinite (value)) {
-            params[i] = value;
-        }
-    }
-
-    float area = params[4]*params[5];
-    fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
-
-    // for PSF fitting, set the shape parameters based on the PSF & source position
-    if (fitMode == PM_SOURCE_FIT_PSF) {
-        source->modelPSF = pmModelFromPSF (model, psf);
-        psFree (model);
-        model = source->modelPSF;
-        params = model->params->data.F32;
-    }
-
-    // list model input shape
-    psEllipseShape shape;
-    shape.sx  = 1.4 / model->params->data.F32[4];
-    shape.sy  = 1.4 / model->params->data.F32[5];
-    shape.sxy = model->params->data.F32[6];
-    axes = psEllipseShapeToAxes (shape, 20.0);
-
-    fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
-
-    fprintf (stderr, "input parameters: \n");
-    for (int i = 0; i < nParams; i++) {
-        fprintf (stderr, "%d : %f\n", i, params[i]);
-    }
-
-    // define the pixels used for the fit
-    psImageKeepCircle (source->maskObj, xObj, yObj, RADIUS, "OR", markVal);
-    psphotSaveImage (NULL, source->maskObj, "mask1.fits");
-
-    char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET");
-    if (status) {
-        status = psphotFitSet (source, model, fitset, fitMode, maskVal);
-        exit (0);
-    }
-
-    if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
-        // build the psf for the object
-        source->modelPSF = pmModelFromPSF (model, psf);
-        source->modelEXT = model;
-
-        // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
-        int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
-        assert (status);
-
-        model = psphotPSFConvModel (readout, source, modelType, maskVal, markVal, psfSize);
-        params = model->params->data.F32;
-    } else {
-        status = pmSourceFitModel (source, model, fitMode, maskVal);
-    }
-
-    // measure the source mags
-    pmSourcePhotometryModel (&fitMag, model);
-    pmSourcePhotometryAper  (NULL, &obsMag, NULL, NULL, model, source->pixels, source->variance, source->maskObj, maskVal);
-    fprintf (stderr, "ap: %f, fit: %f, apmifit: %f, nIter: %d\n", obsMag, fitMag, obsMag - fitMag, model->nIter);
-
-    // write out positive object
-    psphotSaveImage (NULL, source->pixels, "object.fits");
-
-    // subtract object, leave local sky
-    // pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
-    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-
-    fprintf (stderr, "output parameters: \n");
-    for (int i = 0; i < nParams; i++) {
-        fprintf (stderr, "%d : %f\n", i, params[i]);
-    }
-
-    // write out
-    psphotSaveImage (NULL, source->pixels, "resid.fits");
-    psphotSaveImage (NULL, source->maskObj, "mask.fits");
-
-    psLogMsg ("psphot", PS_LOG_INFO, "model test : %f sec\n", psTimerMark ("modelTest"));
-
-    exit (0);
-}
+// all functions which return to this level must raise one of the top-level error codes if they
+// exit with an error.  these error codes are used to specify the program exit status

branches/meh_branches/ppstack_test/psphot/src/psphotOutput.c

@@ -325 +325 @@
     psMetadataItemSupplement (&status, header, analysis, "MSKY_NY");
 
+    psMetadataItemSupplement (&status, header, analysis, "MSKY_DEV");
+    psMetadataItemSupplement (&status, header, analysis, "MSKY_DQ");
+
     psMetadataItemSupplement (&status, header, analysis, "DETEFF.MAGREF");
 
@@ -404 +407 @@
     return true;
 }
+
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule)
+{
+    static int npass = 0;
+    char filename[64];
+
+    // XXX dump tests are disabled unless this is commented out:
+    return true;
+
+    bool status = true;
+
+    int num = psphotFileruleCount(config, filerule);
+
+    snprintf (filename, 64, "testdump.%02d.dat", npass);
+    FILE *f = fopen (filename, "w");
+
+    // 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->newSources ? detections->newSources : detections->allSources;
+        psAssert (sources, "missing sources?");
+
+        if (detections->newSources) {
+            fprintf (f, "## --- from new sources ---\n");
+        } else {
+            fprintf (f, "## --- from all sources ---\n");
+        }
+
+        for (int i = 0; i < sources->n; i++) {
+            pmSource *source = sources->data[i];
+            if (!source) continue;
+
+            pmPeak *peak = source->peak;
+            if (!peak) continue;
+
+            // XXX only dump a given region
+            // if (peak->xf < 20) continue;
+            // if (peak->yf < 20) continue;
+            // if (peak->xf > 40) continue;
+            // if (peak->yf > 70) continue;
+
+            float Msum = source->moments ? source->moments->Sum : NAN;
+            float Mx   = source->moments ? source->moments->Mx : NAN;
+            float My   = source->moments ? source->moments->My : NAN;
+            // float Npix = source->moments ? source->moments->nPixels : NAN;
+            float Io = source->modelPSF ? source->modelPSF->params->data.F32[PM_PAR_I0] : NAN;
+            fprintf (f, "%d %f %f  : %f %f : %f %f\n", source->imageID, peak->xf, peak->yf, Mx, My, Msum, Io);
+        }
+    }
+    fclose (f);
+    npass ++;
+
+    return true;
+}
+
+# if (0)
 bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule, char *filename) {
 
@@ -446 +516 @@
     return true;
 }
-
+# endif

branches/meh_branches/ppstack_test/psphot/src/psphotPSFConvModel.c

@@ -43 +43 @@
     // adjust the pixels based on the footprint
     float radius = psphotSetRadiusEXT (readout, source, markVal);
-    if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false;
+    if (!pmSourceMoments (source, radius, 0.25*radius, 0.0, maskVal)) return false;
 
     // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:

branches/meh_branches/ppstack_test/psphot/src/psphotPetrosianProfile.c

@@ -38 +38 @@
     // convert the isophotal radius vs angle measurements to an elliptical contour
     if (!psphotEllipticalContour (source, petrosian)) {
-        psLogMsg ("psphot", 3, "failed to measure elliptical contour");
+        // psLogMsg ("psphot", 3, "failed to measure elliptical contour");
         psFree (petrosian);
         return false;

branches/meh_branches/ppstack_test/psphot/src/psphotPetrosianStats.c

@@ -143 +143 @@
     // if we failed to reach the PETROSIAN_RATIO, use the lowest significant ratio instead (flag this!)
     if (!anyPetro) {
-        // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
-        if (lowestSignificantRadius == 0) { 
-            // assume Fmax @ R = 0.0
-            petRadius    = InterpolateValues     (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
-            petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[lowestSignificantRadius]);
-
-        } else {
-            int n0 = lowestSignificantRadius-1;
-            int n1 = lowestSignificantRadius;
-            petRadius    = InterpolateValues     (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO);
-            petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO, petRatioErr->data.F32[n0], petRatioErr->data.F32[n1]);
+        petRadius = refRadius->data.F32[lowestSignificantRadius];
+        petRadiusErr = NAN;
+        if (!isfinite(petRadius)) {
+            fprintf (stderr, "nan pet radius\n");
         }
     }
@@ -170 +163 @@
                 petArea    = InterpolateValues     (refRadius->data.F32[i-1], areaSum->data.F32[i-1], refRadius->data.F32[i], areaSum->data.F32[i], apRadius);
                 petApix    = InterpolateValues     (refRadius->data.F32[i-1], apixSum->data.F32[i-1], refRadius->data.F32[i], apixSum->data.F32[i], apRadius);
+                if (!isfinite(petFlux)) {
+                    fprintf (stderr, "nan pet flux\n");
+                }
                 break;
             }
@@ -188 +184 @@
         if (!found50 && (fluxSum->data.F32[i] > flux50)) {
             if (i == 0) {
-                psWarning ("does this case make any sense? (fluxSum[0] %f > flux50 %f)", fluxSum->data.F32[i], flux50);
-                continue;
+                R50    = InterpolateValues     (fluxSum->data.F32[i], refRadius->data.F32[i], fluxSum->data.F32[i+1], refRadius->data.F32[i+1], flux50);
+                R50err = InterpolateValuesErrX (fluxSum->data.F32[i], refRadius->data.F32[i], fluxSum->data.F32[i+1], refRadius->data.F32[i+1], flux50, sqrt(fluxSumErr2->data.F32[i]), sqrt(fluxSumErr2->data.F32[i+1]));
+                found50 = true;
             } else {
                 R50    = InterpolateValues     (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50);
@@ -198 +195 @@
         if (!found90 && (fluxSum->data.F32[i] > flux90)) {
             if (i == 0) {
-                psWarning ("does this case make any sense? (fluxSum[0] %f > flux90 %f)", fluxSum->data.F32[i], flux90);
-                continue;
+                R90    = InterpolateValues     (fluxSum->data.F32[i], refRadius->data.F32[i], fluxSum->data.F32[i+1], refRadius->data.F32[i+1], flux90);
+                R90err = InterpolateValuesErrX (fluxSum->data.F32[i], refRadius->data.F32[i], fluxSum->data.F32[i+1], refRadius->data.F32[i+1], flux90, sqrt(fluxSumErr2->data.F32[i]), sqrt(fluxSumErr2->data.F32[i+1]));
+                found90 = true;
             } else {
                 R90    = InterpolateValues     (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90);

branches/meh_branches/ppstack_test/psphot/src/psphotRadialApertures.c

@@ -3 +3 @@
 bool psphotRadialAperturesSortFlux (psVector *radius, psVector *pixFlux, psVector *pixVar);
 
+// this function measures the radial aperture fluxes for the set of readouts.  this function
+// may be called multiple times, presumably for different versions of PSF-matched or unmatched images.  
+
 // for now, let's store the detections on the readout->analysis for each readout
-bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule)
+bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule, int entry)
 {
     bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Radial Apertures ---");
 
     // select the appropriate recipe information
@@ -20 +26 @@
     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 (!psphotRadialAperturesReadout (config, view, filerule, i, recipe)) {
+        if (i == chisqNum) continue; // skip chisq image
+
+        if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe, entry)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i);
             return false;
@@ -30 +42 @@
 }
 
+// these values are used by all threads repeatedly (and are not modified)
+static psVector *aperRadii = NULL;
+static psVector *aperRadii2 = NULL;
+static float outerRadius = NAN;
+static float SN_LIM = NAN;
+static psImageMaskType maskVal = 0;
+
 // aperture-like measurements for extended sources
 // flux in simple, circular apertures
-bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+// 'entry' tells us which of the matched-PSF images we are working on (0 == unmatched image, also non-stack psphot)
+bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, int entry) {
 
     bool status;
@@ -63 +83 @@
     }
 
-    // radMax stores the upper bounds of the annuli
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+
+    // aperRadii stores the upper bounds of the annuli
     // XXX keep the same name here as for the petrosian / elliptical apertures?
-    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
-    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];
+    aperRadii = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+    psAssert (aperRadii, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
+    psAssert (aperRadii->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
+
+    outerRadius = aperRadii->data.F32[aperRadii->n - 1];
+
+    // save the R^2 values as well for quicker comparison
+    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]);
+    }
 
     // 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
+    maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     assert (maskVal);
 
     // S/N limit to perform full non-linear fits
-    float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
+    SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
 
     // source analysis is done in S/N order (brightest first)
@@ -81 +114 @@
     sources = psArraySort (sources, pmSourceSortByFlux);
 
+    // XXX make this consistent with entry 0 == unmatched
+    int nEntry = 1;
+    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
+    if (fwhmValues) {
+        psAssert (entry < fwhmValues->n, "inconsistent matched-PSF entry");
+        nEntry = fwhmValues->n;
+    }
+    if (entry > 0) {
+        psLogMsg ("psphot", PS_LOG_DETAIL, "Radial Apertures for matched image %s : PSF FWHM = %f pixels\n", file->name, fwhmValues->data.F32[entry]);
+    } else {
+        psLogMsg ("psphot", PS_LOG_DETAIL, "Radial Apertures for unmatched image %s\n", file->name);
+    }
+
     // option to limit analysis to a specific region
     char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
-    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
-    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
+    psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0);
+    *AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined");
+
+    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
+    int Cx = 1, Cy = 1;
+    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
+
+    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
+
+    for (int i = 0; i < cellGroups->n; i++) {
+
+        psArray *cells = cellGroups->data[i];
+
+        for (int j = 0; j < cells->n; j++) {
+
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_RADIAL_APERTURES");
+
+            psArrayAdd(job->args, 1, readout);
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, AnalysisRegion);
+            PS_ARRAY_ADD_SCALAR(job->args, entry,  PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, nEntry, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nradial
+
+            // set this to 0 to run without threading
+# if (1)            
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree(AnalysisRegion);
+                return false;
+            } 
+# else
+            if (!psphotRadialApertures_Threaded(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree(AnalysisRegion);
+                return false;
+            }
+            psScalar *scalar = NULL;
+            scalar = job->args->data[5];
+            Nradial += scalar->data.S32;
+            psFree(job);
+# endif
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false, true)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            psFree(AnalysisRegion);
+            return false;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            } else {
+                psScalar *scalar = NULL;
+                scalar = job->args->data[5];
+                Nradial += scalar->data.S32;
+            }
+            psFree(job);
+        }
+    }
+    psFree (cellGroups);
+    psFree(AnalysisRegion);
+    psFree (aperRadii2);
+
+    psLogMsg ("psphot", PS_LOG_WARN, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
+    return true;
+}
+ 
+bool psphotRadialApertures_Threaded (psThreadJob *job) {
+
+    int Nradial = 0;
+
+    // arguments: readout, sources, models, region, psfSize, maskVal, markVal
+    pmReadout *readout      = job->args->data[0];
+    psArray *sources        = job->args->data[1];
+    psRegion *region        = job->args->data[2];
+    int entry               = PS_SCALAR_VALUE(job->args->data[3],S32); // which psf-matched image are we working on? (0 == unmatched)
+    int nEntry              = PS_SCALAR_VALUE(job->args->data[4],S32); // total number of psf-matched images + 1 unmatched
+
+    // storage for the derived pixel values (these are passed into psphotRadialApertureSource)
+    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
 
     // choose the sources of interest
@@ -90 +223 @@
 
         pmSource *source = sources->data[i];
+
+        // if we have checked the source validity on the basis of the object set, then 
+        // we either skip these tests below or we skip the source completely
+        if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_SKIP) continue;
+        if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_KEEP) goto keepSource;
 
         // skip PSF-like and non-astronomical objects
@@ -104 +242 @@
 
         // limit selection by analysis region
-        if (source->peak->x < AnalysisRegion.x0) continue;
-        if (source->peak->y < AnalysisRegion.y0) continue;
-        if (source->peak->x > AnalysisRegion.x1) continue;
-        if (source->peak->y > AnalysisRegion.y1) continue;
+        if (source->peak->x < region->x0) continue;
+        if (source->peak->y < region->y0) continue;
+        if (source->peak->x > region->x1) continue;
+        if (source->peak->y > region->y1) continue;
+
+    keepSource:
 
         // allocate pmSourceExtendedParameters, if not already defined
-        if (!source->radialAper) {
-            source->radialAper = psArrayAlloc(1);
+        // XXX check that nPSFsizes is consistent with targets
+        if (source->parent) {
+            if (!source->parent->radialAper) {
+                source->parent->radialAper = psArrayAlloc(nEntry);
+            }
+        } else {
+            if (!source->radialAper) {
+                source->radialAper = psArrayAlloc(nEntry);
+            }
         }
 
@@ -119 +266 @@
         }
 
-        // 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 ++;
 
-        // force source image to be a bit larger...
-        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
-
-        if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, 0)) {
+        if (!psphotRadialApertureSource (source, readout, entry, pixRadius2, pixFlux, pixVar)) {
             psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
         } else {
@@ -137 +274 @@
         }
 
-        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 source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
+    psScalar *scalar = job->args->data[5];
+    scalar->data.S32 = Nradial;
+
+    psFree (pixRadius2);
+    psFree (pixFlux);
+    psFree (pixVar);
+
     return true;
 }
 
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
-
+bool psphotRadialApertureSource (pmSource *source, pmReadout *readout, int entry, psVector *pixRadius2, psVector *pixFlux, psVector *pixVar) {
+                                            
     // if we are a child source, save the results to the parent source radial aperture array
     psArray *radialAperSet = source->radialAper;
@@ -163 +300 @@
     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);
+    // find the largest aperture of interest (use only apertures with inner radii <=
+    // source->skyRadius)
+    int lastAp = aperRadii->n;
+    for (int i = 0; i < aperRadii->n; i++) {
+        if (aperRadii->data.F32[i] < source->skyRadius) continue;
+        lastAp = i + 1;
+        break;
+    }
 
     // outer-most radius for initial truncation
-    float Rmax  = aperRadii->data.F32[aperRadii->n - 1];
+    float Rmax  = aperRadii->data.F32[lastAp - 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]);
-    }
+    // in this function, the operatins are relative to the full image (readout->image, etc)
 
     float xCM = NAN, yCM = NAN;
     if (pmSourcePositionUseMoments(source)) {
-        xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
-        yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
+        xCM = source->moments->Mx; // index coord of peak in readout
+        yCM = source->moments->My; // index coord of peak in readout
     } else {
-        xCM = source->peak->xf - 0.5 - source->pixels->col0; // coord of peak in subimage
-        yCM = source->peak->yf - 0.5 - source->pixels->row0; // coord of peak in subimage
-    }
+        xCM = source->peak->xf; // index coord of peak in readout
+        yCM = source->peak->yf; // index coord of peak in readout
+    }
+
+    int Nx = readout->image->numCols;
+    int Ny = readout->image->numRows;
+
+    pixRadius2->n = 0;
+    pixFlux->n = 0;
+    pixVar->n = 0;
 
     // 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;
+    for (int iy = -Rmax; iy < Rmax + 1; iy++) {
+
+        float yDiff = iy + 0.5 + yCM;  // y-coordinate at this offse
+        int yPix = (int) yDiff;
+
+        if (yPix < 0) continue;
+        if (yPix > Ny - 1) continue;
+        if (fabs(iy) > Rmax) continue;
+
+        float *vPix = readout->image->data.F32[yPix];
+        float *vWgt = readout->variance->data.F32[yPix];
+        psImageMaskType  *vMsk = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[yPix];
+
+        for (int ix = -Rmax; ix < Rmax + 1; ix++) {
+
+            float xDiff = ix + 0.5 + xCM;  // x-coordinate at this offse
+            int xPix = (int) xDiff;
+            
+            if (xPix < 0) continue;
+            if (xPix > Nx - 1) continue;
+            if (fabs(ix) > Rmax) continue;
+            
+            if (vMsk[xPix] & maskVal) continue;
+            if (isnan(vPix[xPix])) continue;
 
             // radius is just a function of (xDiff, yDiff)
-            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            float r2  = PS_SQR(ix) + PS_SQR(iy);
             if (r2 > Rmax2) continue;
 
             psVectorAppend(pixRadius2, r2);
-            psVectorAppend(pixFlux, *vPix);
-            psVectorAppend(pixVar, *vWgt);
+            psVectorAppend(pixFlux, vPix[xPix]);
+            psVectorAppend(pixVar, vWgt[xPix]);
         }
     }
@@ -226 +372 @@
     psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
 
+    // init the apertures of interest to 0.0, the rest go to NAN
     psVectorInit (flux,    0.0);
     psVectorInit (fluxStd, 0.0);
     psVectorInit (fluxErr, 0.0);
     psVectorInit (fill,    0.0);
+    for (int i = lastAp; i < flux->n; i++) {
+        flux->data.F32[i] = NAN;
+        fluxStd->data.F32[i] = NAN;
+        fluxErr->data.F32[i] = NAN;
+        fill->data.F32[i] = NAN;
+    }
 
     float *rPix2 = pixRadius2->data.F32;
@@ -236 +389 @@
         int j = 0;
         float *aRad2 = aperRadii2->data.F32;
-        for (; (*aRad2 < *rPix2) && (j < aperRadii2->n); j++, aRad2++);
-        for (; j < aperRadii2->n; j++, aRad2++) {
+        for (; (*aRad2 < *rPix2) && (j < lastAp); j++, aRad2++);
+
+        // XXX I can speed this up by only saving this single aperture
+        for (; j < lastAp; j++, aRad2++) {
             flux->data.F32[j]    += pixFlux->data.F32[i];
             fluxStd->data.F32[j] += PS_SQR(pixFlux->data.F32[i]);
@@ -249 +404 @@
        2) the fractional fill factor (count of valid pixels / effective area of the aperture
        3) the error on the flux within that aperture
-     */
-
-    for (int i = 0; i < flux->n; i++) {
+    */
+
+    for (int i = 0; i < lastAp; i++) {
         // calculate the total flux for bin 'nOut'
         float Area = M_PI*aperRadii2->data.F32[i];
@@ -261 +416 @@
         // XXX report the total flux or the mask-corrected flux?
         // flux->data.F32[i]    = SBmean * Area;
-        // fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
+        // fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / Pinx;
 
         fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
@@ -271 +426 @@
     }
     
+# if (1)
     radialAper->flux = flux;
     radialAper->fluxStdev = fluxStd;
     radialAper->fluxErr = fluxErr;
     radialAper->fill = fill;
-
-    psFree (aperRadii2);
-    psFree (pixRadius2);
-    psFree (pixFlux);
-    psFree (pixVar);
+# else
+    // XXX TEST
+    psFree(flux);
+    psFree(fluxStd);
+    psFree(fluxErr);
+    psFree(fill);
+# endif
 
     return true;
 }
+
+/*** below is a test to use a sort to speed this up, not very successfully ***/
 
 static int nCalls = 0;
@@ -373 +533 @@
 // *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
 # define COMPARE_VECT(A,B) (radius->data.F32[A] < radius->data.F32[B])
-# define SWAP_VECT(TYPE,A,B) { \
-  float tmp; \
-  if (A != B) { \
-    tmp = radius->data.F32[A]; \
-    radius->data.F32[A] = radius->data.F32[B];  \
-    radius->data.F32[B] = tmp; \
-    tmp = pixFlux->data.F32[A]; \
-    pixFlux->data.F32[A] = pixFlux->data.F32[B]; \
-    pixFlux->data.F32[B] = tmp; \
-    tmp = pixVar->data.F32[A]; \
-    pixVar->data.F32[A] = pixVar->data.F32[B]; \
-    pixVar->data.F32[B] = tmp; \
-  } \
-}
+# define SWAP_VECT(TYPE,A,B) {                                  \
+        float tmp;                                              \
+        if (A != B) {                                           \
+            tmp = radius->data.F32[A];                          \
+            radius->data.F32[A] = radius->data.F32[B];          \
+            radius->data.F32[B] = tmp;                          \
+            tmp = pixFlux->data.F32[A];                         \
+            pixFlux->data.F32[A] = pixFlux->data.F32[B];        \
+            pixFlux->data.F32[B] = tmp;                         \
+            tmp = pixVar->data.F32[A];                          \
+            pixVar->data.F32[A] = pixVar->data.F32[B];          \
+            pixVar->data.F32[B] = tmp;                          \
+        }                                                       \
+    }
 
 bool psphotRadialAperturesSortFlux (psVector *radius, psVector *pixFlux, psVector *pixVar) {

branches/meh_branches/ppstack_test/psphot/src/psphotRadialAperturesByObject.c

@@ -162 +162 @@
             // 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);
         }
     }

branches/meh_branches/ppstack_test/psphot/src/psphotRadialProfile.c

@@ -39 +39 @@
     // convert the isophotal radius vs angle measurements to an elliptical contour
     if (!psphotEllipticalContour (source)) {
-        psLogMsg ("psphot", 3, "failed to measure elliptical contour");
+        // psLogMsg ("psphot", 3, "failed to measure elliptical contour");
         return false;
     }

branches/meh_branches/ppstack_test/psphot/src/psphotRadiusChecks.c

@@ -142 +142 @@
 }
 
+# define MIN_WINDOW 5.0
+# define SCALE1 5.0
+# define SCALE2 12.0
+
+// call this function whenever you (re)-define the EXT model
+// XXX this function does not shrink the window
+bool psphotSetRadiusMoments (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal) {
+
+    psAssert (source, "source not defined??");
+    psAssert (source->moments, "moments not defined??");
+
+    *fitRadius = SCALE1 * source->moments->Mrf;
+    *fitRadius = PS_MIN (PS_MAX(*fitRadius, MIN_WINDOW), EXT_FIT_MAX_RADIUS);
+
+    *windowRadius = SCALE2 * source->moments->Mrf;
+    *windowRadius = PS_MIN (PS_MAX(*windowRadius, 2.5*MIN_WINDOW), 2.5*EXT_FIT_MAX_RADIUS);
+
+    // redefine the pixels if needed
+    pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, *windowRadius);
+
+    // set the mask to flag the excluded pixels
+    psImageKeepCircle (source->maskObj, source->peak->xf, source->peak->yf, *fitRadius, "OR", markVal);
+
+    return true;
+}
+# undef SCALE1
+# undef SCALE2
+# undef MIN_WINDOW
+
+# define MIN_WINDOW 5.0
+# define SCALE1 5.0
+# define PAD_WINDOW 3.0
+
+// call this function whenever you (re)-define the EXT model
+// XXX alternate function to set exactly the desired window size
+bool psphotSetRadiusMomentsExact (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal) {
+
+    psRegion newRegion;
+
+    psAssert (source, "source not defined??");
+    psAssert (source->moments, "moments not defined??");
+
+    *fitRadius = SCALE1 * source->moments->Mrf;
+    *fitRadius = PS_MIN (PS_MAX(*fitRadius, MIN_WINDOW), EXT_FIT_MAX_RADIUS);
+
+    *windowRadius = *fitRadius + PAD_WINDOW;
+
+    // check to see if new region is completely contained within old region
+    newRegion = psRegionForSquare (source->peak->xf, source->peak->yf, *windowRadius);
+    newRegion = psRegionForImage (readout->image, newRegion);
+
+    // redefine the pixels to match
+    pmSourceRedefinePixelsByRegion (source, readout, newRegion);
+
+    // set the mask to flag the excluded pixels
+    psImageKeepCircle (source->maskObj, source->peak->xf, source->peak->yf, *fitRadius, "OR", markVal);
+
+    return true;
+}
+
 // call this function whenever you (re)-define the EXT model
 bool psphotSetRadiusFootprint (float *radius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float factor) {
@@ -181 +241 @@
 }
 
+// call this function whenever you (re)-define the EXT model
+bool psphotMaskFootprint (pmReadout *readout, pmSource *source, psImageMaskType markVal) {
+
+    psAssert (source, "source not defined??");
+    psAssert (source->peak, "peak not defined??");
+
+    pmPeak *peak = source->peak;
+
+    // set the radius based on the footprint:
+    if (!peak->footprint) return false;
+    pmFootprint *footprint = peak->footprint;
+    if (!footprint->spans) return false;
+    if (footprint->spans->n < 1) return false;
+
+    int Xo = source->maskObj->col0;
+    int Yo = source->maskObj->row0;
+
+    // mark all pixels 
+    for (int j = 0; j < source->maskObj->numRows; j++) {
+        for (int i = 0; i < source->maskObj->numCols; i++) {
+            source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= markVal;
+        }
+    }
+
+    psImageMaskType clearVal = PS_NOT_IMAGE_MASK(markVal);
+
+    for (int j = 0; j < footprint->spans->n; j++) {
+        pmSpan *span = footprint->spans->data[j];
+
+        // mask the rows before and after each span
+        int minX = span->x0 - Xo - 2;
+        int maxX = span->x1 - Xo + 2;
+        int myY = span->y - Yo;
+
+        // unmark pixels inside the footprint
+        for (int i = minX; i <= maxX; i++) {
+            source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[myY][i] &= clearVal;
+        }
+    }
+    return true;
+}
+
 // alternative EXT radius based on model guess (for use without footprints)
 bool psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep) {

branches/meh_branches/ppstack_test/psphot/src/psphotReadout.c

@@ -9 +9 @@
 }
 
-// for now, let's store the detections on the readout->analysis for each readout
-bool psphotDumpChisqs (pmConfig *config, const pmFPAview *view, const char *filerule)
-{
-    static int npass = 0;
-    char filename[64];
-
-    return true;
-
-    bool status = true;
-
-    int num = psphotFileruleCount(config, filerule);
-
-    snprintf (filename, 64, "chisq.%02d.dat", npass);
-    FILE *f = fopen (filename, "w");
-
-    // 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?");
-
-        for (int i = 0; i < sources->n; i++) {
-            pmSource *source = sources->data[i];
-            if (!source) continue;
-
-            pmModel *model = pmSourceGetModel (NULL, source);
-            if (!model) continue;
-        
-            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
-                fprintf (f, "%f %f %f %d %d %f  1 NONLINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm);
-            } else {
-                fprintf (f, "%f %f %f %d %d %f  0 LINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm);
-            }
-        }
-    }
-    fclose (f);
-    npass ++;
-
-    return true;
-}
-
 bool psphotReadout(pmConfig *config, const pmFPAview *view, const char *filerule) {
 
@@ -102 +51 @@
         return psphotReadoutCleanup (config, view, filerule);
     }
-
-# if (0)
-    // XXX test to mask outliers, not very helpful
-    // mask the high values in the image (with MARK)
-    if (!psphotMaskBackground (config, view, filerule)) {
-        return psphotReadoutCleanup (config, view, filerule);
-    }
-
-    // generate a background model (median, smoothed image)
-    if (!psphotModelBackground (config, view, filerule)) {
-        return psphotReadoutCleanup (config, view, filerule);
-    }
-# endif
-
     if (!psphotSubtractBackground (config, view, filerule)) {
         return psphotReadoutCleanup (config, view, filerule);
@@ -137 +72 @@
     }
 
-    // construct sources and measure basic stats (saved on detections->newSources)
+    // construct sources and measure moments and other basic stats (saved on detections->newSources)
+    // all sources use the auto-scaled window appropriate to a PSF, except for the saturated
+    // stars : these use a larger window (3x the basic window)
     if (!psphotSourceStats (config, view, filerule, true)) { // pass 1
         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
@@ -153 +90 @@
 
     // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
-    if (!psphotBasicDeblend (config, view, filerule)) {
+    // XXX I've deactivated this because it was preventing galaxies close to stars from being
+    // XXX fitted as an extended source.
+    if (false && !psphotBasicDeblend (config, view, filerule)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
         return psphotReadoutCleanup (config, view, filerule);
@@ -188 +127 @@
     // psphotLoadExtSources (config, view, filerule); // pass 1
 
-    // construct an initial model for each object, set the radius to fitRadius, set circular
-    // fit mask (detections->newSources)
-    psphotGuessModels (config, view, filerule); // pass 1
-
     // merge the newly selected sources into the existing list
     // NOTE: merge OLD and NEW
     psphotMergeSources (config, view, filerule);
 
+    // Construct an initial model for each object, set the radius to fitRadius, set circular
+    // fit mask.  NOTE: only applied to sources without guess models
+    // pass 1
+    if (!psphotGuessModels (config, view, filerule)) {
+        psLogMsg ("psphot", 3, "failure to Guess Model - pass 1");
+        return psphotReadoutCleanup (config, view, filerule);
+    }
+
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
     psphotFitSourcesLinear (config, view, filerule, false); // pass 1 (detections->allSources)
-    psphotDumpChisqs (config, view, filerule);
-
-    // XXX re-measure the kron mags with models subtracted
-    psphotKronMasked(config, view, filerule);
+
+    // measure the radial profiles to the sky
+    psphotRadialProfileWings (config, view, filerule);
+
+    // re-measure the kron mags with models subtracted.  this pass uses a circular window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
+    
+    // but this is chosen above to be appropriate for the PSF objects (not galaxies)
+    // psphotKronMasked(config, view, filerule);
+    psphotKronIterate(config, view, filerule);
 
     // identify CRs and extended sources (only unmeasured sources are measured)
@@ -211 +159 @@
     // non-linear PSF and EXT fit to brighter sources
     // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
+    // XXX: can leave faulted job in done queue
     psphotBlendFit (config, view, filerule); // pass 1 (detections->allSources)
-    psphotDumpChisqs (config, view, filerule);
 
     // replace all sources
-    psphotReplaceAllSources (config, view, filerule); // pass 1 (detections->allSources)
+    psphotReplaceAllSources (config, view, filerule, false); // pass 1 (detections->allSources)
 
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
     psphotFitSourcesLinear (config, view, filerule, true); // pass 2 (detections->allSources)
-    psphotDumpChisqs (config, view, filerule);
 
     // if we only do one pass, skip to extended source analysis
@@ -240 +187 @@
         psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources)
 
-        // define new sources based on only the new peaks
+        // define new sources based on only the new peaks & measure moments
         // NOTE: new sources are saved on detections->newSources
         psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources)
@@ -251 +198 @@
         }
 
-        // create full input models, set the radius to fitRadius, set circular fit mask
-        // NOTE: apply only to detections->newSources
-        psphotGuessModels (config, view, filerule); // pass 2 (detections->newSources)
-
         // replace all sources so fit below applies to all at once
         // NOTE: apply only to OLD sources (which have been subtracted)
-        psphotReplaceAllSources (config, view, filerule); // pass 2
+        psphotReplaceAllSources (config, view, filerule, false); // pass 2
 
         // merge the newly selected sources into the existing list
@@ -264 +207 @@
         psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources)
 
+        // Construct an initial model for each object, set the radius to fitRadius, set circular
+        // fit mask.  NOTE: only applied to sources without guess models
+        psphotGuessModels (config, view, filerule); // pass 1
+
         // NOTE: apply to ALL sources
         psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources)
-        psphotDumpChisqs (config, view, filerule);
     }
 
@@ -293 +239 @@
         }
 
-        // create full input models, set the radius to fitRadius, set circular fit mask
-        // NOTE: apply only to detections->newSources
-        psphotGuessModels (config, view, filerule); // pass 2 (detections->newSources)
-
         // replace all sources so fit below applies to all at once
         // NOTE: apply only to OLD sources (which have been subtracted)
-        psphotReplaceAllSources (config, view, filerule); // pass 2
+        psphotReplaceAllSources (config, view, filerule, false); // pass 2
 
         // merge the newly selected sources into the existing list
@@ -306 +248 @@
         psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources)
 
+        // Construct an initial model for each object, set the radius to fitRadius, set circular
+        // fit mask.  NOTE: only applied to sources without guess models
+        psphotGuessModels (config, view, filerule); // pass 1
+
         // NOTE: apply to ALL sources
         psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources)
@@ -312 +258 @@
 pass1finish:
 
-    // XXX re-measure the kron mags with models subtracted
-    psphotKronMasked(config, view, filerule);
+    // measure the radial profiles to the sky (only measures new objects)
+    psphotRadialProfileWings (config, view, filerule);
+
+    // re-measure the kron mags with models subtracted
+    // psphotKronMasked(config, view, filerule);
+    psphotKronIterate(config, view, filerule);
 
     // measure source size for the remaining sources
@@ -321 +271 @@
     psphotExtendedSourceAnalysis (config, view, filerule); // pass 1 (detections->allSources)
     psphotExtendedSourceFits (config, view, filerule); // pass 1 (detections->allSources)
-    psphotRadialApertures(config, view, filerule);
+    psphotRadialApertures(config, view, filerule, 0);
 
 finish:
@@ -359 +309 @@
     }
 
+    psLogMsg ("psphot.readout", PS_LOG_WARN, "complete psphot readout : %f sec\n", psTimerMark ("psphotReadout"));
+
     // create the exported-metadata and free local data
     return psphotReadoutCleanup(config, view, filerule);

branches/meh_branches/ppstack_test/psphot/src/psphotReadoutFindPSF.c

@@ -54 +54 @@
     }
 
+    // merge the newly selected sources into the existing list
+    // NOTE: merge OLD and NEW
+    psphotMergeSources (config, view, filerule); 
+
 # if 0
     // XXX if we want to determine the aperture residual correction here, we either
@@ -61 +65 @@
     psphotGuessModels (config, view, filerule);
 # endif
-
-    // merge the newly selected sources into the existing list
-    // NOTE: merge OLD and NEW
-    psphotMergeSources (config, view, filerule); 
 
 # if 0

branches/meh_branches/ppstack_test/psphot/src/psphotReadoutForcedKnownSources.c

@@ -36 +36 @@
     }
 
-    // construct an initial model for each object
-    psphotGuessModels (config, view, filerule);
-
     // merge the newly selected sources into the existing list
     // NOTE: merge OLD and NEW
     psphotMergeSources (config, view, filerule); 
+
+    // Construct an initial model for each object, set the radius to fitRadius, set circular
+    // fit mask.  NOTE: only applied to sources without guess models
+    psphotGuessModels (config, view, filerule);
 
     // linear PSF fit to source peaks

branches/meh_branches/ppstack_test/psphot/src/psphotReadoutKnownSources.c

@@ -48 +48 @@
     }
 
-    // construct an initial model for each object
-    psphotGuessModels (config, view, filerule);
-
     // merge the newly selected sources into the existing list
     // NOTE: merge OLD and NEW
     psphotMergeSources (config, view, filerule); 
+
+    // Construct an initial model for each object, set the radius to fitRadius, set circular
+    // fit mask.  NOTE: only applied to sources without guess models
+    psphotGuessModels (config, view, filerule);
 
     // linear PSF fit to source peaks

branches/meh_branches/ppstack_test/psphot/src/psphotReadoutMinimal.c

@@ -59 +59 @@
     }
 
-    // construct an initial model for each object
-    psphotGuessModels (config, view, filerule);
-
     // merge the newly selected sources into the existing list
     psphotMergeSources (config, view, filerule);
+
+    // Construct an initial model for each object, set the radius to fitRadius, set circular
+    // fit mask.  NOTE: only applied to sources without guess models
+    psphotGuessModels (config, view, filerule);
 
     // linear PSF fit to source peaks

branches/meh_branches/ppstack_test/psphot/src/psphotReplaceUnfit.c

@@ -9 +9 @@
 
     for (int i = 0; i < sources->n; i++) {
-      source = sources->data[i];
-
-      // replace other sources?
-      if (source->mode & PM_SOURCE_MODE_FAIL) goto replace;
-      continue;
+        source = sources->data[i];
+
+        // replace other sources?
+        if (source->mode & PM_SOURCE_MODE_FAIL) goto replace;
+        continue;
 
     replace:
@@ -23 +23 @@
 
 // for now, let's store the detections on the readout->analysis for each readout
-bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule)
+bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState)
 {
     bool status = true;
@@ -35 +35 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (!psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe)) {
+        if (!psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe, ignoreState)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i);
             return false;
@@ -43 +43 @@
 }
 
-bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
-
-    bool status;
-    pmSource *source;
+bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState) {
+
+    bool status;
 
     psTimerStart ("psphot.replace");
@@ -75 +74 @@
 
     for (int i = 0; i < sources->n; i++) {
-      source = sources->data[i];
-
-      // replace other sources?
-      if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
-
-      pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+        pmSource *source = sources->data[i];
+
+        if (ignoreState) {
+            // rely on the type of source to decide if we subtract it or not
+
+            // skip non-astronomical objects (very likely defects)
+            if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+            if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+      
+            // do not include CRs in the full ensemble fit
+            if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
+        
+            // do not include MOMENTS_FAILURES in the fit
+            if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+        } else {
+            // if we respect the state, do not replace unsubtracted sources
+            if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
+        }
+
+        pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     }
 
@@ -88 +101 @@
 }
 
-bool psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe) {
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotRemoveAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState)
+{
+    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 (!psphotRemoveAllSourcesReadout (config, view, filerule, i, recipe, ignoreState)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i);
+            return false;
+        }
+    }
+    return true;
+}
+
+bool psphotRemoveAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState) {
+
+    bool status;
+
+    psTimerStart ("psphot.replace");
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->allSources;
+    psAssert (sources, "missing sources?");
+
+    // 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?");
+
+    // bit-mask to mark pixels not used in analysis
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
+    assert (markVal);
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+
+        if (ignoreState) {
+            // rely on the type of source to decide if we subtract it or not
+
+            // skip non-astronomical objects (very likely defects)
+            if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+            if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+      
+            // do not include CRs in the full ensemble fit
+            if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
+        
+            // do not include MOMENTS_FAILURES in the fit
+            if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+        } else {
+            // if we respect the state, only remove unsubtracted sources
+            if ((source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
+        }
+
+        pmSourceSub (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;
+}
+
+bool psphotRemoveAllSourcesByArray (const psArray *sources, const psMetadata *recipe) {
 
     bool status;
@@ -100 +192 @@
 
     for (int i = 0; i < sources->n; i++) {
-      source = sources->data[i];
-
-      // replace other sources?
-      if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
-
-      pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+        source = sources->data[i];
+
+        // replace other sources?
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
+
+        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     }
     psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
@@ -158 +250 @@
 
     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;
-      if (!source->modelPSF) continue;
-
-      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);
+        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;
+        if (!source->modelPSF) continue;
+
+        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;
@@ -229 +321 @@
 
     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);
-      if (!model) continue;
-
-      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->rawFlux;
-          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);
+        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);
+        if (!model) continue;
+
+        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->rawFlux;
+            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;
+            model = source->modelPSF;
+        }
+
+        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);
+            // 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 (!)
+            psFree (psfKernel);
+            psFree (rawModelFlux);
+        }
+
+        pmSourceCacheModel (source, maskVal);  // ALLOC x14 (!)
     }
 

branches/meh_branches/ppstack_test/psphot/src/psphotRoughClass.c

@@ -11 +11 @@
 {
     bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Rough Class ---");
 
     // select the appropriate recipe information
@@ -128 +131 @@
     psphotDumpMoments (recipe, sources);
 
-    psLogMsg ("psphot.roughclass", PS_LOG_INFO, "rough classification: %f sec\n", psTimerMark ("psphot.rough"));
+    psLogMsg ("psphot.roughclass", PS_LOG_WARN, "rough classification: %f sec\n", psTimerMark ("psphot.rough"));
 
     psphotVisualPlotMoments (recipe, readout->analysis, sources);

branches/meh_branches/ppstack_test/psphot/src/psphotSetThreads.c

@@ -30 +30 @@
     psFree(task);
 
+    task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 7);
+    task->function = &psphotKronIterate_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
     task = psThreadTaskAlloc("PSPHOT_BLEND_FIT", 10);
     task->function = &psphotBlendFit_Threaded;
@@ -40 +45 @@
     psFree(task);
 
+    task = psThreadTaskAlloc("PSPHOT_EXTENDED_ANALYSIS", 8);
+    task->function = &psphotExtendedSourceAnalysis_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 6);
+    task->function = &psphotRadialApertures_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_RADIAL_PROFILE_WINGS", 3);
+    task->function = &psphotRadialProfileWings_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
     return true;
 }

branches/meh_branches/ppstack_test/psphot/src/psphotSourceFits.c

@@ -211 +211 @@
 bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
 
-    float radius;
+    float fitRadius, windowRadius;
     bool okEXT, okDBL;
     pmModel *ONE = NULL;
@@ -217 +217 @@
     pmModel *EXT = NULL;
     psArray *DBL = NULL;
-    pmMoments psfMoments;
+    // pmMoments psfMoments;
 
     // skip the source if we don't think it is extended
@@ -225 +225 @@
     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
 
+# define TEST_X -420.0
+# define TEST_Y 300.0
+   
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        fprintf (stderr, "test galaxy\n");
+    }
+
+# undef TEST_X
+# undef TEST_Y
+
     // set the radius based on the footprint (also sets the mask pixels)
-    if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false;
-
-    // XXX note that this changes the source moments that are published...
-    // XXX all published moments should use the same measurement
-    // recalculate the source moments using the larger extended-source moments radius
-    // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
-    // this uses the footprint to judge both radius and aperture?
-    // XXX save the psf-based moments for output
-    psfMoments = *source->moments;
-    if (!pmSourceMoments (source, radius, 0.0, 0.5, 0.0, maskVal)) {
-      *source->moments = psfMoments;
-      return false;
-    }
+    if (!psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal)) return false;
+    // fprintf (stderr, "rad: %6.1f %6.1f  | %5.2f %5.2f %5.2f  ", source->peak->xf, source->peak->yf, source->moments->Mrf, fitRadius, windowRadius);
 
     psTrace ("psphot", 5, "trying blob...\n");
@@ -258 +257 @@
         okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
         okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
+        okDBL = false; // XXX this is failing badly...
         // XXX should I keep / save the flags set in the eval functions?
 
@@ -263 +263 @@
         ONE = DBL->data[0];
         if (ONE) {
-            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
+            psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
             chiDBL = ONE->chisqNorm; // save chisq for double-star/galaxy comparison
-            ONE->fitRadius = radius;
+            ONE->fitRadius = fitRadius;
         }
 
@@ -271 +271 @@
         ONE = DBL->data[1];
         if (ONE) {
-            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
-            ONE->fitRadius = radius;
+            psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
+            ONE->fitRadius = fitRadius;
         }
     }
@@ -284 +284 @@
         okEXT = psphotEvalEXT (tmpSrc, EXT);
         chiEXT = EXT ? EXT->chisqNorm : NAN;
+        EXT->fitRadius = fitRadius;
     }
 
@@ -294 +295 @@
 
     if (okEXT && okDBL) {
-        psTrace ("psphot", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL);
         // XXX EAM : a bogus bias: need to examine this better
         if (3*chiEXT > chiDBL) goto keepDBL;
@@ -303 +303 @@
     if (!okEXT && okDBL) goto keepDBL;
 
+    psTrace ("psphot", 4, "both failed: blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
+
     // both models failed; reject them both
     // XXX -- change type flags to psf in this case, and make sure we subtract it?
     // reset the psf moments
-    *source->moments = psfMoments;
+    // XXX *source->moments = psfMoments;
 
     psFree (EXT);
@@ -313 +315 @@
 
 keepEXT:
+    psTrace ("psphot", 4, "goto EXT : blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
     // sub EXT
     psFree (DBL);
@@ -338 +341 @@
 
     // reset the psf moments
-    *source->moments = psfMoments;
+    // XXX *source->moments = psfMoments;
     return true;
 
 keepDBL:
+    psTrace ("psphot", 4, "goto DBL : blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
     // sub DLB
     psFree (EXT);
@@ -354 +358 @@
     pmSource *newSrc = pmSourceCopy (source);
     newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
+    newSrc->peak->footprint = source->peak->footprint; // just a reference; the peak does not own the footprint
+    psArrayAdd(newSrc->peak->footprint->peaks, 1, newSrc->peak); // the footprint owns the peak
 
     // build cached models and subtract
@@ -372 +378 @@
             psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelPSF->params->data.F32[i], source->modelPSF->dparams->data.F32[i]);
         }
-        psphotVisualShowResidualImage (readout, false);
     }
 # endif
-
-    // reset the (original) psf moments
-    *source->moments = psfMoments;
-    *newSrc->moments = psfMoments;
 
     psArrayAdd (newSources, 100, newSrc);
@@ -386 +387 @@
 
 escape:
-    // reset the psf moments
-    *source->moments = psfMoments;
     psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
     psFree (tmpSrc);
@@ -476 +475 @@
     // for sersic models, use a grid search to choose an index, then float the params there
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-        // for the test fits, use a somewhat smaller radius
-        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
         psphotFitSersicIndex (model, readout, source, fitOptions, maskVal, markVal);
-    }
-
-    if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
-        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
     }
 
@@ -494 +487 @@
     pmSourceFitModel (source, model, &options, maskVal);
     // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
-
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+
     return (model);
 }
+
+# define TIMING 0
 
 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
@@ -517 +512 @@
         return NULL;
     }
+
+# define TEST_X -540.0
+# define TEST_Y 540.0
+   
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
+    }
+
+    float t1, t2, t4, t5;
+    if (TIMING) { psTimerStart ("psphotFitPCM"); }
 
     pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
@@ -524 +529 @@
         return model;
     }
-
-    // use the source moments, etc to guess basic model parameters
-    if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
-        psFree (pcm);
-        model->flags |= PM_MODEL_STATUS_BADARGS;
-        return model;
-    }
-
-    // for sersic models, use a grid search to choose an index, then float the params there
+    if (TIMING) { t1 = psTimerMark ("psphotFitPCM"); }
+
+    // NOTE : 65 allocs to here
+    // get the guess for sersic models 
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-        // for the test fits, use a somewhat smaller radius
-        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
-
-        if (!psphotFitSersicIndexPCM (pcm, readout, source, fitOptions, maskVal, markVal, psfSize)) {
+        // use the source moments, etc to guess basic model parameters
+        if (!psphotSersicModelClassGuessPCM (pcm, source)) {
             psFree (pcm);
             model->flags |= PM_MODEL_STATUS_BADARGS;
             return model;
         }
-    }
-
-    if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
-        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
-    }
+    } else {
+        // use the source moments, etc to guess basic model parameters
+        if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
+            psFree (pcm);
+            model->flags |= PM_MODEL_STATUS_BADARGS;
+            return model;
+        }
+    }
+
+    if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); }
 
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
@@ -555 +558 @@
     // update the pcm elements if we have changed the circumstance (options.mode or source->pixels)
     pmPCMupdate(pcm, source, &options, model);
+    if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
-    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
-    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);  // NOTE : 1687 allocs in here
+    if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }
+
+    if (TIMING) {
+        int nPixBig = source->pixels->numCols * source->pixels->numRows;
+        fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t4, t5);
+    }
+
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
+    }
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
@@ -566 +579 @@
 }
 
+# undef TEST_X
+# undef TEST_Y
+
 // note that these should be 1/2n of the standard sersic index
-float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
-# define N_INDEX_GUESS 6
+// float indexGuess[] = {0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0};
+// float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
+float indexGuess[] = {1.0, 2.0, 3.0, 4.0};
+# define N_INDEX_GUESS 4
 
 // A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
@@ -586 +604 @@
     float chiSquare[N_INDEX_GUESS];
 
+# define TEST_X -540.0
+# define TEST_Y 540.0
+   
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
+    }
+
     for (int i = 0; i < N_INDEX_GUESS; i++) {
-        model->params->data.F32[PM_PAR_7] = indexGuess[i];
+        model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[i];
 
         if (!model->modelGuess(model, source)) {
@@ -594 +619 @@
         }
 
-        // each time we change the model guess, we need to adjust the radius
-        // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
-        if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
-            psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
-        }
-        
         pmSourceFitModel (source, model, &options, maskVal);
-        // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+        // fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
 
         chiSquare[i] = model->chisqNorm;
@@ -616 +635 @@
     assert (iMin >= 0);
 
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
+    }
+
     model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
-    model->params->data.F32[PM_PAR_7] = indexGuess[iMin];
+    model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[iMin];
     model->modelGuess(model, source);
-
-    // each time we change the model guess, we need to adjust the radius
-    // if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
-    //  psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal);
-    // }
 
     return true;
@@ -648 +666 @@
     float xMin = NAN;
     float chiSquare[N_INDEX_GUESS];
+
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
+    }
 
     for (int i = 0; i < N_INDEX_GUESS; i++) {
@@ -657 +679 @@
         }
 
-        // each time we change the model guess, we need to adjust the radius
-        // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
-        if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
-            psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
-        }
-        
 # if (0)
+        // this block is to test the relative speed of straight and PCM fits
         pmSourceFitModel (source, model, &options, maskVal);
 # else
@@ -669 +686 @@
         pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
 # endif
+        fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
         // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
 
@@ -684 +702 @@
     assert (iMin >= 0);
     
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
+    }
+
     model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
     model->params->data.F32[PM_PAR_7] = indexGuess[iMin];

branches/meh_branches/ppstack_test/psphot/src/psphotSourceMatch.c

@@ -48 +48 @@
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     psAssert (detections, "missing detections?");
-    psAssert (detections->newSources, "new sources not defined?");
-    psAssert (!detections->allSources, "all sources already defined?");
-
-    // XXX TEST:
-    if (detections->newSources) {
-        psphotMatchSourcesToObjects(objects, detections->newSources, RADIUS);
-    }
+    psAssert (detections->allSources, "all sources not defined?");
+
+    psphotMatchSourcesToObjects(objects, detections->allSources, RADIUS);
 
     return true;
@@ -71 +67 @@
     objects = psArraySort (objects, pmPhotObjSortByX); 
  
-    psVector *found = psVectorAlloc(sources->n, PS_TYPE_U8);
-    psVectorInit (found, 0);
+    psVector *foundSrc = psVectorAlloc(sources->n, PS_TYPE_U8);
+    psVectorInit (foundSrc, 0);
+
+    psVector *foundObj = psVectorAlloc(objects->n, PS_TYPE_U8);
+    psVectorInit (foundObj, 0);
 
     // match sources to existing objects
@@ -97 +96 @@
         if (dx > +1.02*RADIUS) NEXT2; 
  
+        /* this block will match a given detection to the closest object within range of that detection.
+           XXX note that this matches ALL detections within range of the single object to that same object 
+           this is bad, but I cannot just go in linear order (ie, mark off each object as they are
+           used).  I should make a list of all Nobj * Ndet pairs in range and choose the matches
+           based on their separations.  UGH
+        */
+    
         // we are within match range, look for matches: 
         int Jmin = -1;
@@ -102 +108 @@
         for (int J = j; (dx > -1.02*RADIUS) && (J < objects->n); J++) { 
  
+            // skip objects that are already assigned:
+            if (foundObj->data.U8[J]) continue;
             obj = objects->data[J]; 
             
@@ -121 +129 @@
         }
         obj = objects->data[Jmin]; 
+        foundObj->data.U8[Jmin] = 1;
 
         // add to object
         pmPhotObjAddSource (obj, src);
-        found->data.U8[i] = 1;
+        foundSrc->data.U8[i] = 1;
         i++; 
     } 
@@ -132 +141 @@
     for (i = 0; i < sources->n; i++) {
 
-        if (found->data.U8[i]) continue;
+        if (foundSrc->data.U8[i]) continue;
 
         pmSource *src = sources->data[i]; 
@@ -143 +152 @@
     psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n);
 
-    psFree (found);
+    psFree (foundSrc);
+    psFree (foundObj);
     return true;
 } 
@@ -238 +248 @@
             
             // 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);
+            if (footprint) {
+                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
@@ -249 +261 @@
             source->imageID = index;
             source->mode2 |= PM_SOURCE_MODE2_MATCHED; // source is generated based on another image
+            source->type = PM_SOURCE_TYPE_STAR; // until we know more, assume a PSF fit
 
             // add the peak
@@ -261 +274 @@
             peak->assigned = true;
             pmPhotObjAddSource(obj, source);
-            psArrayAdd (detections->newSources, 100, source);
+            psArrayAdd (detections->allSources, 100, source);
             psFree (source);
         }
@@ -272 +285 @@
         pmPhotObj *obj = objects->data[i]; 
         nSources += obj->sources->n;
+        psAssert (obj->sources->n == nImages, "failed to match sources?");
     }
     psLogMsg ("psphot", PS_LOG_DETAIL, "total of %d sources for %d images", nSources, nImages);

branches/meh_branches/ppstack_test/psphot/src/psphotSourceSize.c

@@ -45 +45 @@
     bool status = true;
 
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Source Size ---");
+
     // select the appropriate recipe information
     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
@@ -153 +156 @@
 
     // XXX fix this (was source->n  - first)
-    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.size"));
+    psLogMsg ("psphot.size", PS_LOG_WARN, "measure source sizes for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.size"));
 
     psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
@@ -209 +212 @@
         pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
 
-        float kMag = -2.5*log10(source->moments->KronFlux);
+        float kMag = -2.5*log10(source->moments->KronFluxPSF);
         float dMag = source->psfMag - kMag;
 
@@ -334 +337 @@
         psF32 Mxy = source->moments->Mxy;
 
-        float KronMag = -2.5*log10(source->moments->KronFlux);
+        float KronMag = -2.5*log10(source->moments->KronFluxPSF);
 
         float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
@@ -509 +512 @@
         pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
 
-        float kMag = -2.5*log10(source->moments->KronFlux);
+        float kMag = -2.5*log10(source->moments->KronFluxPSF);
         float dMag = source->psfMag - kMag;
 

branches/meh_branches/ppstack_test/psphot/src/psphotSourceStats.c

@@ -10 +10 @@
     bool status = true;
 
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Source Stats ---");
+
     // select the appropriate recipe information
     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
@@ -98 +101 @@
 
     // generate the array of sources, define the associated pixel
-    bool firstPass = false;
-    if (!detections->newSources) {
-        detections->newSources = psArrayAllocEmpty (peaks->n);
-        firstPass = true;
-    }
-    sources = detections->newSources;
+    psAssert (!detections->newSources, "programming error");
+    sources = detections->newSources = psArrayAllocEmpty (peaks->n);
+
+    bool firstPass = (detections->allSources == NULL);
 
     // if there are no peaks, we save the empty source array and return
@@ -121 +122 @@
         // add the peak
         source->peak = psMemIncrRefCounter(peak);
+
+        // psAssert (source->peak->footprint, "peak without footprint??");
 
         // allocate space for moments
@@ -213 +216 @@
 
         // wait for the threads to finish and manage results
-        if (!psThreadPoolWait (false)) {
+        if (!psThreadPoolWait (false, true)) {
             psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
             psFree(detections->newSources);
@@ -239 +242 @@
     psFree (cellGroups);
 
-    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats"));
+    psLogMsg ("psphot", PS_LOG_WARN, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats"));
 
     psphotVisualShowMoments (sources);
@@ -344 +347 @@
 
         // wait for the threads to finish and manage results
-        if (!psThreadPoolWait (false)) {
+        if (!psThreadPoolWait (false, true)) {
             psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
             return NULL;
@@ -554 +557 @@
         // determine the PSF parameters from the source moment values
         pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX);
-        psLogMsg ("psphot", 3, "radius %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
+        psLogMsg ("psphot", 3, "sigma guess (pix) %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
 
         Rmin[i] = pmSourceMinKronRadius(sources, PSF_SN_LIM);

branches/meh_branches/ppstack_test/psphot/src/psphotStackChisqImage.c

@@ -13 +13 @@
     psAssert (chisqFile, "missing chisq image FPA?");
 
+    // the readout containing the chisq image is generated in the first pass of this loop and
+    // used by the successive passes
     pmReadout *chiReadout = NULL;
 

branches/meh_branches/ppstack_test/psphot/src/psphotStackImageLoop.c

@@ -10 +10 @@
 bool UpdateHeadersForFPA (pmConfig *config, pmFPAview *view);
 bool UpdateHeadersForChip (pmConfig *config, pmFPAview *view);
+bool UpdateHeadersForReadout (pmConfig *config, pmFPAview *view);
 
 bool psphotStackImageLoop (pmConfig *config) {
@@ -52 +53 @@
                 psMemDump("load");
 
-                // PSF matching
+                // Generate the 1st PSF-matched image set (larger target PSFs are generated by smoothing this image)
                 if (!psphotStackMatchPSFs (config, view)) {
                     psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
@@ -66 +67 @@
                     return false;
                 }
+
+                UpdateHeadersForReadout(config, view);
+
                 psMemDump("psphot");
             }
@@ -101 +105 @@
 }
 
+bool UpdateHeadersForReadout (pmConfig *config, pmFPAview *view) {
+
+    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
+    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
+
+    int num = psphotFileruleCount(config, "PSPHOT.INPUT");
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+
+        // find the currently selected readout
+        pmFPAfile *output = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); // File of interest
+        psAssert (output, "missing file?");
+
+        pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
+        pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
+        pmFPAfile *input = useRaw ? inputRaw : inputCnv;
+        psAssert (input, "missing input file");
+
+        // just copy the input headers to the output headers, then update version info
+        pmReadout *inReadout = pmFPAviewThisReadout(view, input->fpa); ///< Chip in the input
+        pmReadout *outReadout = pmFPAviewThisReadout(view, output->fpa); ///< Chip in the output
+
+        psMetadata *header = psMetadataLookupPtr (NULL, inReadout->analysis, "PSPHOT.HEADER");
+        psMetadataAdd (outReadout->analysis, PS_LIST_TAIL, "PSPHOT.HEADER",  PS_DATA_METADATA, "header stats", header);
+    }
+    return true;
+}
+
 bool UpdateHeadersForChip (pmConfig *config, pmFPAview *view) {
 
@@ -131 +164 @@
             outHDU->header = psMetadataCopy(NULL, inHDU->header);
         }
+        psphotVersionHeaderFull(outHDU->header);
         outChip->toFPA = psMemIncrRefCounter(inChip->toFPA);
         outChip->fromFPA = psMemIncrRefCounter(inChip->fromFPA);

branches/meh_branches/ppstack_test/psphot/src/psphotStackMatchPSFs.c

@@ -25 +25 @@
     }
 
-    // loop over the available readouts (ignore chisq image)
+    // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        // set up the PSF-matching parameters describing the input images
         if (!psphotStackMatchPSFsPrepare (config, view, options, i)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to set PSF matching options for entry %d", i);
@@ -33 +34 @@
     }
 
-    // Generate target PSF
+    // XXX convolve == false might not be valid at the moment
     if (options->convolve) {
-        options->psf = psphotStackPSF(config, options->numCols, options->numRows, options->psfs, options->inputMask);
-        if (!options->psf) {
+        // Determine the 1st target PSF (either AUTO or defined by PSPHOT.STACK.TARGET.PSF.FWHM)
+        // NOTE: this also set the full list of target FWHMs (options->targetSeeing)
+        if (!psphotStackPSF(config, options)) {
             psError(psErrorCodeLast(), false, "Unable to determine output PSF.");
             return false;
         }
-        psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf);
-        options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * options->numCols, 0.5 * options->numRows); // FWHM for target
-        psLogMsg("psphotStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing);
-
-        // XXX is this needed to supply the psf to psphot??
-        // pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.TARGET.PSF"); // Output chip
-        // psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "Target PSF", options->psf);
-        // outChip->data_exists = true;
     }
 
@@ -63 +57 @@
 
 // convolve the image to match desired PSF
+// XXX is this code consistent with 'convolve' = false?
 bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) {
 
@@ -106 +101 @@
     // Image Matching (PSFs or just flux)
     if (options->convolve) {
-        matchKernel(config, readoutOut, readoutSrc, options, index);
+        if (!matchKernel(config, readoutOut, readoutSrc, options, index)){
+            psError(psErrorCodeLast(), false, "Unable to match image PSF in readout.");
+            return false;
+        }
         saveMatchData(readoutOut, options, index);
     }
-    rescaleData(readoutOut, config, options, index);
+
+    // renormalize the stack variances to have sigma = 1.0
+    if (!psphotStackRenormaliseVariance(config, readoutOut)) {
+        psError(psErrorCodeLast(), false, "Unable to renormalise variance.");
+        return false;
+    }
 
     // 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
+    // NOTE: fwhmValues as defined here has 1 + nMatched PSF : 0 == unmatched
+    psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
+    psVectorAppend(fwhmValues, NAN); // XXX this corresponds to the unmatched image set
+    for (int i = 0; i < options->targetSeeing->n; i++) {
+        psVectorAppend(fwhmValues, options->targetSeeing->data.F32[i]);
     }
-    psVectorAppend(fwhmValues, options->targetSeeing);
+    psMetadataAddVector(readoutSrc->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
+    psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
+    psFree(fwhmValues); // drops the extra copy
 
     return true;

branches/meh_branches/ppstack_test/psphot/src/psphotStackMatchPSFsNext.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+
+// NOTE : element 0 of fwhmValues if the unmatched image,  
+
+int psphotStackMatchPSFsEntries (pmConfig *config, const pmFPAview *view, const char *filerule) {
+
+    int nRadialEntries = 0;
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest
+    psAssert (file, "missing file?");
+    
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+    
+    bool status = false;
+    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
+    if (!fwhmValues) {
+        return 1;
+    }
+    
+    nRadialEntries = fwhmValues->n;
+    return nRadialEntries;
+}
 
 // smooth the input image to match the next target PSF
@@ -5 +28 @@
 // and that the smoothing can use a 1D Gaussian kernel of width sqrt(TARGET^2 - CURRENT^2)
 // each subsequent call
-bool psphotStackMatchPSFsNext(bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize)
+bool psphotStackMatchPSFsNext(pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize)
 {
-    bool status = true;
-
-    // select the appropriate recipe information
-    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
-    psAssert (recipe, "missing recipe?");
-
-    psVector *fwhmValues = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
-    if (!status) {
-        // must not be a vector, only one value requested
-        *smoothAgain = false;
-        return true;
-    }
-
-    if (lastSize + 1 >= fwhmValues->n) {
-        // all done with target FWHM values 
-        *smoothAgain = false;
-        return true;
-    }
-
     int num = psphotFileruleCount(config, filerule);
 
@@ -33 +37 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, recipe, fwhmValues, lastSize)) {
-            psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF %f", filerule, i, fwhmValues->data.F32[lastSize+1]);
+        if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, lastSize)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF", filerule, i);
             psImageConvolveSetThreads(oldThreads);
             return false;
@@ -44 +48 @@
 }
 
-bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, psVector *fwhmValues, int lastSize) {
+bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, int lastSize) {
 
     bool status = false;
+
+    psTimerStart ("psphot.smooth");
 
     // find the currently selected readout
@@ -58 +64 @@
     psphotVisualShowImage(readout);
 
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
     // 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
@@ -66 +76 @@
         psWarning("PEAKS_MIN_GAUSS is not set in recipe; using default value");
         minGauss = 0.5;
+    }
+
+    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
+    psAssert (fwhmValues, "need target PSFs");
+
+    if (lastSize + 1 >= fwhmValues->n) {
+        return true;
     }
 
@@ -128 +145 @@
     // psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor);
 
-    // 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
-
-    psVector *fwhmValuesOut = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
-    psAssert(fwhmValuesOut, "should already exist..");
-    psVectorAppend(fwhmValuesOut, targetFWHM);
-
     // do not generate a PSF if we already were supplied one
     pmPSF *psfOld = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");

branches/meh_branches/ppstack_test/psphot/src/psphotStackMatchPSFsPrepare.c

@@ -4 +4 @@
 bool determineSeeing (pmPSF *psf, psphotStackOptions *options, int index);
 
-// set up the stacking parameters
+// set up the PSF-matching parameters describing the input images
 bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) {
 

branches/meh_branches/ppstack_test/psphot/src/psphotStackMatchPSFsUtils.c

@@ -2 +2 @@
 # define ARRAY_BUFFER 16                 // Number to add to array at a time
 
-// XXX better name
-bool readImage(psImage **target, // Target for image
-               const char *name, // Name of FITS file
-               const pmConfig *config // Configuration
-    )
-{
-    psString resolved = pmConfigConvertFilename(name, config, false, false); // Resolved filename
-    psFits *fits = psFitsOpen(resolved, "r");
-    psFree(resolved);
-    if (!fits) {
-        psError(PSPHOT_ERR_IO, false, "Unable to open previously produced image: %s", name);
-        return false;
-    }
-    psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest
-    if (!image) {
-        psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image: %s", name);
-        psFitsClose(fits);
-        return false;
-    }
-    psFitsClose(fits);
-
-    psFree(*target);
-    *target = image;
-
-    return true;
-}
+psVector *SetOptWidths (bool *optimum, psMetadata *recipe);
+pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize);
+bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index);
+bool matchKernel(pmConfig *config, pmReadout *cnv, pmReadout *raw, psphotStackOptions *options, int index);
+bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname);
+bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname);
 
 // Get coordinates from a source
@@ -148 +128 @@
 
 // Renormalise a readout's variance map
-bool stackRenormaliseReadout(const pmConfig *config, // Configuration
-                             pmReadout *readout      // Readout to renormalise
+bool psphotStackRenormaliseVariance(const pmConfig *config, // Configuration
+                              pmReadout *readout      // Readout to renormalise
     )
 {
@@ -180 +160 @@
     return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
 }
-
-// This is a hack to use the temporary convolved images and kernel generated previously.
-// This makes the 'matching' operation much faster, allowing debugging of the stack process easier.
-// It implicitly assumes the output root name is the same between invocations.
-
-# if (0)
-bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index) {
-
-    // Read the convolution kernel from the saved file
-    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index);
-    psAssert(file, "Require file");
-
-    pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest
-    view->chip = view->cell = view->readout = 0;
-    psString filename = pmFPAfileNameFromRule(file->filerule, file, view); // Filename of interest
-
-    // Read convolution kernel data
-    psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
-    psFree(filename);
-    psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
-    psFree(resolved);
-    if (!fits || !pmReadoutReadSubtractionKernels(readoutCnv, fits)) {
-        psError(PSPHOT_ERR_IO, false, "Unable to read previously produced kernel");
-        psFitsClose(fits);
-        return false;
-    }
-    psFitsClose(fits);
-
-    // read the convolved pixels (image, mask, variance) -- names are pre-defined
-    if (!readImage(&readoutCnv->image,    options->convImages->data[index],    config) ||
-        !readImage(&readoutCnv->mask,     options->convMasks->data[index],     config) ||
-        !readImage(&readoutCnv->variance, options->convVariances->data[index], config)) {
-        psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image.");
-        return false;
-    }
-
-    // XXX ??? not sure what is happening here -- consult Paul
-    psRegion *region = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
-    pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
-
-    pmSubtractionAnalysis(readoutCnv->analysis, NULL, kernels, region, readoutCnv->image->numCols, readoutCnv->image->numRows);
-
-    psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
-
-    // update the covariance matrix 
-    // XXX why is this needed if we have correctly read the saved data?
-    bool oldThreads = psImageCovarianceSetThreads(true);              // Old thread setting
-    psKernel *covar = psImageCovarianceCalculate(kernel, readoutCnv->covariance); // Covariance matrix
-    psImageCovarianceSetThreads(oldThreads);
-    psFree(readoutCnv->covariance);
-    readoutCnv->covariance = covar;
-    psFree(kernel);
-    return true;
-}
-# endif
 
 bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname) {
@@ -362 +287 @@
         widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths
 
-        // we need to register the FWHM values for use downstream 
-        pmSubtractionSetFWHMs(options->targetSeeing, options->inputSeeing->data.F32[index]);
+        // we need to register the FWHM values for use by pmSubtraction code
+        pmSubtractionSetFWHMs(options->targetSeeing->data.F32[0], options->inputSeeing->data.F32[index]);
 
         pmSubtractionParamScaleOptions(scale, scaleRef, scaleMin, scaleMax);
@@ -469 +394 @@
 }
 
-// Kernel normalisation for convolved readout
-bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index) {
-
-    double sum = 0.0;           // Sum of chi^2
-    int num = 0;                // Number of measurements of chi^2
-    psString regex = NULL;      // Regular expression
-    psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM);
-    psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex);
-    psFree(regex);
-    psMetadataItem *item = NULL;// Item from iteration
-    while ((item = psMetadataGetAndIncrement(iter))) {
-        assert(item->type == PS_TYPE_F32);
-        float norm = item->data.F32; // Normalisation
-        sum += norm;
-        num++;
-    }
-    psFree(iter);
-    float conv = sum/num;       // Mean normalisation from convolution
-    float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars
-    float renorm =  stars / conv; // Renormalisation to apply
-    psLogMsg("psphotStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n", index, renorm, conv, stars);
-
-    psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32));
-    psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32));
-    return true;
-}
-
-// adjust scaling for readout (remove background, ..., determine weighting)
-bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index) {
-
-    psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe
-    psAssert(stackRecipe, "We've thrown an error on this before.");
-
-    // Look up appropriate values from the ppSub recipe
-    psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
-    psAssert(subRecipe, "recipe missing");
-
-    psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in
-    psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad
-
-    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
-    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
-
-    // Ensure the background value is zero
-    psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
-
-    // XXX why is this in config->arguments and not recipe?
-    if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
-        if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
-            psAbort("Can't measure background for image.");
-            // XXX we used to clear error: why is this acceptable? psErrorClear(); 
-        }
-
-        float value = psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN);
-        float stdev = psStatsGetValue(bg, PS_STAT_ROBUST_STDEV);
-
-        psLogMsg("psphotStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", value, stdev);
-        psBinaryOp(readout->image, readout->image, "-", psScalarAlloc(value, PS_TYPE_F32));
-    }
-
-    if (!stackRenormaliseReadout(config, readout)) {
-        psFree(rng);
-        psFree(bg);
-        return false;
-    }
-
-    // Measure the variance level for the weighting
-    if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) {
-        if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
-            psError(PSPHOT_ERR_DATA, false, "Can't measure mean variance for image.");
-            psFree(rng);
-            psFree(bg);
-            return false;
-        }
-        options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance));
-    } else {
-        options->weightings->data.F32[index] = 1.0;
-    }
-    psLogMsg("psphotStack", PS_LOG_INFO, "Weighting for image %d is %f\n", index, options->weightings->data.F32[index]);
-
-    psFree(rng);
-    psFree(bg);
-    return true;
-}
-
 # define NOISE_FRACTION 0.01             // Set minimum flux to this fraction of noise
 # define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources

branches/meh_branches/ppstack_test/psphot/src/psphotStackObjects.c

@@ -47 +47 @@
     return true;
 }
+
+// mark good vs bad objects
+bool psphotStackObjectsSelectForAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
+
+    bool status = false;
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // use the 0-index image to represent the image area
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    // option to limit analysis to a specific region
+    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
+    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
+
+    // S/N limit to perform full non-linear fits
+    float SN_LIM_PETRO  = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
+    float SN_LIM_RADIAL = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
+
+    bool doPetroStars   = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS");
+
+    for (int i = 0; i < objects->n; i++) {
+        pmPhotObj *object = objects->data[i];
+        if (!object) continue;
+        if (!object->sources) continue;
+
+        // we check each source for an object and keep the object if any source is valid
+
+        bool keepObjectRadial = false;
+        bool keepObjectPetro = false;
+        for (int j = 0; j < object->sources->n; j++) {
+
+            pmSource *source = object->sources->data[j];
+            if (!source) continue;
+            if (!source->peak) continue;
+
+            // skip PSF-like and non-astronomical objects
+            if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+            if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+            if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
+            if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
+            
+            // limit selection by analysis region (this automatically apply
+            if (source->peak->x < AnalysisRegion.x0) continue;
+            if (source->peak->y < AnalysisRegion.y0) continue;
+            if (source->peak->x > AnalysisRegion.x1) continue;
+            if (source->peak->y > AnalysisRegion.y1) continue;
+            
+            // SN limit tests for RADIAL APERTURES:
+            bool skipSourceRadial = false;
+            if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+                skipSourceRadial = (source->moments->KronFlux < SN_LIM_RADIAL * source->moments->KronFluxErr);
+            } else {
+                skipSourceRadial = (sqrt(source->peak->detValue) < SN_LIM_RADIAL);
+            }
+            if (!skipSourceRadial) keepObjectRadial = true;
+
+            // SN limit tests for PETRO
+            bool skipSourcePetro = false;
+            if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+                skipSourcePetro = (source->moments->KronFlux < SN_LIM_PETRO * source->moments->KronFluxErr);
+            } else {
+                skipSourcePetro = doPetroStars ? (sqrt(source->peak->detValue) < SN_LIM_PETRO) : true;
+            }
+            if (!skipSourcePetro) keepObjectPetro = true;
+
+            keepObjectPetro = true;
+        }
+
+        for (int j = 0; j < object->sources->n; j++) {
+            pmSource *source = object->sources->data[j];
+            if (!source) continue;
+            if (!source->peak) continue;
+
+            // we have to set a bit in either case to tell psphotExtendedSourceAnalysis to
+            // avoid the single-detection tests
+
+            if (keepObjectPetro) {
+                source->tmpFlags |=  PM_SOURCE_TMPF_PETRO_KEEP;
+                source->tmpFlags &= ~PM_SOURCE_TMPF_PETRO_SKIP;
+            } else {
+                source->tmpFlags |=  PM_SOURCE_TMPF_PETRO_SKIP;
+                source->tmpFlags &= ~PM_SOURCE_TMPF_PETRO_KEEP;
+            }       
+
+            if (keepObjectRadial) {
+                source->tmpFlags |=  PM_SOURCE_TMPF_RADIAL_KEEP;
+                source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_SKIP;
+            } else {
+                source->tmpFlags |=  PM_SOURCE_TMPF_RADIAL_SKIP;
+                source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_KEEP;
+            }       
+        }
+    }
+
+    psLogMsg ("psphot", PS_LOG_INFO, "marked good vs bad objects\n");
+    return true;
+}

branches/meh_branches/ppstack_test/psphot/src/psphotStackOptions.c

@@ -19 +19 @@
     psFree (options->norm);
     psFree (options->matchChi2);
-    psFree (options->weightings);
+    psFree (options->targetSeeing);
 
     return;
@@ -36 +36 @@
     options->convolve = false;
     options->convolveSource = PSPHOT_CNV_SRC_NONE;
-    options->targetSeeing = NAN;
+    options->targetSeeing = NULL;
 
     options->psfs        = psArrayAlloc(num);
@@ -47 +47 @@
     options->norm        = psVectorAlloc(num, PS_TYPE_F32);
     options->matchChi2   = psVectorAlloc(num, PS_TYPE_F32); // chi^2 for stamps when matching
-    options->weightings  = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2)
 
     psVectorInit(options->inputMask,   0);
@@ -53 +52 @@
     psVectorInit(options->norm,        NAN);
     psVectorInit(options->matchChi2,   NAN);
-    psVectorInit(options->weightings,  NAN);
 
     return options;

branches/meh_branches/ppstack_test/psphot/src/psphotStackPSF.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-pmPSF *psphotStackPSF(const pmConfig *config, int numCols, int numRows, const psArray *psfs, const psVector *inputMask)
-{
+// determine the 1st target PSF (either AUTO or defined by PSPHOT.STACK.TARGET.PSF.FWHM)
+bool psphotStackPSF(const pmConfig *config, psphotStackOptions *options)  {
+
     bool mdok = false;
-    pmPSF *psf = NULL;
+
+    int numCols = options->numCols;
+    int numRows = options->numRows;
+    psArray *psfs = options->psfs;
+    psVector *inputMask = options->inputMask;
 
     // Get the recipe values
@@ -17 +22 @@
         
     char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF
+
+    options->targetSeeing = psVectorAllocEmpty(4, PS_TYPE_F32);
 
     if (autoPSF) {
@@ -38 +45 @@
 
         // Solve for the target PSF
-        psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel, psfOrder, psfOrder, maskVal);
-        if (!psf) {
+        options->psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel, psfOrder, psfOrder, maskVal);
+        if (!options->psf) {
             psError(PSPHOT_ERR_PSF, false, "Unable to determine output PSF.");
-            return NULL;
+            return false;
         }
 
-    } else {
-
-        // externally-defined PSF
-        // XXX need to test for compatibility of target with inputs
-
-        float targetFWHM = psMetadataLookupF32 (&mdok, psphotRecipe, "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];
-        }
-
-        // measured scale factors (fwhm = Sxx * 2.35 * scaleFactor / sqrt(2.0))
-        // GAUSS  : 1.000
-        // PGAUSS : 1.006
-        // QGAUSS : 1.151
-        // RGAUSS : 0.883
-        // PS1_V1 : 1.134
-        
-        float scaleFactor = NAN;
-        if (!strcmp(psfModel, "PS_MODEL_GAUSS")) {
-            scaleFactor = 1.000;
-        }
-        if (!strcmp(psfModel, "PS_MODEL_PGAUSS")) {
-            scaleFactor = 1.0006;
-        }
-        if (!strcmp(psfModel, "PS_MODEL_QGAUSS")) {
-            scaleFactor = 1.151;
-        }
-        if (!strcmp(psfModel, "PS_MODEL_RGAUSS")) {
-            scaleFactor = 0.883;
-        }
-        if (!strcmp(psfModel, "PS_MODEL_PS1_V1")) {
-            scaleFactor = 1.134;
-        }
-        psAssert (isfinite(scaleFactor), "invalid model for PSF"); 
-
-        float Sxx = sqrt(2.0)*targetFWHM / 2.35 / scaleFactor;
-
-        // XXX probably should make the model type (and par 7) optional from recipe
-        // psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0);
-        psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 0.2);
-        if (!psf) {
-            psError(PSPHOT_ERR_PSF, false, "Unable to build dummy PSF.");
-            return NULL;
-        }
+        psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf);
+        float targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * options->numCols, 0.5 * options->numRows); // FWHM for target
+        psVectorAppend(options->targetSeeing, targetSeeing);
+        psLogMsg("psphotStack", PS_LOG_INFO, "Target seeing FWHM (auto-scaled): %f\n", targetSeeing);
+        return true;
     }
 
-    return psf;
+    // externally-defined PSF
+    // XXX need to test for compatibility of target with inputs
+
+    // is a single target FWHM specified, or a set of values?  set up the vector options->targetSeeing and the local 1st value
+    float targetSeeing = psMetadataLookupF32 (&mdok, psphotRecipe, "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");
+        for (int i = 0; i < fwhmValues->n; i++) {
+            psVectorAppend(options->targetSeeing, fwhmValues->data.F32[i]);
+        }           
+        targetSeeing = fwhmValues->data.F32[0];
+    } else {
+        psVectorAppend(options->targetSeeing, targetSeeing);
+    }
+
+    // measured scale factors (fwhm = Sxx * 2.35 * scaleFactor / sqrt(2.0))
+    // GAUSS  : 1.000
+    // PGAUSS : 1.006
+    // QGAUSS : 1.151
+    // RGAUSS : 0.883
+    // PS1_V1 : 1.134
+        
+    float scaleFactor = NAN;
+    if (!strcmp(psfModel, "PS_MODEL_GAUSS")) {
+        scaleFactor = 1.000;
+    }
+    if (!strcmp(psfModel, "PS_MODEL_PGAUSS")) {
+        scaleFactor = 1.0006;
+    }
+    if (!strcmp(psfModel, "PS_MODEL_QGAUSS")) {
+        scaleFactor = 1.151;
+    }
+    if (!strcmp(psfModel, "PS_MODEL_RGAUSS")) {
+        scaleFactor = 0.883;
+    }
+    if (!strcmp(psfModel, "PS_MODEL_PS1_V1")) {
+        scaleFactor = 1.134;
+    }
+    psAssert (isfinite(scaleFactor), "invalid model for PSF"); 
+
+    float Sxx = sqrt(2.0)*targetSeeing / 2.35 / scaleFactor;
+
+    // XXX probably should make the model type (and par 7) optional from recipe
+    // psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0);
+    options->psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 0.2);
+    if (!options->psf) {
+        psError(PSPHOT_ERR_PSF, false, "Unable to build dummy PSF.");
+        return false;
+    }
+
+    psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf);
+    psLogMsg("psphotStack", PS_LOG_INFO, "Target seeing FWHM (1 of %ld): %f\n", options->targetSeeing->n, targetSeeing);
+    return true;
 }

branches/meh_branches/ppstack_test/psphot/src/psphotStackReadout.c

@@ -44 +44 @@
 bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
 
+    psArray *objects = NULL; // used below after 'pass1finish' label
+
+    // measure the total elapsed time in psphotReadout.  dtime is the elapsed time used jointly
+    // by the multiple threads, not the total time used by all threads.
     psTimerStart ("psphotReadout");
 
@@ -56 +60 @@
     // optional break-point for processing
     char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
-    PS_ASSERT_PTR_NON_NULL (breakPt, false);
-
-    // XXX or do I set OUT to be a pmFPAfile pointing at the input of interest?
+    psAssert (breakPt, "configuration error: set BREAK_POINT");
+
+    // we have 3 relevant files: RAW (unconvolved), CNV (convolved stack), OUT (psf-matched stack)
+    // select which image (RAW or CNV) is used for analysis (RAW always used for detection)
     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 
+    char *STACK_DET = STACK_RAW;
 
     // set the photcode for each image
@@ -71 +74 @@
     }
 
-    // Generate the mask and weight images
-    // XXX this should be done before we perform the convolutions
+    // Generate the mask and weight images (if not supplied) and set mask bits
     if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
         return psphotReadoutCleanup (config, view, STACK_SRC);
@@ -81 +83 @@
 
     // generate a background model (median, smoothed image)
-    // 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_DET)) {
         return psphotReadoutCleanup (config, view, STACK_SRC);
@@ -89 +89 @@
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
+    if (!psphotSubtractBackground (config, view, STACK_SRC)) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
     if (!strcasecmp (breakPt, "BACKMDL")) {
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
+    // also make the chisq detection image
     if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
@@ -109 +113 @@
     }
 
-    // copy the detections from DET to SRC 
+    // If DET and SRC are different images, copy the detections from DET to SRC.  This 'copy'
+    // is just a copy of the container pointer; the sources on both DET and SRC are the same
+    // memory objects
     if (strcmp(STACK_SRC, STACK_DET)) {
         if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
@@ -122 +128 @@
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
-
-    if (!strcasecmp (breakPt, "TEST1")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
+    if (!strcasecmp (breakPt, "PEAKS")) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+    // psphotDumpTest (config, view, STACK_SRC);
     psMemDump("sourcestats");
-
-    // generate the objects (object unify the sources from the different images)
-    // 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);
-
-    psMemDump("matchsources");
-
-    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_SRC, false)) { // pass 1
-        psFree(objects);
-        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    psMemDump("sourcestats");
-
-    // find blended neighbors of very saturated stars (detections->newSources)
-    // if (!psphotDeblendSatstars (config, view)) {
-    //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
-    //     return psphotReadoutCleanup (config, view, STACK_SRC);
-    // }
-
-    // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
-    // if (!psphotBasicDeblend (config, view)) {
-    //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
-    //     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_SRC)) {
-        psFree(objects);
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
         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_SRC)) { // pass 1
-        psFree(objects);
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "MOMENTS")) {
-        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
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    // use bright stellar objects to measure PSF
     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_SRC);
     }
     if (!strcasecmp (breakPt, "PSFMODEL")) {
-        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_SRC);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
 
     // merge the newly selected sources into the existing list
@@ -200 +164 @@
     psphotMergeSources (config, view, STACK_SRC);
 
+    // Construct an initial model for each object, set the radius to fitRadius, set circular
+    // fit mask.  NOTE: only applied to sources without guess models
+    psphotGuessModels (config, view, STACK_SRC);
+
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
-    psphotFitSourcesLinearStack (config, objects, FALSE);
+    psphotFitSourcesLinear (config, view, STACK_SRC, false);
     psphotStackVisualFilerule(config, view, STACK_SRC);
 
+    // measure the radial profiles to the sky
+    psphotRadialProfileWings (config, view, STACK_SRC);
+
+    // re-measure the kron mags with models subtracted.  this pass starts with a circular
+    // window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
+    // but iterates to an appropriately larger size
+    psphotKronIterate(config, view, STACK_SRC);
+
     // identify CRs and extended sources
-    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_SRC)) {
-        psFree (objects);
-        psLogMsg ("psphot", 3, "failed on psphotApResid");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
+    psphotSourceSize (config, view, STACK_SRC, true);
+
+    // non-linear PSF and EXT fit to brighter sources
+    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
+    psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)
+
+    // replace all sources (do NOT ignore subtraction state)
+    psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
+
+    // if we only do one pass, skip to extended source analysis
+    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
+
+    // linear fit to include all sources (subtract again)
+    // NOTE : apply to ALL sources (extended + psf)
+    // NOTE 2 : this function subtracts the models from the given filerule (SRC), not DET
+    psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources)
+
+    // NOTE: possibly re-measure background model here with objects subtracted / or masked
+
+    // NOTE: this block performs the 2nd pass low-significance PSF detection stage
+    { 
+        // if DET and SRC are different images, generate children sources for all sources in
+        // the SRC image.  This operation replaces the existing DETECTION container on DET
+        // which is currently a view to the one on SRC).  children sources go to
+        // det->allSources
+        if (strcmp(STACK_SRC, STACK_DET)) {
+            psphotSourceChildren (config, view, STACK_DET, STACK_SRC); 
+
+            //  subtract all sources from DET (this will subtract using the psf model for SRC, which
+            //  will somewhat oversubtract the sources -- this is OK
+            psphotRemoveAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
+        }
+
+        // add noise for subtracted objects
+        psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
+
+        // find fainter sources
+        // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
+        psphotFindDetections (config, view, STACK_DET, false); // pass 2 (detections->peaks, detections->footprints)
+
+        // remove noise for subtracted objects (ie, return to normal noise level)
+        // NOTE: this needs to operate only on the OLD sources
+        psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
+
+        // if DET and SRC are different images, copy the detections from DET to SRC 
+        // (this operation just ensures the metadata container has a view on SRC as well
+        if (strcmp(STACK_SRC, STACK_DET)) {
+            // replace all sources in DET
+            psphotReplaceAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
+
+            // copy the newly detected peaks from DET to SRC so SourceStats below can operate on them
+            if (!psphotCopyPeaks (config, view, STACK_SRC, STACK_DET)) {
+                psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
+                return psphotReadoutCleanup (config, view, STACK_SRC);
+            }
+        }
+
+        // define new sources based on only the new peaks & measure moments
+        // NOTE: new sources are saved on detections->newSources
+        psphotSourceStats (config, view, STACK_SRC, false); // pass 2 (detections->newSources)
+
+        // set source type
+        // NOTE: apply only to detections->newSources
+        if (!psphotRoughClass (config, view, STACK_SRC)) { // pass 2 (detections->newSources)
+            psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
+            return psphotReadoutCleanup (config, view, STACK_SRC);
+        }
+
+        // replace all sources so fit below applies to all at once
+        // NOTE: apply only to OLD sources (which have been subtracted)
+        psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 2
+
+        // merge the newly selected sources into the existing list
+        // NOTE: merge OLD and NEW
+        // XXX check on free of sources...
+        psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
+
+        // Construct an initial model for each object, set the radius to fitRadius, set circular
+        // fit mask.  NOTE: only applied to sources without guess models
+        psphotGuessModels (config, view, STACK_SRC);
+    }
+
+pass1finish:
+
+    // generate the objects (objects unify the sources from the different images) NOTE: could
+    // this just match the detections for the chisq image, and not bother measuring the source
+    // stats in that case...?
+    objects = psphotMatchSources (config, view, STACK_SRC);
+    psMemDump("matchsources");
+
+    // Construct an initial model for each object, set the radius to fitRadius, set circular
+    // fit mask.  NOTE: only applied to sources without guess models
+    psphotGuessModels (config, view, STACK_SRC);
+
+    psphotStackObjectsUnifyPosition (objects);
+
+    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
+
+    // NOTE: apply to ALL sources
+    psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources)
+
+    // measure the radial profiles to the sky (only measures new objects)
+    psphotRadialProfileWings (config, view, STACK_SRC);
+
+    // re-measure the kron mags with models subtracted
+    // psphotKronMasked(config, view, STACK_SRC);
+    psphotKronIterate(config, view, STACK_SRC);
+
+    // measure source size for the remaining sources
+    // NOTE: applies only to NEW (unmeasured) sources
+    psphotSourceSize (config, view, STACK_SRC, false); // pass 2 (detections->allSources)
 
     psMemDump("psfstats");
 
-    psphotStackObjectsUnifyPosition (objects);
-
     // measure elliptical apertures, petrosians (objects sorted by S/N)
-    psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
+    // psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotExtendedSourceAnalysis (config, 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_SRC); // pass 1 (detections->allSources)
-
-    // calculate source magnitudes
-    psphotMagnitudes(config, view, STACK_SRC);
 
     // create source children for the OUT filerule (for radial aperture photometry) 
@@ -238 +310 @@
     }
 
+    // measure circular, radial apertures (objects sorted by S/N)
+    // this forces photometry on the undetected sources from other images
+    psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched?
     psMemDump("extmeas");
 
-    bool smoothAgain = true;
-    for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) {
+    int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
+
+    for (int entry = 1; entry < nRadialEntries; entry++) {
+        // NOTE: entry 0 is the unmatched image set
 
         // re-measure the PSF for the smoothed image (using entries in 'allSources')
@@ -253 +330 @@
 
         // measure circular, radial apertures (objects sorted by S/N)
-        psphotRadialAperturesByObject (config, objectsRadial, view, STACK_OUT, nMatchedPSF); 
+        // entry 0 == unmatched? pass entry + 1?
+        psphotRadialApertures (config, view, STACK_OUT, entry); 
 
         // replace the flux in the image so it is returned to its original state
-        psphotReplaceAllSources (config, view, STACK_OUT);
+        psphotReplaceAllSources (config, view, STACK_OUT, false);
 
         // smooth to the next FWHM, or set 'smoothAgain' to false if no more 
-        psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF);
+        psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
         psMemDump("matched");
     }
 
-    if (0 && !psphotEfficiency(config, view, STACK_OUT)) {
+    // measure aperture photometry corrections
+    if (!psphotApResid (config, view, STACK_SRC)) {
+        psFree (objects);
+        psLogMsg ("psphot", 3, "failed on psphotApResid");
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    // calculate source magnitudes
+    psphotMagnitudes(config, view, STACK_SRC);
+
+    // XXX NOTE: this function wants to have the PSF of the image, but we (so far) only measure the 
+    // PSF of the SRC image.  can we fake it by generating the PSF for DET as well (up above)? 
+    if (false && !psphotEfficiency(config, view, STACK_DET)) {
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
         psErrorClear();

branches/meh_branches/ppstack_test/psphot/src/psphotSubtractBackground.c

@@ -103 +103 @@
         psphotSaveImage (NULL, image, name);
     }
-    psLogMsg ("psphot", PS_LOG_INFO, "subtracted background model: %f sec\n", psTimerMark ("psphot.background"));
+    psLogMsg ("psphot", PS_LOG_WARN, "subtracted background model: %f sec\n", psTimerMark ("psphot.background"));
 
     // the pmReadout selected in this function are all view on entries in config->files

branches/meh_branches/ppstack_test/psphot/src/psphotTest.c

@@ -85 +85 @@
 
     // wait for the threads to finish and manage results
-    if (!psThreadPoolWait (true)) {
+    if (!psThreadPoolWait (true, true)) {
         fprintf (stderr, "failure to run FillImage (2)");
         exit (1);

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest