Changeset 38383

Timestamp:

Jun 5, 2015, 12:26:38 PM (11 years ago)

Author:

bills

Message:

Various changes that make psphotStack -updatemode get good residual images for
extended sources and psfs that were modeled well originally

File:

• trunk/psphot/src/psphotStackUpdateReadout.c (modified) (6 diffs)

trunk/psphot/src/psphotStackUpdateReadout.c

@@ -58 +58 @@
     }
 
-    // the intitial run.
+    // Load the psf from the previous run. 
+    // Note this must be done before we load the sources because it is used.
+    // I believe only to set the residuals
     if (!psphotLoadPSF (config, view, STACK_RAW)) {
         // this only happens if we had a programming error in psphotLoadPSF
@@ -78 +80 @@
     }
 
+#ifdef MEASURE_PSF
+    if (!psphotChoosePSF (config, view, STACK_RAW, true)) {
+        psLogMsg ("psphot", 3, "failure to construct a psf model");
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!strcasecmp (breakPt, "PSFMODEL")) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+#endif
+
     if (!psphotDeblendSatstars (config, view, STACK_RAW)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
         return psphotReadoutCleanup (config, view, STACK_RAW);
     }
+
     if (!psphotMergeSources (config, view, STACK_RAW)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to merge sources");
@@ -218 +231 @@
         COPY_PARAM( "MSKY_NX" );
         COPY_PARAM( "MSKY_NY" );
+        COPY_PARAM( "CHIP_SEEING" );
 }
 
@@ -294 +308 @@
     psphotInitRadiusEXT (recipe, readout);
 
+    // not really fitting but need to get some things set up to instantiate the pcm models
+
+    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
+    fitOptions->mode = PM_SOURCE_FIT_EXT_AND_SKY; 
+    fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+    // XXX: change all of this status checking to asserts. All should be defined with the current recipe.
+    float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolution
+    if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) {
+            fitNsigmaConv = 5.0;
+    }
+    fitOptions->nsigma         = fitNsigmaConv;
+    // Poisson or Constant weight for chisq tests?
+    fitOptions->poissonErrors = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_POISSON");
+    if (!status) fitOptions->poissonErrors = true;
+    int psfSize  = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
+    assert (status);
+
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
+        source->imageID = index;     // XXX is this necessary?
         pmModel *modelPSF = source->modelPSF;
-        if (!modelPSF) {
-            continue;
-        }
         // free any previous radial aperture measurements
         psFree(source->radialAper);
 
         // Guess Models ususally does this
-        psphotCheckRadiusPSF (readout, source,  modelPSF, markVal);
-        source->modelPSF->residuals = psf->residuals;
+        if (modelPSF) {
+            psphotCheckRadiusPSF (readout, source,  modelPSF, markVal);
+            source->modelPSF->residuals = psf->residuals;
+        }
+        if (!(source->mode & PM_SOURCE_MODE_PSFMODEL)) {
+            // The cmf loaders unconditionallly create psf models even if the source did not have one.
+            fprintf(stderr, "source %d %5d should not have a psf model 0x%08X%08X\n", index, source->seq, source->mode2, source->mode);
+            psFree(source->modelPSF) ;
+        }
 
         // make sure that the window radius is large enough. Should we do this regardless of whether
@@ -318 +354 @@
         }
         
+        // If we find a good model we will cache it below.
         bool goodModel = false;
-        // the cmf readers do not set PM_PAR_I0 for the models. 
-        // Calculate it here.
+
+        // the cmf readers do not set PM_PAR_I0 for the models so we calculate it here.
         // We may want to put this code in the cmf readers but for now I don't want to change psphot operation
-        // for anything except psphotStack -updatemode
-        if (isfinite(source->psfFlux)) {
+        // for anything except psphotStack -updatemode.
+        if (modelPSF && isfinite(source->psfFlux)) {
+            // The cmf reader is incorrectly setting the sky parameter for psf models
+            modelPSF->params->data.F32[PM_PAR_SKY] = 0.0;
+            modelPSF->dparams->data.F32[PM_PAR_SKY] = 0.0;
+            // Calculate flux for normalized model
+            float saveI0 = modelPSF->params->data.F32[PM_PAR_I0];
             modelPSF->params->data.F32[PM_PAR_I0] = 1.0;
             float normFlux = modelPSF->class->modelFlux(modelPSF->params);
-            modelPSF->params->data.F32[PM_PAR_I0] =  source->psfFlux / normFlux;
-            goodModel = true;
-        }
-        if (source->modelFits && source->modelFits->n) {
-            for (int m = 0; m < source->modelFits->n; m++) {
-                pmModel *model = source->modelFits->data[m];
-                if (isfinite(model->mag)) {
-                    // calculate flux from model magnitude
-                    float modelFlux = pow(10., -0.4 * model->mag);
-                    model->params->data.F32[PM_PAR_I0] = 1.0;
-                    float normFlux = modelPSF->class->modelFlux(modelPSF->params);
-                    model->params->data.F32[PM_PAR_I0] =  modelFlux / normFlux;
-                    if (model == source->modelEXT) {
-                        goodModel = true;
-                    }
-                } else {
-                    // nan model magnitude no way to calculate the flux. 
-                    // The flux was probably negative.
-                    // If this is the extended model for this source we can't use the model.
-                    // Set source type to star to prevent source subtraction from crashing. 
-                    // This is somewhat rare .2% of extened sources. Is there another way that we can estimate 
-                    // the flux, or more to the point the value for PM_PAR_I0?
-                    if (model == source->modelEXT && source->type == PM_SOURCE_TYPE_EXTENDED) {
-                        source->type = PM_SOURCE_TYPE_STAR;
-                        goodModel = false;
-                    }
+            float psfFlux = source->psfFlux;
+#define APPLY_APCORR
+#ifdef APPLY_APCORR
+            double apTrend = pmTrend2DEval (psf->ApTrend, (float)source->peak->x, (float)source->peak->y);
+            double adjustment = pow(10.0, -0.4*apTrend);
+            if (isfinite(adjustment) && adjustment != 0.0) {
+                psfFlux = psfFlux / adjustment;
+            } else {
+                fprintf(stderr, "apTrend adjustment for source %5d invalid: %f\n", source->seq, adjustment);
+            }
+#endif
+            float newI0 =  psfFlux / normFlux;
+            if (isfinite (newI0) ) {
+                // psf model looks to be good. 
+                modelPSF->params->data.F32[PM_PAR_I0] = newI0;
+#ifdef USE_SAVED_I0
+                // I am testing with a modified cmf that saves PM_PAR_I0
+                modelPSF->params->data.F32[PM_PAR_I0] = saveI0;
+#endif
+                // fprintf (stderr, "%5d %12f %12f %12f\n", source->seq, saveI0, newI0, saveI0 - newI0);
+                goodModel = true;
+            } else {
+                psAssert(isfinite (modelPSF->params->data.F32[PM_PAR_I0]), "got infinite I0 for psf model");
+            }
+        }
+
+        // It is rather expensive initializing all of the extended models. Lets only initialize
+        // the extended model that will be used if any.
+
+        bool isPSF = false;
+        pmModel *model = pmSourceGetModel (&isPSF, source);
+        if (model && !isPSF) {
+            // selected model is extended
+            goodModel = false;
+            // calculate flux from the previously measured model magnitude. It will be nan if the
+            // flux was negative unfortunately.
+            if (isfinite(model->mag)) {
+                float modelFlux = pow(10., -0.4 * model->mag);
+
+                // XXX: We are assuming here that we only use PCM models. Get from recipe.
+                pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
+                if (!pcm) {
+                    psError (PSPHOT_ERR_UNKNOWN, false, "failed to to initialize PCM for model");
+                    return false;
                 }
+                model->params->data.F32[PM_PAR_I0] = 1.0;
+                pmPCMMakeModel (source, model, pcm->nsigma, maskVal, psfSize);
+                float normFlux = model->class->modelFlux(model->params);
+                float I0 = modelFlux / normFlux;
+                if (isfinite(I0)) {
+                    model->params->data.F32[PM_PAR_I0] = I0;
+                    goodModel = true;
+                    // Usually this it is set when doing the fit. Since we are instantiating the model we
+                    // set it explicitly. 
+                    model->isPCM = true;
+                }
+                psFree(pcm);
+            }
+            if (!goodModel) {
+                // We do not have a usable extended model for this source.
+                // The original flux was probably negative.
+                // This is somewhat rare. I saw it in about .2% of extended sources in my first test.
+                // Set source type to star to prevent source subtraction from crashing. 
+                // this may cause model psf to be used in places.
+                // XXX However we aren't going to cache it below.
+                source->type = PM_SOURCE_TYPE_STAR;
+                model = NULL;
             }
         }
@@ -373 +455 @@
         }
 
-        // Should we use pmPCMCacheModel for extended sources?
-        // XXX: are the extended models ready?
-        if (goodModel) {
-            pmSourceCacheModel (source, maskVal); // ALLOC x14 (!)
-        }
-    }
-
-    // as a test drop psf so we can try and remeasure it
-//    psMetadataRemoveKey(readout->analysis, "PSPHOT.PSF");
+        if (goodModel && model) {
+            if (model->isPCM) {
+                pmPCMCacheModel (source, maskVal, psfSize, fitNsigmaConv);
+            } else {
+                pmSourceCacheModel (source, maskVal);
+            }
+        }
+    }
+
+    psFree(fitOptions);
+
+#ifdef MEASURE_PSF
+    // need to drop the psf model in order for it to be re-measured
+    // Hmm, will this invalidate the results from above where we used the loaded one?
+    psMetadataRemoveKey(readout->analysis, "PSPHOT.PSF");
+#endif
 
     return true;