Changeset 27657

Timestamp:

Apr 11, 2010, 5:08:29 PM (16 years ago)

Author:

eugene

Message:

updates to support psphotStack

Location:

trunk

Files:

• ippconfig (modified) (1 prop)
• ippconfig/recipes/filerules-simple.mdc (modified) (1 diff)
• ppSim (modified) (1 prop)
• ppSim/src/ppSimInsertGalaxies.c (modified) (2 diffs)
• ppSim/src/ppSimInsertStars.c (modified) (1 diff)
• ppSim/src/ppSimPhotomReadout.c (modified) (4 diffs)
• ppSim/src/ppSimPhotomReadoutFake.c (modified) (1 diff)
• ppSim/src/ppSimPhotomReadoutForce.c (modified) (2 diffs)
• ppSim/src/ppSimUtils.c (modified) (1 diff)
• psModules (modified) (1 prop)
• psModules/src/camera/pmFPAfile.c (modified) (3 diffs)
• psModules/src/camera/pmFPAfile.h (modified) (1 diff)
• psModules/src/camera/pmFPAfileDefine.c (modified) (2 diffs)
• psModules/src/camera/pmFPAfileDefine.h (modified) (1 diff)
• psModules/src/objects/Makefile.am (modified) (2 diffs)
• psModules/src/objects/pmPeaks.h (modified) (1 diff)
• psModules/src/objects/pmPhotObj.c (modified) (2 diffs)
• psModules/src/objects/pmPhotObj.h (modified) (1 diff)
• psModules/src/objects/pmSource.c (modified) (2 diffs)
• psModules/src/objects/pmSource.h (modified) (3 diffs)
• psModules/src/psmodules.h (modified) (1 diff)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/Makefile.am (modified) (3 diffs)
• psphot/src/psphot.h (modified) (1 diff)
• psphot/src/psphotApResid.c (modified) (1 diff)
• psphot/src/psphotChoosePSF.c (modified) (1 diff)
• psphot/src/psphotEfficiency.c (modified) (1 diff)
• psphot/src/psphotFitSourcesLinearStack.c (modified) (11 diffs)
• psphot/src/psphotGuessModels.c (modified) (1 diff)
• psphot/src/psphotImageLoop.c (modified) (1 diff)
• psphot/src/psphotImageQuality.c (modified) (1 diff)
• psphot/src/psphotMagnitudes.c (modified) (1 diff)
• psphot/src/psphotMergeSources.c (modified) (3 diffs)
• psphot/src/psphotModelBackground.c (modified) (1 diff)
• psphot/src/psphotReadout.c (modified) (1 diff)
• psphot/src/psphotReadoutCleanup.c (modified) (2 diffs)
• psphot/src/psphotRoughClass.c (modified) (1 diff)
• psphot/src/psphotSkyReplace.c (modified) (1 diff)
• psphot/src/psphotSourceMatch.c (modified) (1 diff)
• psphot/src/psphotSourceMerge.c (deleted)
• psphot/src/psphotSourceSize.c (modified) (1 diff)
• psphot/src/psphotSourceStats.c (modified) (10 diffs)
• psphot/src/psphotStackArguments.c (modified) (4 diffs)
• psphot/src/psphotStackChisqImage.c (modified) (6 diffs)
• psphot/src/psphotStackImageLoop.c (modified) (4 diffs)
• psphot/src/psphotStackParseCamera.c (modified) (4 diffs)
• psphot/src/psphotStackReadout.c (modified) (3 diffs)
• psphot/src/psphotSubtractBackground.c (modified) (1 diff)

trunk/ippconfig/recipes/filerules-simple.mdc

@@ -167 +167 @@
 PSPHOT.PSF.RAW.SAVE     OUTPUT {OUTPUT}.psf          PSF       NONE       FPA        TRUE      NONE
 PSPHOT.PSF.SKY.SAVE     OUTPUT {OUTPUT}.psf          PSF       NONE       FPA        TRUE      NONE
-                                                     
+
 # outputs for psphotStack:
 PSPHOT.CHISQ.IMAGE      OUTPUT {OUTPUT}.chisq.im.fits IMAGE    NONE       FPA        TRUE      SIMPLE

trunk/ppSim/src/ppSimInsertGalaxies.c

@@ -74 +74 @@
     int dY = PM_CELL_TO_CHIP (0.0, y0Cell, yParityCell, binning);
 
-    // psMetadataLookupPtr (readout->analysis, "PSPHOT.SOURCES", 0);
-
-    psArray *sources = psArrayAllocEmpty (galaxies->n);
-
+    pmDetections *detections = psMetadataLookupPtr (&mdok, readout->analysis, "PSPHOT.DETECTIONS");
+    if (!detections) {
+        detections = pmDetectionsAlloc();
+        detections->allSources = psArrayAllocEmpty (galaxies->n);
+        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY | PS_META_REPLACE, "psphot detections", detections);
+    } else {
+        psMemIncrRefCounter (detections);
+    }
+    psArray *sources = sources = detections->allSources;
 
     // add sources to the readout image & weight
@@ -181 +186 @@
     }
 
-    // NOTE: readout must be part of the pmFPAfile named "PPSIM.OUTPUT"
-    // psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE, "psphot sources", sources);
+    // XXX many leaks in here, i think 
+    psFree (detections);
 
-    // XXX many leaks in here, i think 
     return true;
 }

trunk/ppSim/src/ppSimInsertStars.c

@@ -145 +145 @@
     fclose (outfile);
 
-    // NOTE: the pmFPAfile "PPSIM.OUTPUT" points at these sources
-    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "fake sources", sources);
-    psFree(sources);
+    pmDetections *detections = pmDetectionsAlloc();
+    detections->allSources = sources;
+
+    // save detections on the readout->analysis
+    if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "fake sources", detections)) {
+        psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout");
+        return false;
+    }
+    psFree(detections);
 
     // XXX many leaks in here, i think

trunk/ppSim/src/ppSimPhotomReadout.c

@@ -6 +6 @@
     PS_ASSERT_PTR_NON_NULL (readout, NULL);
 
-    psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->allSources;
     return sources;
 }
@@ -151 +154 @@
 
     // create the exported-metadata and free local data
-    // XXX this places the sources on readout->analysis as PSPHOT.SOURCES.  modify?
-    // (or don't supply the sources, and do this with a different function)
     psphotReadoutCleanup(config, readout, recipe, NULL, psf, NULL);
 
@@ -167 +168 @@
     psAssert (forceReadout, "no forceReadout?");
     pmChipSetDataStatus (forceChip, true);
-    psMetadataAddArray (forceReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "forced photometry ", forceSources);
-    psFree (forceSources);
+
+    pmDetections *detectionsForce = pmDetectionsAlloc();
+    detectionForce->allSources = forceSources;
+
+    // save detections on the readout->analysis
+    if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "force sources", detectionsForce)) {
+        psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout");
+        return false;
+    }
+    psFree(detectionsForce);
 
     pmCell    *fakeCell    = pmFPAfileThisCell (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeCell, "no cell?");
@@ -179 +188 @@
     psAssert (fakeReadout, "no fakeReadout?");
     pmChipSetDataStatus (fakeChip, true);
-    psMetadataAddArray (fakeReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "fake photometry ", fakeSources);
+
+    pmDetections *detectionsFake = pmDetectionsAlloc();
+    detectionFake->allSources = fakeSources;
+
+    // save detections on the readout->analysis
+    if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "fake sources", detectionsFake)) {
+        psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout");
+        return false;
+    }
+    psFree(detectionsFake);
 
     return true;

trunk/ppSim/src/ppSimPhotomReadoutFake.c

@@ -111 +111 @@
     psAssert (fakeReadout, "no fakeReadout?");
     pmChipSetDataStatus (fakeChip, true);
-    psMetadataAddArray (fakeReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "fake photometry ", fakeSources);
+
+    pmDetections *detections = pmDetectionsAlloc();
+    detection->allSources = fakeSources;
+
+    psMetadataAddArray (fakeReadout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE, "fake photometry ", detections);
+    psFree(detections);
 
     return true;

trunk/ppSim/src/ppSimPhotomReadoutForce.c

@@ -98 +98 @@
 
     // create the exported-metadata and free local data
-    // XXX this places the sources on readout->analysis as PSPHOT.SOURCES.  modify?
-    // (or don't supply the sources, and do this with a different function)
     psphotReadoutCleanup(config, readout, recipe, NULL, psf, NULL);
 
@@ -111 +109 @@
     psAssert (forceReadout, "no forceReadout?");
     pmChipSetDataStatus (forceChip, true);
-    psMetadataAddArray (forceReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "force photometry ", forceSources);
+
+    pmDetections *detections = pmDetectionsAlloc();
+    detection->allSources = forceSources;
+    psMetadataAddArray (forceReadout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE, "force photometry ", detections);
+    psFree(detections);
 
     return true;

trunk/ppSim/src/ppSimUtils.c

@@ -270 +270 @@
 psArray *ppSimSelectSources (pmConfig *config, const pmFPAview *view, const char *filename) {
 
+    bool status;
+
     pmReadout *readout = pmFPAfileThisReadout (config->files, view, filename);
     PS_ASSERT_PTR_NON_NULL (readout, NULL);
 
-    psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->allSources;
     return sources;
 }

trunk/psModules/src/camera/pmFPAfile.c

@@ -111 +111 @@
     file->save = false;
 
-    file->index = fileNum++;
+    file->fileIndex = fileNum++;
+    file->fileID = 0;
 
     file->imageId = 0;
@@ -372 +373 @@
         // Number of the file in list
         psString num = NULL;            // Number to use
-        psStringAppend(&num, "%d", file->index);
+        psStringAppend(&num, "%d", file->fileIndex);
         psStringSubstitute(&newRule, num, "{FILE.INDEX}");
+        psFree(num);
+    }
+
+    if (strstr(newRule, "{FILE.ID}")) {
+        // Number of the file in list
+        psString num = NULL;            // Number to use
+        psStringAppend(&num, "%03d", file->fileID);
+        psStringSubstitute(&newRule, num, "{FILE.ID}");
         psFree(num);
     }
@@ -638 +647 @@
     psFree(iter);
 
-    psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find instance %d of file %s", num, name);
+    psLogMsg("psModules.camera", PS_LOG_MINUTIA, "Unable to find instance %d of file %s", num, name);
     return NULL;
 }

trunk/psModules/src/camera/pmFPAfile.h

@@ -111 +111 @@
     psString formatName;                // name of the camera format
 
-    int index;                          // Index of file
+    int fileIndex;                      // Index of file
+    int fileID;                         // internal sequence number
 
     psS64 imageId, sourceId;            // Image and source identifiers

trunk/psModules/src/camera/pmFPAfileDefine.c

@@ -1271 +1271 @@
     file->name = psStringCopy (name);
 
+    // free a previously existing readout
+    psFree(file->readout);
     file->readout = readout;
-    psMetadataAddPtr(files, PS_LIST_TAIL, name, PS_DATA_UNKNOWN, "", file);
+
+    // allow for multiple entries
+    // XXX handle replace vs multiple?
+    psMetadataAddPtr(files, PS_LIST_TAIL, name, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", file);
     psFree(file);
     // we free this copy of file, but 'files' still has a copy
@@ -1313 +1318 @@
                                 const char *name, // name of internal/external file
                                 const pmFPA *fpa, // use this fpa to generate
-                                const psImageBinning *binning) {
+                                const psImageBinning *binning,
+                                int index) {
   pmReadout *readout = NULL;
 
-  bool status = true;
-  pmFPAfile *file = psMetadataLookupPtr(&status, config->files, name);
+  pmFPAfile *file = pmFPAfileSelectSingle(config->files, name, index);
 
   // if the file does not exist, it is not being used as an I/O file: define an internal version
   if (file == NULL) {
-    readout = pmFPAfileDefineInternal (config->files, name, binning->nXruff, binning->nYruff, PS_TYPE_F32);
-    return readout;
+      // XXX currently, we do not guarantee that the defined file lands on entry 'index'
+      psAssert (binning, "internal files must be supplied a psImageBinning for the output images size"); 
+      readout = pmFPAfileDefineInternal (config->files, name, binning->nXruff, binning->nYruff, PS_TYPE_F32);
+      return readout;
   }
 

trunk/psModules/src/camera/pmFPAfileDefine.h

@@ -172 +172 @@
                                 const char *name, // name of internal/external file
                                 const pmFPA *fpa, // use this fpa to generate
-                                const psImageBinning *binning);
+                                const psImageBinning *binning,
+                                int index
+    );
 
 /// @}

trunk/psModules/src/objects/Makefile.am

@@ -20 +20 @@
         pmModelUtils.c \
         pmSource.c \
+        pmPhotObj.c \
         pmSourceMasks.c \
         pmSourceMoments.c \
@@ -80 +81 @@
         pmModelUtils.h \
         pmSource.h \
+        pmPhotObj.h \
         pmSourceMasks.h \
         pmSourceDiffStats.h \

trunk/psModules/src/objects/pmPeaks.h

@@ -63 +63 @@
     bool assigned;                      ///< is peak assigned to a source?
     pmPeakType type;                    ///< Description of peak.
-    pmFootprint *footprint;     ///< reference to containing footprint
+    pmFootprint *footprint;             ///< reference to containing footprint
 }
 pmPeak;

trunk/psModules/src/objects/pmPhotObj.c

@@ -17 +17 @@
 #include <pslib.h>
 #include "pmPhotObj.h"
+#include "pmSource.h"
 
 static void pmPhotObjFree (pmPhotObj *tmp)
@@ -38 +39 @@
 }
 
+bool pmPhotObjAddSource(pmPhotObj *object, pmSource *source) {
+
+    psAssert (source, "programming error: NULL source");
+    if (!source->peak) {
+        psError(PS_ERR_UNKNOWN, true, "source missing peak");
+        return false; 
+    }
+    if (!finite(source->peak->xf)) {
+        psError(PS_ERR_UNKNOWN, true, "NAN peak coordinate");
+        return false; 
+    }
+    if (!finite(source->peak->yf)) {
+        psError(PS_ERR_UNKNOWN, true, "NAN peak coordinate");
+        return false; 
+    }
+
+    // XXX we should probably use the fitted position if it exists
+    if (!object->sources) {
+        object->sources = psArrayAllocEmpty(1);
+        object->x = source->peak->xf;
+        object->y = source->peak->yf;
+    }
+    psArrayAdd (object->sources, 1, source);
+    return true;
+}

trunk/psModules/src/objects/pmPhotObj.h

@@ -34 +34 @@
  */
 typedef struct {
-  int seq;                            ///< ID for output (generated on write OR set on read)
-  psArray *sources;
+    int id;                            ///< ID for output (generated on write OR set on read)
+    psArray *sources;
+    int flags;
+    float x;
+    float y;
 } pmPhotObj;
+
+bool pmPhotObjAddSource(pmPhotObj *object, pmSource *source);
+pmPhotObj *pmPhotObjAlloc(void);
 
 /// @}
 # endif /* PM_PHOT_OBJ_H */
+

trunk/psModules/src/objects/pmSource.c

@@ -180 +180 @@
     source->type = in->type;
     source->mode = in->mode;
+    source->imageID = in->imageID;
 
     return(source);
@@ -1058 +1059 @@
     psF32 fA = (A->peak == NULL) ? 0 : A->peak->y;
     psF32 fB = (B->peak == NULL) ? 0 : B->peak->y;
+
+    psF32 diff = fA - fB;
+    if (diff > FLT_EPSILON) return (+1);
+    if (diff < FLT_EPSILON) return (-1);
+    return (0);
+}
+
+// sort by X (ascending)
+int pmSourceSortByX (const void **a, const void **b)
+{
+    pmSource *A = *(pmSource **)a;
+    pmSource *B = *(pmSource **)b;
+
+    psF32 fA = (A->peak == NULL) ? 0 : A->peak->x;
+    psF32 fB = (B->peak == NULL) ? 0 : B->peak->x;
 
     psF32 diff = fA - fB;

trunk/psModules/src/objects/pmSource.h

@@ -43 +43 @@
     PM_SOURCE_TMPF_SIZE_MEASURED     = 0x0004,
     PM_SOURCE_TMPF_SIZE_CR_CANDIDATE = 0x0008,
+    PM_SOURCE_TMPF_MOMENTS_MEASURED  = 0x0010,
 } pmSourceTmpF;
 
@@ -90 +91 @@
     pmSourceExtendedPars *extpars;      ///< extended source parameters
     pmSourceDiffStats *diffStats;       ///< extra parameters for difference detections
+    int imageID;
 };
 
@@ -245 +247 @@
 int  pmSourceSortBySN (const void **a, const void **b);
 int  pmSourceSortByY (const void **a, const void **b);
+int  pmSourceSortByX (const void **a, const void **b);
 int  pmSourceSortBySeq (const void **a, const void **b);
 

trunk/psModules/src/psmodules.h

@@ -124 +124 @@
 #include <pmSourceMasks.h>
 #include <pmSource.h>
+#include <pmPhotObj.h>
 #include <pmSourceUtils.h>
 #include <pmSourceIO.h>

trunk/psphot/src

@@ -21 +21 @@
 psphotForced
 psphotMakePSF
+psphotStack

@@ -21 +21 @@
 psphotForced
 psphotMakePSF
+psphotStack

trunk/psphot/src/Makefile.am

@@ -17 +17 @@
 # Force recompilation of psphotVersion.c, since it gets the version information
 psphotVersion.c: psphotVersionDefinitions.h
-psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE
-        -$(RM) psphotVersionDefinitions.h
-        $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h
-FORCE: ;
+# psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE
+#       -$(RM) psphotVersionDefinitions.h
+#       $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h
+# FORCE: ;
 
 libpsphot_la_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
 libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 
-bin_PROGRAMS = psphot psphotForced psphotMakePSF
+bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotStack
 # bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy 
 
@@ -39 +39 @@
 psphotMakePSF_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 psphotMakePSF_LDADD = libpsphot.la
+
+psphotStack_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+psphotStack_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
+psphotStack_LDADD = libpsphot.la
 
 # psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
@@ -80 +84 @@
         psphotMosaicChip.c         \
         psphotCleanup.c
+
+# a psphot-variant for stack photometry
+psphotStack_SOURCES = \
+        psphotStack.c                 \
+        psphotStackArguments.c        \
+        psphotStackParseCamera.c      \
+        psphotStackImageLoop.c        \
+        psphotStackReadout.c          \
+        psphotStackChisqImage.c       \
+        psphotFitSourcesLinearStack.c \
+        psphotSourceMatch.c           \
+        psphotCleanup.c
+
+
 
 # # psphot analysis of the detectability of specified positions

trunk/psphot/src/psphot.h

@@ -341 +341 @@
 bool psphotMakePSFReadout(pmConfig *config, const pmFPAview *view);
 
+/**** psphotStack prototypes ****/
+
+pmConfig *psphotStackArguments(int argc, char **argv);
+bool psphotStackParseCamera (pmConfig *config);
+bool psphotStackImageLoop (pmConfig *config);
+bool psphotStackReadout (pmConfig *config, const pmFPAview *view);
+bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view);
+bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration
+                                     const pmFPAview *view,
+                                     pmReadout **chiReadout,
+                                     char *filename, 
+                                     int index);
+
+bool psphotStackRemoveChisqFromInputs (pmConfig *config);
+bool pmFPAfileRemoveSingle(psMetadata *files, const char *name, int num);
+
+psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view);
+bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index);
+bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS);
+
+bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final);
+int pmPhotObjSortBySN (const void **a, const void **b);
+int pmPhotObjSortByX (const void **a, const void **b);
+
 #endif

trunk/psphot/src/psphotApResid.c

@@ -16 +16 @@
     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
         if (!psphotApResidReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for PSPHOT.INPUT entry %d", i);

trunk/psphot/src/psphotChoosePSF.c

@@ -13 +13 @@
     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
         if (!psphotChoosePSFReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for PSPHOT.INPUT entry %d", i);

trunk/psphot/src/psphotEfficiency.c

@@ -160 +160 @@
     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
         if (!psphotEfficiencyReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure detection efficiency for PSPHOT.INPUT entry %d", i);

trunk/psphot/src/psphotFitSourcesLinearStack.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+// XXX need array of covar factors for each image
+// XXX define the 'good' / 'bad' flags?
 
-// for now, let's store the detections on the readout->analysis for each readout
-bool psphotFitSourcesLinearStack (pmConfig *config, const pmFPAview *view, bool final)
-{
-    bool status = true;
+# define COVAR_FACTOR 1.0
 
-    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
-    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
-
-    // loop over the available readouts
-    for (int i = 0; i < num; i++) {
-        if (!psphotFitSourcesLinearReadoutStack (config, view, "PSPHOT.INPUT", i, final)) {
-            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for PSPHOT.INPUT entry %d", i);
-            return false;
-        }
-    }
-    return true;
-}
-
-bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final) {
+bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final) {
 
     bool status;
-    float x;
-    float y;
     float f;
-    // float r;
 
     // select the appropriate recipe information
     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     assert (recipe);
-
-    // 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?");
-
-    if (!sources->n) {
-        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit");
-        return true;
-    }
-
-    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
-    psAssert (sources, "missing psf?");
 
     psTimerStart ("psphot.linear");
@@ -65 +27 @@
     maskVal |= markVal;
 
-    // source analysis is done in spatial order
-    sources = psArraySort (sources, pmSourceSortByY);
+    // analysis is done in spatial order (to speed up overlap search)
+    // sort by first element in each source list
+    objects = psArraySort (objects, pmPhotObjSortByX);
 
     // storage array for fitSources
-    psArray *fitSources = psArrayAllocEmpty (sources->n);
+    psArray *fitSources = psArrayAllocEmpty (objects->n);
 
-    // option to limit analysis to a specific region
-    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
-    psRegion AnalysisRegion = psRegionFromString (region);
-    AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion);
-    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
+    bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
+    psAssert (status, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
 
-    bool CONSTANT_PHOTOMETRIC_WEIGHTS =
-        psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
-    if (!status) {
-        psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
-    }
-    int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
-    if (!status) {
-        SKY_FIT_ORDER = 0;
-    }
-    bool SKY_FIT_LINEAR = psMetadataLookupBool(&status, recipe, "SKY_FIT_LINEAR");
-    if (!status) {
-        SKY_FIT_LINEAR = false;
-    }
-
-    // XXX test: choose a larger-than expected radius:
-    float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
-    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor);
-
-    // XXX do not apply covarFactor for the moment...
-    // covarFactor = 1.0;
+    // XXX store a local static array of covar factors for each of the images (by image ID)
+    // float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+    // psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor);
 
     // select the sources which will be used for the fitting analysis
@@ -104 +47 @@
         if (!object->sources) continue;
 
-        // check an element of the group to see if we should use it
-        if (!object->flags & PM_PHOT_OBJ_BAD) continue;
+        // XXX check an element of the group to see if we should use it
+        // if (!object->flags & PM_PHOT_OBJ_BAD) continue;
 
         for (int j = 0; j < object->sources->n; j++) {
@@ -123 +66 @@
         }
     }
-    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), sources->n);
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), objects->n);
 
     if (fitSources->n == 0) {
@@ -143 +86 @@
 
         // diagonal elements of the sparse matrix (auto-cross-product)
-        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
+        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
         psSparseMatrixElement (sparse, i, i, f);
 
         // the formal error depends on the weighting scheme
         if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
-            float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor);
+            float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR);
             errors->data.F32[i] = 1.0 / sqrt(var);
         } else {
@@ -155 +98 @@
 
         // find the image x model value
-        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
+        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
         psSparseVectorElement (sparse, i, f);
 
@@ -162 +105 @@
             pmSource *SRCj = fitSources->data[j];
 
-            // XXX I need to know if this source is on the same image as SRCi --
-            if (!sameImge) { continue; }
+            // we only need to generate dot terms for source on the same image
+            if (SRCj->imageID != SRCi->imageID) { continue; }
 
             // skip over disjoint source images, break after last possible overlap
-            if (SRCi->pixels->row0 + SRCi->pixels->numRows < SRCj->pixels->row0) break;
-            if (SRCj->pixels->row0 + SRCj->pixels->numRows < SRCi->pixels->row0) continue;
-            if (SRCi->pixels->col0 + SRCi->pixels->numCols < SRCj->pixels->col0) continue;
-            if (SRCj->pixels->col0 + SRCj->pixels->numCols < SRCi->pixels->col0) continue;
+            if (SRCj->pixels->row0 + SRCj->pixels->numRows < SRCi->pixels->row0) continue;  // source(i) is above source(j)
+            if (SRCi->pixels->row0 + SRCi->pixels->numRows < SRCj->pixels->row0) continue;  // source(i) is below source(j)
+            if (SRCj->pixels->col0 + SRCj->pixels->numCols < SRCi->pixels->col0) continue;  // source(i) is right of source(j)
+            if (SRCi->pixels->col0 + SRCi->pixels->numCols < SRCj->pixels->col0) break;     // source(i) is left of source(j) [no other source(j) can overlap source(i)]
 
             // got an overlap; calculate cross-product and add to output array
-            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
+            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
             psSparseMatrixElement (sparse, j, i, f);
         }
@@ -190 +133 @@
     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
-
     // adjust I0 for fitSources and subtract
     for (int i = 0; i < fitSources->n; i++) {
@@ -208 +145 @@
         model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
         model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
-        // XXX is the value of 'errors' modified by the sky fit?
 
         // subtract object
@@ -221 +157 @@
         if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue;
         pmModel *model = pmSourceGetModel (NULL, source);
-        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, covarFactor);
+        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, COVAR_FACTOR);
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
@@ -229 +165 @@
     psFree (fitSources);
     psFree (norm);
-    psFree (skyfit);
     psFree (errors);
-    psFree (border);
 
     psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
-
-    psphotVisualShowResidualImage (readout);
-    psphotVisualPlotChisq (sources);
-    // psphotVisualShowFlags (sources);
-
-    // We have to place this visualization here because the models are not realized until
-    // psphotGuessModels or fitted until psphotFitSourcesLinear.
-    psphotVisualShowPSFStars (recipe, psf, sources);
 
     return true;
 }
 
-// XXX do we need this?
-// XXX disallow the simultaneous sky fit and remove this code...
+// sort by X (ascending)
+int pmPhotObjSortByX (const void **a, const void **b)
+{
+    pmPhotObj *objA = *(pmPhotObj **)a;
+    pmPhotObj *objB = *(pmPhotObj **)b;
 
-// Calculate the weight terms for the sky fit component of the matrix.  This function operates
-// on the pixels which correspond to all of the sources of interest.  These elements fill in
-// the border matrix components in the sparse matrix equation.
-static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal) {
+    psF32 fA = objA->x;
+    psF32 fB = objB->x;
 
-    // generate the image-wide weight terms
-    // turn on MARK for all image pixels
-    psRegion fullArray = psRegionSet (0, 0, 0, 0);
-    fullArray = psRegionForImage (readout->mask, fullArray);
-    psImageMaskRegion (readout->mask, fullArray, "OR", markVal);
-
-    // turn off MARK for all object pixels
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-        pmModel *model = pmSourceGetModel (NULL, source);
-        if (model == NULL) continue;
-        float x = model->params->data.F32[PM_PAR_XPOS];
-        float y = model->params->data.F32[PM_PAR_YPOS];
-        psImageMaskCircle (source->maskView, x, y, model->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal));
-    }
-
-    // accumulate the image statistics from the masked regions
-    psF32 **image  = readout->image->data.F32;
-    psF32 **variance = readout->variance->data.F32;
-    psImageMaskType  **mask   = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA;
-
-    double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy;
-    w = x = y = x2 = xy = y2 = fo = fx = fy = 0;
-
-    int col0 = readout->image->col0;
-    int row0 = readout->image->row0;
-
-    for (int j = 0; j < readout->image->numRows; j++) {
-        for (int i = 0; i < readout->image->numCols; i++) {
-            if (mask[j][i]) continue;
-            if (constant_weights) {
-                wt = 1.0;
-            } else {
-                wt = variance[j][i];
-            }
-            f = image[j][i];
-            w   += 1/wt;
-            fo  += f/wt;
-            if (SKY_FIT_ORDER == 0) continue;
-
-            xc  = i + col0;
-            yc  = j + row0;
-            x  +=    xc/wt;
-            y  +=    yc/wt;
-            x2 += xc*xc/wt;
-            xy += xc*yc/wt;
-            y2 += yc*yc/wt;
-            fx +=  f*xc/wt;
-            fy +=  f*yc/wt;
-        }
-    }
-
-    // turn off MARK for all image pixels
-    psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_IMAGE_MASK(markVal));
-
-    // set the Border T elements
-    psSparseBorderElementG (border, 0, fo);
-    psSparseBorderElementT (border, 0, 0, w);
-    if (SKY_FIT_ORDER > 0) {
-        psSparseBorderElementG (border, 0, fx);
-        psSparseBorderElementG (border, 0, fy);
-        psSparseBorderElementT (border, 1, 0, x);
-        psSparseBorderElementT (border, 2, 0, y);
-        psSparseBorderElementT (border, 0, 1, x);
-        psSparseBorderElementT (border, 1, 1, x2);
-        psSparseBorderElementT (border, 2, 1, xy);
-        psSparseBorderElementT (border, 0, 2, y);
-        psSparseBorderElementT (border, 1, 2, xy);
-        psSparseBorderElementT (border, 2, 2, y2);
-    }
-
-    return true;
+    psF32 diff = fA - fB;
+    if (diff > FLT_EPSILON) return (+1);
+    if (diff < FLT_EPSILON) return (-1);
+    return (0);
 }

trunk/psphot/src/psphotGuessModels.c

@@ -15 +15 @@
     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
         if (!psphotGuessModelsReadout (config, view, "PSPHOT.INPUT", i)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on to guess models for PSPHOT.INPUT entry %d", i);

trunk/psphot/src/psphotImageLoop.c

@@ -91 +91 @@
             }
 
-            status = true;
-            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");
-            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV");
-            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND");
-            if (!status) {
-                psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files");
-                psFree (view);
-                return false;
-            }
+            // drop all versions of the internal files
+            status = true;
+            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");
+            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV");
+            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND");
+            if (!status) {
+                psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files");
+                psFree (view);
+                return false;
+            }
         }
-
         // save output which is saved at the chip level
         if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot.");

trunk/psphot/src/psphotImageQuality.c

@@ -16 +16 @@
     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
         if (!psphotImageQualityReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on to measure image quality for PSPHOT.INPUT entry %d", i);

trunk/psphot/src/psphotMagnitudes.c

@@ -12 +12 @@
     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
 
+    // 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

trunk/psphot/src/psphotMergeSources.c

@@ -76 +76 @@
     pmDetections *extCMF = NULL;
     psArray *extSourcesTXT = NULL;
-
-    // find the currently selected readout
-    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest
+    int index = 0;
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest
     psAssert (file, "missing file?");
 
@@ -118 +119 @@
                 psFree (source->modelPSF);
                 source->modelPSF = NULL;
+                source->imageID = index;
 
                 psArrayAdd (detections->newSources, 100, source);
@@ -140 +142 @@
                 psFree (source->modelPSF);
                 source->modelPSF = NULL;
+                source->imageID = index;
 
                 psArrayAdd (detections->newSources, 100, source);

trunk/psphot/src/psphotModelBackground.c

@@ -371 +371 @@
 
     psImageBinning *binning = psphotBackgroundBinning(readout->image, config); // Image binning parameters
-    pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning);
-    pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning);
+    pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning, index);
+    pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning, index);
 
     if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config)) {

trunk/psphot/src/psphotReadout.c

@@ -34 +34 @@
 
     // Generate the mask and weight images, including the user-defined analysis region of interest
-    psphotSetMaskAndVariance (config, view);
+    if (!psphotSetMaskAndVariance (config, view)) {
+        return psphotReadoutCleanup(config, view);
+    }
     if (!strcasecmp (breakPt, "NOTHING")) {
         return psphotReadoutCleanup(config, view);

trunk/psphot/src/psphotReadoutCleanup.c

@@ -53 +53 @@
     pmPSF        *psf        = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
-    psArray      *sources    = detections->allSources;
+    psArray      *sources    = detections ? detections->allSources : NULL;
     // XXX where do we free these, in here (psMetadataRemove?)
 
@@ -73 +73 @@
     // Check to see if the image quality was measured
     // XXX not sure we want / need this test
-    if (!psf) {
+    if (0 && !psf) {
         bool mdok;                      // Status of MD lookup
         int nIQ = psMetadataLookupS32(&mdok, recipe, "IQ_NSTAR"); // Number of stars for IQ measurement

trunk/psphot/src/psphotRoughClass.c

@@ -19 +19 @@
     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
         if (!psphotRoughClassReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for PSPHOT.INPUT entry %d", i);

trunk/psphot/src/psphotSkyReplace.c

@@ -8 +8 @@
     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
         if (!psphotSkyReplaceReadout (config, view, "PSPHOT.INPUT", i)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to replace sky for PSPHOT.INPUT entry %d", i);

trunk/psphot/src/psphotSourceMatch.c

@@ -1 @@
-# include "psphotInternal.h"
-                                                                                                                                                                               
-// we need a couple of functions to distinguish coincident sources:                                                                                                              
-// XXX these are similar (identical?) to the goals of pmSourceMatch.c                                                                                                            
-                                                                                                                                                                                 
-# define NEXT1 { i++; continue; }                                                                                                                                                
-# define NEXT2 { j++; continue; }                                                                                                                                                
-                                                                                                                                                                                 
-psphotSourceMerge () {                                                                                                                                                           
-                                                                                                                                                                                 
-    float dx, dy;                                                                                                                                                                
-                                                                                                                                                                                 
-    // sort the source list by X                                                                                                                                                 
-    pmSourceSortByX (sources1);                                                                                                                                                  
-    pmSourceSortByX (sources2);                                                                                                                                                  
-                                                                                                                                                                                 
-    int i, j;                                                                                                                                                                    
-    for (i = j = 0; (i < sources1->n) && (j < sources2->n); ) {                                                                                                                  
-                                                                                                                                                                                 
-        pmSource *src1 = sources1->data[i];                                                                                                                                      
-        pmSource *src2 = sources2->data[j];                                                                                                                                      
-                                                                                                                                                                                 
-        if (!src1) NEXT1;                                                                                                                                                        
-        if (!src1->peak) NEXT1;                                                                                                                                                  
-        if (!finite(src1->peak->xf)) NEXT1;                                                                                                                                      
-        if (!finite(src1->peak->yf)) NEXT1;                                                                                                                                      
-                                                                                                                                                                                 
-        if (!src2) NEXT2;                                                                                                                                                        
-        if (!src2->peak) NEXT2;                                                                                                                                                  
-        if (!finite(src2->peak->xf)) NEXT2;                                                                                                                                      
-        if (!finite(src2->peak->yf)) NEXT2;                                                                                                                                      
-                                                                                                                                                                                 
-        dx = src1->peak->xf - src2->peak->xf;                                                                                                                                    
-        if (dx < -1.02*RADIUS) NEXT1;                                                                                                                                            
-        if (dx > +1.02*RADIUS) NEXT2;                                                                                                                                            
-                                                                                                                                                                                 
-        // we are within match range, look for matches:                                                                                                                          
-        for (int J = j; (dx > -1.02*radius) && (J < sources2->n); J++) {                                                                                                         
-                                                                                                                                                                                 
-            dx = src1->peak->xf - src2->peak->xf;                                                                                                                                
-            dy = src1->peak->yf - src2->peak->yf;                                                                                                                                
-                                                                                                                                                                                 
-            dr = dx*dx + dy*dy;                                                                                                                                                  
-            if (dr > RADIUS2) continue;                                                                                                                                          
-                                                                                                                                                                                 
-            // add to group?                                                                                                                                                     
-        }                                                                                                                                                                        
-        i++;                                                                                                                                                                     
-    }                                                                                                                                                                            
-}      
+# include "psphotInternal.h" 
+
+bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects);
+ 
+psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view) 
+{
+    bool status = true;
+
+    psArray *objects = psArrayAllocEmpty(100);
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (!psphotMatchSourcesReadout (objects, config, view, "PSPHOT.INPUT", i)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for PSPHOT.INPUT entry %d", i);
+            psFree (objects);
+            return NULL;
+        }
+    }
+
+    psphotMatchSourcesGenerate (config, view, objects);
+
+    return objects;
+}
+
+bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index) { 
+ 
+    bool status = false;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int RADIUS = psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.MATCH.RADIUS");
+    psAssert (status, "programming error: must define PSPHOT.STACK.MATCH.RADIUS");
+
+    // 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?");
+    psAssert (detections->newSources, "new sources not defined?");
+    psAssert (!detections->allSources, "all sources already defined?");
+
+    // XXX TEST:
+    if (detections->newSources) {
+        psphotMatchSourcesToObjects(objects, detections->newSources, RADIUS);
+    }
+
+    return true;
+}
+
+# define NEXT1 { i++; continue; } 
+# define NEXT2 { j++; continue; } 
+bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS) { 
+ 
+    float dx, dy; 
+ 
+    float RADIUS2 = RADIUS*RADIUS;
+
+    // sort the source list by X 
+    sources = psArraySort (sources, pmSourceSortByX); 
+    objects = psArraySort (objects, pmPhotObjSortByX); 
+ 
+    psVector *found = psVectorAlloc(sources->n, PS_TYPE_U8);
+    psVectorInit (found, 0);
+
+    // match sources to existing objects
+
+    psLogMsg ("psphot", PS_LOG_DETAIL, "attempt to match sources (%ld vs %ld)", sources->n, objects->n);
+
+    int i, j; 
+    for (i = j = 0; (i < sources->n) && (j < objects->n); ) { 
+ 
+        pmSource  *src = sources->data[i]; 
+        pmPhotObj *obj = objects->data[j]; 
+ 
+        if (!src) NEXT1; 
+        if (!src->peak) NEXT1; 
+        if (!finite(src->peak->xf)) NEXT1; 
+        if (!finite(src->peak->yf)) NEXT1; 
+ 
+        if (!obj) NEXT2; 
+        if (!finite(obj->x)) NEXT2; 
+        if (!finite(obj->y)) NEXT2; 
+ 
+        dx = src->peak->xf - obj->x; 
+        if (dx < -1.02*RADIUS) NEXT1; 
+        if (dx > +1.02*RADIUS) NEXT2; 
+ 
+        // we are within match range, look for matches: 
+        int Jmin = -1;
+        float Rmin = RADIUS2;
+        for (int J = j; (dx > -1.02*RADIUS) && (J < objects->n); J++) { 
+ 
+            obj = objects->data[J]; 
+            
+            dx = src->peak->xf - obj->x; 
+            dy = src->peak->yf - obj->y; 
+ 
+            float dr = dx*dx + dy*dy; 
+            if (dr > RADIUS2) continue; 
+            if (dr < Rmin) {
+                Rmin = dr;
+                Jmin  = J;
+            }
+        }
+
+        // no match, try next source
+        if (Jmin == -1) {
+            i++;
+            continue;
+        }
+        obj = objects->data[Jmin]; 
+
+        // add to object
+        pmPhotObjAddSource (obj, src);
+        found->data.U8[i] = 1;
+        i++; 
+    } 
+
+    // add missed sources to new objects
+
+    for (i = 0; i < sources->n; i++) {
+
+        if (found->data.U8[i]) continue;
+
+        pmSource *src = sources->data[i]; 
+
+        pmPhotObj *obj = pmPhotObjAlloc();
+        pmPhotObjAddSource(obj, src);
+        psArrayAdd (objects, 100, obj);
+        psFree (obj);
+    }
+    psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n);
+
+    psFree (found);
+    return true;
+} 
+
+bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects) {
+
+    bool status = false;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    // determine properties (sky, moments) of initial sources
+    float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
+    psAssert (status, "missing SKY_OUTER_RADIUS in recipe?");
+
+    int nImages = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // generate look-up arrays for detections and readouts
+    psArray *detArrays = psArrayAlloc(nImages);
+    psArray *readouts = psArrayAlloc(nImages);
+
+    for (int i = 0; i < nImages; i++) {
+
+        // find the currently selected readout
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 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?");
+
+        detArrays->data[i] = psMemIncrRefCounter(detections);
+        readouts->data[i] = psMemIncrRefCounter(readout);
+    }
+
+    // vector to track if source for an image is found
+    psVector *found = psVectorAlloc(nImages, PS_TYPE_U8);
+
+    for (int i = 0; i < objects->n; i++) { 
+        pmPhotObj *obj = objects->data[i]; 
+
+        // mark the images for which sources have been found
+        psVectorInit (found, 0);
+        for (int j = 0; j < obj->sources->n; j++) {
+
+            pmSource *src = obj->sources->data[j]; 
+            int index = src->imageID;
+            psAssert (index >= 0, "invalid index");
+            psAssert (index < found->n, "invalid index");
+
+            found->data.U8[index] = 1;
+        }
+
+        // generate new sources for the image that are missing
+        for (int index = 0; index < found->n; index++) {
+            if (found->data.U8[index]) continue;
+
+            pmDetections *detections = detArrays->data[index];
+            pmReadout *readout = readouts->data[index];
+            int row0 = readout->image->row0;
+            int col0 = readout->image->col0;
+
+            // XXX the peak type is not really used in psphot
+            // PM_PEAK_LONE is certainly not true, but irrelevant
+            float peakFlux = readout->image->data.F32[(int)(obj->y-row0-0.5)][(int)(obj->x-col0-0.5)];
+            pmPeak *peak = pmPeakAlloc(obj->x, obj->y, peakFlux, PM_PEAK_LONE);
+            peak->flux = peakFlux;
+            peak->SN = 1.0;
+            peak->xf = obj->x;
+            peak->yf = obj->y;
+            peak->dx = NAN;
+            peak->dy = NAN;
+            
+            // XXX assign to a footprint?
+
+            // create a new source
+            pmSource *source = pmSourceAlloc();
+            source->imageID = index;
+
+            // add the peak
+            source->peak = psMemIncrRefCounter(peak);
+
+            // allocate space for moments
+            source->moments = pmMomentsAlloc();
+
+            // allocate image, weight, mask arrays for each peak (square of radius OUTER)
+            pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
+
+            peak->assigned = true;
+            pmPhotObjAddSource(obj, source);
+            psArrayAdd (detections->newSources, 100, source);
+            psFree (source);
+            psFree (peak);
+        }
+    }
+
+    // how many sources do we have now?
+    int nSources = 0;
+    for (int i = 0; i < objects->n; i++) { 
+        pmPhotObj *obj = objects->data[i]; 
+        nSources += obj->sources->n;
+    }
+    psLogMsg ("psphot", PS_LOG_DETAIL, "total of %d sources for %d images", nSources, nImages);
+
+
+    psFree (found);
+    psFree (detArrays);
+    psFree (readouts);
+    return true;
+}

trunk/psphot/src/psphotSourceSize.c

@@ -44 +44 @@
     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
         if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe, getPSFsize)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on source size analysis for PSPHOT.INPUT entry %d", i);

trunk/psphot/src/psphotSourceStats.c

@@ -2 +2 @@
 
 // convert detections to sources and measure their basic properties (moments, local sky, sky
-// variance) Note: this function only generates sources for the new peaks (peak->assigned)
+// variance) Note: this function only generates sources for the new peaks (peak->assigned).
+// The new sources are added to any existing sources on detections->newSources.  The sources 
+// on detections->allSources are ignored.
 bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow)
 {
@@ -40 +42 @@
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     psAssert (detections, "missing detections?");
-    psAssert (!detections->newSources, "new sources already defined?");
 
     // XXX TEST:
@@ -84 +85 @@
 
     // generate the array of sources, define the associated pixel
-    sources = psArrayAllocEmpty (peaks->n);
+    if (!detections->newSources) {
+        detections->newSources = psArrayAllocEmpty (peaks->n);
+    }
+    sources = detections->newSources;
 
     // if there are no peaks, we save the empty source array and return
     if (!peaks->n) {
-        // save the new sources on the detection structure:
-        detections->newSources = sources;
         return true;
     }
@@ -100 +102 @@
         // create a new source
         pmSource *source = pmSourceAlloc();
+        source->imageID = index;
 
         // add the peak
@@ -119 +122 @@
         psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
         psphotVisualShowMoments (sources);
-        detections->newSources = sources;
         return true;
     }
@@ -126 +128 @@
         if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) {
             psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
-            psFree(sources);
+            psFree(detections->newSources);
             return false;
         }
@@ -178 +180 @@
                 psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
                 psFree (job);
-                psFree(sources);
+                psFree(detections->newSources);
                 return false;
             }
@@ -187 +189 @@
         if (!psThreadPoolWait (false)) {
             psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
-            psFree(sources);
+            psFree(detections->newSources);
             return false;
         }
@@ -215 +217 @@
     psphotVisualShowMoments (sources);
 
-    // save the new sources on the detection structure:
-    detections->newSources = sources;
-
-
     if (detections->allSources) {
         psphotMaskCosmicRayFootprintCheck(detections->allSources);
@@ -377 +375 @@
         pmSource *source = sources->data[i];
 
+        if (source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue;
+        source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED;
+
         // skip faint sources for moments measurement
         if (source->peak->SN < MIN_SN) {

trunk/psphot/src/psphotStackArguments.c

@@ -2 +2 @@
 
 static void usage(pmConfig *config, int exitCode);
-static void writeHelpInfo(pmConfig* config, FILE* ofile);
+static void writeHelpInfo(FILE* ofile);
 
 pmConfig *psphotStackArguments(int argc, char **argv) {
 
     int N;
-    bool status;
 
     // print help info
-    if (psArgumentGet(argc, argv, "-help")) writeHelpInfo(argv[0], config, stdout);
-    if (psArgumentGet(argc, argv, "-h")) writeHelpInfo(argv[0], config, stdout);
+    if (psArgumentGet(argc, argv, "-help")) writeHelpInfo(stdout);
+    if (psArgumentGet(argc, argv, "-h")) writeHelpInfo(stdout);
 
     // load config data from default locations
@@ -29 +28 @@
     // Number of threads is handled
     PSARGUMENTS_INSTANTIATE_THREADSARG( psphot, config, argc, argv )
-
-    // photcode : used in output to supplement header data (argument or recipe?)
-    if ((N = psArgumentGet (argc, argv, "-photcode"))) {
-        if (argc <= N+1) {
-          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
-          exit(PS_EXIT_CONFIG_ERROR);
-        }
-        psArgumentRemove (N, &argc, argv);
-        psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", argv[N]);
-        psArgumentRemove (N, &argc, argv);
-    }
 
     // visual : interactive display mode
@@ -90 +78 @@
 }
 
-static void writeHelpInfo(pmConfig* config, FILE* ofile)
+static void writeHelpInfo(FILE* ofile)
 {
   fprintf(ofile,
-          "Usage: psphotStack -input (input.mdc) outroot\n"
+          "Usage: psphotStack -input (INPUTS.mdc) (OUTROOT)\n"
           "\n"
-          "where:\n"
-          "  FileNameList is a text file containing filenames, one per line\n"
-          "  MaskFileNameList is a text file of mask filenames, one per line\n"
-          "  VarFileNameList is a text file of variance filenames, one per line\n"
-          "  OutFileBaseName is the 'root name' for output files\n"
+          "where INPUTS.mdc contains various METADATAs, each with:\n"
+          "\tIMAGE(STR):     Image filename\n"
+          "\tMASK(STR):      Mask filename\n"
+          "\tVARIANCE(STR):  Variance map filename\n"
+          "OUTROOT is the 'root name' for output files\n"
           "\n"
           "additional options:\n"
@@ -136 +124 @@
           "  -trace-levels\n"
           "     print current trace levels\n");
-    psFree(config);
-    pmConfigDone();
     psLibFinalize();
     exit(PS_EXIT_SUCCESS);

trunk/psphot/src/psphotStackChisqImage.c

@@ -12 +12 @@
     psTimerStart ("psphot.chisq.image");
 
-    // create a holder for the image
-    pmReadout *chiImage = pmReadoutAlloc();
+    pmFPAfile *chisqFile = pmFPAfileSelectSingle(config->files, "PSPHOT.CHISQ.IMAGE", 0);
+    psAssert (chisqFile, "missing chisq image FPA?");
+
+    pmReadout *chiReadout = NULL;
 
     int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
@@ -20 +22 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (!psphotStackChisqImageAddReadout(config, view, chiImage, "PSPHOT.INPUT", i)) {
+        if (!psphotStackChisqImageAddReadout(config, view, &chiReadout, "PSPHOT.INPUT", i)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to model background for PSPHOT.INPUT entry %d", i);
             return false;
@@ -26 +28 @@
     }
 
+    num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.CHISQ.NUM", PS_META_REPLACE, "", num);
+    num++;
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", num);
+
     // need to save the resulting image somewhere (PSPHOT.INPUT?)
+    if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, "PSPHOT.INPUT", PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) {
+        psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files");
+        return false;
+    }
 
     psLogMsg ("psphot", PS_LOG_INFO, "built chisq image: %f sec\n", psTimerMark ("psphot.chisq.image"));
@@ -34 +47 @@
 
 bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration
-                                     pmFPAview *view,
-                                     pmReadout *chiReadout,
+                                     const pmFPAview *view,
+                                     pmReadout **chiReadout,
                                      char *filename, 
                                      int index) 
@@ -49 +62 @@
 
     psImage *inImage     = inReadout->image;
+    psAssert (inImage, "missing image?");
+
     psImage *inVariance  = inReadout->variance;
+    psAssert (inVariance, "missing variance?");
+
     psImage *inMask      = inReadout->mask;
+    psAssert (inMask, "missing mask?");
 
-    psImage *chiImage    = chiReadout->image;
-    psImage *chiVariance = chiReadout->variance;
-    psImage *chiMask     = chiReadout->mask;
+    if (*chiReadout == NULL) {
+        *chiReadout = pmFPAGenerateReadout(config, view, "PSPHOT.CHISQ.IMAGE", input->fpa, NULL, 0);
+    }
+
+    psImage *chiImage = (*chiReadout)->image;
+    psAssert (chiImage, "missing chi image");
+
+    psImage *chiVariance = (*chiReadout)->variance;
+    psAssert (chiVariance, "missing chi variance");
+
+    psImage *chiMask = (*chiReadout)->mask;
+    psAssert (chiMask, "missing chi mask");
 
     // select the appropriate recipe information
@@ -83 +110 @@
     return true;
 }
+
+bool psphotStackRemoveChisqFromInputs (pmConfig *config) {
+
+    bool status = false;
+
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    psAssert (status, "programming error: must define PSPHOT.CHISQ.NUM");
+
+    int inputNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    pmFPAfileRemoveSingle (config->files, "PSPHOT.INPUT", chisqNum);
+
+    inputNum --;
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", inputNum);
+
+    return true;
+}
+
+bool pmFPAfileRemoveSingle(psMetadata *files, const char *name, int num)
+{
+    PS_ASSERT_PTR_NON_NULL(files, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(num, NULL);
+
+    psList* mdList = files->list;
+    psHash* mdHash = files->hash;
+
+    // Generate a REGEX to select only items that match 'name'
+    psString regex = NULL;              // Regular expression
+    if (name) {
+        if (!psMetadataLookup(files, name)) {
+            // No files match the requested name
+            return false;
+        }
+        psStringAppend(&regex, "^%s$", name);
+    }
+
+    psMetadataIterator *iter = psMetadataIteratorAlloc(files, PS_LIST_HEAD, regex); // Iterator
+    psFree(regex);
+    psMetadataItem *item;               // Item from iteration
+
+    bool found = false;
+    int i = 0;                          // Counter
+    for (i = 0; !found && (item = psMetadataGetAndIncrement(iter)); i++) {
+        if (i == num) found = true;
+    }
+    psFree(iter);
+    if (!found) {
+        return false;
+    }
+
+    char *key = item->name;
+
+    // look up the name via hash to see if we have a multi or not
+    psMetadataItem* hashItem = psHashLookup(mdHash, name);
+    if (hashItem == NULL) {
+        psError(PS_ERR_UNKNOWN, false, _("Failed to remove metadata item, %s, from metadata table."), key);
+        return false;
+    }
+
+    if (hashItem->type == PS_DATA_METADATA_MULTI) {
+        // multiple entries with same key, remove just the specified one
+        psListRemoveData(hashItem->data.list, item);
+    } else {
+        psHashRemove(mdHash, key);
+    }
+    psListRemoveData(mdList, item);
+
+    return true;
+}

trunk/psphot/src/psphotStackImageLoop.c

@@ -7 +7 @@
 }
 
-bool psphotImageLoop (pmConfig *config) {
+bool psphotStackImageLoop (pmConfig *config) {
 
     bool status;
@@ -14 +14 @@
     pmReadout *readout;
 
+    pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
+    if (!status) {
+        psError(PSPHOT_ERR_PROG, false, "Can't find input data!");
+        return false;
+    }
+
     pmFPAview *view = pmFPAviewAlloc (0);
 
@@ -19 +25 @@
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
 
-    // XXX for now, we assume there is only a single chip in the PHU:
-    psphotStackReadout (config, view);
+    // for psphot, we force data to be read at the chip level
+    while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
+        psLogMsg ("psphot", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (! chip->process || ! chip->file_exists) { continue; }
+        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphotStack.");
+
+        // there is now only a single chip (multiple readouts?). loop over it and process
+        while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
+            psLogMsg ("psphot", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+            if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Cell in psphotStack.");
+
+            // process each of the readouts
+            while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
+                psLogMsg ("psphot", 6, "Readout %d: %x %x\n", view->readout, cell->file_exists, cell->process);
+                if (! readout->data_exists) { continue; }
+
+                // XXX for now, we assume there is only a single chip in the PHU:
+                if (!psphotStackReadout (config, view)) {
+                    psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
+                    psFree (view);
+                    return false;
+                }
+
+            }
+            // drop all versions of the internal files
+            status = true;
+            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");
+            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV");
+            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND");
+            if (!status) {
+                psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files");
+                psFree (view);
+                return false;
+            }
+        }
+        // save output which is saved at the chip level
+        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot.");
+    }
+    // save output which is saved at the fpa level
+    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot.");
 
     // fail if we failed to handle an error
@@ -32 +76 @@
    the easiest way to implement this is to assume we can pre-load the full set of images up front.
    with 5 filters and 6000^2 (image, mask, var = 10 byte per pixel), we need 1.8GB, which is not too bad.
- */
+*/
 

trunk/psphot/src/psphotStackParseCamera.c

@@ -15 +15 @@
     }
 
-    psMetadata *item == NULL;
-    int nInputs = 0;
-    while ((item = psMetadataGetIndex(inputs, nInputs)) != NULL) {
+    int nInputs = inputs->list->n;
+    for (int i = 0; i < nInputs; i++) {
+        psMetadataItem *item = psMetadataGet(inputs, i);
         if (item->type != PS_DATA_METADATA) {
             psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Component %s of the input metadata is not of type METADATA", item->name);
@@ -40 +40 @@
         psString mask = psMetadataLookupStr(&status, input, "MASK"); // Name of mask
         if (mask && strlen(mask) > 0) {
-            if (!defineFile(config, imageFile, "PSPHOT.INPUT.MASK", mask, PM_FPA_FILE_MASK)) {
+            if (!defineFile(config, imageFile, "PSPHOT.MASK", mask, PM_FPA_FILE_MASK)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask);
                 return false;
@@ -46 +46 @@
         }
 
-        psString variance = psMetadataLookupStr(&mdok, input, "VARIANCE"); // Name of variance map
+        psString variance = psMetadataLookupStr(&status, input, "VARIANCE"); // Name of variance map
         if (variance && strlen(variance) > 0) {
-            if (!defineFile(config, imageFile, "PSPHOT.INPUT.VARIANCE", variance, PM_FPA_FILE_VARIANCE)) {
+            if (!defineFile(config, imageFile, "PSPHOT.VARIANCE", variance, PM_FPA_FILE_VARIANCE)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance);
                 return false;
             }
         }
-        nInputs ++;
+        // the output sources are carried on the input->fpa structures
+        pmFPAfile *outsources = pmFPAfileDefineOutputFromFile (config, imageFile, "PSPHOT.STACK.OUTPUT");
+        if (!outsources) {
+            psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.OUTPUT");
+            return false;
+        }
+        outsources->save = true;
+        outsources->fileID = i;         // this is used to generate output names
     }
     psMetadataRemoveKey(config->arguments, "FILENAMES");
     psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", nInputs);
 
+    if (!psphotSetMaskBits (config)) {
+        psError (PS_ERR_UNKNOWN, false, "failed to set mask bit values");
+        return NULL;
+    }
+
+    // generate an pmFPAimage for the chisqImage
+    pmFPAfile *chisqImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.CHISQ.IMAGE");
+    if (!chisqImage) {
+        psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.CHISQ.IMAGE");
+        return false;
+    }
+    chisqImage->save = true;
+
+    pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.MASK");
+    if (!chisqMask) {
+        psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.MASK"));
+        return NULL;
+    }
+    if (chisqMask->type != PM_FPA_FILE_MASK) {
+        psError(PS_ERR_IO, true, "PSPHOT.CHISQ.MASK is not of type MASK");
+        return NULL;
+    }
+    chisqMask->save = true;
+
+    pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.VARIANCE");
+    if (!chisqVariance) {
+        psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.VARIANCE"));
+        return NULL;
+    }
+    if (chisqVariance->type != PM_FPA_FILE_VARIANCE) {
+        psError(PS_ERR_IO, true, "PSPHOT.CHISQ.VARIANCE is not of type VARIANCE");
+        return NULL;
+    }
+    chisqVariance->save = true;
+
+# if (0)    
     // define the additional input/output files associated with psphot
     // XXX figure out which files are needed by psphotStack
@@ -64 +107 @@
         return false;
     }
+# endif
 
     psTrace("psphot", 1, "Done with psphotStackParseCamera...\n");

trunk/psphot/src/psphotStackReadout.c

@@ -13 +13 @@
         return false;
     }
+    // optional break-point for processing
+    char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
+    PS_ASSERT_PTR_NON_NULL (breakPt, false);
 
-    // set the photcode for this image (XXX currently goes into recipe, should it go into analysis?)
+    // set the photcode for each image
     if (!psphotAddPhotcode (config, view)) {
         psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
@@ -20 +23 @@
     }
 
-    // generate a background model (median, smoothed image)
-    if (!psphotSetMaskAndVariance (config, view, recipe)) {
-        return psphotReadoutCleanup (config, NULL, recipe, NULL, NULL, NULL);
+    // Generate the mask and weight images
+    if (!psphotSetMaskAndVariance (config, view)) {
+        return psphotReadoutCleanup (config, view);
     }
     if (!strcasecmp (breakPt, "NOTHING")) {
-        return psphotReadoutCleanup(config, NULL, recipe, NULL, NULL, NULL);
+        return psphotReadoutCleanup (config, view);
     }
 
-    // optional break-point for processing
-    char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
-    PS_ASSERT_PTR_NON_NULL (breakPt, false);
-
     // generate a background model (median, smoothed image)
+    // XXX I think this is not defined correctly for an array of images.
     if (!psphotModelBackground (config, view)) {
-        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
+        return psphotReadoutCleanup (config, view);
     }
     if (!psphotSubtractBackground (config, view)) {
-        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
+        return psphotReadoutCleanup (config, view);
     }
     if (!strcasecmp (breakPt, "BACKMDL")) {
-        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
+        return psphotReadoutCleanup (config, view);
+    }
+
+    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are determined and saved on
+    // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT
+    if (!psphotLoadPSF (config, view)) {
+        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
+        return psphotReadoutCleanup (config, view);
+    }
+
+    if (!psphotStackChisqImage(config, view)) {
+        psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
+        return psphotReadoutCleanup (config, view);
+    }
+    if (!strcasecmp (breakPt, "CHISQ")) {
+        return psphotReadoutCleanup (config, view);
     }
 
     // find the detections (by peak and/or footprint) in the image.
-    pmDetections *detections = psphotFindDetections (NULL, readout, recipe);
-    if (!detections) {
-        psLogMsg ("psphot", 3, "unable to find detections in this image");
-        return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL);
+    // This finds the detections on Chisq image as well as the individuals
+    if (!psphotFindDetections (config, view, true)) { // pass 1
+        // this only happens if we had an error in psphotFindDetections
+        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
+        return psphotReadoutCleanup (config, view);
     }
 
-    // construct sources and measure basic stats
-    psArray *sources = psphotSourceStats (config, readout, detections, true);
-    if (!sources) return false;
-    if (!strcasecmp (breakPt, "PEAKS")) {
-        return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
+    // construct sources and measure basic stats (saved on detections->newSources)
+    // only run this on detections from the input images, not chisq image
+    if (!psphotSourceStats (config, view, true)) { // pass 1
+        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
+        return psphotReadoutCleanup (config, view);
     }
 
-    // find blended neighbors of very saturated stars
-    // XXX merge this with Basic Deblend?
-    psphotDeblendSatstars (readout, sources, recipe);
+    // *** generate the objects (which unify the sources from the different images)
+    psArray *objects = psphotMatchSources (config, view);
 
-    // mark blended peaks PS_SOURCE_BLEND
-    if (!psphotBasicDeblend (sources, recipe)) {
-        psLogMsg ("psphot", 3, "failed on deblend analysis");
-        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
+    // construct sources for the newly-generated sources (from other images)
+    if (!psphotSourceStats (config, view, false)) { // pass 1
+        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
+        return psphotReadoutCleanup (config, view);
     }
 
+    // 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);
+    // }
+
+    // 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);
+    // }
+
     // classify sources based on moments, brightness
-    if (!psphotRoughClass (readout, sources, recipe, havePSF)) {
-        psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
-        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
+    // only run this on detections from the input images, not chisq image
+    if (!psphotRoughClass (config, view)) {
+        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
+        return psphotReadoutCleanup (config, view);
+    }
+    // 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)) { // pass 1
+        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
+        return psphotReadoutCleanup(config, view);
     }
     if (!strcasecmp (breakPt, "MOMENTS")) {
-        return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
+        return psphotReadoutCleanup (config, view);
     }
 
-    if (!havePSF && !psphotImageQuality (recipe, sources)) {
-        psLogMsg("psphot", 3, "failed to measure image quality");
-        return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
+    // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
+    // this step is skipped
+    if (!psphotChoosePSF (config, view)) { // pass 1
+        psLogMsg ("psphot", 3, "failure to construct a psf model");
+        return psphotReadoutCleanup (config, view);
+    }
+    if (!strcasecmp (breakPt, "PSFMODEL")) {
+        return psphotReadoutCleanup (config, view);
     }
 
-    // if we were not supplied a PSF, choose one here
-    if (psf == NULL) {
-        // use bright stellar objects to measure PSF
-        // XXX if we do not have enough stars to generate the PSF, build one
-        // from the SEEING guess and model class
-        psf = psphotChoosePSF (readout, sources, recipe);
-        if (psf == NULL) {
-            psLogMsg ("psphot", 3, "failure to construct a psf model");
-            return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
-        }
-        havePSF = true;
-    }
-    if (!strcasecmp (breakPt, "PSFMODEL")) {
-        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
-    }
-    psphotVisualShowPSFModel (readout, psf);
-
     // include externally-supplied sources
-    psphotLoadExtSources (config, view, sources);
+    // XXX fix this in the new multi-input context
+    // psphotLoadExtSources (config, view); // pass 1
 
     // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
-    psphotGuessModels (config, readout, sources, psf);
+    psphotGuessModels (config, view);
+
+    // merge the newly selected sources into the existing list
+    // NOTE: merge OLD and NEW
+    psphotMergeSources (config, view);
 
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
-    psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
-
-    // We have to place these visualizations here because the models are not realized until
-    // psphotGuessModels or fitted until psphotFitSourcesLinear.
-    psphotVisualShowPSFStars (recipe, psf, sources);
+    psphotFitSourcesLinearStack (config, objects, FALSE);
+    psFree (objects);
 
     // identify CRs and extended sources
-    psphotSourceSize (config, readout, sources, recipe, psf, 0);
-    if (!strcasecmp (breakPt, "ENSEMBLE")) {
-        goto finish;
-    }
-    psphotVisualShowSatStars (recipe, psf, sources);
-
-    // non-linear PSF and EXT fit to brighter sources
-    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
-    psphotBlendFit (config, readout, sources, psf);
-
-    // replace all sources
-    psphotReplaceAllSources (sources, recipe);
-
-    // linear fit to include all sources (subtract again)
-    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
-
-    // if we only do one pass, skip to extended source analysis
-    if (!strcasecmp (breakPt, "PASS1")) {
-        goto pass1finish;
-    }
-    // NOTE: possibly re-measure background model here with objects subtracted
-
-    // add noise for subtracted objects
-    psphotAddNoise (readout, sources, recipe);
-
-    // find fainter sources (pass 2)
-    detections = psphotFindDetections (detections, readout, recipe);
-
-    // remove noise for subtracted objects (ie, return to normal noise level)
-    psphotSubNoise (readout, sources, recipe);
-
-    // define new sources based on only the new peaks
-    psArray *newSources = psphotSourceStats (config, readout, detections, false);
-
-    // set source type
-    if (!psphotRoughClass (readout, newSources, recipe, havePSF)) {
-        psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
-        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
-    }
-
-    // create full input models, set the radius to fitRadius, set circular fit mask
-    psphotGuessModels (config, readout, newSources, psf);
-
-    // replace all sources so fit below applies to all at once
-    psphotReplaceAllSources (sources, recipe);
-
-    // merge the newly selected sources into the existing list
-    psphotMergeSources (sources, newSources);
-    psFree (newSources);
-
-    // linear fit to all sources
-    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
-
-pass1finish:
-
-    // measure source size for the remaining sources
-    psphotSourceSize (config, readout, sources, recipe, psf, 0);
-
-    psphotExtendedSourceAnalysis (readout, sources, recipe);
-
-    psphotExtendedSourceFits (readout, sources, recipe);
-
-finish:
-
-    // plot positive sources
-    // psphotSourcePlots (readout, sources, recipe);
+    psphotSourceSize (config, view, TRUE);
 
     // measure aperture photometry corrections
-    if (!psphotApResid (config, readout, sources, psf)) {
+    if (!psphotApResid (config, view)) {
         psLogMsg ("psphot", 3, "failed on psphotApResid");
-        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
+        return psphotReadoutCleanup (config, view);
     }
 
     // calculate source magnitudes
-    psphotMagnitudes(config, readout, view, sources, psf);
+    psphotMagnitudes(config, view);
 
-    if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) {
+    if (!psphotEfficiency(config, view)) {
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
         psErrorClear();
@@ -200 +159 @@
 
     // drop the references to the image pixels held by each source
-    psphotSourceFreePixels (sources);
+    psphotSourceFreePixels (config, view);
+
+    // remove chisq image from config->file:PSPHOT.INPUT (why?)
+    psphotStackRemoveChisqFromInputs(config);
 
     // create the exported-metadata and free local data
-    return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
+    return psphotReadoutCleanup (config, view);
 }
 

trunk/psphot/src/psphotSubtractBackground.c

@@ -23 +23 @@
     pmFPAfile *modelFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest
     assert (modelFile);
-    pmReadout *model = pmFPAviewThisReadout (view, modelFile->fpa);
+
+    pmReadout *model = NULL;
+    if (modelFile->mode == PM_FPA_MODE_INTERNAL) {
+        model = modelFile->readout;
+    } else {
+        model = pmFPAviewThisReadout (view, modelFile->fpa);
+    }
     assert (model);
 

@@ -21 +21 @@
 psphotForced
 psphotMakePSF
+psphotStack