Changeset 33098 for branches/meh_branches/ppsub_test/psphot/src/psphotSourceFits.c

Timestamp:

Jan 11, 2012, 11:19:21 PM (15 years ago)

Author:

mhuber

Message:

merged latest r33093 trunk changes to meh_branches/ppsub_test

Location:

branches/meh_branches/ppsub_test

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (2 props)
• psphot/src/psphotSourceFits.c (modified) (27 diffs)

branches/meh_branches/ppsub_test/psphot/src

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

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

branches/meh_branches/ppsub_test/psphot/src/psphotSourceFits.c

@@ -211 +211 @@
 bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
 
-    float radius;
+    float fitRadius, windowRadius;
     bool okEXT, okDBL;
     pmModel *ONE = NULL;
@@ -217 +217 @@
     pmModel *EXT = NULL;
     psArray *DBL = NULL;
-    pmMoments psfMoments;
+    // pmMoments psfMoments;
 
     // skip the source if we don't think it is extended
@@ -225 +225 @@
     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
 
+# define TEST_X -420.0
+# define TEST_Y 300.0
+   
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        fprintf (stderr, "test galaxy\n");
+    }
+
+# undef TEST_X
+# undef TEST_Y
+
     // set the radius based on the footprint (also sets the mask pixels)
-    if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false;
-
-    // XXX note that this changes the source moments that are published...
-    // XXX all published moments should use the same measurement
-    // 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
-    // this uses the footprint to judge both radius and aperture?
-    // XXX save the psf-based moments for output
-    psfMoments = *source->moments;
-    if (!pmSourceMoments (source, radius, 0.0, 0.5, 0.0, maskVal)) {
-      *source->moments = psfMoments;
-      return false;
-    }
+    if (!psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal)) return false;
+    // fprintf (stderr, "rad: %6.1f %6.1f  | %5.2f %5.2f %5.2f  ", source->peak->xf, source->peak->yf, source->moments->Mrf, fitRadius, windowRadius);
 
     psTrace ("psphot", 5, "trying blob...\n");
@@ -258 +257 @@
         okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
         okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
+        okDBL = false; // XXX this is failing badly...
         // XXX should I keep / save the flags set in the eval functions?
 
@@ -263 +263 @@
         ONE = DBL->data[0];
         if (ONE) {
-            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
+            psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
             chiDBL = ONE->chisqNorm; // save chisq for double-star/galaxy comparison
-            ONE->fitRadius = radius;
+            ONE->fitRadius = fitRadius;
         }
 
@@ -271 +271 @@
         ONE = DBL->data[1];
         if (ONE) {
-            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
-            ONE->fitRadius = radius;
+            psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
+            ONE->fitRadius = fitRadius;
         }
     }
@@ -284 +284 @@
         okEXT = psphotEvalEXT (tmpSrc, EXT);
         chiEXT = EXT ? EXT->chisqNorm : NAN;
+        EXT->fitRadius = fitRadius;
     }
 
@@ -294 +295 @@
 
     if (okEXT && okDBL) {
-        psTrace ("psphot", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL);
         // XXX EAM : a bogus bias: need to examine this better
         if (3*chiEXT > chiDBL) goto keepDBL;
@@ -303 +303 @@
     if (!okEXT && okDBL) goto keepDBL;
 
+    psTrace ("psphot", 4, "both failed: blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
+
     // both models failed; reject them both
     // XXX -- change type flags to psf in this case, and make sure we subtract it?
     // reset the psf moments
-    *source->moments = psfMoments;
+    // XXX *source->moments = psfMoments;
 
     psFree (EXT);
@@ -313 +315 @@
 
 keepEXT:
+    psTrace ("psphot", 4, "goto EXT : blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
     // sub EXT
     psFree (DBL);
@@ -338 +341 @@
 
     // reset the psf moments
-    *source->moments = psfMoments;
+    // XXX *source->moments = psfMoments;
     return true;
 
 keepDBL:
+    psTrace ("psphot", 4, "goto DBL : blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
     // sub DLB
     psFree (EXT);
@@ -354 +358 @@
     pmSource *newSrc = pmSourceCopy (source);
     newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
+    newSrc->peak->footprint = source->peak->footprint; // just a reference; the peak does not own the footprint
+    psArrayAdd(newSrc->peak->footprint->peaks, 1, newSrc->peak); // the footprint owns the peak
 
     // build cached models and subtract
@@ -372 +378 @@
             psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelPSF->params->data.F32[i], source->modelPSF->dparams->data.F32[i]);
         }
-        psphotVisualShowResidualImage (readout, false);
     }
 # endif
-
-    // reset the (original) psf moments
-    *source->moments = psfMoments;
-    *newSrc->moments = psfMoments;
 
     psArrayAdd (newSources, 100, newSrc);
@@ -386 +387 @@
 
 escape:
-    // reset the psf moments
-    *source->moments = psfMoments;
     psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
     psFree (tmpSrc);
@@ -476 +475 @@
     // for sersic models, use a grid search to choose an index, then float the params there
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-        // for the test fits, use a somewhat smaller radius
-        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
         psphotFitSersicIndex (model, readout, source, fitOptions, maskVal, markVal);
-    }
-
-    if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
-        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
     }
 
@@ -494 +487 @@
     pmSourceFitModel (source, model, &options, maskVal);
     // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
-
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+
     return (model);
 }
+
+# define TIMING 0
 
 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
@@ -517 +512 @@
         return NULL;
     }
+
+# define TEST_X -540.0
+# define TEST_Y 540.0
+   
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
+    }
+
+    float t1, t2, t4, t5;
+    if (TIMING) { psTimerStart ("psphotFitPCM"); }
 
     pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
@@ -524 +529 @@
         return model;
     }
-
-    // use the source moments, etc to guess basic model parameters
-    if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
-        psFree (pcm);
-        model->flags |= PM_MODEL_STATUS_BADARGS;
-        return model;
-    }
-
-    // for sersic models, use a grid search to choose an index, then float the params there
+    if (TIMING) { t1 = psTimerMark ("psphotFitPCM"); }
+
+    // NOTE : 65 allocs to here
+    // get the guess for sersic models 
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-        // for the test fits, use a somewhat smaller radius
-        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
-
-        if (!psphotFitSersicIndexPCM (pcm, readout, source, fitOptions, maskVal, markVal, psfSize)) {
+        // use the source moments, etc to guess basic model parameters
+        if (!psphotSersicModelClassGuessPCM (pcm, source)) {
             psFree (pcm);
             model->flags |= PM_MODEL_STATUS_BADARGS;
             return model;
         }
-    }
-
-    if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
-        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
-    }
+    } else {
+        // use the source moments, etc to guess basic model parameters
+        if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
+            psFree (pcm);
+            model->flags |= PM_MODEL_STATUS_BADARGS;
+            return model;
+        }
+    }
+
+    if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); }
 
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
@@ -555 +558 @@
     // update the pcm elements if we have changed the circumstance (options.mode or source->pixels)
     pmPCMupdate(pcm, source, &options, model);
+    if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
-    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
-    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);  // NOTE : 1687 allocs in here
+    if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }
+
+    if (TIMING) {
+        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);
+    }
+
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
+    }
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
@@ -566 +579 @@
 }
 
+# undef TEST_X
+# undef TEST_Y
+
 // note that these should be 1/2n of the standard sersic index
-float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
-# define N_INDEX_GUESS 6
+// float indexGuess[] = {0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0};
+// float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
+float indexGuess[] = {1.0, 2.0, 3.0, 4.0};
+# define N_INDEX_GUESS 4
 
 // A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
@@ -586 +604 @@
     float chiSquare[N_INDEX_GUESS];
 
+# define TEST_X -540.0
+# define TEST_Y 540.0
+   
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
+    }
+
     for (int i = 0; i < N_INDEX_GUESS; i++) {
-        model->params->data.F32[PM_PAR_7] = indexGuess[i];
+        model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[i];
 
         if (!model->modelGuess(model, source)) {
@@ -594 +619 @@
         }
 
-        // each time we change the model guess, we need to adjust the radius
-        // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
-        if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
-            psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
-        }
-        
         pmSourceFitModel (source, model, &options, maskVal);
-        // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+        // fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
 
         chiSquare[i] = model->chisqNorm;
@@ -616 +635 @@
     assert (iMin >= 0);
 
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
+    }
+
     model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
-    model->params->data.F32[PM_PAR_7] = indexGuess[iMin];
+    model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[iMin];
     model->modelGuess(model, source);
-
-    // each time we change the model guess, we need to adjust the radius
-    // if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
-    //  psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal);
-    // }
 
     return true;
@@ -648 +666 @@
     float xMin = NAN;
     float chiSquare[N_INDEX_GUESS];
+
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
+    }
 
     for (int i = 0; i < N_INDEX_GUESS; i++) {
@@ -657 +679 @@
         }
 
-        // each time we change the model guess, we need to adjust the radius
-        // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
-        if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
-            psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
-        }
-        
 # if (0)
+        // this block is to test the relative speed of straight and PCM fits
         pmSourceFitModel (source, model, &options, maskVal);
 # else
@@ -669 +686 @@
         pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
 # endif
+        fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
         // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
 
@@ -684 +702 @@
     assert (iMin >= 0);
     
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
+    }
+
     model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
     model->params->data.F32[PM_PAR_7] = indexGuess[iMin];

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