Changeset 34399

Timestamp:

Sep 5, 2012, 4:04:34 PM (14 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ipp-20120805/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/Makefile.am (modified) (2 diffs)
• src/psphot.h (modified) (4 diffs)
• src/psphotApResid.c (modified) (1 diff)
• src/psphotBlendFit.c (modified) (2 diffs)
• src/psphotEfficiency.c (modified) (6 diffs)
• src/psphotKronIterate.c (modified) (13 diffs)
• src/psphotMaskReadout.c (modified) (1 diff)
• src/psphotReadout.c (modified) (2 diffs)
• src/psphotReadoutMinimal.c (modified) (1 diff)
• src/psphotSetThreads.c (modified) (1 diff)
• src/psphotSourceMatch.c (modified) (13 diffs)
• src/psphotSourceMemory.c (copied) (copied from trunk/psphot/src/psphotSourceMemory.c )
• src/psphotSourceStats.c (modified) (2 diffs)
• src/psphotStackAllocateOutput.c (copied) (copied from trunk/psphot/src/psphotStackAllocateOutput.c )
• src/psphotStackImageLoop.c (modified) (3 diffs, 1 prop)
• src/psphotStackMatchPSFs.c (modified) (5 diffs)
• src/psphotStackMatchPSFsNext.c (modified) (2 diffs)
• src/psphotStackMatchPSFsUtils.c (modified) (4 diffs)
• src/psphotStackReadout.c (modified) (18 diffs)

branches/eam_branches/ipp-20120805/psphot/src/Makefile.am

@@ -117 +117 @@
         psphotStackObjects.c          \
         psphotStackPSF.c              \
+        psphotStackAllocateOutput.c   \
         psphotCleanup.c
 
@@ -219 +220 @@
         psphotPetrosianVisual.c        \
         psphotEfficiency.c             \
-        psphotSetNFrames.c
+        psphotSetNFrames.c             \
+        psphotSourceMemory.c
 
 # not currently used

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

@@ -364 +364 @@
                                      pmReadout **chiReadout,
                                      int index);
+bool psphotStackAllocateOutput( const pmConfig *config, pmFPAview *view, psMetadata *recipe);
 
 bool psphotStackRemoveChisqFromInputs (pmConfig *config, const char *filerule);
@@ -371 +372 @@
 bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index);
 bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS);
+bool psphotDropBadMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
 
 bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final);
@@ -487 +489 @@
 bool psphotMaskFootprint (pmReadout *readout, pmSource *source, psImageMaskType markVal);
 
-bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule);
-bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index);
+bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule, int pass);
+bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass);
 bool psphotKronIterate_Threaded (psThreadJob *job);
 
@@ -520 +522 @@
 bool psphotSatstarPhotometry (pmSource *source);
 
+bool psphotSourceMemory(pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotSourceMemoryReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index);
+
 #endif

branches/eam_branches/ipp-20120805/psphot/src/psphotApResid.c

@@ -207 +207 @@
         if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY");
         if (source->mode &  PM_SOURCE_MODE_DEFECT) SKIPSTAR ("DEFECT");
+        if (source->mode2 & PM_SOURCE_MODE2_MATCHED) SKIPSTAR ("MATCHED");
 
         if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {

branches/eam_branches/ipp-20120805/psphot/src/psphotBlendFit.c

@@ -94 +94 @@
     assert (status && isfinite(skySig) && skySig > 0);
 
+# if (0)
     // **** test block : generate an ID image where pixels are set based the source models (flux > 0.1 peak && flux > 2.0 skySig)
     psImage *IDimage = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
@@ -104 +105 @@
     }
     psphotSaveImage (NULL, IDimage, "idimage.fits");
+# endif
 
     // Define source fitting parameters for extended source fits

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

@@ -57 +57 @@
     psFree(stats);
 
+#ifndef USE_SINGLE_POINT_FOR_MODEL
     // Need to normalise out difference between Gaussian and real PSF
+    int sizeX = ro->variance->numCols / 16;
+    int sizeY = ro->variance->numRows / 16;
+    int numPoints = 0;
+    float sum =  0;
+    for (int y = sizeY * 0.5 ; y < ro->variance->numRows; y += sizeY) {
+        for (int x = sizeX * 0.5 ; x < ro->variance->numCols; x += sizeX) {
+            pmModel *normModel = pmModelFromPSFforXY(psf, (float) x, (float) y, 1.0); // model for normalization
+
+            if (!normModel || (normModel->flags & MODEL_MASK)) {
+                psFree(normModel);
+                continue;
+            }
+            float flux = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0
+            psFree(normModel);
+            if (!isfinite(flux)) {
+                continue;
+            }
+            numPoints++;
+            sum += flux;
+        }
+    }
+    if (!isfinite(sum) || numPoints == 0) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate PSF model for any of %d points.", sizeX * sizeY);
+        return false;
+    }
+    *norm = sum / numPoints;
+#else
+    // This can fail for readout that has few good pixels. It's better to sample many points.
     pmModel *normModel = pmModelFromPSFforXY(psf, 0.5 * ro->variance->numCols,
-                                             0.5 * ro->variance->numRows, 1.0); // Model for normalisation
+                                         0.5 * ro->variance->numRows, 1.0); // Model for normalisation
+    psFree(normModel);
     if (!normModel || (normModel->flags & MODEL_MASK)) {
         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate PSF model.");
@@ -67 +97 @@
     *norm = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0
     psFree(normModel);
+#endif
 
     // The signal-to-noise of the smoothed image is: S/N ~ I_smooth / sqrt(variance * factor)
@@ -188 +219 @@
         if (!psphotEfficiencyReadout (config, view, filerule, i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure detection efficiency for %s entry %d", filerule, i);
-            return false;
+            psErrorStackPrint(stderr, " ");
+            psErrorClear();
+            // send message to log as well
+            psLogMsg ("psphot", PS_LOG_WARN, "failed to measure detection efficiency for %s entry %d", filerule, i);
         }
     }
@@ -376 +410 @@
         float mag = magLim + magOffsets->data.F32[i]; // Magnitude for bin
         float peak = powf(10.0, -0.4 * mag) / norm;   // Peak flux
+#ifndef USE_SINGLE_POINT_FOR_MODEL
+        int sizeX = numCols / 16;
+        int sizeY = numRows / 16;
+        int numPoints = 0;
+        float sum =  0;
+        for (int y = sizeY * 0.5 ; y < numRows; y += sizeY) {
+            for (int x = sizeX * 0.5 ; x < numCols; x += sizeX) {
+                // Need to normalise out difference between Gaussian and real PSF
+                pmModel *model = pmModelFromPSFforXY(psf, (float) x, (float) y, peak); 
+
+                if (!model || (model->flags & MODEL_MASK)) {
+                    psFree(model);
+                    continue;
+                }
+                float sourceRadius = PS_MAX(radius, model->modelRadius(model->params, minFlux)); // Radius for source
+                psFree(model);
+                if (!isfinite(sourceRadius)) {
+                    continue;
+                }
+                numPoints++;
+                sum += sourceRadius;
+            }
+        }
+        if (!isfinite(sum) || numPoints == 0) {
+            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate PSF model for any of %d points.", sizeX * sizeY);
+            return false;
+        }
+        float sourceRadius = sum / numPoints;
+#else
+        // This old way can fail for heavily masked images
         pmModel *model = pmModelFromPSFforXY(psf, numCols / 2.0, numRows / 2.0, peak); // Model for source
         if (!model || (model->flags & MODEL_MASK)) {
@@ -384 +448 @@
         float sourceRadius = PS_MAX(radius, model->modelRadius(model->params, minFlux)); // Radius for source
         psFree(model);
+#endif
 
         psArray *sources = psArrayAllocEmpty(numSources); // Sources in this bin
@@ -401 +466 @@
                 source->peak->yf = y;
                 source->modelPSF = pmModelFromPSFforXY(psf, x, y, 2.0 * sqrtf(sig));
-                source->type = PM_SOURCE_TYPE_STAR;
-
-                source->modelPSF->fitRadius = sourceRadius;
-                source->apRadius = sourceRadius;
-
-                numFound++;
-                psArrayAdd(sources, sources->n, source);
-                psArrayAdd(fakeSourcesAll, fakeSourcesAll->n, source);
+                if (source->modelPSF) {
+                    source->type = PM_SOURCE_TYPE_STAR;
+
+                    source->modelPSF->fitRadius = sourceRadius;
+                    source->apRadius = sourceRadius;
+
+                    numFound++;
+                    psArrayAdd(sources, sources->n, source);
+                    psArrayAdd(fakeSourcesAll, fakeSourcesAll->n, source);
+                } else {
+                    psLogMsg("psphot", PS_LOG_WARN, "pmModelFromPSFforXY failed for sig: %f x: %6.1f y: %6.1f\n", 
+                        sig, x, y);
+                }
                 psFree(source);
             }

branches/eam_branches/ipp-20120805/psphot/src/psphotKronIterate.c

@@ -5 +5 @@
 
 
-bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule)
+bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule, int pass)
 {
     bool status = true;
@@ -42 +42 @@
         // psAssert (psf, "missing psf?");
 
-        if (!psphotKronIterateReadout (config, recipe, view, filerule, readout, sources, psf, i)) {
+        if (!psphotKronIterateReadout (config, recipe, view, filerule, readout, sources, psf, i, pass)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
             return false;
@@ -54 +54 @@
 bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
 
-bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index) {
+#ifdef DUMP_KRS
+FILE *dumpFile = NULL;
+#endif
+
+bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) {
 
     bool status = false;
+
+#ifdef DUMP_KRS
+    if (!dumpFile) {
+        dumpFile = fopen("kr.txt", "w");
+        psAssert (dumpFile, "failed to open kr.txt");
+    }
+    fprintf(dumpFile, "\n\n Input %d\n", index);
+#endif
 
     if (!sources->n) {
@@ -64 +76 @@
 
     psTimerStart ("psphot.kron");
-
 
     // determine the number of allowed threads
@@ -160 +171 @@
 
     // generate the window image: multiply the flux by this to downweight neighbors
+    // XXX: we don't need this image if we aren't going to apply the window
     psImage *kronWindow = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
     psImageInit (kronWindow, 1.0);
@@ -235 +247 @@
             PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA,  PS_TYPE_F32);
             PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, pass,               PS_TYPE_S32);
 
 // set this to 0 to run without threading
@@ -297 +310 @@
     float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[11],F32);
     float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[12],F32);
+    int pass                        = PS_SCALAR_VALUE(job->args->data[13],S32);
+#ifndef REVERT_ON_BAD_MEASUREMENT
+    (void) pass;
+#endif
 
     for (int j = 0; j < KRON_ITERATIONS; j++) {
@@ -308 +325 @@
             if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
             if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+            if (!isfinite(source->moments->Mrf)) {
+                // Once we save a bad Mrf measurement we give up on this source
+                // checking here allows us to avoid adding and subtracting the model
+                continue;
+            }
 
             // skip saturated stars modeled with a radial profile 
@@ -324 +346 @@
             }
 
-            // On first iteration set window radius to sky radius (if valid) on second iteration
-            // use a factor times the previous radial moment value up to a maximum value that
-            // depends on the surface brightness of the source
-            float maxWindow;
-            if (j == 0) {
-                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
-            } else {
+            // On first iteration set window radius to sky radius (if valid). We also use this on subsequent
+            // iterations if we cannot find a better limit
+            float maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
+            if (j > 0) {
+                // on subsequent iterations we use a factor times the previous radial moment value
+                // limited to a maximum value that depends on the surface brightness of the source
                 if (KRON_SB_MIN_DIVISOR) {
-                    if (isfinite(source->moments->KronFlux) && (source->moments->KronFlux > 0)) {
-                        // Limit window radius based on surface brightness
+                    // Limit window radius based on surface brightness if we have a good measurement of kron flux
+                    if (isfinite(source->moments->KronFlux)  && (source->moments->KronFlux > 0)) {
                         float Rmax = sqrt(source->moments->KronFlux) / KRON_SB_MIN_DIVISOR;
 
-                        if (source->moments->Mrf > 0) {
+                        if (isfinite(source->moments->Mrf) && source->moments->Mrf > 0) {
                             maxWindow = PS_MIN(6.0*source->moments->Mrf, Rmax);
                         } else {
-                            maxWindow = Rmax;
+                            maxWindow = PS_MIN(Rmax, maxWindow);
                         }
-                    } else {
-                        maxWindow = RADIUS;
                     }
                 } else {
-                    // old code
+                    // old recipe, no surface brightness cut
                     maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS;
                 }
             }
             float windowRadius = PS_MAX(RADIUS, maxWindow);
+
+#ifdef REVERT_ON_BAD_MEASURMENT
+            // save previous measurements. We might revert back to them if this round fails
+            float MrfPrior = source->moments->Mrf;
+            float KronFluxPrior = source->moments->KronFlux;
+            float KronFluxErrPrior = source->moments->KronFluxErr;
+#endif
 
             // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
@@ -376 +402 @@
             }
 
+#ifdef REVERT_ON_BAD_MEASUREMENT
+            // on pass 2 if we get an invalid measurement on a pass 1 source where we had a good one previously
+            // in pass 1 keep that measurement
+            bool reverted = false;
+            if (pass > 1 && !isfinite(source->moments->Mrf) && (source->mode2 & PM_SOURCE_MODE2_PASS1_SRC)) {
+                source->moments->Mrf = MrfPrior;    // This is finite otherwise we wouldn't have gotten here
+                source->moments->KronFlux = KronFluxPrior;
+                source->moments->KronFluxErr = KronFluxErrPrior;
+                reverted = true;
+            }
+#endif
+#ifdef DUMP_KRS
+#ifndef REVERT_ON_BAD_MEASUREMENT
+            bool reverted = false;
+#endif
+            fprintf(dumpFile, "%7d %1d %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %2d\n", source->id, reverted, source->moments->Mrf, MrfPrior, maxWindow, windowRadius, source->peak->xf, source->peak->yf, source->imageID);
+#endif
+
             // if we subtracted it above, re-subtract the object, leave local sky
             if (reSubtract) {
@@ -494 +538 @@
     }
 
+    float MrfTry = RF/RS;
+    if (RF <= 0. || RS <= 0 || !isfinite(MrfTry)) {
+        // We did not get a good measurement
+        source->moments->Mrf = NAN;
+        source->moments->KronFlux  = NAN;
+        source->moments->KronFluxErr  = NAN;
+        return false;
+    }
+
+    float Mrf = MAX(minKronRadius, MrfTry);
     // Saturate the 1st radial moment
-    float Mrf = MAX(minKronRadius, RF/RS);
     if (sqrt(source->peak->detValue) < 10.0) {
         Mrf = MIN (radius, Mrf);
@@ -544 +597 @@
     }
 
-    source->moments->Mrf = Mrf;
+    source->moments->Mrf         = Mrf;
     source->moments->KronFlux    = Sum;
     source->moments->KronFluxErr = sqrt(Var);
@@ -562 +615 @@
     psAssert(kronWindow, "need a window");
 
+    // XXX: If we are not applying the window then we don't need to check for valid Mrf here.
+    // We should give the this module a chance to measure a good value. 
+    // However experiments show that it hardly ever succeeds in getting a better value
     if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;
 

branches/eam_branches/ipp-20120805/psphot/src/psphotMaskReadout.c

@@ -58 +58 @@
             return false;
         }
+    }
+
+    // insure than any non-finite pixels in image or variance are masked
+    // get the PSPHOT.MASK value from the config
+    psImageMaskType maskValue;
+    if (!pmConfigMaskSetBits (&maskValue, NULL, config)) {
+        psError (PS_ERR_UNKNOWN, false, "Unable to define the mask bit values");
+        return false;
+    }
+    if (!pmReadoutMaskInvalid(readout, maskValue, maskSat)) {
+        psError (PS_ERR_UNKNOWN, false,  "Unable to mask invalid pixels in readout.");
+        return false;
     }
 

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

@@ -194 +194 @@
     // but this is chosen above to be appropriate for the PSF objects (not galaxies)
     // psphotKronMasked(config, view, filerule);
-    psphotKronIterate(config, view, filerule);
+    psphotKronIterate(config, view, filerule, 1);
 
     // identify CRs and extended sources (only unmeasured sources are measured)
@@ -310 +310 @@
     // re-measure the kron mags with models subtracted
     // psphotKronMasked(config, view, filerule);
-    psphotKronIterate(config, view, filerule);
+    psphotKronIterate(config, view, filerule, 2);
 
     // measure source size for the remaining sources

branches/eam_branches/ipp-20120805/psphot/src/psphotReadoutMinimal.c

@@ -88 +88 @@
 
     // re-measure the kron mags with models subtracted and more appropriate windows
-    psphotKronIterate(config, view, filerule);
+    psphotKronIterate(config, view, filerule, 1);
 
     // measure source size for the remaining sources

branches/eam_branches/ipp-20120805/psphot/src/psphotSetThreads.c

@@ -30 +30 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 13);
+    task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 14);
     task->function = &psphotKronIterate_Threaded;
     psThreadTaskAdd(task);

branches/eam_branches/ipp-20120805/psphot/src/psphotSourceMatch.c

@@ -1 +1 @@
 # include "psphotInternal.h" 
 
-bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
+// Structure for storing the contents of the PSPHOT.STACK.MATCH.FILTERS list from recipe
+#define MATCH_INFO_FILTER_STR_LEN 16
+typedef struct {
+    int     inputNum;
+    int     order;
+    bool    matchAll;
+    char    filterID[MATCH_INFO_FILTER_STR_LEN];
+} psphotStackMatchInfo;
+
+static psphotStackMatchInfo * psphotStackGetMatchInfo (pmConfig *config, const pmFPAview *view, const char *filerule);
+
+// functions for sorting the match info array
+static int compareMatchInfoByInputNum(const void *a, const void *b);
+static int compareMatchInfoByOrder(const void *a, const void *b);
+
+bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects, psphotStackMatchInfo *matchInfo);
 bool psphotMatchSourcesSetIDs (psArray *objects);
 
@@ -14 +29 @@
     int num = psphotFileruleCount(config, filerule);
 
-    // loop over the available readouts
-    for (int i = 0; i < num; i++) {
+    psphotStackMatchInfo *matchInfo = psphotStackGetMatchInfo(config, view, filerule);
+    
+    // loop over the available readouts matching sources found to objects. The inputs are processed
+    // in the order specified by the recipe
+    for (int j = 0; j < num; j++) {
+        int i = matchInfo[j].inputNum;
         if (!psphotMatchSourcesReadout (objects, config, view, filerule, i)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for %s entry %d", filerule, i);
@@ -23 +42 @@
     }
 
+    // Now re-order the matchInfo array by input number so we can find each inputs entry easily
+    qsort(matchInfo, num, sizeof(psphotStackMatchInfo), compareMatchInfoByInputNum);
+
     // create sources for images where an object has been detected in the other images
-    psphotMatchSourcesAddMissing (config, view, filerule, objects);
+    psphotMatchSourcesAddMissing (config, view, filerule, objects, matchInfo);
 
     // choose a consistent position; set common sequence values
     psphotMatchSourcesSetIDs (objects);
 
+    psFree(matchInfo);
+
     return objects;
 }
 
-bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index) { 
+bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index) {
  
     bool status = false;
@@ -61 +85 @@
 # define NEXT1 { i++; continue; } 
 # define NEXT2 { j++; continue; } 
-bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS) { 
+bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS) {
  
     float dx, dy; 
@@ -70 +94 @@
     sources = psArraySort (sources, pmSourceSortByX); 
     objects = psArraySort (objects, pmPhotObjSortByX); 
- 
+
     psVector *foundSrc = psVectorAlloc(sources->n, PS_TYPE_U8);
     psVectorInit (foundSrc, 0);
@@ -87 +111 @@
         pmPhotObj *obj = objects->data[j]; 
  
-        if (!src) NEXT1; 
+        if (!src) NEXT1;
         if (!src->peak) NEXT1; 
         if (!isfinite(src->peak->xf)) NEXT1; 
         if (!isfinite(src->peak->yf)) NEXT1; 
- 
+
         if (!obj) NEXT2; 
         if (!isfinite(obj->x)) NEXT2; 
@@ -141 +165 @@
     } 
 
-    // add missed sources to new objects
-
+    // create new objects for unmatched sources
     for (i = 0; i < sources->n; i++) {
 
-        if (foundSrc->data.U8[i]) continue;
+        if (foundSrc->data.U8[i]) continue;
 
         pmSource *src = sources->data[i]; 
 
-        pmPhotObj *obj = pmPhotObjAlloc();
-        pmPhotObjAddSource(obj, src);
-        psArrayAdd (objects, 100, obj);
-        psFree (obj);
+        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);
@@ -161 +184 @@
 } 
 
-bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
+bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects, psphotStackMatchInfo *matchInfo) {
 
     bool status = false;
@@ -169 +192 @@
     psAssert (recipe, "missing recipe?");
 
+    int minDetectionsForForced = psMetadataLookupS32 (&status, recipe, "PSPHOT.STACK.MIN.DETECT.FOR.FORCED");
+    if (!status) {
+        minDetectionsForForced = 2;
+    }
+
     // determine properties (sky, moments) of initial sources
     float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
@@ -205 +233 @@
         int nSpansMax = 0;
         int iSpansMax = -1;
+
+        bool matchAll = false;
 
         // mark the images for which sources have been found
@@ -219 +249 @@
                 iSpansMax = j;
             }
+            // If this detection was on an input that has the matchAll flag set we create matched sources
+            // irespective of the number of inputs in which the object was found.
+            if (matchInfo[index].matchAll) {
+                matchAll = true;
+            }
 
             found->data.U8[index] = 1;
         }
+
+        // skip adding matched sources for this object if the number of detections for less than
+        // the supplied mininum unless one of the detections was in a band for which the matchAll flag is set
+        if (!matchAll && obj->sources->n < minDetectionsForForced) {
+            continue;
+        }
 
         // we make a copy of the largest footprint; this will be used for all new sources associated with this object
@@ -299 +340 @@
         pmPhotObj *obj = objects->data[i]; 
         nSources += obj->sources->n;
-        psAssert (obj->sources->n == nImages, "failed to match sources?");
+        if (minDetectionsForForced <= 1) {
+            psAssert (obj->sources->n == nImages, "failed to match sources?");
+        }
     }
     psLogMsg ("psphot", PS_LOG_DETAIL, "total of %d sources for %d images", nSources, nImages);
@@ -375 +418 @@
     return copy;
 }
+
+// qsort function to sort match info by order
+static int compareMatchInfoByOrder(const void *a, const void *b) {
+    psphotStackMatchInfo *infoA = (psphotStackMatchInfo *) a;
+    psphotStackMatchInfo *infoB = (psphotStackMatchInfo *) b;
+
+    int     result;
+    if (infoA->order < infoB->order) {
+        result = -1;
+    } else if (infoA->order == infoB->order) {
+        result = 0;
+    } else {
+        result = 1;
+    }
+    return result;
+}
+
+// qsort function to sort match info by input number
+static int compareMatchInfoByInputNum(const void *a, const void *b) {
+    psphotStackMatchInfo *infoA = (psphotStackMatchInfo *) a;
+    psphotStackMatchInfo *infoB = (psphotStackMatchInfo *) b;
+
+    int     result;
+    if (infoA->inputNum < infoB->inputNum) {
+        result = -1;
+    } else if (infoA->inputNum == infoB->inputNum) {
+        result = 0;
+    } else {
+        result = 1;
+    }
+    return result;
+}
+
+// psphotStackGetMatchInfo
+// process the recipe PSPHOT.STACK.MATCH.FILTERS which controls the order that the inputs
+// are processed for source matching and the rules for creating matched sources for sources that
+// are only detected in a single band.
+static psphotStackMatchInfo *psphotStackGetMatchInfo (pmConfig *config, const pmFPAview *view, const char *filerule) {
+
+    bool status = false;
+    psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+
+    psMetadata *filterInfo = psMetadataLookupPtr (&status, recipe, "PSPHOT.STACK.MATCH.FILTERS");
+    if (!status || !filterInfo) {   
+        psLogMsg ("psphot", PS_LOG_WARN, "PSPHOT.MATCH.STACK.FILTERS not found in the recipe. Will process inputs in order supplied.");
+        filterInfo = NULL;
+    }
+
+    int num = psphotFileruleCount(config, filerule);
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // Find the filter for each of the inputs in fpa concepts
+    psArray *inputFilters = psArrayAlloc(num);
+    for (int i = 0 ; i < num; i++) {
+        if (i != chisqNum) {
+            pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i);
+            psAssert (file, "missing file?");
+
+            psString filterid = psMetadataLookupStr(&status, file->fpa->concepts, "FPA.FILTERID");
+            psAssert (filterid, "missing FPA.FILTERID?");
+            psArraySet(inputFilters, i, filterid);
+        } else {
+            // The chisq image inherits its filterid from the first input so we shouldn't use that.
+            psString tmp  = psStringCopy("chisq");
+            psArraySet(inputFilters, i, tmp);
+            psFree(tmp);
+        }
+    }
+    // Allocate the match info array
+    psphotStackMatchInfo *matchInfo = psAlloc(num *sizeof(psphotStackMatchInfo));
+    memset(matchInfo, 0, num*sizeof(psphotStackMatchInfo));
+    int highest_order = -1;
+    if (filterInfo) {
+        // PSPHOT.STACK.MATCH.FILTERS is a metadata which contains a list of metadata objects each
+        // entry pertaining to a filter.
+        // The objects contained in each filter's metadata are strings.
+      
+        // Loop over the entries in the recipe and find the entry for each our our input readouts.
+        // XXX: Will this work if more than one input with a given filter is supplied?
+        // I think so.
+        psMetadataIterator *iter = psMetadataIteratorAlloc(filterInfo, PS_LIST_HEAD, NULL);
+        psMetadataItem *item = NULL;
+        while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
+            if (item->type != PS_DATA_METADATA) {
+                psAbort ("Invalid type for PSPHOT.STACK.MATCH.FILTER: %s, not a metadata folder", item->name);
+            }
+            psString thisFilter = psMetadataLookupStr (&status, item->data.md, "FILTER.ID");
+            psAssert(thisFilter, "missing FILTER.ID");
+
+            psString orderStr = psMetadataLookupStr (&status, item->data.md, "ORDER");
+            psAssert(orderStr, "missing ORDER");
+            psS32 order = atoi(orderStr);
+
+            psString matchAllStr = psMetadataLookupStr (&status, item->data.md, "MATCH.ALL");
+            psAssert(matchAllStr, "missing MATCH.ALL");
+            bool matchAll = (strcasecmp("true", matchAllStr) == 0);
+
+            // look for this filter in the inputs
+            for (int i = 0; i < num; i++) {
+                if (!strcmp((char *)inputFilters->data[i], thisFilter)) {
+                    // We have an input for this one
+                    matchInfo[i].inputNum = i;
+                    strncpy(matchInfo[i].filterID, thisFilter, MATCH_INFO_FILTER_STR_LEN - 1 );
+                    matchInfo[i].order = order;
+                    matchInfo[i].matchAll = matchAll;
+                    psLogMsg ("psphot", PS_LOG_DETAIL, "input: %d %s match order: %d match all: %d\n",
+                                                               i, thisFilter, order, matchAll);
+                    if (order > highest_order) {
+                        highest_order = order;
+                    }
+                    break;
+                }
+            }
+        }
+        psFree(iter);
+    }
+
+    // Make sure we have a matchInfo for all of the inputs. Fill in an entry for any that were not found.
+    for (int i = 0; i < num; i++) {
+        if (!*matchInfo[i].filterID) {
+            matchInfo[i].inputNum = i;
+            strncpy(matchInfo[i].filterID, inputFilters->data[i], MATCH_INFO_FILTER_STR_LEN - 1);
+            matchInfo[i].order = ++highest_order;
+            matchInfo[i].matchAll = false;
+            psLogMsg ("psphot", PS_LOG_WARN, "Entry in PSPHOT.MATCH.STACK.FILTERS not found for input: %d filter: %s using order %d",
+                i, (char *) inputFilters->data[i], highest_order);
+        }
+    }
+    psFree(inputFilters);
+
+    if (filterInfo) {
+        // Sort the array by ORDER.
+        qsort(matchInfo, num, sizeof(psphotStackMatchInfo), compareMatchInfoByOrder);
+    } else {
+        // no need to sort we just built the list in the order that the inputs were supplied
+    }
+
+    return matchInfo;
+}
+
+bool psphotDropBadMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
+
+    bool status = false;
+
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Drop Bad Matched Sources ---");
+
+    psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+
+    // select the appropriate recipe information
+    bool keepBadMatches = psMetadataLookupBool (&status, recipe, "PSPHOT.STACK.KEEP.BAD.MATCHES");
+    if (!status) {
+        keepBadMatches = true;
+    }
+
+    if (keepBadMatches) {
+        psLogMsg ("psphot", PS_LOG_INFO, "keeping bad matches");
+        return true;
+    }
+
+    int numImages = psphotFileruleCount(config, filerule);
+    psVector *dropped = psVectorAlloc(numImages, PS_TYPE_U32);
+    psVectorInit(dropped, 0);
+
+    int nDropped = 0;
+    for (int i = 0; i < objects->n; i++) { 
+        pmPhotObj *obj = objects->data[i]; 
+
+        // traverse the array from the end so that sources don't move until after we've processed them
+        for (int j = obj->sources->n - 1; j >= 0; j--) {
+
+            pmSource *source = obj->sources->data[j]; 
+            // This applies only to matched sources
+            if (!(source->mode2 & PM_SOURCE_MODE2_MATCHED)) continue;
+
+            if (isfinite(source->apFlux)) continue;
+
+            psTrace ("psphot", 7, "Dropping matched source from image %d at (%d, %d) no valid flux",
+                source->imageID, source->peak->x, source->peak->y);
+
+            psAssert(source->imageID >= 0 && source->imageID < numImages, "bad imageID %d", source->imageID);
+
+            dropped->data.U32[source->imageID]++;
+            
+            nDropped++;
+
+            psArrayRemoveIndex(obj->sources, j);
+        }
+    }
+
+    psLogMsg ("psphot", PS_LOG_DETAIL, "Dropped %d matched sources with no valid flux", nDropped);
+    psLogMsg ("psphot", PS_LOG_DETAIL, "       Input  Num Dropped");
+    for (int i=0; i<numImages; i++) {
+        psLogMsg ("psphot", PS_LOG_DETAIL, "    %8d     %8d", i, dropped->data.U32[i]);
+    }
+    psFree(dropped);
+
+    return true;
+}
+
+

branches/eam_branches/ipp-20120805/psphot/src/psphotSourceStats.c

@@ -524 +524 @@
     float PSF_SN_LIM = 2.0*psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
 
-    // XXX move this to a config file?
+    // XXX this will cause an error in the vector length is > 8
     # define NSIGMA 8
-    float sigma[NSIGMA] = {1.0, 2.0, 3.0, 4.5, 6.0, 9.0, 12.0, 18.0};
+    // moved to config file
+    psVector  *sigmavec = psMetadataLookupPtr (&status, recipe, "PSF.SIGMA.VALUES"); 
+
+    float sigma[NSIGMA];
     float Sout[NSIGMA];
     float Rmin[NSIGMA];
     int   Nout[NSIGMA]; // number of stars found in clump : use this to control the number of regions measured by psphotRoughClass
 
+    int nsigma;
+    if (sigmavec) {
+        psAssert(sigmavec->n <= NSIGMA, "too many sigma values in recipe %ld maximum is %d", sigmavec->n, NSIGMA);
+        // copy the data from vector to local array to keep the code below easier to read
+        for (int i = 0 ; i < sigmavec->n; i++) {
+            sigma[i] = sigmavec->data.F32[i];
+        }
+        assert (sigmavec->n <= 8);
+        nsigma = sigmavec->n;
+    } else {
+        // requiring this causes updates to pop an assertion
+        //  psAssert(status, "missing PSF.SIGMA.VALUES");
+        float defaultsigma[NSIGMA]  = {1.0, 2.0, 3.0, 4.5, 6.0, 9.0, 12.0, 18.0};
+        for (int i = 0 ; i < NSIGMA; i++) {
+            sigma[i] = defaultsigma[i];
+        }
+        nsigma = NSIGMA;
+    }
+
     // this sorts by peak->rawFlux
     sources = psArraySort (sources, pmSourceSortByFlux);
 
     // loop over radii:
-    for (int i = 0; i < NSIGMA; i++) {
+    for (int i = 0; i < nsigma; i++) {
 
         // XXX move max source number to config
@@ -592 +614 @@
     float minS = Sout[0];
     float maxS = Sout[0];
-    for (int i = 0; i < NSIGMA; i++) {
+    for (int i = 0; i < nsigma; i++) {
         minS = PS_MIN(Sout[i], minS);
         maxS = PS_MAX(Sout[i], maxS);
     }
-    if (minS > 0.65) Sigma = sigma[NSIGMA-1];
+    if (minS > 0.65) Sigma = sigma[nsigma-1];
     if (maxS < 0.65) Sigma = sigma[0];
 
-    for (int i = 0; i < NSIGMA - 1 && isnan(Sigma); i++) {
+    for (int i = 0; i < nsigma - 1 && isnan(Sigma); i++) {
         if (!isfinite(Sout[i]) || !isfinite(Sout[i+1])) continue;
         if ((Sout[i] > 0.65) && (Sout[i+1] > 0.65)) continue;

branches/eam_branches/ipp-20120805/psphot/src/psphotStackImageLoop.c

@@ -21 +21 @@
     psMemDump("startloop");
 
+    pmFPAview *view = pmFPAviewAlloc (0);
     pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
     pmFPAfile *inputCnv = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.CNV");
@@ -37 +38 @@
     }
 
-    pmFPAview *view = pmFPAviewAlloc (0);
-
-
-    // XXX for now, just load the full set of images up front except for EXPNUM which we defer
+    bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES");
+    bool match_psfs = radial_apertures;
+    bool needConvolved = radial_apertures || !useRaw;
+    if (!needConvolved) {
+        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.INPUT.CNV");
+    }
+
+    // just load the full set of images up front except for EXPNUM which we defer
     pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.RAW");
     pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.CNV");
@@ -63 +68 @@
                 psMemDump("load");
 
-                // Generate the 1st PSF-matched image set (larger target PSFs are generated by smoothing this image)
-                if (!psphotStackMatchPSFs (config, view)) {
-                    psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
-                    psFree (view);
-                    return false;
-                }
+                if (match_psfs) {
+                    // Generate the 1st PSF-matched image set (larger target PSFs are generated by smoothing this image)
+                    if (!psphotStackMatchPSFs (config, view)) {
+                        psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
+                        psFree (view);
+                        return false;
+                    }
+                } else {
+                    if (!psphotStackAllocateOutput (config, view, recipe)) {
+                        psError(psErrorCodeLast(), false, "failure in psphotStackAllocateOutput for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
+                        psFree (view);
+                        return false;
+                    }
+                }
                 psMemDump("stackmatch");
 

branches/eam_branches/ipp-20120805/psphot/src/psphotStackMatchPSFs.c

@@ -4 +4 @@
 {
     bool status = true;
+
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphotStack Match PSFs ---");
 
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PSPHOT"); 
@@ -57 +59 @@
 
 // convolve the image to match desired PSF
-// XXX is this code consistent with 'convolve' = false?
+// XXX is this code consistent with 'convolve' = false?  XXX: No
 bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) {
 
     psImageMaskType maskValue;
     psImageMaskType markValue;
+
+    psLogMsg("psphot", PS_LOG_DETAIL, "-- starting PSF matching for readout %d --", index);
 
     // get the PSPHOT.MASK value from the config
@@ -102 +106 @@
         saveMatchData(readoutOut, options, index);
     }
+    bool notMatched = psMetadataLookupBool(NULL, readoutOut->analysis, "NOT.PSF.MATCHED");
+    if (notMatched) {
+        if (!readoutOut->image) {
+            // our images got freed in the matching process. Recreate them so that they
+            // are available when sources are copied to this readout.
+            readoutOut->image = psImageCopy(NULL, readoutSrc->image, PS_TYPE_F32);
+            if (readoutSrc->variance) {
+                psFree(readoutOut->variance);
+                readoutOut->variance = psImageCopy(NULL, readoutSrc->variance, PS_TYPE_F32);
+            }
+            if (readoutSrc->mask) {
+                psFree(readoutOut->mask);
+                readoutOut->mask = psImageCopy(NULL, readoutSrc->mask, PS_TYPE_IMAGE_MASK);
+            }
+        }
+    }
 
     // renormalize the stack variances to have sigma = 1.0
@@ -111 +131 @@
     // save the output fwhm values in the readout->analysis.  we may have / will have multiple output PSF sizes,
     // so we save this in a vector.  if the vector is not yet defined, create it
-    // Skip this if the readout deconvolution fraction was over the limit.
+    // Skip this if psf matching failed for this readout. For example if deconvolution fraction was over the limit.
     // NOTE: fwhmValues as defined here has 1 + nMatched PSF : 0 == unmatched
     psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
     psVectorAppend(fwhmValues, NAN); // XXX this corresponds to the unmatched image set
 
-    bool overLimit = psMetadataLookupBool(NULL, readoutOut->analysis, "DECONV.OVERLIMIT");
-    if (!overLimit) {
+    if (!notMatched) {
         for (int i = 0; i < options->targetSeeing->n; i++) {
             psVectorAppend(fwhmValues, options->targetSeeing->data.F32[i]);
@@ -124 +143 @@
     psMetadataAddVector(readoutSrc->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
     psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
-    psFree(fwhmValues); // drops the extra copy
+    psFree(fwhmValues); // drops this function's reference
 
     return true;

branches/eam_branches/ipp-20120805/psphot/src/psphotStackMatchPSFsNext.c

@@ -81 +81 @@
     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     psAssert (maskVal, "missing mask value?");
+    psImageMaskType maskSat = psMetadataLookupImageMask(&status, recipe, "MASK.SAT"); // Mask value for bad pixels
+    psAssert (maskSat, "missing mask value?");
+
 
     float minGauss = psMetadataLookupF32(NULL, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel
@@ -123 +126 @@
     psImageSmoothMask_Threaded(readout->variance, readout->variance, readout->mask, maskVal, SIGMA_SMTH * M_SQRT1_2, NSIGMA_SMTH, minGauss);
     psLogMsg("psphot", PS_LOG_MINUTIA, "smooth variance: %f sec\n", psTimerMark("psphot.smooth"));
+
+    // Insure that invalid pixels are masked
+    // XXX: the smoothing seems to generate nan pixels in the variance image
+    // XXX: We may need to loop over the cached sources and redefine the maskObj images...
+    pmReadoutMaskInvalid(readout, maskVal, maskSat);
+    {
+        // Now go rebuild the sources' copies of the mask
+        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+        psAssert (detections, "missing detections?");
+
+        psArray *sources = detections->allSources;
+
+        if (sources && sources->n) {
+            for (int i = 0 ; i < sources->n; i++) {
+                // XXX: move this to a function in pmSource.c
+                pmSource *source = sources->data[i];
+                if (source->maskObj && source->maskView) {
+                    psFree(source->maskObj);
+                    source->maskObj = psImageCopy (source->maskObj, source->maskView, PS_TYPE_IMAGE_MASK);
+                }
+            }
+        }
+    }
 
     psLogMsg("psphot", PS_LOG_INFO, "smoothed");

branches/eam_branches/ipp-20120805/psphot/src/psphotStackMatchPSFsUtils.c

@@ -280 +280 @@
 
     // Do the image matching
+    bool rejectReadout = false; 
     pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readoutSrc->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel
     if (kernel) {
@@ -301 +302 @@
 
         if (!pmSubtractionMatch(NULL, readoutOut, fake, readoutSrc, footprint, stride, regionSize, spacing, threshold, stampSources, stampsName, type, size, order, widthsCopy, orders, inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) {
-            psError(psErrorCodeLast(), false, "Unable to match images.");
-            goto escape;
+            int errorCode = psErrorCodeLast();
+            if (errorCode == PM_ERR_SMALL_AREA) {
+                // failed to match but, don't fault. Just drop this input from measurements that need the
+                // matched readout.
+                rejectReadout = true;
+                psLogMsg("psphot", PS_LOG_WARN, "PSF match failed for input %d. Too few pixels.", index);
+            } else {
+                psError(psErrorCodeLast(), false, "Unable to match images.");
+                goto escape;
+            }
         }
     }
 
-    // Reject image completely if the maximum deconvolution fraction exceeds the limit
+    // If the maximum deconvolution fraction exceeds the limit, reject this input for analyses that require PSF matching
     float deconvLimit = psMetadataLookupF32(NULL, stackRecipe, "DECONV.LIMIT"); // Limit on deconvolution fraction
     float deconv = psMetadataLookupF32(NULL, readoutOut->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
     if (deconv > deconvLimit) {
-#if (1)
-        // XXX: don't reject the image set 
-        psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image for matched psf analysis%d\n", deconv, deconvLimit, index);
-        psMetadataAddBool(readoutOut->analysis, PS_LIST_TAIL, "DECONV.OVERLIMIT", PS_META_REPLACE, "", true);
-#else
-        psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index);
-        goto escape;
-#endif
+        rejectReadout = true;
+        psLogMsg("psphot", PS_LOG_WARN, "Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image for matched psf analysis%d\n", deconv, deconvLimit, index);
+    }
+    
+    if (rejectReadout) {
+        psMetadataAddBool(readoutOut->analysis, PS_LIST_TAIL, "NOT.PSF.MATCHED", PS_META_REPLACE, "", true);
     } else {
-        psMetadataAddBool(readoutOut->analysis, PS_LIST_TAIL, "DECONV.OVERLIMIT", PS_META_REPLACE, "", false);
-    }
-
+        psMetadataAddBool(readoutOut->analysis, PS_LIST_TAIL, "NOT.PSF.MATCHED", PS_META_REPLACE, "", false);
+    }
     // save the PSF on the new readout->analysis:
     // if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot psf model", options->psf)) {
@@ -470 +476 @@
 }
 
-// Set input to be skipped if the decovolution fraction was overlimit. Use for radial apertures
-// This interface can be potentiall be extended for other uses
+// Set input to be skipped if the analysis reports that psf matching has failed for some reason.
+// For exmaple if decovolution fraction was overlimit or the bad fraction was too high.
+//
 bool psphotStackSetInputsToSkip(pmConfig *config, const pmFPAview *view, const char *filerule, bool set) {
     int num = psphotFileruleCount(config, filerule);
@@ -487 +494 @@
         psAssert (readout, "missing readout?");
         if (set) {
-            bool overLimit = psMetadataLookupBool(&status, readout->analysis, "DECONV.OVERLIMIT");
-            psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSPHOT.SKIP.INPUT", PS_META_REPLACE, "Skip analysis", overLimit);
+            bool notMatched = psMetadataLookupBool(&status, readout->analysis, "NOT.PSF.MATCHED");
+            psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSPHOT.SKIP.INPUT", PS_META_REPLACE, "Skip analysis", notMatched);
         } else {
             psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSPHOT.SKIP.INPUT", PS_META_REPLACE, "Skip analysis", false);

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

@@ -2 +2 @@
 
 static bool psphotStackLoadWCS(pmConfig *config, const pmFPAview *view, const char *filerule);
+static void logMemStats(const char *heading);
 
 // we have 3 possible real filesets:
@@ -51 +52 @@
     // by the multiple threads, not the total time used by all threads.
     psTimerStart ("psphotReadout");
+
+    logMemStats("Start");
 
     pmModelClassSetLimits(PM_MODEL_LIMITS_LAX); // allow models to have ugly fits (eg, central cusp)
@@ -82 +85 @@
     }
 
-    // Generate the mask and weight images (if not supplied) and set mask bits
+    // Generate the mask and weight images (if not supplied) and set mask bits. 
+    // This also insures that all invalid pixels are masked (this is done for STACK_CNV in psphotStackMatchPSFs)
     if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+    if (!psphotSetMaskAndVariance (config, view, STACK_OUT)) {
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
@@ -114 +121 @@
     }
 
+#ifdef MAKE_CHISQ_IMAGE
     // also make the chisq detection image
     if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
@@ -119 +127 @@
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
+#endif
     if (!strcasecmp (breakPt, "CHISQ")) {
         return psphotReadoutCleanup (config, view, STACK_SRC);
@@ -151 +160 @@
     // psphotDumpTest (config, view, STACK_SRC);
     psMemDump("sourcestats");
+    logMemStats("sourcestats");
 
     // classify sources based on moments, brightness
@@ -196 +206 @@
     // window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
     // but iterates to an appropriately larger size
-    psphotKronIterate(config, view, STACK_SRC);
+        logMemStats("before.kron.1");
+    psphotKronIterate(config, view, STACK_SRC, 1);
+        logMemStats("after.kron.1");
 
     // identify CRs and extended sources
@@ -208 +220 @@
     psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
 
+    logMemStats("pass1");
 
     // if we only do one pass, skip to extended source analysis
@@ -299 +312 @@
     }
 
+    logMemStats("prematch");
+
     // generate the objects (objects unify the sources from the different images) NOTE: could
     // this just match the detections for the chisq image, and not bother measuring the source
@@ -331 +346 @@
     // re-measure the kron mags with models subtracted
     // psphotKronMasked(config, view, STACK_SRC);
-    psphotKronIterate(config, view, STACK_SRC);
+    logMemStats("before.kron.2");
+    psphotKronIterate(config, view, STACK_SRC, 2);
+    logMemStats("after.kron.2");
 
     // measure source size for the remaining sources
@@ -338 +355 @@
 
     psMemDump("psfstats");
+
+    // drop matched sources without any useful measurements
+    psphotDropBadMatchedSources (config, view, STACK_SRC, objects);
 
     // measure elliptical apertures, petrosians (objects sorted by S/N)
@@ -354 +374 @@
     }
 
-
     bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES");
     if (radial_apertures) {
@@ -361 +380 @@
 
         // NOTE: we always do the radial apertures analysis on the convolved image since
-        // those are the ones that are psf matched and the source of STACK_OUT's pixels
+        // those are the ones that are psf matched and are the source of STACK_OUT's pixels
         // XXX: Actually if PSPHOT.STACK.MATCH.PSF.SOURCE were set to RAW this wouldn't be true.
         // but in that case we don't get past the psf matching step because there is no
@@ -419 +438 @@
         // psphotEfficiency wants to have the PSF of the image, but since we are measuring on
         // the convolved images we need to generate PSFs for the DET images
-        if (!psphotChoosePSF (config, view, STACK_DET, false)) { // pass 1
+        if (!psphotChoosePSF (config, view, STACK_DET, false)) {
             psLogMsg ("psphot", 3, "failure to construct a psf model for raw input");
             return psphotReadoutCleanup (config, view, STACK_DET);
@@ -430 +449 @@
     psphotCopyEfficiency (config, view, STACK_OUT, STACK_DET);
 
+    logMemStats("final");
+#if (1)
+    psphotSourceMemory(config, view, STACK_SRC);
+    psphotSourceMemory(config, view, STACK_OUT);
+#endif
+
     // replace failed sources?
     // psphotReplaceUnfitSources (sources);
@@ -440 +465 @@
     psphotSourceFreePixels (config, view, STACK_SRC);
 
+#ifdef MAKE_CHISQ_IMAGE
     // remove chisq image from config->file:PSPHOT.INPUT
     psphotStackRemoveChisqFromInputs(config, STACK_DET);
@@ -445 +471 @@
         psphotStackRemoveChisqFromInputs(config, STACK_SRC);
     }
+#endif
 
     psFree (objects);
@@ -474 +501 @@
     return true;
 }
+
+// read the memory usage data from /proc and log them out
+// This will only work on a system that has /proc (not MacOS for instance) but since this function just
+// tries to open and read from a file it is safe
+// XXX: refine this and move it to psLib
+static void logMemStats(const char *heading) {
+
+    // file containing memory statistics for this process proc. See proc(5)
+    const char* statm_path = "/proc/self/statm";
+
+    FILE *f = fopen(statm_path,"r");
+    if (!f) {
+        psLogMsg ("psphot", PS_LOG_WARN, "failed to open %s", statm_path);
+        return;
+    }
+
+    unsigned long vmSize, resident, share, text, lib, data;
+
+    int nread;
+    nread = fscanf(f,"%ld %ld %ld %ld %ld %ld", &vmSize, &resident, &share, &text, &lib, &data);
+    fclose(f);
+    if (nread != 6) {
+        psLogMsg ("psphot", PS_LOG_WARN, "failed to read 6 items from %s", statm_path);
+        return;
+    }
+  
+    // XXX: assuming 4 KB page size here
+#   define PAGES_TO_MB(_v) (_v * 4.096 / 1024.)
+    psLogMsg ("psphot", PS_LOG_INFO, "Memory usage at %20s: Total VmSize: %8.2f MB   Resident %8.2f MB   Data: %8.2f MB\n",
+        heading, PAGES_TO_MB(vmSize), PAGES_TO_MB(resident), PAGES_TO_MB(data));
+}
+
+
 
 /* here is the process: