Changeset 25755 for trunk/psphot/src/psphotSourceFits.c

Timestamp:

Oct 2, 2009, 3:12:47 PM (17 years ago)

Author:

eugene

Message:

check in changes from gene@development branch : extensive changes to moments calculation, psf model generation, aperture residuals; improvements to extended source analysis

Location:

trunk/psphot/src

Files:

• . (modified) (1 prop)
• psphotSourceFits.c (modified) (19 diffs)

trunk/psphot/src

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy

trunk/psphot/src/psphotSourceFits.c

@@ -90 +90 @@
     psphotCheckRadiusPSFBlend (readout, source, PSF, markVal, dR);
 
-    // fit PSF model (set/unset the pixel mask)
+    // fit PSF model
     pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal);
+
+    // clear the circular mask
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
+
+    if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
 
     // correct model chisq for flux trend
@@ -101 +106 @@
         pmSource *blend = sourceSet->data[i];
         pmModel *model  = modelSet->data[i];
+
+        if (!isfinite(model->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
 
         // correct model chisq for flux trend
@@ -120 +127 @@
         pmSourceCacheModel (blend, maskVal);
         pmSourceSub (blend, PM_MODEL_OP_FULL, maskVal);
-        blend->tmpFlags |=  PM_SOURCE_TMPF_SUBTRACTED;
         blend->mode |=  PM_SOURCE_MODE_BLEND_FIT;
     }
@@ -144 +150 @@
     pmSourceCacheModel (source, maskVal);
     pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-    source->tmpFlags |=  PM_SOURCE_TMPF_SUBTRACTED;
     source->mode |=  PM_SOURCE_MODE_BLEND_FIT;
     return true;
@@ -167 +172 @@
     // fit PSF model (set/unset the pixel mask)
     pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF, maskVal);
+
+    if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
+
+    // clear the circular mask
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
 
     // correct model chisq for flux trend
@@ -186 +196 @@
     pmSourceCacheModel (source, maskVal);
     pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-    source->tmpFlags |=  PM_SOURCE_TMPF_SUBTRACTED;
-    return true;
-}
-
-static float EXT_MIN_SN;
-static float EXT_MOMENTS_RAD;
+    return true;
+}
+
 static pmModelType modelTypeEXT;
 
@@ -197 +204 @@
 
     bool status;
-
-    // extended source model parameters
-    EXT_MIN_SN       = psMetadataLookupF32 (&status, recipe, "EXT_MIN_SN");
-    EXT_MOMENTS_RAD  = psMetadataLookupF32 (&status, recipe, "EXT_MOMENTS_RADIUS");
 
     // extended source model descriptions
@@ -221 +224 @@
     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
-    if (source->peak->SN < EXT_MIN_SN) return false;
-
+
+    // set the radius based on the footprint (also sets the mask pixels)
+    float radius = psphotSetRadiusEXT (readout, source, markVal);
+
+    // XXX note that this changes the source moments that are published...
     // recalculate the source moments using the larger extended-source moments radius
     // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
-    if (!pmSourceMoments (source, EXT_MOMENTS_RAD, 0.0, 0.0)) return false;
+    // this uses the footprint to judge both radius and aperture?
+    if (!pmSourceMoments (source, radius, 0.0, 0.0)) return false;
 
     psTrace ("psphot", 5, "trying blob...\n");
@@ -237 +244 @@
     // XXX need to handle failures better here
     pmModel *EXT = psphotFitEXT (readout, source, modelTypeEXT, maskVal, markVal);
+    if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
+
     okEXT = psphotEvalEXT (tmpSrc, EXT);
     chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN;
@@ -246 +255 @@
     // XXX should I keep / save the flags set in the eval functions?
 
+    // clear the circular mask
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
+
     // correct first model chisqs for flux trend
     chiDBL = NAN;
     ONE = DBL->data[0];
     if (ONE) {
+        if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
       chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
       ONE->chisqNorm = ONE->chisq / chiTrend;
@@ -258 +271 @@
     ONE = DBL->data[1];
     if (ONE) {
+        if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
       chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
       ONE->chisqNorm = ONE->chisq / chiTrend;
@@ -277 +291 @@
 
     // both models failed; reject them both
+    // XXX -- change type flags to psf in this case and keep original moments?
     psFree (EXT);
     psFree (DBL);
@@ -287 +302 @@
     // save new model
     source->modelEXT = EXT;
+    source->modelEXT->fitRadius = radius;
     source->type = PM_SOURCE_TYPE_EXTENDED;
     source->mode |= PM_SOURCE_MODE_EXTMODEL;
@@ -293 +309 @@
     pmSourceCacheModel (source, maskVal);
     pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-    source->tmpFlags |=  PM_SOURCE_TMPF_SUBTRACTED;
+
+# if (PS_TRACE_ON)   
     psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]);
+    if (psTraceGetLevel("psphot") >= 6) {
+        psLogMsg ("psphot", 1, "source 2:\n");
+        for (int i = 0; i < source->modelEXT->params->n; i++) {
+            psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelEXT->params->data.F32[i], source->modelEXT->dparams->data.F32[i]);
+        }
+    }
+# endif
+
     return true;
 
@@ -304 +329 @@
     psFree (source->modelPSF);
     source->modelPSF = psMemIncrRefCounter (DBL->data[0]);
-    source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     source->mode     |= PM_SOURCE_MODE_PAIR;
+    source->modelPSF->fitRadius = radius;
 
     // copy most data from the primary source (modelEXT, blends stay NULL)
-    // XXX use pmSourceCopy?
     pmSource *newSrc = pmSourceCopy (source);
     newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
+    newSrc->modelPSF->fitRadius = radius;
 
     // build cached models and subtract
@@ -317 +342 @@
     pmSourceCacheModel (newSrc, maskVal);
     pmSourceSub (newSrc, PM_MODEL_OP_FULL, maskVal);
+
+# if (PS_TRACE_ON)   
     psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]);
+    if (psTraceGetLevel("psphot") >= 6) {
+        psLogMsg ("psphot", 1, "source 1:\n");
+        for (int i = 0; i < newSrc->modelPSF->params->n; i++) {
+            psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, newSrc->modelPSF->params->data.F32[i], newSrc->modelPSF->dparams->data.F32[i]);
+        }
+        psLogMsg ("psphot", 1, "source 2:\n");
+        for (int i = 0; i < source->modelPSF->params->n; i++) {
+            psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelPSF->params->data.F32[i], source->modelPSF->dparams->data.F32[i]);
+        }
+        psphotVisualShowResidualImage (readout);
+    }
+# endif
 
     psArrayAdd (newSources, 100, newSrc);
@@ -356 +395 @@
     // save the PSF model from the Ensemble fit
     PSF = source->modelPSF;
-    psphotCheckRadiusPSFBlend (readout, source, PSF, markVal, 8.0);
     if (isnan(PSF->params->data.F32[1])) psAbort("nan in dbl fit");
 
@@ -389 +427 @@
     PS_ASSERT (EXT, NULL);
 
-    // if (isnan(EXT->params->data.F32[1])) psAbort("nan in ext fit");
-
-    psphotCheckRadiusEXT (readout, source, EXT, markVal);
-
     if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
         psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
@@ -401 +435 @@
     return (EXT);
 }
-

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy