Changeset 32685 for branches/eam_branches/ipp-20111110/psphot/src

Timestamp:

Nov 17, 2011, 5:41:38 AM (15 years ago)

Author:

eugene

Message:

psphotStack is generating duplicate detections because the detection removal before the second pass only occurs on the SRC, not DET image; working on fixing this by modifying the removal process

Location:

branches/eam_branches/ipp-20111110/psphot/src

Files:

• psphot.h (modified) (2 diffs)
• psphotEfficiency.c (modified) (1 diff)
• psphotFake.c (modified) (1 diff)
• psphotOutput.c (modified) (2 diffs)
• psphotRadialApertures.c (modified) (1 diff)
• psphotReadout.c (modified) (9 diffs)
• psphotReplaceUnfit.c (modified) (10 diffs)
• psphotStackReadout.c (modified) (13 diffs)

branches/eam_branches/ipp-20111110/psphot/src/psphot.h

@@ -108 +108 @@
 bool            psphotBlendFit_Threaded (psThreadJob *job);
 
-bool            psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule);
-bool            psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+bool            psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState);
+bool            psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool ignoreState);
 
 bool            psphotAddNoise (pmConfig *config, const pmFPAview *view, const char *filerule);
@@ -176 +176 @@
 
 // in psphotReplaceUnfit.c:
-bool            psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe);
+bool            psphotRemoveAllSourcesByArray (const psArray *sources, const psMetadata *recipe);
+bool            psphotRemoveAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState);
+bool            psphotRemoveAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState);
 bool            psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal);
 

branches/eam_branches/ipp-20111110/psphot/src/psphotEfficiency.c

@@ -255 +255 @@
 
     // remove all sources, adding noise for subtracted sources
-    psphotRemoveAllSources(realSources, recipe);
+    psphotRemoveAllSourcesByArray(realSources, recipe);
 
 #if TESTING

branches/eam_branches/ipp-20111110/psphot/src/psphotFake.c

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

branches/eam_branches/ipp-20111110/psphot/src/psphotOutput.c

@@ -407 +407 @@
     return true;
 }
+
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule)
+{
+    static int npass = 0;
+    char filename[64];
+
+    // XXX uncomment to disreturn true;
+
+    bool status = true;
+
+    int num = psphotFileruleCount(config, filerule);
+
+    snprintf (filename, 64, "testdump.%02d.dat", npass);
+    FILE *f = fopen (filename, "w");
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+
+        // find the currently selected readout
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
+        psAssert (file, "missing file?");
+
+        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+        psAssert (readout, "missing readout?");
+
+        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+        psAssert (detections, "missing detections?");
+
+        psArray *sources = detections->newSources ? detections->newSources : detections->allSources;
+        psAssert (sources, "missing sources?");
+
+        if (detections->newSources) {
+            fprintf (f, "## --- from new sources ---\n");
+        } else {
+            fprintf (f, "## --- from all sources ---\n");
+        }
+
+        for (int i = 0; i < sources->n; i++) {
+            pmSource *source = sources->data[i];
+            if (!source) continue;
+
+            pmPeak *peak = source->peak;
+            if (!peak) continue;
+
+            // XXX only dump a given region
+            if (peak->xf < 20) continue;
+            if (peak->yf < 20) continue;
+            if (peak->xf > 40) continue;
+            if (peak->yf > 70) continue;
+
+            float Msum = source->moments ? source->moments->Sum : NAN;
+            float Mx   = source->moments ? source->moments->Mx : NAN;
+            float My   = source->moments ? source->moments->My : NAN;
+            float Npix = source->moments ? source->moments->nPixels : NAN;
+            fprintf (f, "%f %f  : %f %f : %f %f\n", peak->xf, peak->yf, Mx, My, Msum, Npix);
+        }
+    }
+    fclose (f);
+    npass ++;
+
+    return true;
+}
+
+# if (0)
 bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule, char *filename) {
 
@@ -449 +514 @@
     return true;
 }
-
+# endif

branches/eam_branches/ipp-20111110/psphot/src/psphotRadialApertures.c

@@ -156 +156 @@
 
             // set this to 0 to run without threading
-# if (1)            
+# if (0)            
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");

branches/eam_branches/ipp-20111110/psphot/src/psphotReadout.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule);
 
 // this should be called by every program that links against libpsphot
@@ -9 +10 @@
 }
 
-// for now, let's store the detections on the readout->analysis for each readout
-bool psphotDumpChisqs (pmConfig *config, const pmFPAview *view, const char *filerule)
-{
-    static int npass = 0;
-    char filename[64];
-
-    return true;
-
-    bool status = true;
-
-    int num = psphotFileruleCount(config, filerule);
-
-    snprintf (filename, 64, "chisq.%02d.dat", npass);
-    FILE *f = fopen (filename, "w");
-
-    // loop over the available readouts
-    for (int i = 0; i < num; i++) {
-
-        // find the currently selected readout
-        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
-        psAssert (file, "missing file?");
-
-        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
-        psAssert (readout, "missing readout?");
-
-        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
-        psAssert (detections, "missing detections?");
-
-        psArray *sources = detections->allSources;
-        psAssert (sources, "missing sources?");
-
-        for (int i = 0; i < sources->n; i++) {
-            pmSource *source = sources->data[i];
-            if (!source) continue;
-
-            pmModel *model = pmSourceGetModel (NULL, source);
-            if (!model) continue;
-        
-            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
-                fprintf (f, "%f %f %f %d %d %f  1 NONLINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm);
-            } else {
-                fprintf (f, "%f %f %f %d %d %f  0 LINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm);
-            }
-        }
-    }
-    fclose (f);
-    npass ++;
-
-    return true;
-}
-
 bool psphotReadout(pmConfig *config, const pmFPAview *view, const char *filerule) {
 
@@ -133 +83 @@
         return psphotReadoutCleanup (config, view, filerule);
     }
+    psphotDumpTest (config, view, filerule);
 
     // find blended neighbors of very saturated stars (detections->newSources)
@@ -188 +139 @@
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
     psphotFitSourcesLinear (config, view, filerule, false); // pass 1 (detections->allSources)
-    psphotDumpChisqs (config, view, filerule);
+    psphotDumpTest (config, view, filerule);
 
     // measure the radial profiles to the sky
@@ -208 +159 @@
     // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
     psphotBlendFit (config, view, filerule); // pass 1 (detections->allSources)
-    psphotDumpChisqs (config, view, filerule);
 
     // replace all sources
-    psphotReplaceAllSources (config, view, filerule); // pass 1 (detections->allSources)
+    psphotReplaceAllSources (config, view, filerule, false); // pass 1 (detections->allSources)
 
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
     psphotFitSourcesLinear (config, view, filerule, true); // pass 2 (detections->allSources)
-    psphotDumpChisqs (config, view, filerule);
 
     // if we only do one pass, skip to extended source analysis
@@ -239 +188 @@
         // NOTE: new sources are saved on detections->newSources
         psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources)
+        psphotDumpTest (config, view, filerule);
 
         // set source type
@@ -253 +203 @@
         // replace all sources so fit below applies to all at once
         // NOTE: apply only to OLD sources (which have been subtracted)
-        psphotReplaceAllSources (config, view, filerule); // pass 2
+        psphotReplaceAllSources (config, view, filerule, false); // pass 2
 
         // merge the newly selected sources into the existing list
@@ -259 +209 @@
         // XXX check on free of sources...
         psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources)
+        psphotDumpTest (config, view, filerule);
 
         // NOTE: apply to ALL sources
         psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources)
-        psphotDumpChisqs (config, view, filerule);
     }
 
@@ -295 +245 @@
         // replace all sources so fit below applies to all at once
         // NOTE: apply only to OLD sources (which have been subtracted)
-        psphotReplaceAllSources (config, view, filerule); // pass 2
+        psphotReplaceAllSources (config, view, filerule, false); // pass 2
 
         // merge the newly selected sources into the existing list

branches/eam_branches/ipp-20111110/psphot/src/psphotReplaceUnfit.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+
+static int replace_pass = 0;
+static int remove_pass = 0;
 
 // replace the flux for sources which failed
@@ -9 +12 @@
 
     for (int i = 0; i < sources->n; i++) {
-      source = sources->data[i];
-
-      // replace other sources?
-      if (source->mode & PM_SOURCE_MODE_FAIL) goto replace;
-      continue;
+        source = sources->data[i];
+
+        // replace other sources?
+        if (source->mode & PM_SOURCE_MODE_FAIL) goto replace;
+        continue;
 
     replace:
@@ -23 +26 @@
 
 // for now, let's store the detections on the readout->analysis for each readout
-bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule)
+bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState)
 {
     bool status = true;
@@ -35 +38 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (!psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe)) {
+        if (!psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe, ignoreState)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i);
             return false;
         }
     }
-    return true;
-}
-
-bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
+    replace_pass ++;
+    return true;
+}
+
+bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState) {
 
     bool status;
-    pmSource *source;
 
     psTimerStart ("psphot.replace");
@@ -57 +60 @@
     psAssert (readout, "missing readout?");
 
+    char name[128];
+    snprintf (name, 128, "testadd.neg.%02d.%02d.fits", replace_pass, index);
+    psphotSaveImage (NULL, readout->image, name);
+
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     psAssert (detections, "missing detections?");
@@ -75 +82 @@
 
     for (int i = 0; i < sources->n; i++) {
-      source = sources->data[i];
-
-      // replace other sources?
-      if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
-
-      pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-    }
+        pmSource *source = sources->data[i];
+
+        if (ignoreState) {
+            // rely on the type of source to decide if we subtract it or not
+
+            // skip non-astronomical objects (very likely defects)
+            if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+            if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+      
+            // do not include CRs in the full ensemble fit
+            if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
+        
+            // do not include MOMENTS_FAILURES in the fit
+            if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+        } else {
+            // if we respect the state, only remove unsubtracted sources
+            if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
+        }
+
+        pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+    }
+
+    snprintf (name, 128, "testadd.pos.%02d.%02d.fits", replace_pass, index);
+    psphotSaveImage (NULL, readout->image, name);
 
     psphotVisualShowImage(readout);
@@ -88 +112 @@
 }
 
-bool psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe) {
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotRemoveAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (!psphotRemoveAllSourcesReadout (config, view, filerule, i, recipe, ignoreState)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i);
+            return false;
+        }
+    }
+    remove_pass ++;
+    return true;
+}
+
+bool psphotRemoveAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState) {
+
+    bool status;
+
+    psTimerStart ("psphot.replace");
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    char name[128];
+    snprintf (name, 128, "testsub.pos.%02d.%02d.fits", remove_pass, index);
+    psphotSaveImage (NULL, readout->image, name);
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->allSources;
+    psAssert (sources, "missing sources?");
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    psAssert (maskVal, "missing mask value?");
+
+    // bit-mask to mark pixels not used in analysis
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
+    assert (markVal);
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+
+        // skip non-astronomical objects (very likely defects)
+        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+        
+        // do not include CRs in the full ensemble fit
+        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
+        
+        // do not include MOMENTS_FAILURES in the fit
+        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+        
+        // if we respect the state, only remove unsubtracted sources
+        if (!ignoreState && (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
+        
+        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+    }
+
+    snprintf (name, 128, "testsub.neg.%02d.%02d.fits", remove_pass, index);
+    psphotSaveImage (NULL, readout->image, name);
+
+    psphotVisualShowImage(readout);
+    psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
+    return true;
+}
+
+bool psphotRemoveAllSourcesByArray (const psArray *sources, const psMetadata *recipe) {
 
     bool status;
@@ -100 +207 @@
 
     for (int i = 0; i < sources->n; i++) {
-      source = sources->data[i];
-
-      // replace other sources?
-      if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
-
-      pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+        source = sources->data[i];
+
+        // replace other sources?
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
+
+        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     }
     psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
@@ -158 +265 @@
 
     for (int i = 0; i < sources->n; i++) {
-      source = sources->data[i];
-
-      // sources have not yet been subtracted in this image (but this flag may be raised)
-      source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
-      if (!source->modelPSF) continue;
-
-      float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
-      float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
-      float radius = source->modelPSF->fitRadius;
-
-      // force a redefine to this image
-      pmSourceFreePixels(source);
-      pmSourceRedefinePixels (source, readout, Xo, Yo, radius);
+        source = sources->data[i];
+
+        // sources have not yet been subtracted in this image (but this flag may be raised)
+        source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+        if (!source->modelPSF) continue;
+
+        float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+        float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+        float radius = source->modelPSF->fitRadius;
+
+        // force a redefine to this image
+        pmSourceFreePixels(source);
+        pmSourceRedefinePixels (source, readout, Xo, Yo, radius);
     }
     return true;
@@ -229 +336 @@
 
     for (int i = 0; i < sources->n; i++) {
-      source = sources->data[i];
-
-      // *** we need to cache the 'best' model, and we have 3 cases:
-      // 1) model is the psf model --> generate from the new psf
-      // 2) model is an unconvolved extended model --> just cache the copy (not perfect)
-      // 3) model is a convolved extended model --> re-generate
-
-      // use the 'best' model to cache the model (PSF or EXT : EXT may point at one of modelFits
-      bool isPSF = false;
-      pmModel *model = pmSourceGetModel(&isPSF, source);
-      if (!model) continue;
-
-      float radius = model->fitRadius; // save for future use below
-
-      // regenerate the PSF if the model is a PSF, or if we need the PSF for a PCM
-      if (isPSF || model->isPCM) {
-          // the guess central intensity comes from the peak:
-          float Io = source->peak->rawFlux;
-          float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
-          float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
-
-          // generate a model for this object with Io = 1.0
-          pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io);
-          if (modelPSF == NULL) {
-              psWarning ("Failed to determine PSF model for source at (%f,%f), skipping", Xo, Yo);
-              continue;
-          }
-
-          // set the source PSF model
-          psFree (source->modelPSF);
-          source->modelPSF = modelPSF;
-          source->modelPSF->fitRadius = radius;
-      }
-
-      if (model->isPCM) {
-          psAssert(false, "this section is not complete");
-
-          pmSourceCachePSF (source, maskVal);
-
-          psKernel *psfKernel = pmPCMkernelFromPSF(source, psfSize);
-          if (!psfKernel) { 
-              psWarning ("no psf kernel");
-          }
-
-          // generate an image of the right size
-          psImage *rawModelFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
-          psImageInit (rawModelFlux, 0.0);
+        source = sources->data[i];
+
+        // *** we need to cache the 'best' model, and we have 3 cases:
+        // 1) model is the psf model --> generate from the new psf
+        // 2) model is an unconvolved extended model --> just cache the copy (not perfect)
+        // 3) model is a convolved extended model --> re-generate
+
+        // use the 'best' model to cache the model (PSF or EXT : EXT may point at one of modelFits
+        bool isPSF = false;
+        pmModel *model = pmSourceGetModel(&isPSF, source);
+        if (!model) continue;
+
+        float radius = model->fitRadius; // save for future use below
+
+        // regenerate the PSF if the model is a PSF, or if we need the PSF for a PCM
+        if (isPSF || model->isPCM) {
+            // the guess central intensity comes from the peak:
+            float Io = source->peak->rawFlux;
+            float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+            float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+
+            // generate a model for this object with Io = 1.0
+            pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io);
+            if (modelPSF == NULL) {
+                psWarning ("Failed to determine PSF model for source at (%f,%f), skipping", Xo, Yo);
+                continue;
+            }
+
+            // set the source PSF model
+            psFree (source->modelPSF);
+            source->modelPSF = modelPSF;
+            source->modelPSF->fitRadius = radius;
+        }
+
+        if (model->isPCM) {
+            psAssert(false, "this section is not complete");
+
+            pmSourceCachePSF (source, maskVal);
+
+            psKernel *psfKernel = pmPCMkernelFromPSF(source, psfSize);
+            if (!psfKernel) { 
+                psWarning ("no psf kernel");
+            }
+
+            // generate an image of the right size
+            psImage *rawModelFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
+            psImageInit (rawModelFlux, 0.0);
           
-          // insert the model image normalized to 1.0
-          pmModelAdd (rawModelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
-
-          psImageConvolveFFT (source->modelFlux, rawModelFlux, NULL, 0, psfKernel);
+            // insert the model image normalized to 1.0
+            pmModelAdd (rawModelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
+
+            psImageConvolveFFT (source->modelFlux, rawModelFlux, NULL, 0, psfKernel);
           
-          psFree (psfKernel);
-          psFree (rawModelFlux);
-      }
-
-      pmSourceCacheModel (source, maskVal);  // ALLOC x14 (!)
+            psFree (psfKernel);
+            psFree (rawModelFlux);
+        }
+
+        pmSourceCacheModel (source, maskVal);  // ALLOC x14 (!)
     }
 

branches/eam_branches/ipp-20111110/psphot/src/psphotStackReadout.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule);
 
 // we have 3 possible real filesets:
@@ -130 +131 @@
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
+    psphotDumpTest (config, view, STACK_SRC);
     psMemDump("sourcestats");
 
@@ -170 +172 @@
     psphotFitSourcesLinear (config, view, STACK_SRC, false);
     psphotStackVisualFilerule(config, view, STACK_SRC);
+    psphotDumpTest (config, view, STACK_SRC);
 
     // measure the radial profiles to the sky
@@ -186 +189 @@
     psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)
 
-    // replace all sources
-    psphotReplaceAllSources (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    // replace all sources (do NOT ignore subtraction state)
+    psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
 
     // if we only do one pass, skip to extended source analysis
@@ -194 +197 @@
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
+    // NOTE 2 : this function subtracts the models from the given filerule (SRC)
     psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources)
 
@@ -200 +204 @@
     // NOTE: this block performs the 2nd pass low-significance PSF detection stage
     { 
+        //  subtract all sources from DET (this will subtract using the psf model for SRC, which
+        //  will somewhat oversubtract the sources -- this is OK
+        psphotRemoveAllSources (config, view, STACK_DET, true); // ignore subtraction state for sources
+
         // add noise for subtracted objects
         psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
@@ -212 +220 @@
 
         // if DET and SRC are different images, copy the detections from DET to SRC 
+        // (this operation just ensures the metadata container has a view on SRC as well
         if (strcmp(STACK_SRC, STACK_DET)) {
-            // XXX how does this handle 1st vs 2nd pass sources?
             if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
                 psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
@@ -220 +228 @@
         }
 
+        // replace all sources in DET
+        psphotReplaceAllSources (config, view, STACK_DET, true); // ignore subtraction state for sources
+
         // define new sources based on only the new peaks & measure moments
         // NOTE: new sources are saved on detections->newSources
         psphotSourceStats (config, view, STACK_SRC, false); // pass 2 (detections->newSources)
+        psphotDumpTest (config, view, STACK_SRC);
 
         // set source type
@@ -237 +249 @@
         // replace all sources so fit below applies to all at once
         // NOTE: apply only to OLD sources (which have been subtracted)
-        psphotReplaceAllSources (config, view, STACK_SRC); // pass 2
+        psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 2
 
         // merge the newly selected sources into the existing list
@@ -243 +255 @@
         // XXX check on free of sources...
         psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
+        psphotDumpTest (config, view, STACK_SRC);
     }
 
@@ -256 +269 @@
 
     psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
+    psphotDumpTest (config, view, STACK_SRC);
 
     // NOTE: apply to ALL sources
@@ -295 +309 @@
     int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
 
+    psphotDumpTest (config, view, STACK_SRC);
+
     for (int entry = 1; entry < nRadialEntries; entry++) {
         // NOTE: entry 0 is the unmatched image set
@@ -312 +328 @@
 
         // replace the flux in the image so it is returned to its original state
-        psphotReplaceAllSources (config, view, STACK_OUT);
+        psphotReplaceAllSources (config, view, STACK_OUT, false);
 
         // smooth to the next FWHM, or set 'smoothAgain' to false if no more