Changeset 14652 for trunk/psModules/src/objects/models/pmModel_PGAUSS.c

Timestamp:

Aug 23, 2007, 2:11:02 PM (19 years ago)

Author:

magnier

Message:

• added function pointers to pmModel to provide class-dependent functions
• dropped pmModel*_GetFunction functions (use pmModel->func functions instead)
• added modelParamsFromPSF functions to model classes
• changed pmModelGroup to pmModelClass
• added pmSourceFitSet.[ch]
• updated pmSourceFitSet to allow variable model classes
• added functions to add/sub and eval models & sources with an offset between image and chip
• added function to set a pmModel flux
• added function to instatiate a pmModel from a pmPSF at a coordinate
• moved pmPSF I/O to chip->analysis from readout->analysis
• changed pmPSF I/O function names from *ForPSFmodel to pmPSFmodel*
• changed pmModel.params_NEW back to pmModel.params
• changed pmFind*Peaks to pmPeaksIn* (* = Image,Vector)
• dropped pmCullPeaks (deprecated)
• changed pmModelSetType to pmModelClassGetType
• changed pmModelGetType to pmModelClassGetName
• fixed PGAUSS implementation of modelRadius function

File:

• trunk/psModules/src/objects/models/pmModel_PGAUSS.c (modified) (10 diffs)

trunk/psModules/src/objects/models/pmModel_PGAUSS.c

@@ -19 +19 @@
  *****************************************************************************/
 
-# define PM_MODEL_FUNC       pmModelFunc_PGAUSS
-# define PM_MODEL_FLUX       pmModelFlux_PGAUSS
-# define PM_MODEL_GUESS      pmModelGuess_PGAUSS
-# define PM_MODEL_LIMITS     pmModelLimits_PGAUSS
-# define PM_MODEL_RADIUS     pmModelRadius_PGAUSS
-# define PM_MODEL_FROM_PSF   pmModelFromPSF_PGAUSS
-# define PM_MODEL_FIT_STATUS pmModelFitStatus_PGAUSS
+# define PM_MODEL_FUNC            pmModelFunc_PGAUSS
+# define PM_MODEL_FLUX            pmModelFlux_PGAUSS
+# define PM_MODEL_GUESS           pmModelGuess_PGAUSS
+# define PM_MODEL_LIMITS          pmModelLimits_PGAUSS
+# define PM_MODEL_RADIUS          pmModelRadius_PGAUSS
+# define PM_MODEL_FROM_PSF        pmModelFromPSF_PGAUSS
+# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_PGAUSS
+# define PM_MODEL_FIT_STATUS      pmModelFitStatus_PGAUSS
 
 // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
@@ -66 +67 @@
 # define AR_RATIO 0.99
 
-// we constraint limits by the min valid minor axis:
-# define MIN_MINOR_AXIS 0.5
-
-// f3 = (s_b^-2 - s_a^-2); F3_SQ_MAX is MIN_MINOR_AXIS^-4
-# define F3_SQ_MAX 16.0 
-
-static float saveParams[8];
-
 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
 {
@@ -150 +143 @@
             psAbort("invalid parameter %d for param min test", nParam);
         }
-        saveParams[nParam] = params[nParam];
         if (params[nParam] < params_min) {
             params[nParam] = params_min;
@@ -184 +176 @@
             psAbort("invalid parameter %d for param max test", nParam);
         }
-        saveParams[nParam] = params[nParam];
         if (params[nParam] > params_max) {
             params[nParam] = params_max;
@@ -263 +254 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
+    psF64 z, f;
+    int Nstep = 0;
     psEllipseShape shape;
 
@@ -280 +273 @@
     // this estimates the radius assuming f(z) is roughly exp(-z)
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux));
+    psF64 sigma = axes.major;
+
+    psF64 limit = flux / PAR[PM_PAR_I0];
+
+    // use the fact that f is monotonically decreasing
+    z = 0;
+    Nstep = 0;
+
+    // choose a z value guaranteed to be beyond our limit
+    float z0 = pow((1.0 / limit), (1.0 / 3.0));
+    float z1 = (1.0 / limit);
+    z1 = PS_MAX (z0, z1);
+    z0 = 0.0;
+
+    // perform a type of bisection to find the value
+    float f0 = 1.0 / (1 + z0 + z0*z0/2.0 + z0*z0*z0/6.0);
+    float f1 = 1.0 / (1 + z1 + z1*z1/2.0 + z1*z1*z1/6.0);
+    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
+        z = 0.5*(z0 + z1);
+        f = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
+        if (f > limit) {
+            z0 = z;
+            f0 = f;
+        } else {
+            z1 = z;
+            f1 = f;
+        }
+        Nstep ++;
+    }
+    psF64 radius = sigma * sqrt (2.0 * z);
+
+    // psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux));
 
     if (isnan(radius))
@@ -297 +321 @@
 
     // we require these two parameters to exist
-    assert (psf->params_NEW->n > PM_PAR_YPOS);
-    assert (psf->params_NEW->n > PM_PAR_XPOS);
-
-    for (int i = 0; i < psf->params_NEW->n; i++) {
-        if (psf->params_NEW->data[i] == NULL) {
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
+
+    for (int i = 0; i < psf->params->n; i++) {
+        if (psf->params->data[i] == NULL) {
             out[i] = in[i];
         } else {
-            psPolynomial2D *poly = psf->params_NEW->data[i];
+            psPolynomial2D *poly = psf->params->data[i];
             out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
         }
@@ -322 +346 @@
     // apply the model limits here: this truncates excessive extrapolation
     // XXX do we need to do this still?  should we put in asserts to test?
-    for (int i = 0; i < psf->params_NEW->n; i++) {
+    for (int i = 0; i < psf->params->n; i++) {
         // apply the limits to all components or just the psf-model parameters?
-        if (psf->params_NEW->data[i] == NULL)
+        if (psf->params->data[i] == NULL)
             continue;
 
@@ -339 +363 @@
 }
 
+// construct the PSF model from the FLT model and the psf
+// XXX is this sufficiently general do be a global function, not a pmModelClass function?
+bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, pmPSF *psf, float Xo, float Yo, float Io)
+{
+    psF32 *PAR = model->params->data.F32;
+
+    // we require these two parameters to exist
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
+
+    PAR[PM_PAR_SKY]  = 0.0;
+    PAR[PM_PAR_I0]   = Io;
+    PAR[PM_PAR_XPOS] = Xo;
+    PAR[PM_PAR_YPOS] = Yo;
+    
+    // supply the model-fitted parameters, or copy from the input
+    for (int i = 0; i < psf->params->n; i++) {
+        if (i == PM_PAR_SKY) continue;
+        if (i == PM_PAR_I0) continue;
+        if (i == PM_PAR_XPOS) continue;
+        if (i == PM_PAR_YPOS) continue;
+        psPolynomial2D *poly = psf->params->data[i];
+        assert (poly);
+        PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
+    }
+
+    // the 2D PSF model fits polarization terms (E0,E1,E2)
+    // convert to shape terms (SXX,SYY,SXY)
+    // XXX user-defined value for limit?
+    if (!pmPSF_FitToModel (PAR, 0.1)) {
+        psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
+        return false;
+    }
+
+    // apply the model limits here: this truncates excessive extrapolation
+    // XXX do we need to do this still?  should we put in asserts to test?
+    for (int i = 0; i < psf->params->n; i++) {
+        // apply the limits to all components or just the psf-model parameters?
+        if (psf->params->data[i] == NULL)
+            continue;
+
+        bool status = true;
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL);
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL);
+        if (!status) {
+            psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo);
+            model->flags |= PM_MODEL_STATUS_LIMITS;
+        }
+    }
+    return(true);
+}
+
 bool PM_MODEL_FIT_STATUS (pmModel *model)
 {
@@ -367 +443 @@
 # undef PM_MODEL_RADIUS
 # undef PM_MODEL_FROM_PSF
+# undef PM_MODEL_PARAMS_FROM_PSF
 # undef PM_MODEL_FIT_STATUS