Changeset 31990 for branches/eam_branches/ipp-20110710/psphot/src/psphotSourceFits.c

Timestamp:

Aug 3, 2011, 10:30:02 AM (15 years ago)

Author:

eugene

Message:

add new window function using first radial moment; use new sersic guess function

File:

• branches/eam_branches/ipp-20110710/psphot/src/psphotSourceFits.c (modified) (19 diffs)

branches/eam_branches/ipp-20110710/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;
@@ -225 +225 @@
     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
 
+# 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)) {
+        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;
+    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);
 
     // XXX note that this changes the source moments that are published...
@@ -234 +245 @@
     // this uses the footprint to judge both radius and aperture?
     // XXX save the psf-based moments for output
+    // XXX do not limit to pixels in the marked circle
     psfMoments = *source->moments;
-    if (!pmSourceMoments (source, radius, 0.0, 0.5, 0.0, maskVal)) {
+    if (!pmSourceMoments (source, windowRadius, 0.0, 0.5, 0.0, maskVal & PS_NOT_IMAGE_MASK(markVal))) {
       *source->moments = psfMoments;
+      fprintf (stderr, "skip (bad moments)\n");
       return false;
     }
+    // fprintf (stderr, "r1: %f\n", source->moments->Mrf);
+
+    // XXX note that we do not re-try to measure the moments if the resulting Mrf is large compared 
+    // to the initial value
 
     psTrace ("psphot", 5, "trying blob...\n");
@@ -263 +280 @@
         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 +288 @@
         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 +301 @@
         okEXT = psphotEvalEXT (tmpSrc, EXT);
         chiEXT = EXT ? EXT->chisqNorm : NAN;
+        EXT->fitRadius = fitRadius;
     }
 
@@ -476 +494 @@
     // 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);
     }
 
@@ -499 +511 @@
 }
 
+# define TIMING 0
+
 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
 
@@ -517 +531 @@
         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 +548 @@
         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"); }
+
+    // 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 +576 @@
     // 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);
+    if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }
+
+    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);
+    }
+
+    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 +596 @@
 }
 
+# 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 +621 @@
     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 +636 @@
         }
 
-        // 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 +652 @@
     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 +683 @@
     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 +696 @@
         }
 
-        // 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 +703 @@
         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 +719 @@
     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];