Changeset 14612 for branches/eam_branch_20070817/psModules/src/objects/models/pmModel_PGAUSS.c

Timestamp:

Aug 22, 2007, 2:55:49 PM (19 years ago)

Author:

magnier

Message:

adjusted pmModel.h to allow pmModel to use functions which use pmModel
as an argument.

adjusted the order of the pmPSF.h entries to allow that as an argument
as well

changed pmPSF I/O functions to load/save the psf on the
chip->analysis, not the readout->analysis

fixed the modelRadius function for PGAUSS

File:

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

branches/eam_branch_20070817/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