Changeset 32157

Timestamp:

Aug 21, 2011, 10:44:17 AM (15 years ago)

Author:

eugene

Message:

distinguish window used for moment calculation and window used for model fitting more clearly; re-generate the cached model after PCM fitting using new function (instead of caching the fit models); add function to set window which shrinks as well as expands as needed

Location:

branches/eam_branches/ipp-20110710/psphot/src

Files:

• . (modified) (1 prop)
• psphot.h (modified) (3 diffs)
• psphotExtendedSourceFits.c (modified) (7 diffs)
• psphotKronMasked.c (modified) (3 diffs)
• psphotRadiusChecks.c (modified) (2 diffs)
• psphotReadout.c (modified) (5 diffs)
• psphotSourceFits.c (modified) (2 diffs)

branches/eam_branches/ipp-20110710/psphot/src

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

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

@@ -19 +19 @@
     PSPHOT_FORCED,
     PSPHOT_MAKE_PSF,
+    PSPHOT_MODEL_TEST,
 } psphotImageLoopMode;
 
@@ -328 +329 @@
 bool psphotMakePSFReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
 
+pmConfig *psphotModelTestArguments(int argc, char **argv);
+bool psphotModelTestReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
+
 int psphotFileruleCount(const pmConfig *config, const char *filerule);
 bool psphotFileruleCountSet(const pmConfig *config, const char *filerule, int num);
@@ -476 +480 @@
 
 bool psphotSetRadiusMoments (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal);
+bool psphotSetRadiusMomentsExact (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal);
 
 #endif

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

@@ -316 +316 @@
         Next ++;
 
-        // set the radius based on the first radial moment (also sets the mask pixels)
-        if (!psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal)) {
-            fprintf (stderr, "skipping (1) %f, %f\n", source->peak->xf, source->peak->yf);
-            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-            // XXX raise an error of some kind or set a flag bit
-            continue;
-        }
+        // set the fit radius based on the first radial moment (also sets the mask pixels)
+        psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal);
 
         // Recalculate the source moments using the larger extended-source moments radius.  At
@@ -332 +327 @@
             fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf);
             *source->moments = psfMoments;
+
+            // probably need to bump up the size for subtraction
             pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
             // XXX raise an error flag of some kind
             continue;
         }
-
-        // save the PSF-based modelFlux here in case we need to subtract it (for failure).  If
-        // it does not exist, attempt to generate it.  Give up if we cannot even do that.
-        psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux);
-        if (!modelFluxStart) {
-            pmSourceCacheModel (source, maskVal);
-            modelFluxStart = psMemIncrRefCounter (source->modelFlux);
-            if (!modelFluxStart) {
-                fprintf (stderr, "skipping (3) %f, %f\n", source->peak->xf, source->peak->yf);
-                *source->moments = psfMoments;
-                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-                // XXX raise an error of some kind?
-                continue;
-            }
-        }
-        
-        // array to store the pointers to the model flux images while the models are being fitted
-        psArray *modelFluxes = psArrayAllocEmpty (models->list->n);
 
         // allocate the array to store the model fits
@@ -417 +396 @@
                   continue;
               }
-              pmSourceCacheModel (source, maskVal);
               psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", 
                        source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
@@ -426 +404 @@
               }
           }
-
-          // save each of the model flux images for now, and (below) store the best
-          psArrayAdd (modelFluxes, 4, source->modelFlux);
 
           pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
@@ -477 +452 @@
           psTrace ("psphot", 5, "failed to fit extended source model to object at %f, %f", source->moments->Mx, source->moments->My);
 
-          // replace original model, subtract it
-          psFree (source->modelFlux);
-          source->modelFlux = modelFluxStart;
+          // ensure the modelEXT is cached
+          pmSourceCacheModel (source, maskVal);
 
           *source->moments = psfMoments;
           pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-
-          psFree (modelFluxes);
 
           continue;
@@ -493 +465 @@
         source->modelEXT = psMemIncrRefCounter (source->modelFits->data[minModel]);
 
-        // save the modelFlux for the best model
-        psFree (source->modelFlux);
-        source->modelFlux = psMemIncrRefCounter (modelFluxes->data[minModel]);
+        // adjust the window so the subtraction covers the faint wings
+        psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal);
+
+        // cache the model flux
+        if (source->modelEXT->isPCM) {
+            pmPCMCacheModel (source, maskVal, psfSize);
+        } else {
+            pmSourceCacheModel (source, maskVal);
+        }
 
         // replace the original moments
@@ -502 +480 @@
         // subtract the best fit from the object, leave local sky
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-
-        // the initial model flux is no longer needed
-        psFree (modelFluxStart);
-        psFree (modelFluxes);
 
         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);

branches/eam_branches/ipp-20110710/psphot/src/psphotKronMasked.c

@@ -96 +96 @@
         if (!source->moments) continue;
 
+        float windowRadius = RADIUS;
+        if (source->moments->KronFlux / source->moments->KronFluxErr > 10.0) {
+            windowRadius = PS_MAX (RADIUS, 10*source->moments->Mrf);
+        }
+
+        // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
+        pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
+
+        pmMoments psfMoments = *source->moments;
+
         // replace object in image
         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
@@ -101 +111 @@
         }
 
-        psphotKronMag (source, RADIUS, MIN_KRON_RADIUS, maskVal);
+        psphotKronMag (source, windowRadius, MIN_KRON_RADIUS, maskVal);
 
         // re-subtract the object, leave local sky
@@ -136 +146 @@
         }
     }           
+
     // psphotSaveImage (NULL, kronMask, "kronmask.fits");
     psLogMsg ("psphot.kron", PS_LOG_DETAIL, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);

branches/eam_branches/ipp-20110710/psphot/src/psphotRadiusChecks.c

@@ -144 +144 @@
 # define MIN_WINDOW 5.0
 # define SCALE1 5.0
-# define SCALE2 12.5
+# define SCALE2 12.0
 
 // call this function whenever you (re)-define the EXT model
+// XXX this function does not shrink the window
 bool psphotSetRadiusMoments (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal) {
 
@@ -160 +161 @@
     // redefine the pixels if needed
     pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, *windowRadius);
+
+    // set the mask to flag the excluded pixels
+    psImageKeepCircle (source->maskObj, source->peak->xf, source->peak->yf, *fitRadius, "OR", markVal);
+
+    return true;
+}
+# undef SCALE1
+# undef SCALE2
+# undef MIN_WINDOW
+
+# define MIN_WINDOW 5.0
+# define SCALE1 5.0
+# define PAD_WINDOW 3.0
+
+// call this function whenever you (re)-define the EXT model
+// XXX alternate function to set exactly the desired window size
+bool psphotSetRadiusMomentsExact (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal) {
+
+    psRegion newRegion;
+
+    psAssert (source, "source not defined??");
+    psAssert (source->moments, "moments not defined??");
+
+    *fitRadius = SCALE1 * source->moments->Mrf;
+    *fitRadius = PS_MIN (PS_MAX(*fitRadius, MIN_WINDOW), EXT_FIT_MAX_RADIUS);
+
+    *windowRadius = *fitRadius + PAD_WINDOW;
+
+    // check to see if new region is completely contained within old region
+    newRegion = psRegionForSquare (source->peak->xf, source->peak->yf, *windowRadius);
+    newRegion = psRegionForImage (readout->image, newRegion);
+
+    // redefine the pixels to match
+    pmSourceRedefinePixelsByRegion (source, readout, newRegion);
 
     // set the mask to flag the excluded pixels

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

@@ -137 +137 @@
     }
 
-    // construct sources and measure basic stats (saved on detections->newSources)
+    // construct sources and measure moments and other basic stats (saved on detections->newSources)
+    // all sources use the auto-scaled window appropriate to a PSF, except for the saturated
+    // stars : these use a larger window (3x the basic window)
     if (!psphotSourceStats (config, view, filerule, true)) { // pass 1
         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
@@ -200 +202 @@
     psphotDumpChisqs (config, view, filerule);
 
-    // XXX re-measure the kron mags with models subtracted
+    // re-measure the kron mags with models subtracted.  this pass uses a circular window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
+    
+    // but this is chosen above to be appropriate for the PSF objects (not galaxies)
     psphotKronMasked(config, view, filerule);
 
@@ -214 +218 @@
     psphotDumpChisqs (config, view, filerule);
 
+    // XXX note above the extended sources have their moments remeasured with a larger window,
+    // but the result is replaced with the original PSF-scaled window
+
     // replace all sources
     psphotReplaceAllSources (config, view, filerule); // pass 1 (detections->allSources)
@@ -240 +247 @@
         psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources)
 
-        // define new sources based on only the new peaks
+        // define new sources based on only the new peaks & measure moments
         // NOTE: new sources are saved on detections->newSources
         psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources)
@@ -313 +320 @@
 
     // XXX re-measure the kron mags with models subtracted
+    // Note that this uses the PSF_MOMENTS_RADIUS as a window
     psphotKronMasked(config, view, filerule);
 

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

@@ -511 +511 @@
 }
 
-# define TIMING 0
+# define TIMING 1
 
 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
@@ -583 +583 @@
 
     if (TIMING) {
-        fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, t1, t2, t4, t5);
+        int nPixBig = source->pixels->numCols * source->pixels->numRows;
+        fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t4, t5);
     }
 

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest