Changeset 36025

Timestamp:

Aug 25, 2013, 12:11:59 PM (13 years ago)

Author:

eugene

Message:

grid search on index

File:

• branches/eam_branches/ipp-20130711/psphot/src/psphotSourceFits.c (modified) (5 diffs)

branches/eam_branches/ipp-20130711/psphot/src/psphotSourceFits.c

@@ -558 +558 @@
 }
 
-# define TIMING 0
+# define TIMING 1
+
+bool psphotSersicModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize);
 
 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
@@ -580 +582 @@
 
     float t1, t2, t4, t5;
+    t1 = t2 = t4 = t5 = 0.0;
     if (TIMING) { psTimerStart ("psphotFitPCM"); }
+
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        options.mode = PM_SOURCE_FIT_NO_INDEX;
+    } else {
+        options.mode = PM_SOURCE_FIT_EXT_AND_SKY; // XXX note that there may be a conflict with psphotExtendedSourceFits.c:133
+    }
 
     pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
@@ -594 +603 @@
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
         // use the source moments, etc to guess basic model parameters
-        if (!psphotSersicModelClassGuessPCM (pcm, source)) {
+        if (!psphotSersicModelGuessPCM (pcm, source, maskVal, psfSize)) {
             psFree (pcm);
             model->flags |= PM_MODEL_STATUS_BADARGS;
@@ -613 +622 @@
     if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); }
 
-    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-        options.mode = PM_SOURCE_FIT_NO_INDEX;
-    } else {
-      options.mode = PM_SOURCE_FIT_EXT_AND_SKY; // XXX note that there may be a conflict with psphotExtendedSourceFits.c:133
-    }
     // 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"); }
+    // pmPCMupdate(pcm, source, &options, model);
+    // if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
@@ -751 +755 @@
     return true;
 }
+
+// float indexGuessInv[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
+float indexGuessInv[] = {0.5, 0.25, 0.167, 0.125};
+# define N_INDEX_GUESS_INV 4
+
+// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
+// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
+bool psphotSersicModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize) {
+
+    // we get a reasonable guess from:
+    // * Reff = Kron R1
+    // * Rmajor / Rminor & Theta from moments
+    // * Io from total Kron flux
+
+    // the guesses are used to fill in PAR:
+    psF32 *PAR = pcm->modelConv->params->data.F32;
+
+    // convert the moments to Major,Minor,Theta
+    psEllipseMoments moments;
+
+    if (!isfinite(source->moments->Mrf)) return false;
+    if (!isfinite(source->moments->Mxx)) return false;
+    if (!isfinite(source->moments->Mxy)) return false;
+    if (!isfinite(source->moments->Myy)) return false;
+
+    moments.x2 = source->moments->Mxx;
+    moments.y2 = source->moments->Myy;
+    moments.xy = source->moments->Mxy;
+    
+    // limit axis ratio < 20.0
+    psEllipseAxes momentAxes = psEllipseMomentsToAxes (moments, 20.0);
+
+    psEllipseAxes guessAxes;
+    guessAxes.major = source->moments->Mrf;
+    guessAxes.minor = (momentAxes.minor / momentAxes.minor) * guessAxes.major;
+    guessAxes.theta = momentAxes.theta;
+
+    if (!isfinite(guessAxes.major)) return false;
+    if (!isfinite(guessAxes.minor)) return false;
+    if (!isfinite(guessAxes.theta)) return false;
+
+    // convert the major,minor,theta to shape parameters for an Reff-like model
+    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
+
+    // set the model position
+    if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
+      return false;
+    }
+
+    // sky is zero (no longer fitted, but not yet deprecated)
+    PAR[PM_PAR_SKY]  = 0.0;
+
+    // for the index loop, use Io = 1.0, use fitted values to determine Io
+    PAR[PM_PAR_I0] = 1.0;
+
+    float xMin = NAN;
+    float iMin = NAN;
+    float sMin = NAN;
+
+    // loop over index and Reff, keeping the ARatio and Theta constant?
+    // loop over index guesses and find the best fit
+    for (int i = 0; i < N_INDEX_GUESS_INV; i++) {
+        PAR[PM_PAR_7] = indexGuessInv[i];
+
+        // generated the modelFlux
+        pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
+        
+        float YY = 0.0;
+        float YM = 0.0;
+        float MM = 0.0;
+        bool usePoisson = false;
+
+        for (int iy = 0; iy < source->pixels->numRows; iy++) {
+            for (int ix = 0; ix < source->pixels->numCols; ix++) {
+                // skip masked points
+                if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
+                    continue;
+                }
+                // skip zero-variance points
+                if (source->variance->data.F32[iy][ix] == 0) {
+                    continue;
+                }
+                // skip nan value points
+                if (!isfinite(source->pixels->data.F32[iy][ix])) {
+                    continue;
+                }
+
+                float fy = source->pixels->data.F32[iy][ix];
+                float fm = source->modelFlux->data.F32[iy][ix];
+                float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0;
+
+                YY += PS_SQR(fy) * wt;
+                YM += fm * fy * wt;
+                MM += PS_SQR(fm) * wt;
+            }
+        }
+
+        float Io = YM / MM;
+        float Chisq = YY - 2 * Io * YM + Io * Io * MM;
+        if ((i == 0) || (Chisq < xMin)) {
+            xMin = Chisq;
+            iMin = Io;
+            sMin = indexGuessInv[i];
+        }
+        fprintf (stderr, "%d | %f %f %f | %f %f %f", i, indexGuessInv[i], Io, Chisq, sMin, iMin, xMin);
+        fprintf (stderr, "\n");
+    }
+
+    PAR[PM_PAR_I0] = iMin;
+    PAR[PM_PAR_7] = sMin;
+
+    return true;
+}