Changeset 27649

Timestamp:

Apr 9, 2010, 4:41:16 PM (16 years ago)

Author:

eugene

Message:

working on psphotStack (nearly done)

Location:

branches/eam_branches/stackphot.20100406/psphot/src

Files:

• Makefile.am (modified) (1 diff)
• psphot.h (modified) (1 diff)
• psphotApResid.c (modified) (1 diff)
• psphotChoosePSF.c (modified) (1 diff)
• psphotEfficiency.c (modified) (1 diff)
• psphotFitSourcesLinearStack.c (modified) (9 diffs)
• psphotGuessModels.c (modified) (1 diff)
• psphotImageLoop.c (modified) (1 diff)
• psphotImageQuality.c (modified) (1 diff)
• psphotMagnitudes.c (modified) (1 diff)
• psphotModelBackground.c (modified) (1 diff)
• psphotReadoutCleanup.c (modified) (2 diffs)
• psphotRoughClass.c (modified) (1 diff)
• psphotSkyReplace.c (modified) (1 diff)
• psphotSourceMatch.c (modified) (1 diff)
• psphotSourceSize.c (modified) (1 diff)
• psphotStackChisqImage.c (modified) (6 diffs)
• psphotStackImageLoop.c (modified) (3 diffs)
• psphotStackParseCamera.c (modified) (5 diffs)
• psphotStackReadout.c (modified) (5 diffs)
• psphotSubtractBackground.c (modified) (1 diff)

branches/eam_branches/stackphot.20100406/psphot/src/Makefile.am

@@ -87 +87 @@
 # a psphot-variant for stack photometry
 psphotStack_SOURCES = \
-        psphotStack.c              \
-        psphotStackArguments.c     \
-        psphotStackParseCamera.c   \
-        psphotStackImageLoop.c     \
-        psphotStackReadout.c       \
-        psphotStackChisqImage.c    \
-        psphotCleanup.c
-
-# psphotFitSourceLinearStack.c    
+        psphotStack.c                 \
+        psphotStackArguments.c        \
+        psphotStackParseCamera.c      \
+        psphotStackImageLoop.c        \
+        psphotStackReadout.c          \
+        psphotStackChisqImage.c       \
+        psphotFitSourcesLinearStack.c \
+        psphotSourceMatch.c           \
+        psphotCleanup.c
+
 
 

branches/eam_branches/stackphot.20100406/psphot/src/psphot.h

@@ -350 +350 @@
 bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration
                                      const pmFPAview *view,
-                                     pmReadout *chiReadout,
+                                     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

branches/eam_branches/stackphot.20100406/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);

branches/eam_branches/stackphot.20100406/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);

branches/eam_branches/stackphot.20100406/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);

branches/eam_branches/stackphot.20100406/psphot/src/psphotFitSourcesLinearStack.c

@@ -3 +3 @@
 // XXX define the 'good' / 'bad' flags?
 
+# define COVAR_FACTOR 1.0
+
 bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final) {
 
     bool status;
-    float x;
-    float y;
     float f;
 
@@ -29 +29 @@
     // analysis is done in spatial order (to speed up overlap search)
     // sort by first element in each source list
-    objects = psArraySort (objects, pmPhotObjSortByY);
+    objects = psArraySort (objects, pmPhotObjSortByX);
 
     // storage array for fitSources
@@ -66 +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) {
@@ -86 +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 {
@@ -98 +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);
 
@@ -109 +109 @@
 
             // 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);
         }
@@ -157 +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);
@@ -172 +172 @@
 }
 
-// sort by Y (ascending)
-int pmPhotObjSortBySN (const void **a, const void **b)
+// sort by X (ascending)
+int pmPhotObjSortByX (const void **a, const void **b)
 {
     pmPhotObj *objA = *(pmPhotObj **)a;
     pmPhotObj *objB = *(pmPhotObj **)b;
 
-    psAssert (objA->sources, "missing sources?");
-    psAssert (objB->sources, "missing sources?");
-
-    psAssert (objA->sources->n, "missing sources");
-    psAssert (objB->sources->n, "missing sources");
-
-    psAssert (objA->sources->data[0], "missing sources");
-    psAssert (objB->sources->data[0], "missing sources");
-
-    pmSource *A = objA->sources->data[0];
-    pmSource *B = objB->sources->data[0];
-
-    psF32 fA = (A->peak == NULL) ? 0 : A->peak->y;
-    psF32 fB = (B->peak == NULL) ? 0 : B->peak->y;
+    psF32 fA = objA->x;
+    psF32 fB = objB->x;
 
     psF32 diff = fA - fB;
@@ -198 +186 @@
     return (0);
 }
-

branches/eam_branches/stackphot.20100406/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);

branches/eam_branches/stackphot.20100406/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.");

branches/eam_branches/stackphot.20100406/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);

branches/eam_branches/stackphot.20100406/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

branches/eam_branches/stackphot.20100406/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)) {

branches/eam_branches/stackphot.20100406/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

branches/eam_branches/stackphot.20100406/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);

branches/eam_branches/stackphot.20100406/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);

branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceMatch.c

@@ -14 +14 @@
         if (!psphotMatchSourcesReadout (objects, config, view, "PSPHOT.INPUT", i)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for PSPHOT.INPUT entry %d", i);
-            return false;
+            psFree (objects);
+            return NULL;
         }
     }
+
+    // psphotMatchSourcesGenerate (objects, config, view, "PSPHOT.INPUT", i);
+
+    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;
 }
 
-// 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; } 
- 
-bool psphotSourceMerge (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index) { 
+bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS) { 
  
     float dx, dy; 
  
+    float RADIUS2 = RADIUS*RADIUS;
+
     // sort the source list by X 
-    pmSourceSortByX (sources1); 
-    pmSourceSortByX (sources2); 
+    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 < sources1->n) && (j < sources2->n); ) { 
+    for (i = j = 0; (i < sources->n) && (j < objects->n); ) { 
  
-        pmSource *src1 = sources1->data[i]; 
-        pmSource *src2 = sources2->data[j]; 
+        pmSource  *src = sources->data[i]; 
+        pmPhotObj *obj = objects->data[j]; 
  
-        if (!src1) NEXT1; 
-        if (!src1->peak) NEXT1; 
-        if (!finite(src1->peak->xf)) NEXT1; 
-        if (!finite(src1->peak->yf)) NEXT1; 
+        if (!src) NEXT1; 
+        if (!src->peak) NEXT1; 
+        if (!finite(src->peak->xf)) NEXT1; 
+        if (!finite(src->peak->yf)) NEXT1; 
  
-        if (!src2) NEXT2; 
-        if (!src2->peak) NEXT2; 
-        if (!finite(src2->peak->xf)) NEXT2; 
-        if (!finite(src2->peak->yf)) NEXT2; 
+        if (!obj) NEXT2; 
+        if (!finite(obj->x)) NEXT2; 
+        if (!finite(obj->y)) NEXT2; 
  
-        dx = src1->peak->xf - src2->peak->xf; 
+        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: 
-        for (int J = j; (dx > -1.02*radius) && (J < sources2->n); J++) { 
+        int Jmin = -1;
+        float Rmin = RADIUS2;
+        for (int J = j; (dx > -1.02*RADIUS) && (J < objects->n); J++) { 
  
-            dx = src1->peak->xf - src2->peak->xf; 
-            dy = src1->peak->yf - src2->peak->yf; 
+            obj = objects->data[J]; 
+            
+            dx = src->peak->xf - obj->x; 
+            dy = src->peak->yf - obj->y; 
  
-            dr = dx*dx + dy*dy; 
+            float dr = dx*dx + dy*dy; 
             if (dr > RADIUS2) continue; 
- 
-            // add to group? 
-        } 
+            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);
+    }
+    psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n);
+
+    return true;
 } 

branches/eam_branches/stackphot.20100406/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);

branches/eam_branches/stackphot.20100406/psphot/src/psphotStackChisqImage.c

@@ -12 +12 @@
     psTimerStart ("psphot.chisq.image");
 
-    pmFPAfile *chisqFile = pmFPAfileSelectSingle(config->files, "PSPHOT.CHISQ.OUTPUT", 0);
+    pmFPAfile *chisqFile = pmFPAfileSelectSingle(config->files, "PSPHOT.CHISQ.IMAGE", 0);
     psAssert (chisqFile, "missing chisq image FPA?");
 
-    pmCell *chisqCell = pmFPAviewThisCell(view, chisqFile->fpa);
-
-    // create a holder for the image
-    pmReadout *chiReadout = pmFPAviewThisReadout(view, chisqFile->fpa);
-    if (!chiReadout) {
-      chiReadout = pmReadoutAlloc(chisqCell);
-    } else {
-      psMemIncrRefCounter(chiReadout);
-    }
+    pmReadout *chiReadout = NULL;
 
     int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
@@ -30 +22 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (!psphotStackChisqImageAddReadout(config, view, chiReadout, "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;
@@ -51 +43 @@
     psLogMsg ("psphot", PS_LOG_INFO, "built chisq image: %f sec\n", psTimerMark ("psphot.chisq.image"));
 
-    psFree (chiReadout);
     return true;
 }
@@ -57 +48 @@
 bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration
                                      const pmFPAview *view,
-                                     pmReadout *chiReadout,
+                                     pmReadout **chiReadout,
                                      char *filename, 
                                      int index) 
@@ -71 +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
@@ -105 +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;
+}

branches/eam_branches/stackphot.20100406/psphot/src/psphotStackImageLoop.c

@@ -9 +9 @@
 bool psphotStackImageLoop (pmConfig *config) {
 
+    bool status;
+    pmChip *chip;
+    pmCell *cell;
+    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);
 
@@ -14 +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
@@ -27 +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.
- */
+*/
 

branches/eam_branches/stackphot.20100406/psphot/src/psphotStackParseCamera.c

@@ -15 +15 @@
     }
 
-    psMetadataItem *item = NULL;
-    int nInputs = 0;
-    while ((item = psMetadataGet(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);
@@ -34 +34 @@
         pmFPAfile *imageFile = defineFile(config, NULL, "PSPHOT.INPUT", image, PM_FPA_FILE_IMAGE); // File for image
         if (!imageFile) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", nInputs, image);
+            psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, image);
             return false;
         }
@@ -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)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", nInputs, 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;
             }
@@ -48 +48 @@
         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)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", nInputs, 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
@@ -64 +76 @@
         return false;
     }
-    pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PPIMAGE.CHISQ.MASK");
+    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 PPIMAGE.CHISQ.MASK"));
+        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, "PPIMAGE.CHISQ.MASK is not of type MASK");
+        psError(PS_ERR_IO, true, "PSPHOT.CHISQ.MASK is not of type MASK");
         return NULL;
     }
-    pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PPIMAGE.CHISQ.VARIANCE");
+    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 PPIMAGE.CHISQ.VARIANCE"));
+        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, "PPIMAGE.CHISQ.VARIANCE is not of type VARIANCE");
+        psError(PS_ERR_IO, true, "PSPHOT.CHISQ.VARIANCE is not of type VARIANCE");
         return NULL;
     }
+    chisqVariance->save = true;
 
 # if (0)    

branches/eam_branches/stackphot.20100406/psphot/src/psphotStackReadout.c

@@ -59 +59 @@
 
     // find the detections (by peak and/or footprint) in the image.
+    // This finds the detections on Chisq image as well as the individuals
     if (!psphotFindDetections (config, view, true)) { // pass 1
         // this only happens if we had an error in psphotFindDetections
@@ -66 +67 @@
 
     // 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);
     }
-    if (!strcasecmp (breakPt, "PEAKS")) {
-        return psphotReadoutCleanup(config, view);
+
+    // *** generate the objects (which unify the sources from the different images)
+    psArray *objects = psphotMatchSources (config, view);
+
+    // 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);
-    }
+    // 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);
-    }
+    // if (!psphotBasicDeblend (config, view)) {
+    //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
+    //     return psphotReadoutCleanup (config, view);
+    // }
 
     // classify sources based on moments, brightness
+    // 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");
@@ -92 +101 @@
     }
     // 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");
@@ -121 +131 @@
     psphotMergeSources (config, view);
 
-    // *** generate the objects (which unify the sources from the different images)
-    // pmArray *objects = psphotMatchSources (config, view);
-    
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
-    // psphotFitSourcesLinearStack (config, objects, FALSE);
+    psphotFitSourcesLinearStack (config, objects, FALSE);
+    psFree (objects);
 
     // identify CRs and extended sources
@@ -153 +161 @@
     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, view);

branches/eam_branches/stackphot.20100406/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);