Changeset 36198

Timestamp:

Oct 9, 2013, 4:14:19 PM (13 years ago)

Author:

eugene

Message:

merge from trunk

Location:

branches/eam_branches/ipp-20130904/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/psphot.h (modified) (1 diff)
• src/psphotExtendedSourceAnalysis.c (modified) (11 diffs)
• src/psphotExtendedSourceFits.c (modified) (3 diffs)
• src/psphotMergeSources.c (modified) (2 diffs)
• src/psphotPetrosianRadialBins.c (modified) (2 diffs)
• src/psphotPetrosianStats.c (modified) (4 diffs)
• src/psphotRadialProfile.c (modified) (7 diffs)
• src/psphotSourceFits.c (modified) (1 diff)
• src/psphotStackImageLoop.c (modified) (9 diffs, 1 prop)
• src/psphotStackMatchPSFsNext.c (modified) (2 diffs)
• src/psphotStackReadout.c (modified) (16 diffs)

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

@@ -479 +479 @@
 bool psphotSourceChildren (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc);
 bool psphotSourceChildrenReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
-psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc);
+psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc, bool sourcesSubtracted);
 
 bool psphotSourceParents (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc);

branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceAnalysis.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+void psphotRadialProfileShowSkips ();
 
 // measure the elliptical radial profile and use this to measure the petrosian parameters for the sources
@@ -33 +34 @@
 }
 
+/*** for the moment, this test code : it is not thread safe ***/
+int    Nall = 0;
+int  Nskip1 = 0;
+int  Nskip2 = 0;
+int  Nskip3 = 0;
+int  Nskip4 = 0;
+int  Nskip5 = 0;
+int  Nskip6 = 0;
+int  Nskip7 = 0;
+int  Nskip8 = 0;
+int  Nskip9 = 0;
+int Nskip10 = 0;
+int Nskip11 = 0;
+int Nskip12 = 0;
+int Nskip13 = 0;
+int Nskip14 = 0;
+
+# define SKIP(VALUE) { VALUE++; continue; }
+
 // aperture-like measurements for extended sources
 bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
 
     bool status;
+    int NfaintEXT = 0;
+    int NfaintPSF = 0;
+
     int Next = 0;
     int Npetro = 0;
@@ -108 +131 @@
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NfaintEXT
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NfaintPSF
 
 // set this to 0 to run without threading
-# if (1)            
+# if (0)            
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
@@ -129 +154 @@
             scalar = job->args->data[7];
             Nannuli += scalar->data.S32;
+            scalar = job->args->data[8];
+            NfaintEXT += scalar->data.S32;
+            scalar = job->args->data[9];
+            NfaintPSF += scalar->data.S32;
             psFree(job);
 # endif
@@ -153 +182 @@
                 scalar = job->args->data[7];
                 Nannuli += scalar->data.S32;
+                scalar = job->args->data[8];
+                NfaintEXT += scalar->data.S32;
+                scalar = job->args->data[9];
+                NfaintPSF += scalar->data.S32;
             }
             psFree(job);
@@ -163 +196 @@
     psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
+    psLogMsg ("psphot", PS_LOG_INFO, "  skipped: %d EXT, %d PSF\n", NfaintEXT, NfaintPSF);
+
+    fprintf (stderr, "ext analysis skipped @ 1  : %d\n", Nskip1);
+    fprintf (stderr, "ext analysis skipped @ 2  : %d\n", Nskip2);
+    fprintf (stderr, "ext analysis skipped @ 3  : %d\n", Nskip3);
+    fprintf (stderr, "ext analysis skipped @ 4  : %d\n", Nskip4);
+    fprintf (stderr, "ext analysis skipped @ 5  : %d\n", Nskip5);
+    fprintf (stderr, "ext analysis skipped @ 6  : %d\n", Nskip6);
+    fprintf (stderr, "ext analysis skipped @ 7  : %d\n", Nskip7);
+    fprintf (stderr, "ext analysis skipped @ 8  : %d\n", Nskip8);
+    fprintf (stderr, "ext analysis skipped @ 9  : %d\n", Nskip9);
+    fprintf (stderr, "ext analysis skipped @ 10 : %d\n", Nskip10);
+    fprintf (stderr, "ext analysis skipped @ 11 : %d\n", Nskip11);
+    fprintf (stderr, "ext analysis skipped @ 12 : %d\n", Nskip12);
+    fprintf (stderr, "ext analysis skipped @ 13 : %d\n", Nskip13);
+    fprintf (stderr, "ext analysis skipped @ 14 : %d\n", Nskip14);
+
+    psphotRadialProfileShowSkips ();
 
     psphotVisualShowResidualImage (readout, false);
@@ -177 +228 @@
 
     bool status;
+
+    int NfaintEXT = 0;
+    int NfaintPSF = 0;
 
     int Next = 0;
@@ -207 +261 @@
         pmSource *source = sources->data[i];
 
+        Nall ++;
+
         // if we have checked the source validity on the basis of the object set, then 
         // we either skip these tests below or we skip the source completely
-        if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) continue;
+        if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) SKIP (Nskip1);
         if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_KEEP) goto keepSource;
 
         // skip PSF-like and non-astronomical objects
-        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
-        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
-        if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
+        if (source->type == PM_SOURCE_TYPE_DEFECT) SKIP (Nskip2);
+        if (source->type == PM_SOURCE_TYPE_SATURATED) SKIP (Nskip3);
+        if (source->mode & PM_SOURCE_MODE_DEFECT) SKIP (Nskip4);
+        if (source->mode & PM_SOURCE_MODE_SATSTAR) SKIP (Nskip5);
 
         // skip saturated stars modeled with a radial profile 
-        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) SKIP (Nskip6);
 
         // optionally allow non-extended objects to get petrosians as well
         if (!doPetroStars) {
-            if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
-            if (source->type == PM_SOURCE_TYPE_STAR) continue;
+            if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) SKIP (Nskip7);
+            // if (source->type == PM_SOURCE_TYPE_STAR) SKIP (Nskip8); -- XXX this might be set for input lists (instead of EXT_LIMIT)
         }
 
@@ -233 +289 @@
         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
             skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr);
+            NfaintEXT ++;
         } else {
             skipSource = (sqrt(source->peak->detValue) < SN_LIM);
-        }
-        if (skipSource) continue;
+            NfaintPSF ++;
+        }
+        if (skipSource) SKIP (Nskip9);
 
         // limit selection by analysis region (this automatically apply
-        if (source->peak->x < region->x0) continue;
-        if (source->peak->y < region->y0) continue;
-        if (source->peak->x > region->x1) continue;
-        if (source->peak->y > region->y1) continue;
+        if (source->peak->x < region->x0) SKIP (Nskip10);
+        if (source->peak->y < region->y0) SKIP (Nskip11);
+        if (source->peak->x > region->x1) SKIP (Nskip12);
+        if (source->peak->y > region->y1) SKIP (Nskip13);
 
     keepSource:
@@ -266 +324 @@
                 psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
                 pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-                continue;
+                SKIP (Nskip14);
             }
             Nannuli ++;
@@ -304 +362 @@
     scalar->data.S32 = Nannuli;
 
+    scalar = job->args->data[8];
+    scalar->data.S32 = NfaintEXT;
+
+    scalar = job->args->data[9];
+    scalar->data.S32 = NfaintPSF;
+
     return true;
 }

branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c

@@ -475 +475 @@
         }
 
+        bool testObject = false;
+        // testObject |= ((fabs(source->peak->xf -  179) < 5) && (fabs(source->peak->yf - 1138) < 5));
+        // testObject |= ((fabs(source->peak->xf - 5047) < 5) && (fabs(source->peak->yf -  151) < 5));
+        // testObject |= ((fabs(source->peak->xf - 3929) < 5) && (fabs(source->peak->yf - 4109) < 5));
+        // testObject |= ((fabs(source->peak->xf -  915) < 5) && (fabs(source->peak->yf - 5998) < 5));
+        // testObject |= ((fabs(source->peak->xf - 5406) < 5) && (fabs(source->peak->yf -  326) < 5));
+        if (testObject) {
+            fprintf (stderr, "test object @ %f, %f\n", source->peak->xf, source->peak->yf);
+            psTraceSetLevel ("psModules.objects.pmPCM_MinimizeChisq", 5);
+            psTraceSetLevel ("psphot.psphotExtendedSourceFits_Threaded", 5);
+        }
+
         // loop here over the models chosen for each source (exclude by S/N)
         // Reset the iterator
@@ -531 +543 @@
                   psFree (source->modelFlux);
                   source->modelFlux = NULL;
-                  // if ((fabs(source->peak->xf - 2572) < 20) && (fabs(source->peak->yf - 5874) < 20)) {
-                  //   fprintf (stderr, "test object\n");
-                  // }
-
                   modelFit = psphotFitEXT (modelFit, readout, source, fitOptions, modelType, maskVal, markVal);
                   if (!modelFit) {
@@ -646 +654 @@
         psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq);
         psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+
+        if (testObject) {
+            psTraceSetLevel ("psModules.objects.pmPCM_MinimizeChisq", 0);
+            psTraceSetLevel ("psphot.psphotExtendedSourceFits_Threaded", 0);
+        }
     }
 

branches/eam_branches/ipp-20130904/psphot/src/psphotMergeSources.c

@@ -844 +844 @@
 // array containing the child sources.  XXX currently, this is only used by psphotStackReadout
 // (sources go on allSources so that psphotChoosePSF can be called repeatedly)
-psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc) {
+psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc, bool sourcesSubtracted) {
 
     bool status;
@@ -951 +951 @@
             // child sources have not been subtracted in this image, but this flag may be raised if
             // they were subtracted in the parent's image
-            sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+            // XXX NOTE : in the pre-20130914 version of psphotStack, we carried a copy of the pixels 
+            // generated before the subtraction took place (and then we smoothed to match the desired PSF).  
+            // in the new version, we copy the image after subtraction; we need to distinguish these cases
+            if (!sourcesSubtracted) {
+                sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+            }
 
             // set the output detections:

branches/eam_branches/ipp-20130904/psphot/src/psphotPetrosianRadialBins.c

@@ -36 +36 @@
     pmSourceRadialProfile *profile = source->extpars->petProfile;
 
-    float skyModelErrorSQ = PS_SQR(skynoise);
+    // float skyModelErrorSQ = PS_SQR(skynoise);
 
     int nMax = radiusMax;
@@ -133 +133 @@
         if (radius->data.F32[i] > Rmax) {
             // calculate the value for the nOut bin
-            float value, dvalue;
+            float value; // , dvalue;
             if (values->n > 0) {
                 psVectorStats (stats, values, NULL, NULL, 0);
                 value = stats->robustMedian;
-                dvalue = stats->robustStdev;
+                // dvalue = stats->robustStdev;
             } else {
                 value = NAN;
-                dvalue = NAN;
+                // dvalue = NAN;
             }
 
             binSB->data.F32[nOut] = value;
-            binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ);
+            // binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ);
+            binSBstdev->data.F32[nOut] = skynoise / sqrt(values->n);
             binFill->data.F32[nOut] = values->n / binArea->data.F32[nOut];
 

branches/eam_branches/ipp-20130904/psphot/src/psphotPetrosianStats.c

@@ -6 +6 @@
 // generate the Petrosian radius and flux from the mean surface brightness (r_i)
 
+float InterpolateValuesQuadratic (float *Xin, float *Yin, float X);
 float InterpolateValues     (float X0, float Y0, float X1, float Y1, float X);
 float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1);
@@ -53 +54 @@
     float dFsum2 = 0.0;
 
-    float nSigma = 3.0;
+    float nSigma = 2.0;
     int lowestSignificantRadius = 0;
     float lowestSignificantRatio = 1.0;
@@ -118 +119 @@
                 petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[nOut]);
             } else {
-                petRadius    = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
-                petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]);
+              // petRadius    = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
+              if (nOut > 1) {
+                petRadius    = InterpolateValuesQuadratic (&petRatio->data.F32[nOut-2],   &refRadius->data.F32[nOut-2],   PETROSIAN_RATIO);
+              } else {
+                petRadius    = InterpolateValuesQuadratic (&petRatio->data.F32[nOut-3], &refRadius->data.F32[nOut-3], PETROSIAN_RATIO);
+              }
+              float petRadiusLinear = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
+              if (fabs(petRadius - petRadiusLinear) > fabs(refRadius->data.F32[nOut] - refRadius->data.F32[nOut-1])) { 
+                fprintf (stderr, "big difference : %f vs %f\n", petRadius, petRadiusLinear); 
+              }
+              petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]);
             }
             above = false;
@@ -235 +245 @@
 }
 
+// Lagrange's form of the interpolating polynomial...
+float InterpolateValuesQuadratic (float *Xin, float *Yin, float X) {
+
+  float dx01 = Xin[0] - Xin[1];
+  float dx02 = Xin[0] - Xin[2];
+  float dx12 = Xin[1] - Xin[2];
+
+  float dx0  = X - Xin[0];
+  float dx1  = X - Xin[1];
+  float dx2  = X - Xin[2];
+
+  float y0 = Yin[0]*dx1*dx2/(dx01*dx02);
+  float y1 = Yin[1]*dx0*dx2/(dx01*dx12); // need - sign
+  float y2 = Yin[2]*dx0*dx1/(dx02*dx12); 
+
+  float Y = y0 - y1 + y2;
+  return Y;
+}
+
 float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) {
     float dydx = (Y1 - Y0) / (X1 - X0);

branches/eam_branches/ipp-20130904/psphot/src/psphotRadialProfile.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+
+static int Nskip1 = 0;
+static int Nskip2 = 0;
+static int Nskip3 = 0;
+static int Nskip4 = 0;
+static int Nskip5 = 0;
+
+# define SKIP(VALUE) { VALUE++; return false; }
 
 bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) {
@@ -23 +31 @@
     if (!psphotRadialProfilesByAngles (source, Nsec, Rmax)) {
         psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian");
-        return false;
+        SKIP (Nskip1);
     }
     // allocate: extpars->radFlux->radii,fluxes,theta
@@ -32 +40 @@
     if (!psphotRadiiFromProfiles (source, fluxMin, fluxMax)) {
         psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles");
-        return false;
+        SKIP (Nskip2);
     }
     // allocate : extpars->radFlux->isophotalRadii (use profile->radii,fluxes)
@@ -40 +48 @@
     if (!psphotEllipticalContour (source)) {
         // psLogMsg ("psphot", 3, "failed to measure elliptical contour");
-        return false;
+        SKIP (Nskip3);
     }
     // use extpars->radFlux->isophotalRadii,theta (result in extpars->axes)
@@ -48 +56 @@
     if (!psphotEllipticalProfile (source, RAW_RADIUS)) {
         psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
-        return false;
+        SKIP (Nskip4);
     }
     // allocate extpars->ellipticalFlux->radiusElliptical,fluxElliptical (use axes to scale raw pixels)
@@ -55 +63 @@
     if (!psphotRadialBins (recipe, source, Rmax, skynoise)) {
         psError (PS_ERR_UNKNOWN, false, "failed to generate radial bins");
-        return false;
+        SKIP (Nskip5);
     }
     // allocate extpars->radProfile->binSB, binSBstdv, binSum, binFill, radialBins, area (small lengths)
@@ -62 +70 @@
     return true;
 }
+
+void psphotRadialProfileShowSkips () {
+  
+  fprintf (stderr, "radial profile skipped @ 1  : %d\n", Nskip1);
+  fprintf (stderr, "radial profile skipped @ 2  : %d\n", Nskip2);
+  fprintf (stderr, "radial profile skipped @ 3  : %d\n", Nskip3);
+  fprintf (stderr, "radial profile skipped @ 4  : %d\n", Nskip4);
+  fprintf (stderr, "radial profile skipped @ 5  : %d\n", Nskip5);
+}

branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c

@@ -589 +589 @@
     if (TIMING) { psTimerStart ("psphotFitPCM"); }
 
+    // if we are ever (in a given psphot implementation) going to fit a parameter, we must set the options here to include
+    // that parameter (otherwise pmPCMupdate will fail to allocate the dmodelFlux image)
+    // thus, if the sersic analysis below uses an index fit, need to use this EXT_AND_SKY mode for init
+
+    options.mode = PM_SOURCE_FIT_EXT;
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-        options.mode = PM_SOURCE_FIT_NO_INDEX; // XXX note that there may be a conflict with psphotExtendedSourceFits.c:133
         options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
-        // if we are ever (in a given psphot implementation) going to fit a parameter, we must set the options here to include
-        // that parameter (otherwise pmPCMupdate will fail to allocate the dmodelFlux image)
-        // thus, if the sersic analysis below uses an index fit, need to use this EXT_AND_SKY mode for init
-    } else {
+    }
+    if (modelType == pmModelClassGetType("PS_MODEL_DEV")) {
+        options.mode = PM_SOURCE_FIT_SHAPE;
+        options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
+    }
+    if (modelType == pmModelClassGetType("PS_MODEL_EXP")) {
         options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
     }

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

@@ -22 +22 @@
 
     pmFPAview *view = pmFPAviewAlloc (0);
-    pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
-    pmFPAfile *inputCnv = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.CNV");
-    pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+    pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
 
     if (!input) {
@@ -32 +30 @@
 
     psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
-    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
-    if (useRaw && inputRaw == NULL) {
-        psLogMsg ("psphot", 1, "PSPHOT.STACK.USE.RAW set but no raw input.");
-        useRaw = false;
-    }
-
-    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");
-        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.MASK.CNV");
-        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.VARIANCE.CNV");
-        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.PSF.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");
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
 
@@ -71 +53 @@
                 psMemDump("load");
 
-                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;
-                    }
+                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");
@@ -126 +99 @@
 
     // Load the appropriate EXPNUM image
-    pmFPAfileActivate (config->files, true, useRaw ? "PSPHOT.STACK.EXPNUM.RAW" : "PSPHOT.STACK.EXPNUM.CNV");
+    pmFPAfileActivate (config->files, true, "PSPHOT.STACK.EXPNUM.RAW");
 
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa EXPNUM in psphot.");
@@ -146 +119 @@
                 if (! readout->data_exists) { continue; }
 
-                if (!psphotSetNFrames (config, view, useRaw ? inputRaw->name : inputCnv->name)) ESCAPE ("failed to setNFrames.");
+                if (!psphotSetNFrames (config, view, input->name)) ESCAPE ("failed to setNFrames.");
             }
         }
@@ -167 +140 @@
 bool UpdateHeadersForReadout (pmConfig *config, pmFPAview *view) {
 
-    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
-    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
-
     int num = psphotFileruleCount(config, "PSPHOT.INPUT");
 
@@ -179 +149 @@
         psAssert (output, "missing file?");
 
-        pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
-        pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
-        pmFPAfile *input = useRaw ? inputRaw : inputCnv;
+        pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
         psAssert (input, "missing input file");
 
@@ -205 +173 @@
         psAssert (output, "missing file?");
 
-        pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
-        pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
-        pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+        pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
         psAssert (input, "missing input file");
 
@@ -242 +208 @@
         psAssert (output, "missing file?");
 
-        pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
-        pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
-        pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+        pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
         psAssert (input, "missing input file");
 

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

@@ -43 +43 @@
     for (int i = 0; i < num; i++) {
         if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, lastSize)) {
-            psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF", filerule, i);
+            psLogMsg ("psphot", PS_LOG_INFO, "failed to smooth image %s (%d) to target PSF", filerule, i);
             psImageConvolveSetThreads(oldThreads);
             return false;
@@ -102 +102 @@
 
     if (targetFWHM <= currentFWHM) {
-        psError (PSPHOT_ERR_CONFIG, true, "target FWHM cannot be smaller than current FWHM");
+        // psError (PSPHOT_ERR_CONFIG, true, "target FWHM cannot be smaller than current FWHM");
+        psLogMsg ("psphot", PS_LOG_INFO, "target FWHM (%f) is smaller than current FWHM (%f), not smoothing\n", targetFWHM, currentFWHM);
+        fwhmValues->data.F32[lastSize + 1] = currentFWHM;
         return false;
     }

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

@@ -1 +1 @@
 # include "psphotInternal.h"
 
+static bool psphotStackMatchPSFsetup (pmConfig *config, const pmFPAview *view, const char *filerule, const char *fPSF);
+static bool psphotStackMatchPSFsetupReadout (pmConfig *config, const pmFPAview *view, const char *filerule, const char *fPSF, int index);
 static bool psphotStackLoadWCS(pmConfig *config, const pmFPAview *view, const char *filerule);
 static void logMemStats(const char *heading);
 
-// we have 3 possible real filesets:
+// relevant filesets:
 # define STACK_RAW "PSPHOT.STACK.INPUT.RAW"
-# define STACK_CNV "PSPHOT.STACK.INPUT.CNV"
-# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"  /* the psf-matched image */
-
-// we have 3 files on which we operate:
-// DET (detection image)       : nominally RAW (optionally CNV?)
-// SRC (source analysis image) : nominally CNV (optionally RAW)
-// OUT (psf-matched images)    : always OUT
+# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"
+
+// XXX STACK_OUT currently is a copy of STACK_RAW, but should be a pointer to is as in psphot (single)
+
+// TEST CODE, can be removed
+bool psphotDumpImages (pmConfig *config, const pmFPAview *view, const char *filerule, char *base) {
+
+    // XXX do nothing
+    return true;
+
+    int num = psphotFileruleCount(config, "PSPHOT.INPUT");
+
+    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?");
+
+        char line[256];
+        snprintf (line, 256, "%s.%d.im.fits", base, i);
+        psphotSaveImage (NULL, readout->image, line);
+
+        snprintf (line, 256, "%s.%d.wt.fits", base, i);
+        psphotSaveImage (NULL, readout->variance, line);
+
+        snprintf (line, 256, "%s.%d.mk.fits", base, i);
+        psphotSaveImage (NULL, readout->mask, line);
+    }
+    return true;
+}
 
 bool psphotStackVisualFilerule(pmConfig *config, const pmFPAview *view, const char *filerule) {
@@ -67 +94 @@
     psAssert (breakPt, "configuration error: set BREAK_POINT");
 
-    // we have 3 relevant files: RAW (unconvolved), CNV (convolved stack), OUT (psf-matched stack)
-    // select which image (RAW or CNV) is used for analysis (RAW always used for detection)
-    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
-    char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV;
-    char *STACK_DET = STACK_RAW;
-
     // load WCS 
-    if (!psphotStackLoadWCS(config, view, STACK_SRC)) {
-        psError (PSPHOT_ERR_CONFIG, false, "trouble loading WCS for %s", STACK_SRC);
+    if (!psphotStackLoadWCS(config, view, STACK_RAW)) {
+        psError (PSPHOT_ERR_CONFIG, false, "trouble loading WCS for %s", STACK_RAW);
         return false;
     }
 
     // set the photcode for each image
-    if (!psphotAddPhotcode (config, view, STACK_SRC)) {
+    if (!psphotAddPhotcode (config, view, STACK_RAW)) {
         psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
         return false;
@@ -86 +107 @@
 
     // 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);
+    if (!psphotSetMaskAndVariance (config, view, STACK_RAW)) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "NOTHING")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // generate a background model (median, smoothed image)
-    if (!psphotModelBackground (config, view, STACK_DET)) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    if (!psphotSubtractBackground (config, view, STACK_DET)) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    if (strcmp(STACK_SRC, STACK_DET)) {
-#define MODEL_BACKGROUND_SRC 1
-#ifdef MODEL_BACKGROUND_SRC
-        // work around the fact that the background levels on the convolved
-        // and unconvolved stacks can be different
-        if (!psphotModelBackground (config, view, STACK_SRC)) {
-            return psphotReadoutCleanup (config, view, STACK_SRC);
-        }
-#endif
-        if (!psphotSubtractBackground (config, view, STACK_SRC)) {
-            return psphotReadoutCleanup (config, view, STACK_SRC);
-        }
+    if (!psphotModelBackground (config, view, STACK_RAW)) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "BACKMDL")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+
+// XXX TEST for background:
+    if (!psphotModelBackground (config, view, STACK_RAW)) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!psphotModelBackground (config, view, STACK_RAW)) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+// XXX TEST END 
 
 #ifdef MAKE_CHISQ_IMAGE
     // also make the chisq detection image
-    if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
+    if (!psphotStackChisqImage(config, view, STACK_RAW, STACK_RAW)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!strcasecmp (breakPt, "CHISQ")) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 #endif
-    if (!strcasecmp (breakPt, "CHISQ")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
 
     // 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, STACK_DET, true)) { // pass 1
-        // this only happens if we had an error in psphotFindDetections
+    if (!psphotFindDetections (config, view, STACK_RAW, true)) { // pass 1
         psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    // If DET and SRC are different images, copy the detections from DET to SRC.  This 'copy'
-    // is just a copy of the container pointer; the sources on both DET and SRC are the same
-    // memory objects
-    if (strcmp(STACK_SRC, STACK_DET)) {
-        if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
-            psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-            return psphotReadoutCleanup (config, view, STACK_SRC);
-        }
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // construct sources and measure basic stats (saved on detections->newSources)
-    if (!psphotSourceStats (config, view, STACK_SRC, true)) { // pass 1
+    if (!psphotSourceStats (config, view, STACK_RAW, true)) { // pass 1
         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "PEAKS")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    // psphotDumpTest (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    // psphotDumpTest (config, view, STACK_RAW);
     psMemDump("sourcestats");
     logMemStats("sourcestats");
@@ -164 +172 @@
     // classify sources based on moments, brightness
     // only run this on detections from the input images, not chisq image
-    if (!psphotRoughClass (config, view, STACK_SRC)) {
+    if (!psphotRoughClass (config, view, STACK_RAW)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    // If DET and SRC are different images, subtract radial profiles for the convolved
-    // image first.  The profiles found on the convolved image will be replaced by those
-    // found for the unconvolved image.  Downstream, in psphotFindDetections, we will
-    // replace the the profiles (and re-subtract them) for the detection image, so we want
-    // to keep those versions of the profiles on the sources.
-    if (strcmp(STACK_SRC, STACK_DET)) {
-      // find and subtract radial profile models for saturated stars (XXX change name eventually)
-      if (!psphotDeblendSatstars (config, view, STACK_SRC)) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+
+    // find and subtract radial profile models for saturated stars (XXX change name eventually)
+    if (!psphotDeblendSatstars (config, view, STACK_RAW)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-      }
-    }
-    // find and subtract radial profile models for saturated stars (XXX change name eventually)
-    if (!psphotDeblendSatstars (config, view, STACK_DET)) {
-        psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     // only run this on detections from the input images, not chisq image
-    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
+    if (!psphotImageQuality (config, view, STACK_RAW)) { // pass 1
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "MOMENTS")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // use bright stellar objects to measure PSF
-    if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
+    if (!psphotChoosePSF (config, view, STACK_RAW, true)) { // pass 1
         psLogMsg ("psphot", 3, "failure to construct a psf model");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "PSFMODEL")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // merge the newly selected sources into the existing list
     // NOTE: merge OLD and NEW
-    psphotMergeSources (config, view, STACK_SRC);
+    psphotMergeSources (config, view, STACK_RAW);
 
     // Construct an initial model for each object, set the radius to fitRadius, set circular
     // fit mask.  NOTE: only applied to sources without guess models
-    psphotGuessModels (config, view, STACK_SRC);
+    psphotGuessModels (config, view, STACK_RAW);
 
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
-    psphotFitSourcesLinear (config, view, STACK_SRC, false, false);
-    psphotStackVisualFilerule(config, view, STACK_SRC);
+    psphotFitSourcesLinear (config, view, STACK_RAW, false, false);
+    psphotStackVisualFilerule(config, view, STACK_RAW);
 
     // measure the radial profiles to the sky
-    psphotRadialProfileWings (config, view, STACK_SRC);
+    psphotRadialProfileWings (config, view, STACK_RAW);
 
     // re-measure the kron mags with models subtracted.  this pass starts with a circular
@@ -225 +221 @@
     // but iterates to an appropriately larger size
     logMemStats("before.kron.1");
-    psphotKronIterate(config, view, STACK_SRC, 1);
+    psphotKronIterate(config, view, STACK_RAW, 1);
     logMemStats("after.kron.1");
         
     // identify CRs and extended sources
-    psphotSourceSize (config, view, STACK_SRC, true);
+    psphotSourceSize (config, view, STACK_RAW, true);
 
     // non-linear PSF and EXT fit to brighter sources
     // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
-    psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotBlendFit (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
     // replace all sources (do NOT ignore subtraction state)
-    psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
+    psphotReplaceAllSources (config, view, STACK_RAW, false); // pass 1 (detections->allSources)
 
     logMemStats("pass1");
@@ -245 +241 @@
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
-    // NOTE 2 : this function subtracts the models from the given filerule (SRC), not DET
-    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 2 (detections->allSources)
+    // NOTE 2 : this function subtracts the models from the given filerule
+    psphotFitSourcesLinear (config, view, STACK_RAW, true, false); // pass 2 (detections->allSources)
 
     // NOTE: possibly re-measure background model here with objects subtracted / or masked
@@ -252 +248 @@
     // NOTE: this block performs the 2nd pass low-significance PSF detection stage
     { 
-        // if DET and SRC are different images, generate children sources for all sources in
-        // the SRC image.  This operation replaces the existing DETECTION container on DET
-        // which is currently a view to the one on SRC).  children sources go to
-        // det->allSources
-        if (strcmp(STACK_SRC, STACK_DET)) {
-            psphotSourceChildren (config, view, STACK_DET, STACK_SRC); 
-
-            //  subtract all sources from DET (this will subtract using the psf model for SRC, which
-            //  will somewhat oversubtract the sources -- this is OK
-            psphotRemoveAllSources (config, view, STACK_DET, false); // do not ignore subtraction state for sources
-        }
-
         // add noise for subtracted objects
-        psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
+        psphotAddNoise (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
         // find fainter sources
         // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
-        psphotFindDetections (config, view, STACK_DET, false); // pass 2 (detections->peaks, detections->footprints)
+        psphotFindDetections (config, view, STACK_RAW, false); // pass 2 (detections->peaks, detections->footprints)
 
         // remove noise for subtracted objects (ie, return to normal noise level)
@@ -276 +260 @@
         bool footprintsUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED");
         if (!footprintsUseUnsubtracted) {
-            psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
+            psphotSubNoise (config, view, STACK_RAW); // pass 1 (detections->allSources)
         }
-
-        // if DET and SRC are different images, copy the detections from DET to SRC 
-        // (this operation just ensures the metadata container has a view on SRC as well
-        if (strcmp(STACK_SRC, STACK_DET)) {
-            // replace all sources in DET
-            psphotReplaceAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
-
-            // copy the newly detected peaks from DET to SRC so SourceStats below can operate on them
-            if (!psphotCopyPeaks (config, view, STACK_SRC, STACK_DET)) {
-                psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-                return psphotReadoutCleanup (config, view, STACK_SRC);
-            }
-        }
 
         // define new sources based on only the new peaks & measure moments
         // NOTE: new sources are saved on detections->newSources
-        psphotSourceStats (config, view, STACK_SRC, false); // pass 2 (detections->newSources)
+        psphotSourceStats (config, view, STACK_RAW, false); // pass 2 (detections->newSources)
 
         // set source type
         // NOTE: apply only to detections->newSources
-        if (!psphotRoughClass (config, view, STACK_SRC)) { // pass 2 (detections->newSources)
+        if (!psphotRoughClass (config, view, STACK_RAW)) { // pass 2 (detections->newSources)
             psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
-            return psphotReadoutCleanup (config, view, STACK_SRC);
+            return psphotReadoutCleanup (config, view, STACK_RAW);
         }
 
         // replace all sources so fit below applies to all at once
         // NOTE: apply only to OLD sources (which have been subtracted)
-        psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 2
+        psphotReplaceAllSources (config, view, STACK_RAW, false); // pass 2
 
         // merge the newly selected sources into the existing list
         // NOTE: merge OLD and NEW
         // XXX check on free of sources...
-        psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
+        psphotMergeSources (config, view, STACK_RAW); // (detections->newSources + detections->allSources -> detections->allSources)
 
         // Construct an initial model for each object, set the radius to fitRadius, set circular
         // fit mask.  NOTE: only applied to sources without guess models
-        psphotGuessModels (config, view, STACK_SRC);
+        psphotGuessModels (config, view, STACK_RAW);
     }
 
@@ -325 +296 @@
     if (splitLinearFit) {
         psLogMsg ("psphot", 3, "splitting fit of detected and matched soures\n");
-        // Fit the detected sources separately from matched that wea are about to create.
+        // Fit the detected sources separately from matched ones that we are about to create.
         // NOTE: apply to ALL sources but only include sources with postitive flux in the fit
-        psphotFitSourcesLinear (config, view, STACK_SRC, true, true); // pass 3 (detections->allSources)
+        psphotFitSourcesLinear (config, view, STACK_RAW, true, true); // pass 3 (detections->allSources)
     }
 
@@ -335 +306 @@
     // this just match the detections for the chisq image, and not bother measuring the source
     // stats in that case...?
-    objects = psphotMatchSources (config, view, STACK_SRC);
+    objects = psphotMatchSources (config, view, STACK_RAW);
     psMemDump("matchsources");
 
@@ -344 +315 @@
     // Construct an initial model for each object, set the radius to fitRadius, set circular
     // fit mask.  NOTE: only applied to sources without guess models
-    psphotGuessModels (config, view, STACK_SRC);
+    psphotGuessModels (config, view, STACK_RAW);
 
     psphotStackObjectsUnifyPosition (objects);
 
-    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
+    psphotStackObjectsSelectForAnalysis (config, view, STACK_RAW, objects);
 
     // final linear fit. NOTE: if splitLinearFit is true above, this pass will only fit
     // the unsubtracted (matched) sources (the sources that we fit above are subtracted)
-    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
+    psphotFitSourcesLinear (config, view, STACK_RAW, true, false); // pass 4 (detections->allSources)
 
     // measure the radial profiles to the sky (only measures new objects)
-    psphotRadialProfileWings (config, view, STACK_SRC);
+    psphotRadialProfileWings (config, view, STACK_RAW);
 
     // re-measure the kron mags with models subtracted
     // psphotKronMasked(config, view, STACK_SRC);
     logMemStats("before.kron.2");
-    psphotKronIterate(config, view, STACK_SRC, 2);
+    psphotKronIterate(config, view, STACK_RAW, 2);
     logMemStats("after.kron.2");
 
     // measure source size for the remaining sources
     // NOTE: applies only to NEW (unmeasured) sources
-    psphotSourceSize (config, view, STACK_SRC, false); // pass 2 (detections->allSources)
+    psphotSourceSize (config, view, STACK_RAW, false); // pass 2 (detections->allSources)
 
     psMemDump("psfstats");
@@ -371 +342 @@
     // drop matched sources without any useful measurements and set kron radii for the ones
     // we decide to keep
-    psphotFilterMatchedSources (config, view, STACK_SRC, objects);
+    psphotFilterMatchedSources (config, view, STACK_RAW, objects);
 
     // measure kron fluxes for the matched sources only
-    psphotKronIterate(config, view, STACK_SRC, 3);
+    psphotKronIterate(config, view, STACK_RAW, 3);
 
     // measure elliptical apertures, petrosians (objects sorted by S/N)
     // psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
-    psphotExtendedSourceAnalysis (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotExtendedSourceAnalysis (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
     // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
-    psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotExtendedSourceFits (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
     // create source children for the OUT filerule (for radial aperture photometry and output) 
-    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
+    // NOTE: The new source children have image arrays pointing to the readout associated with 
+    // STACK_OUT.  in psphotStackMatchPSFsetup, we copy the current pixel values from RAW to OUT, 
+    // but keep the pointers the same so we do not break these source image references
+    // XXX NOTE : if we use the pre-20130914 psphotStackReadout code, we need to use 'false' for the
+    // sourcesSubtracted argument
+    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects, true);
     if (!objectsOut) {
         psFree(objects);
         psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
@@ -396 +372 @@
         // this forces photometry on the undetected sources from other images
 
-        // NOTE: we always do the radial apertures analysis on the convolved image since
-        // 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
-        // target psf for the RAW inputs
-
-        // If useRaw copy the sources to the convolved readout
-        if (strcmp(STACK_SRC, STACK_CNV)) {
-            if (!psphotCopySources (config, view, STACK_CNV, STACK_SRC)) {
-                psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-                return psphotReadoutCleanup (config, view, STACK_SRC);
-            }
-        }
-        // mark any inputs that we want to skip the matched apertures for
-        psphotStackSetInputsToSkip(config, view, STACK_CNV, true);
-        psphotStackSetInputsToSkip(config, view, STACK_OUT, true);
-        psphotRadialApertures (config, view, STACK_CNV, 0); // entry 0 == unmatched
-        psMemDump("extmeas");
-
+        // set up the FWHM vector
+        psphotStackMatchPSFsetup (config, view, STACK_OUT, STACK_RAW);
+        psphotDumpImages (config, view, STACK_RAW, "raw.t0");
+        psphotDumpImages (config, view, STACK_OUT, "out.t0");
 
         int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
-        for (int entry = 1; entry < nRadialEntries; entry++) {
+
+        for (int entry = 0; entry < nRadialEntries; entry++) {
             // NOTE: entry 0 is the unmatched image set
 
-            // re-measure the PSF for the smoothed image (using entries in 'allSources')
-            psphotChoosePSF (config, view, STACK_OUT, false);
-
-            // this is necessary to update the models based on the new PSF
-            psphotResetModels (config, view, STACK_OUT);
-
-            // this is necessary to get the right normalization for the new models
-            psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
+            char line[256];
 
             // measure circular, radial apertures (objects sorted by S/N)
             psphotRadialApertures (config, view, STACK_OUT, entry); 
+            snprintf (line, 256, "%s.%d", "out.t1", entry);
+            psphotDumpImages (config, view, STACK_OUT, line);
 
             // replace the flux in the image so it is returned to its original state
             psphotReplaceAllSources (config, view, STACK_OUT, false);
-
-            // smooth to the next FWHM, or set 'smoothAgain' to false if no more 
-            psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
-            psMemDump("matched");
+            snprintf (line, 256, "%s.%d", "out.t2", entry);
+            psphotDumpImages (config, view, STACK_OUT, line);
+
+            if (entry < nRadialEntries - 1) {
+                // smooth to the next FWHM
+                // this function does nothing if the targetFWHM is smaller than the currentFWHM
+                psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
+                snprintf (line, 256, "%s.%d", "out.t3", entry);
+                psphotDumpImages (config, view, STACK_OUT, line);
+                psMemDump("matched");
+
+                // re-measure the PSF for the smoothed image (using entries in 'allSources')
+                psphotChoosePSF (config, view, STACK_OUT, false);
+
+                // this is necessary to update the models based on the new PSF
+                psphotResetModels (config, view, STACK_OUT);
+
+                // this is necessary to get the right normalization for the new models
+                // and to subtract the sources
+                psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
+                snprintf (line, 256, "%s.%d", "out.t4", entry);
+                psphotDumpImages (config, view, STACK_OUT, line);
+            }
         }
     }
-    psphotStackSetInputsToSkip(config, view, STACK_CNV, false);
-    psphotStackSetInputsToSkip(config, view, STACK_OUT, false);
 
     // measure aperture photometry corrections
-    if (!psphotApResid (config, view, STACK_SRC)) {
+    if (!psphotApResid (config, view, STACK_RAW)) {
         psFree (objects);
         psFree (objectsOut);
         psLogMsg ("psphot", 3, "failed on psphotApResid");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // calculate source magnitudes
-    psphotMagnitudes(config, view, STACK_SRC);
-
-    if (!useRaw) {
-        // 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)) {
-            psLogMsg ("psphot", 3, "failure to construct a psf model for raw input");
-            return psphotReadoutCleanup (config, view, STACK_DET);
-        }
-    }
-    if (!psphotEfficiency(config, view, STACK_DET)) {
+    psphotMagnitudes(config, view, STACK_RAW);
+
+    if (!psphotEfficiency(config, view, STACK_RAW)) {
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
         psErrorClear();
     }
-    psphotCopyEfficiency (config, view, STACK_OUT, STACK_DET);
+    psphotCopyEfficiency (config, view, STACK_OUT, STACK_RAW);
 
     logMemStats("final");
 #if (1)
-    psphotSourceMemory(config, view, STACK_SRC);
+    psphotSourceMemory(config, view, STACK_RAW);
     psphotSourceMemory(config, view, STACK_OUT);
 #endif
 
-    // replace failed sources?
-    // psphotReplaceUnfitSources (sources);
-
     // replace background in residual image
-    psphotSkyReplace (config, view, STACK_DET);
+    psphotSkyReplace (config, view, STACK_RAW);
 
     // drop the references to the image pixels held by each source
+    psphotSourceFreePixels (config, view, STACK_RAW);
     psphotSourceFreePixels (config, view, STACK_OUT);
-    psphotSourceFreePixels (config, view, STACK_SRC);
 
 #ifdef MAKE_CHISQ_IMAGE
     // remove chisq image from config->file:PSPHOT.INPUT
-    psphotStackRemoveChisqFromInputs(config, STACK_DET);
-    if (strcmp(STACK_SRC, STACK_DET)) {
-        psphotStackRemoveChisqFromInputs(config, STACK_SRC);
-    }
+    psphotStackRemoveChisqFromInputs(config, STACK_RAW);
 #endif
 
@@ -496 +456 @@
 
     // create the exported-metadata and free local data
-    return psphotReadoutCleanup (config, view, STACK_SRC);
+    return psphotReadoutCleanup (config, view, STACK_RAW);
 }
 
@@ -552 +512 @@
 }
 
-
-
 /* here is the process:
 
- * we have three(*) images:
- * RAW : unconvolved image stack 
- * CNV : input convolved image
+ * we have two image sets:
+ * RAW : unconvolved image stacks
 
  * OUT : psf-matched output image (there may be more than one of
@@ -649 +606 @@
 
    */
+
+
+// generate a vector fwhmValues where the first has the fwhm of the raw image and the
+// successive entries have the target values
+
+bool psphotStackMatchPSFsetup (pmConfig *config, const pmFPAview *view, const char *filerule, const char *filerulePSF) {
+
+    bool status;
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+
+        if (!psphotStackMatchPSFsetupReadout (config, view, filerule, filerulePSF, i)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to define target PSF sizes");
+            return false;
+        }
+    }
+
+    return true;
+}
+
+float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index);
+
+// copy the pixels from RAW to OUT (
+
+bool psphotStackMatchPSFsetupReadout (pmConfig *config, const pmFPAview *view, const char *fileruleOut, const char *fileruleRaw, int index) {
+
+    bool status;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+
+    // find the currently selected readout
+    pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, fileruleOut, index); // File of interest
+    psAssert (fileOut, "missing file?");
+
+    pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa);
+    psAssert (readoutOut, "missing readout?");
+
+    // find the currently selected readout
+    pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, fileruleRaw, index); // File of interest
+    psAssert (fileRaw, "missing file?");
+
+    pmReadout *readoutRaw = pmFPAviewThisReadout(view, fileRaw->fpa);
+    psAssert (readoutRaw, "missing readout?");
+
+    readoutOut->image = psImageCopy(readoutOut->image, readoutRaw->image, PS_TYPE_F32);
+    if (readoutRaw->variance) {
+        readoutOut->variance = psImageCopy(readoutOut->variance, readoutRaw->variance, PS_TYPE_F32);
+    }
+    if (readoutRaw->mask) {
+        readoutOut->mask = psImageCopy(readoutOut->mask, readoutRaw->mask, PS_TYPE_IMAGE_MASK);
+    }
+
+    // pmChip *chipRaw = pmFPAviewThisChip(view, fileRaw->fpa); // The chip holds the PSF
+    // psAssert (chipRaw, "missing chip");
+
+    pmPSF *psf = psMetadataLookupPtr(&status, readoutRaw->analysis, "PSPHOT.PSF"); // PSF
+    if (!psf) {
+        // we should have a PSF by this point in psphot
+        psError(PSPHOT_ERR_PROG, true, "Unable to find PSF.");
+        return false;
+    }
+
+    float fwhmRaw = psphotPSFseeing (psf, readoutRaw, index);
+
+    psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
+    psVectorAppend(fwhmValues, fwhmRaw);
+
+    // is a single target FWHM specified, or a set of values?  set up the vector options->targetSeeing and the local 1st value
+    float targetSeeing = psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
+    if (!status) {
+        psVector *targetSeeing = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
+        psAssert (status, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
+        for (int i = 0; i < targetSeeing->n; i++) {
+            psVectorAppend(fwhmValues, targetSeeing->data.F32[i]);
+        }           
+    } else {
+        psVectorAppend(fwhmValues, targetSeeing);
+    }
+
+    psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
+    psFree (fwhmValues);
+
+    return true;
+}
+
+float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index) {
+
+    psImage *image = readout->image;
+
+    int Nx = image->numCols;
+    int Ny = image->numRows;
+
+    float sumFWHM = 0.0;                  // FWHM for image
+    int numFWHM = 0;                      // Number of FWHM measurements
+    for (float x = 0; x < Nx; x += 0.25*Nx) {
+        for (float y = 0; y < Ny; y += 0.25*Ny) {
+            float fwhm = pmPSFtoFWHM(psf, x, y);
+            if (isfinite(fwhm)) {
+                sumFWHM += fwhm;
+                numFWHM++;
+            }
+        }
+    }
+    if (numFWHM == 0) {
+        psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", index);
+        return NAN;
+    } 
+
+    float fwhm = sumFWHM / (float) numFWHM;
+    psLogMsg ("psphotStack", PS_LOG_INFO, "Input Seeing for %d: %f\n", index, fwhm);
+
+    return fwhm;
+}