Changeset 38515

Timestamp:

Jun 22, 2015, 3:23:53 PM (11 years ago)

Author:

bills

Message:

Handle pcm extended source models correctly in psphotSourceChildrenByReadout and psphotResetModels.
Avoid extra model convolutions.

Location:

trunk/psphot/src

Files:

• psphot.h (modified) (1 diff)
• psphotExtendedSourceFits.c (modified) (1 diff)
• psphotMergeSources.c (modified) (6 diffs)
• psphotReplaceUnfit.c (modified) (3 diffs)
• psphotStackUpdateReadout.c (modified) (8 diffs)

trunk/psphot/src/psphot.h

@@ -513 +513 @@
 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, bool sourcesSubtracted);
+psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *fileruleOut, const char *fileruleSrc, psArray *objectsSrc, bool sourcesSubtracted);
 
 bool psphotSourceParents (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc);

trunk/psphot/src/psphotExtendedSourceFits.c

@@ -149 +149 @@
     fitOptions->poissonErrors = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_POISSON");
     if (!status) fitOptions->poissonErrors = true;
+
+    // Save pointer to fitOptions for possible use later reinstantiating the pcm models on the output readout
+    psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PCM_FIT_OPTIONS", PS_DATA_UNKNOWN | PS_META_REPLACE, "pcm fit options", 
+        fitOptions);
 
     // maskVal is used to test for rejected pixels, and must include markVal

trunk/psphot/src/psphotMergeSources.c

@@ -853 +853 @@
 // 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, bool sourcesSubtracted) {
+psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *fileruleOut, const char *fileruleSrc, psArray *objectsSrc, bool sourcesSubtracted) {
 
     bool status;
 
-    int nImages = psphotFileruleCount(config, filerule);
+    int nImages = psphotFileruleCount(config, fileruleOut);
 
     // generate look-up arrays for detections and readouts
     psArray *detArrays = psArrayAlloc(nImages);
     psArray *readouts = psArrayAlloc(nImages);
+    psArray *fitOptionsArray = psArrayAlloc(nImages);
+
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    assert (recipe);
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
+    assert (maskVal);
+    int psfSize  = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
+    assert (status);
 
     for (int i = 0; i < nImages; i++) {
 
         // find the currently selected readout
-        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, fileruleOut, i); // File of interest
         psAssert (file, "missing file?");
 
         pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
         psAssert (readout, "missing readout?");
+
+        pmFPAfile *fileSrc = pmFPAfileSelectSingle(config->files, fileruleSrc, i); // File of interest
+        psAssert (file, "missing file?");
+
+        pmReadout *readoutSrc = pmFPAviewThisReadout(view, fileSrc->fpa);
+        psAssert (readoutSrc, "missing readout?");
+
 
         // create DETECTIONS containers for each image, in case one lacks it
@@ -886 +901 @@
             psAssert (detections, "missing detections?");
         }
+        pmSourceFitOptions *fitOptions = psMetadataLookupPtr (&status, readoutSrc->analysis, "PCM_FIT_OPTIONS");
+        psAssert (fitOptions, "missing pcm fit options");
+        psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PCM_FIT_OPTIONS", PS_DATA_UNKNOWN | PS_META_REPLACE, "pcm fit options", fitOptions);
 
         // we need to save the new sources on the detection arrays of the appropriate image
         detArrays->data[i] = psMemIncrRefCounter(detections);
         readouts->data[i] = psMemIncrRefCounter(readout);
+        fitOptionsArray->data[i] = psMemIncrRefCounter(fitOptions);
     }
 
@@ -932 +951 @@
             // does this copy all model data? (NO)
             sourceOut->modelPSF = pmModelCopy(sourceSrc->modelPSF);
-            sourceOut->modelEXT = pmModelCopy(sourceSrc->modelEXT);
-
+
+            bool foundModelEXT = false;
             if (sourceSrc->modelFits) {
                 sourceOut->modelFits = psArrayAlloc(sourceSrc->modelFits->n);
                 for (int j = 0; j < sourceSrc->modelFits->n; j++) {
-                    sourceOut->modelFits->data[j] = pmModelCopy(sourceSrc->modelFits->data[j]);
-                }
+                    pmModel *modelSrc = sourceSrc->modelFits->data[j];
+                    pmModel *modelOut = sourceOut->modelFits->data[j] = pmModelCopy(modelSrc);
+                    if (modelSrc == sourceSrc->modelEXT) {
+                        foundModelEXT = true;
+                        sourceOut->modelEXT = psMemIncrRefCounter (modelOut);
+                    }
+                    modelOut->isPCM = modelSrc->isPCM;
+                }
             }
+            if (!foundModelEXT && sourceSrc->modelEXT) {
+                // Will this ever happen?
+                sourceOut->modelEXT = pmModelCopy(sourceSrc->modelEXT);
+            }
 
             // drop the references to the original image pixels:
@@ -949 +978 @@
             pmReadout *readout = readouts->data[index];
 
+            pmSourceFitOptions *fitOptions = fitOptionsArray->data[index];
+
             // allocate image, weight, mask for the new image for each peak
             if (sourceOut->modelPSF) {
                 pmSourceRedefinePixels (sourceOut, readout, sourceOut->peak->x, sourceOut->peak->y, 
-                                                                        sourceOut->modelPSF->fitRadius);
+                                                                        sourceSrc->windowRadius);
             } else {
                 // if we have no pixels we can't use it to determine the psf so make sure this bit is off
@@ -965 +996 @@
             if (!sourcesSubtracted) {
                 sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
-            }
+            } else {
+                if (sourceSrc->modelFlux) {
+                    bool isPSF = false;
+                    pmModel *model = pmSourceGetModel (&isPSF, sourceOut);
+                    if (model->isPCM) {
+                        pmPCMdata *pcm = pmPCMinit (sourceOut, fitOptions, model, maskVal, psfSize);
+                        if (pcm) {
+                            // pmPCMMakeModel (sourceOut, model, pcm->nsigma, maskVal, psfSize);
+                            pmPCMCacheModel (sourceOut, maskVal, psfSize, pcm->nsigma);
+                            psFree(pcm);
+                        } else {
+                            // What to do here? 
+                            psAssert (pcm, "pmPCMinit failed!");
+                        }
+                    } else {
+                        pmSourceCacheModel (sourceOut, maskVal);
+                    }
+                }
+            }
 
             // set the output detections:
@@ -981 +1030 @@
     psFree (detArrays);
     psFree (readouts);
+    psFree (fitOptionsArray);
 
     return objectsOut;

trunk/psphot/src/psphotReplaceUnfit.c

@@ -316 +316 @@
     }
 
+    pmSourceFitOptions *fitOptions = psMetadataLookupPtr (&status, readout->analysis, "PCM_FIT_OPTIONS");
+    psAssert (fitOptions, "missing pcm fit options");
 
     // XXX the sources have already been copied (merge into here?)
@@ -368 +370 @@
             source->modelPSF = modelPSF;
             source->modelPSF->fitRadius = radius;
-            model = source->modelPSF;
+            //  model = source->modelPSF;
         }
 
         if (model->isPCM) {
+#if 0
             psAssert(false, "this section is not complete");
 
@@ -389 +392 @@
 
             psImageConvolveFFT (source->modelFlux, rawModelFlux, NULL, 0, psfKernel);
-          
             psFree (psfKernel);
             psFree (rawModelFlux);
-        }
-
-        pmSourceCacheModel (source, maskVal);  // ALLOC x14 (!)
+#endif
+           pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
+            if (pcm) {
+                // pmPCMMakeModel (source, model, pcm->nsigma, maskVal, psfSize);
+                pmPCMCacheModel (source, maskVal, psfSize, pcm->nsigma);
+                psFree(pcm);
+            } else {
+                psAssert (pcm, "pmPCMinit failed!");
+            }
+          
+        } else {
+            pmSourceCacheModel (source, maskVal);  // ALLOC x14 (!)
+        }
     }
 

trunk/psphot/src/psphotStackUpdateReadout.c

@@ -116 +116 @@
     // 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);
+    // XXX: do this here so that the variance is availble to rebuild the models
+    psphotStackMatchPSFsetup (config, view, STACK_OUT, STACK_RAW);
+    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, STACK_RAW, objects, true);
     if (!objectsOut) {
         psFree(objects);
@@ -135 +137 @@
         // this forces photometry on the undetected sources from other images
 
-        // set up the FWHM vector
-        psphotStackMatchPSFsetup (config, view, STACK_OUT, STACK_RAW);
+//      Moved this up above
+//      psphotStackMatchPSFsetup (config, view, STACK_OUT, STACK_RAW);
         psphotDumpImages (config, view, STACK_RAW, "raw.t0");
         psphotDumpImages (config, view, STACK_OUT, "out.t0");
@@ -177 +179 @@
                 // and to subtract the sources
                 psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
+
+#ifdef notyet
+            psphotRemoveAllSources (config, view, STACK_OUT, false);
+#endif
+
                 snprintf (line, 256, "%s.%d", "out.t4", entry);
                 psphotDumpImages (config, view, STACK_OUT, line);
@@ -221 +228 @@
         COPY_PARAM( "MSKY_NX" );
         COPY_PARAM( "MSKY_NY" );
-        COPY_PARAM( "CHIP_SEEING" );
+        // COPY_PARAM( "CHIP_SEEING" );
 }
 
@@ -312 +319 @@
     fitOptions->poissonErrors = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_POISSON");
     if (!status) fitOptions->poissonErrors = true;
+
     int psfSize  = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
     assert (status);
+    
+    // Remember these data for later
+    psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PCM_FIT_OPTIONS", PS_DATA_UNKNOWN | PS_META_REPLACE, "pcm fit options", fitOptions);
 
     for (int i = 0; i < sources->n; i++) {
@@ -331 +342 @@
             // 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) ;
+            modelPSF = NULL;
         }
 
         // make sure that the window radius is large enough. Should we do this regardless of whether
         // or not the source is extended?
-        if (source->modelEXT) {
-            float fitRadius = NAN;
-            float windowRadius = NAN;
+        float fitRadius = NAN;
+        float windowRadius = NAN;
+        if (1 || source->modelEXT) {
             if (!psphotSetRadiusMoments (&fitRadius, &windowRadius, readout, source, markVal)) {
                 psError (PSPHOT_ERR_UNKNOWN, false, "failed to set radius for ext source");
@@ -387 +399 @@
                 float modelFlux = pow(10., -0.4 * model->mag);
 
-                // XXX: We are assuming here that we only use PCM models. Should get that information from the recipe.
+                model->fitRadius = fitRadius;
+
+                // XXX: We are assuming here that we only use PCM models. 
+                // I Should be getting that information from the recipe.
                 pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
                 if (!pcm) {
@@ -394 +409 @@
                 }
                 model->params->data.F32[PM_PAR_I0] = 1.0;
-                pmPCMMakeModel (source, model, pcm->nsigma, maskVal, psfSize);
+                // XXX: this duplicates work that pmPCMCacheModel does
+                // pmPCMMakeModel (source, model, pcm->nsigma, maskVal, psfSize);
                 float normFlux = model->class->modelFlux(model->params);
                 float I0 = modelFlux / normFlux;