Changeset 14652

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

Location:

trunk/psModules/src/objects

Files:

• Makefile.am (modified) (3 diffs)
• models/pmModel_GAUSS.c (modified) (9 diffs)
• models/pmModel_PGAUSS.c (modified) (10 diffs)
• models/pmModel_QGAUSS.c (modified) (5 diffs)
• models/pmModel_RGAUSS.c (modified) (5 diffs)
• models/pmModel_SERSIC.c (modified) (5 diffs)
• pmGrowthCurve.c (modified) (2 diffs)
• pmModel.c (modified) (14 diffs)
• pmModel.h (modified) (7 diffs)
• pmModelClass.c (added)
• pmModelClass.h (added)
• pmModelGroup.c (modified) (3 diffs)
• pmModelGroup.h (modified) (4 diffs)
• pmModelUtils.c (modified) (3 diffs)
• pmModelUtils.h (added)
• pmPSF.c (modified) (10 diffs)
• pmPSF.h (modified) (4 diffs)
• pmPSF_IO.c (modified) (29 diffs)
• pmPSF_IO.h (modified) (2 diffs)
• pmPSFtry.c (modified) (7 diffs)
• pmPeaks.c (modified) (14 diffs)
• pmPeaks.h (modified) (5 diffs)
• pmSource.c (modified) (7 diffs)
• pmSource.h (modified) (5 diffs)
• pmSourceContour.c (modified) (2 diffs)
• pmSourceFitModel.c (modified) (7 diffs)
• pmSourceFitSet.c (added)
• pmSourceFitSet.h (added)
• pmSourceIO.c (modified) (2 diffs)
• pmSourceIO_CMP.c (modified) (3 diffs)
• pmSourceIO_OBJ.c (modified) (2 diffs)
• pmSourceIO_PS1_DEV_0.c (modified) (3 diffs)
• pmSourceIO_RAW.c (modified) (2 diffs)
• pmSourceIO_SMPDATA.c (modified) (3 diffs)
• pmSourceIO_SX.c (modified) (2 diffs)
• pmSourcePhotometry.c (modified) (5 diffs)
• pmSourcePlotApResid.c (modified) (2 diffs)
• pmSourcePlotMoments.c (modified) (2 diffs)
• pmSourcePlotPSFModel.c (modified) (2 diffs)
• pmSourceSky.c (modified) (2 diffs)
• pmSourceUtils.c (added)
• pmSourceUtils.h (added)

trunk/psModules/src/objects/Makefile.am

@@ -7 +7 @@
      pmMoments.c \
      pmModel.c \
-     pmModelGroup.c \
+     pmModelClass.c \
+     pmModelUtils.c \
      pmSource.c \
+     pmSourceUtils.c \
      pmSourceSky.c \
      pmSourceContour.c \
      pmSourceFitModel.c \
+     pmSourceFitSet.c \
      pmSourcePhotometry.c \
      pmSourceIO.c \
@@ -31 +34 @@
 
 EXTRA_DIST = \
-        models/pmModel_GAUSS.c \
-        models/pmModel_PGAUSS.c \
-        models/pmModel_QGAUSS.c \
-        models/pmModel_SGAUSS.c \
-        models/pmModel_RGAUSS.c \
-        models/pmModel_SERSIC.c
+     models/pmModel_GAUSS.c \
+     models/pmModel_PGAUSS.c \
+     models/pmModel_QGAUSS.c \
+     models/pmModel_SGAUSS.c \
+     models/pmModel_RGAUSS.c \
+     models/pmModel_SERSIC.c
 
 pkginclude_HEADERS = \
@@ -42 +45 @@
      pmMoments.h \
      pmModel.h \
-     pmModelGroup.h \
+     pmModelClass.h \
+     pmModelUtils.h \
      pmSource.h \
+     pmSourceUtils.h \
      pmSourceSky.h \
      pmSourceContour.h \
      pmSourceFitModel.h \
+     pmSourceFitSet.h \
      pmSourcePhotometry.h \
      pmSourceIO.h \

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

@@ -19 +19 @@
  *****************************************************************************/
 
-# define PM_MODEL_FUNC       pmModelFunc_GAUSS
-# define PM_MODEL_FLUX       pmModelFlux_GAUSS
-# define PM_MODEL_GUESS      pmModelGuess_GAUSS
-# define PM_MODEL_LIMITS     pmModelLimits_GAUSS
-# define PM_MODEL_RADIUS     pmModelRadius_GAUSS
-# define PM_MODEL_FROM_PSF   pmModelFromPSF_GAUSS
-# define PM_MODEL_FIT_STATUS pmModelFitStatus_GAUSS
-
+# define PM_MODEL_FUNC            pmModelFunc_GAUSS
+# define PM_MODEL_FLUX            pmModelFlux_GAUSS
+# define PM_MODEL_GUESS           pmModelGuess_GAUSS
+# define PM_MODEL_LIMITS          pmModelLimits_GAUSS
+# define PM_MODEL_RADIUS          pmModelRadius_GAUSS
+# define PM_MODEL_FROM_PSF        pmModelFromPSF_GAUSS
+# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_GAUSS
+# define PM_MODEL_FIT_STATUS      pmModelFitStatus_GAUSS
+
+// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
 psF32 PM_MODEL_FUNC(psVector *deriv,
                     const psVector *params,
@@ -38 +40 @@
     psF32 py = Y / PAR[PM_PAR_SYY];
     psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
+    assert (z >= 0.0);
+
     psF32 r  = exp(-z);
     psF32 q  = PAR[PM_PAR_I0]*r;
@@ -61 +65 @@
 # define AR_MAX 20.0
 # define AR_RATIO 0.99
+
 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
 {
@@ -93 +98 @@
             break;
         case PM_PAR_SXX:
-            beta_lim = 0.5;
+            beta_lim = 2.0;
             break;
         case PM_PAR_SYY:
-            beta_lim = 0.5;
+            beta_lim = 2.0;
             break;
         case PM_PAR_SXY:
@@ -257 +262 @@
 bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)
 {
-
     psF32 *out = modelPSF->params->data.F32;
     psF32 *in  = modelFLT->params->data.F32;
 
     // we require these two parameters to exist
-    assert (psf->params_NEW->n > PM_PAR_YPOS);
-    assert (psf->params_NEW->n > PM_PAR_XPOS);
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
 
     // supply the model-fitted parameters, or copy from the input
-    for (int i = 0; i < psf->params_NEW->n; i++) {
-        if (psf->params_NEW->data[i] == NULL) {
+    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]);
         }
@@ -289 +293 @@
     // 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;
 
@@ -306 +310 @@
 }
 
+// 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);
+}
+
 // check the status of the fitted model
 // this test is invalid if the parameters are derived
@@ -311 +367 @@
 bool PM_MODEL_FIT_STATUS (pmModel *model)
 {
-
     psF32 dP;
     bool  status;
@@ -339 +394 @@
 # undef PM_MODEL_RADIUS
 # undef PM_MODEL_FROM_PSF
+# undef PM_MODEL_PARAMS_FROM_PSF
 # undef PM_MODEL_FIT_STATUS

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

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

@@ -20 +20 @@
    *****************************************************************************/
 
-# define PM_MODEL_FUNC       pmModelFunc_QGAUSS
-# define PM_MODEL_FLUX       pmModelFlux_QGAUSS
-# define PM_MODEL_GUESS      pmModelGuess_QGAUSS
-# define PM_MODEL_LIMITS     pmModelLimits_QGAUSS
-# define PM_MODEL_RADIUS     pmModelRadius_QGAUSS
-# define PM_MODEL_FROM_PSF   pmModelFromPSF_QGAUSS
-# define PM_MODEL_FIT_STATUS pmModelFitStatus_QGAUSS
+# define PM_MODEL_FUNC            pmModelFunc_QGAUSS
+# define PM_MODEL_FLUX            pmModelFlux_QGAUSS
+# define PM_MODEL_GUESS           pmModelGuess_QGAUSS
+# define PM_MODEL_LIMITS          pmModelLimits_QGAUSS
+# define PM_MODEL_RADIUS          pmModelRadius_QGAUSS
+# define PM_MODEL_FROM_PSF        pmModelFromPSF_QGAUSS
+# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_QGAUSS
+# define PM_MODEL_FIT_STATUS      pmModelFitStatus_QGAUSS
 
 psF32 PM_MODEL_FUNC (psVector *deriv,
@@ -350 +351 @@
 
     // 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]);
         }
@@ -372 +373 @@
     // 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;
 
@@ -390 +391 @@
 }
 
+// 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)
 {
@@ -418 +471 @@
 # undef PM_MODEL_RADIUS
 # undef PM_MODEL_FROM_PSF
+# undef PM_MODEL_PARAMS_FROM_PSF
 # undef PM_MODEL_FIT_STATUS

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

@@ -20 +20 @@
  *****************************************************************************/
 
-# define PM_MODEL_FUNC       pmModelFunc_RGAUSS
-# define PM_MODEL_FLUX       pmModelFlux_RGAUSS
-# define PM_MODEL_GUESS      pmModelGuess_RGAUSS
-# define PM_MODEL_LIMITS     pmModelLimits_RGAUSS
-# define PM_MODEL_RADIUS     pmModelRadius_RGAUSS
-# define PM_MODEL_FROM_PSF   pmModelFromPSF_RGAUSS
-# define PM_MODEL_FIT_STATUS pmModelFitStatus_RGAUSS
+# define PM_MODEL_FUNC            pmModelFunc_RGAUSS
+# define PM_MODEL_FLUX            pmModelFlux_RGAUSS
+# define PM_MODEL_GUESS           pmModelGuess_RGAUSS
+# define PM_MODEL_LIMITS          pmModelLimits_RGAUSS
+# define PM_MODEL_RADIUS          pmModelRadius_RGAUSS
+# define PM_MODEL_FROM_PSF        pmModelFromPSF_RGAUSS
+# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_RGAUSS
+# define PM_MODEL_FIT_STATUS      pmModelFitStatus_RGAUSS
 
 psF32 PM_MODEL_FUNC (psVector *deriv,
@@ -343 +344 @@
 
     // 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]);
         }
@@ -365 +366 @@
     // 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;
 
@@ -383 +384 @@
 }
 
+// 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)
 {
@@ -411 +464 @@
 # undef PM_MODEL_RADIUS
 # undef PM_MODEL_FROM_PSF
+# undef PM_MODEL_PARAMS_FROM_PSF
 # undef PM_MODEL_FIT_STATUS

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

@@ -23 +23 @@
    *****************************************************************************/
 
-# define PM_MODEL_FUNC       pmModelFunc_SERSIC
-# define PM_MODEL_FLUX       pmModelFlux_SERSIC
-# define PM_MODEL_GUESS      pmModelGuess_SERSIC
-# define PM_MODEL_LIMITS     pmModelLimits_SERSIC
-# define PM_MODEL_RADIUS     pmModelRadius_SERSIC
-# define PM_MODEL_FROM_PSF   pmModelFromPSF_SERSIC
-# define PM_MODEL_FIT_STATUS pmModelFitStatus_SERSIC
+# define PM_MODEL_FUNC            pmModelFunc_SERSIC
+# define PM_MODEL_FLUX            pmModelFlux_SERSIC
+# define PM_MODEL_GUESS           pmModelGuess_SERSIC
+# define PM_MODEL_LIMITS          pmModelLimits_SERSIC
+# define PM_MODEL_RADIUS          pmModelRadius_SERSIC
+# define PM_MODEL_FROM_PSF        pmModelFromPSF_SERSIC
+# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_SERSIC
+# define PM_MODEL_FIT_STATUS      pmModelFitStatus_SERSIC
 
 psF32 PM_MODEL_FUNC (psVector *deriv,
@@ -360 +361 @@
 
     // 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]);
         }
@@ -382 +383 @@
     // 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;
 
@@ -400 +401 @@
 }
 
+// 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)
 {
@@ -428 +481 @@
 # undef PM_MODEL_RADIUS
 # undef PM_MODEL_FROM_PSF
+# undef PM_MODEL_PARAMS_FROM_PSF
 # undef PM_MODEL_FIT_STATUS

trunk/psModules/src/objects/pmGrowthCurve.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -23 +23 @@
 #include "pmMoments.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
 #include "pmSource.h"
-#include "pmGrowthCurve.h"
-#include "pmPSF.h"
 #include "psVectorBracket.h"
 

trunk/psModules/src/objects/pmModel.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-20 02:22:26 $
+ *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -21 +21 @@
 #include <string.h>
 #include <pslib.h>
+#include "pmHDU.h"
+#include "pmFPA.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-
-/* The following types and functions need to be known by pmModel.c and are defined in
-   pmModelGroup.h.  The use of pmSource and other types in pmModelGroup.h prevent us from
-   directly including it here  ***/
-
-typedef psMinimizeLMChi2Func pmModelFunc;
-typedef bool (*pmModelFitStatusFunc)(pmModel *model);
-pmModelFunc pmModelFunc_GetFunction (pmModelType type);
-pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type);
-int pmModelParameterCount(pmModelType type);
+#include "pmModelClass.h"
 
 static void modelFree(pmModel *tmp)
@@ -45 +40 @@
 pmModelAlloc(): Allocate the pmModel structure, along with its parameters,
 and initialize the type member.  Initialize the params to 0.0.
-XXX EAM: simplifying code with pmModelParameterCount
 *****************************************************************************/
 pmModel *pmModelAlloc(pmModelType type)
@@ -53 +47 @@
     pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel));
     psMemSetDeallocator(tmp, (psFreeFunc) modelFree);
+
+    pmModelClass *class = pmModelClassSelect (type);
+    if (class == NULL) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return(NULL);
+    }
 
     tmp->type = type;
@@ -62 +62 @@
     tmp->residuals = NULL;              // XXX should the model own this memory?
 
-    psS32 Nparams = pmModelParameterCount(type);
-    if (Nparams == 0) {
-        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
-        return(NULL);
-    }
+    psS32 Nparams = pmModelClassParameterCount(type);
+    assert (Nparams);
 
     tmp->params  = psVectorAlloc(Nparams, PS_TYPE_F32);
     tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32);
+    assert (tmp->params);
+    assert (tmp->dparams);
 
     for (psS32 i = 0; i < tmp->params->n; i++) {
@@ -75 +74 @@
         tmp->dparams->data.F32[i] = 0.0;
     }
+
+    tmp->modelFunc          = class->modelFunc;
+    tmp->modelFlux          = class->modelFlux;
+    tmp->modelRadius        = class->modelRadius;
+    tmp->modelLimits        = class->modelLimits;
+    tmp->modelGuess         = class->modelGuess;
+    tmp->modelFromPSF       = class->modelFromPSF;
+    tmp->modelParamsFromPSF = class->modelParamsFromPSF;
+    tmp->modelFitStatus     = class->modelFitStatus;
 
     psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
@@ -92 +100 @@
     new->radiusFit = model->radiusFit;
 
-    // XXX note that model->residuals is just a reference
-    new->residuals = model->residuals;
-
     for (int i = 0; i < new->params->n; i++) {
         new->params->data.F32[i]  = model->params->data.F32[i];
         new->dparams->data.F32[i] = model->dparams->data.F32[i];
     }
+
+    // note that model->residuals is just a reference
+    new->residuals = model->residuals;
 
     return (new);
@@ -123 +131 @@
     x->data.F32[1] = (psF32) (row + image->row0);
     psF32 tmpF;
-    pmModelFunc modelFunc;
-
-    modelFunc = pmModelFunc_GetFunction (model->type);
-    tmpF = modelFunc (NULL, model->params, x);
+
+    tmpF = model->modelFunc (NULL, model->params, x);
+    psFree(x);
+    psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
+    return(tmpF);
+}
+
+psF32 pmModelEvalWithOffset(pmModel *model, psImage *image, psS32 col, psS32 row, int dx, int dy)
+{
+    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(image, false);
+    PS_ASSERT_PTR_NON_NULL(model, false);
+    PS_ASSERT_PTR_NON_NULL(model->params, false);
+
+    // Allocate the x coordinate structure and convert row/col to image space.
+    //
+    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
+    x->data.F32[0] = (psF32) (col + image->col0 + dx);
+    x->data.F32[1] = (psF32) (row + image->row0 + dy);
+    psF32 tmpF;
+
+    tmpF = model->modelFunc (NULL, model->params, x);
     psFree(x);
     psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
@@ -137 +163 @@
                           pmModelOpMode mode,
                           bool add,
-                          psMaskType maskVal
+                          psMaskType maskVal,
+                          int dx,
+                          int dy
     )
 {
@@ -148 +176 @@
     psVector *x = psVectorAlloc(2, PS_TYPE_F32);
     psVector *params = model->params;
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+
     psS32 imageCol;
     psS32 imageRow;
@@ -209 +237 @@
             // Convert i/j to image coord space:
             // XXX should we use using 0.5 pixel offset?
-            imageCol = ix + image->col0;
-            imageRow = iy + image->row0;
+            imageCol = ix + image->col0 + dx;
+            imageRow = iy + image->row0 + dy;
 
             x->data.F32[0] = (float) imageCol;
@@ -219 +247 @@
             // add in the desired components for this coordinate
             if (mode & PM_MODEL_OP_FUNC) {
-                pixelValue += modelFunc (NULL, params, x);
+                pixelValue += model->modelFunc (NULL, params, x);
             }
 
@@ -280 +308 @@
 {
     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
-    psBool rc = AddOrSubModel(image, mask, model, mode, true, maskVal);
+    psBool rc = AddOrSubModel(image, mask, model, mode, true, maskVal, 0.0, 0.0);
     psTrace("psModules.objects", 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);
@@ -294 +322 @@
 {
     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
-    psBool rc = AddOrSubModel(image, mask, model, mode, false, maskVal);
+    psBool rc = AddOrSubModel(image, mask, model, mode, false, maskVal, 0.0, 0.0);
     psTrace("psModules.objects", 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);
 }
 
-// check for success of model fit
-bool pmModelFitStatus (pmModel *model)
-{
-
-    bool status;
-
-    pmModelFitStatusFunc statusFunc = pmModelFitStatusFunc_GetFunction (model->type);
-    status = statusFunc (model);
-
-    return (status);
-}
-
+/******************************************************************************
+ *****************************************************************************/
+bool pmModelAddWithOffset(psImage *image,
+                          psImage *mask,
+                          pmModel *model,
+                          pmModelOpMode mode,
+                          psMaskType maskVal,
+                          int dx,
+                          int dy)
+{
+    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
+    psBool rc = AddOrSubModel(image, mask, model, mode, true, maskVal, dx, dy);
+    psTrace("psModules.objects", 3, "---- %s(%d) end ----\n", __func__, rc);
+    return(rc);
+}
+
+/******************************************************************************
+ *****************************************************************************/
+bool pmModelSubWithOffset(psImage *image,
+                          psImage *mask,
+                          pmModel *model,
+                          pmModelOpMode mode,
+                          psMaskType maskVal,
+                          int dx,
+                          int dy)
+{
+    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
+    psBool rc = AddOrSubModel(image, mask, model, mode, false, maskVal, dx, dy);
+    psTrace("psModules.objects", 3, "---- %s(%d) end ----\n", __func__, rc);
+    return(rc);
+}

trunk/psModules/src/objects/pmModel.h

@@ -1 @@
-/* @file  pmObjects.h
+/* @file  pmModel.h
  * @brief Functions to define and manipulate object models
  *
@@ -5 +5 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-06-20 02:22:26 $
+ * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
  *
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -17 +17 @@
 /// @{
 
-// type of model carried by the pmModel structure
-typedef int pmModelType;
+/* pointers for the functions types below are supplied to each pmModel, and can be used by
+   the programmer without needing to know the model class */
 
 typedef enum {
@@ -41 +41 @@
 } pmModelOpMode;
 
+typedef struct pmModel pmModel;
+typedef struct pmSource pmSource;
+
+//  This function is the model chi-square minimization function for this model.
+typedef psMinimizeLMChi2Func pmModelFunc;
+
+//  This function sets the parameter limits for this model.
+typedef psMinimizeLMLimitFunc pmModelLimits;
+
+// This function returns the integrated flux for the given model parameters.
+typedef psF64 (*pmModelFlux)(const psVector *params);
+
+// This function returns the radius at which the given model and parameters
+// achieves the given flux.
+typedef psF64 (*pmModelRadius)(const psVector *params, double flux);
+
+//  This function provides the model guess parameters based on the details of
+//  the given source.
+typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source);
+
+//  This function constructs the PSF model for the given source based on the
+//  supplied psf and the EXT model for the object.
+typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, pmPSF *psf);
+
+//  This function sets the model parameters based on the PSF for a given coordinate and central
+//  intensity
+typedef bool (*pmModelParamsFromPSF)(pmModel *model, pmPSF *psf, float Xo, float Yo, float Io);
+
+//  This function returns the success / failure status of the given model fit
+typedef bool (*pmModelFitStatusFunc)(pmModel *model);
+
 /** pmModel data structure
  *
@@ -50 +81 @@
  *
  */
-typedef struct
-{
+struct pmModel {
     pmModelType type;                   ///< Model to be used.
     psVector *params;                   ///< Paramater values.
@@ -62 +92 @@
     float radiusFit;                    ///< fit radius actually used
     pmResiduals *residuals;             ///< normalized PSF residuals
-}
-pmModel;
 
-/* XXX we are currently saving the residuals with the pmModel.  It might be better to save this
- * in the pmSource.  we may want an API to Add/Sub a pmModel (analytical model only) or a
- * pmSource (analytical + residuals).
- */
+    // functions for this model which depend on the model class
+    pmModelFunc          modelFunc;
+    pmModelFlux          modelFlux;
+    pmModelRadius        modelRadius;
+    pmModelLimits        modelLimits;
+    pmModelGuessFunc     modelGuess;
+    pmModelFromPSFFunc   modelFromPSF;
+    pmModelParamsFromPSF modelParamsFromPSF;
+    pmModelFitStatusFunc modelFitStatus;
+};
 
 /** Symbolic names for the elements of [d]params
@@ -130 +164 @@
 );
 
+bool pmModelAddWithOffset(psImage *image,
+                          psImage *mask,
+                          pmModel *model,
+                          pmModelOpMode mode,
+                          psMaskType maskVal,
+                          int dx,
+                          int dy);
+
+bool pmModelSubWithOffset(psImage *image,
+                          psImage *mask,
+                          pmModel *model,
+                          pmModelOpMode mode,
+                          psMaskType maskVal,
+                          int dx,
+                          int dy);
+
 /** pmModelFitStatus()
  *

trunk/psModules/src/objects/pmModelGroup.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-20 00:31:43 $
+ *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -28 +28 @@
 #include "pmResiduals.h"
 #include "pmModel.h"
-#include "pmResiduals.h"
-#include "pmPSF.h"
-#include "pmSource.h"
 #include "pmModelGroup.h"
 #include "pmErrorCodes.h"
@@ -198 +195 @@
     return (models[type].name);
 }
-
-/******************************************************************************
-    pmSourceModelGuess(source, model): This function allocates a new
-    pmModel structure based on the given modelType specified in the argument list.
-    The corresponding pmModelGuess function is returned, and used to
-    supply the values of the params array in the pmModel structure.
- 
-    XXX: Many parameters are based on the src->moments structure, which is in
-    image, not subImage coords.  Therefore, the calls to the model evaluation
-    functions will be in image, not subImage coords.  Remember this.
-*****************************************************************************/
-pmModel *pmSourceModelGuess(pmSource *source,
-                            pmModelType modelType)
-{
-    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(source, NULL);
-    PS_ASSERT_PTR_NON_NULL(source->moments, NULL);
-    PS_ASSERT_PTR_NON_NULL(source->peak, NULL);
-
-    pmModel *model = pmModelAlloc(modelType);
-
-    pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType);
-    if (!modelGuessFunc(model, source)) {
-        psFree (model);
-        return NULL;
-    }
-
-    psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
-    return(model);
-}
-

trunk/psModules/src/objects/pmModelGroup.h

@@ -1 +1 @@
 /* @file  pmModelGroup.h
  *
- * The object model function types are defined to allow for the flexible addition
- * of new object models. Every object model, with parameters represented by
- * pmModel, has an associated set of functions which provide necessary support
- * operations. A set of abstract functions allow the programmer to select the
- * approriate function or property for a specific named object model.
+ * The object model function types are desined to allow for the flexible addition of new object
+ * models. Every object model, with parameters represented by pmModel, has an associated set of
+ * functions which provide necessary support operations. A set of abstract functions allow the
+ * programmer to select the approriate function or property for a specific named object model.
+ *
+ * Every model instance belongs to a class of models, defined by the value of
+ * the pmModelType type entry. Various functions need access to information about
+ * each of the models. Some of this information varies from model to model, and
+ * may depend on the current parameter values or other data quantities. In order
+ * to keep the code from requiring the information about each model to be coded
+ * into the low-level fitting routines, we define a collection of functions which
+ * allow us to abstract this type of model-dependent information. These generic
+ * functions take the model type and return the corresponding function pointer
+ * for the specified model. Each model is defined by creating this collection of
+ * specific functions, and placing them in a single file for each model. We
+ * define the following structure to carry the collection of information about
+ * the models.
  *
  * @author EAM, IfA
  *
- * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-01-24 02:54:15 $
+ * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -23 +35 @@
 typedef psMinimizeLMChi2Func pmModelFunc;
 
-//  This function is the model chi-square minimization function for this model.
+//  This function sets the parameter limits for this model.
 typedef psMinimizeLMLimitFunc pmModelLimits;
 
@@ -29 +41 @@
 typedef psF64 (*pmModelFlux)(const psVector *params);
 
-
 // This function returns the radius at which the given model and parameters
 // achieves the given flux.
 typedef psF64 (*pmModelRadius)(const psVector *params, double flux);
 
-/*  This function sets the model parameter limits vectors for the given model
- */
-// typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max);
-
-/*  This function provides the model guess parameters based on the details of
- *   the given source.
- */
+//  This function provides the model guess parameters based on the details of
+//  the given source.
 typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source);
 
-
-/*  This function constructs the PSF model for the given source based on the
- *  supplied psf and the EXT model for the object.
- */
+//  This function constructs the PSF model for the given source based on the
+//  supplied psf and the EXT model for the object.
 typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, pmPSF *psf);
 
-/*  This function returns the success / failure status of the given model fit
- */
+//  This function sets the model parameters based on the PSF for a given coordinate and central
+//  intensity
+typedef bool (*pmModelParamsFromPSF)(pmModel *model, pmPSF *psf, float Xo, float Yo, float Io);
+
+//  This function returns the success / failure status of the given model fit
 typedef bool (*pmModelFitStatusFunc)(pmModel *model);
 
-/* Every model instance belongs to a class of models, defined by the value of
- * the pmModelType type entry. Various functions need access to information about
- * each of the models. Some of this information varies from model to model, and
- * may depend on the current parameter values or other data quantities. In order
- * to keep the code from requiring the information about each model to be coded
- * into the low-level fitting routines, we define a collection of functions which
- * allow us to abstract this type of model-dependent information. These generic
- * functions take the model type and return the corresponding function pointer
- * for the specified model. Each model is defined by creating this collection of
- * specific functions, and placing them in a single file for each model. We
- * define the following structure to carry the collection of information about
- * the models.
- */
 typedef struct
 {
@@ -76 +70 @@
     pmModelGuessFunc     modelGuessFunc;
     pmModelFromPSFFunc   modelFromPSFFunc;
+    pmModelParamsFromPSF modelParamsFromPSF;
     pmModelFitStatusFunc modelFitStatusFunc;
 }

trunk/psModules/src/objects/pmModelUtils.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-08-16 18:33:37 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -20 +20 @@
 #include <string.h>
 #include <pslib.h>
+#include "pmHDU.h"
+#include "pmFPA.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
+#include "pmErrorCodes.h"
+#include "pmModelUtils.h"
 
-// instantiate a model for the PSF at this location (normalized or not?)
-// NOTE: psf and (x,y) are defined wrt chip coordinates
-pmModel *pmModelFromPSFforXY (pmPSF *psf, float x, float y, float flux) {
-
-    assert (psf);
-
-    // allocate a new pmModel to hold the PSF version
-    pmModel *modelEXT = pmModelAlloc (psf->type);
-
-    modelEXT->params->data.F32[PM_PAR_SKY]  = 0.0;
-    modelEXT->params->data.F32[PM_PAR_I0]   = 1.0;
-    modelEXT->params->data.F32[PM_PAR_XPOS] = x;
-    modelEXT->params->data.F32[PM_PAR_YPOS] = y;
-
-    // find function used to set the model parameters
-    pmModelFromPSFFunc modelFromPSFFunc = pmModelFromPSFFunc_GetFunction (psf->type);
-
+/*****************************************************************************
+pmModelFromPSF (*modelEXT, *psf):  use the model position parameters to
+construct a realization of the PSF model at the object coordinates
+ *****************************************************************************/
+pmModel *pmModelFromPSF (pmModel *modelEXT, pmPSF *psf)
+{
     // allocate a new pmModel to hold the PSF version
     pmModel *modelPSF = pmModelAlloc (psf->type);
 
-    // adjust the normalization: 
-    pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
-    normFlux = modelFluxFunc (model->params);
-    assert (isfinite(normFlux));
-    assert (normFlux > 0);
-
-    // set the correct normalization
-    modelEXT->params->data.F32[PM_PAR_I0] = flux / normFlux;
-
     // set model parameters for this source based on PSF information
-    if (!modelFromPSFFunc (modelPSF, modelEXT, psf)) {
+    if (!modelEXT->modelFromPSF (modelPSF, modelEXT, psf)) {
         psError(PM_ERR_PSF, false, "Failed to set model params from PSF");
         psFree(modelPSF);
@@ -60 +46 @@
     // XXX note that model->residuals is just a reference
     modelPSF->residuals = psf->residuals;
-    psFree (modelEXT);
 
     return (modelPSF);
 }
 
-pmSource *pmSourceFromModel (pmModel *model, pmReadout *readout, pmSourceType type) {
+// instantiate a model for the PSF at this location with peak flux
+// NOTE: psf and (Xo,Yo) are defined wrt chip coordinates
+pmModel *pmModelFromPSFforXY (pmPSF *psf, float Xo, float Yo, float Io) {
 
-    pmSource *source = pmSourceAlloc ();
+    assert (psf);
 
-    // use the model centroid for the peak
-    switch (type) {
-      case PM_SOURCE_TYPE_STAR:
-        source->modelPSF = model;
-        break;
-      case PM_SOURCE_TYPE_EXTENDED:
-        source->modelEXT = model;
-        break;
-      default:
-        psAbort ("psphot", "error");
+    // allocate a new pmModel to hold the PSF version
+    pmModel *modelPSF = pmModelAlloc (psf->type);
+
+    // set model parameters for this source based on PSF information
+    if (!modelPSF->modelParamsFromPSF (modelPSF, psf, Xo, Yo, Io)) {
+        psError(PM_ERR_PSF, false, "Failed to set model params from PSF");
+        psFree(modelPSF);
+        return NULL;
     }
 
-    source->peak = pmPeakAlloc ();
+    // XXX note that model->residuals is just a reference
+    modelPSF->residuals = psf->residuals;
 
-    float x = model->params->data.F32[PM_PAR_XPOS];
-    float y = model->params->data.F32[PM_PAR_YPOS];
+    return (modelPSF);
+}
 
-    // XXX need to define the radius in some rational way
-    // XXX x,y are defined wrt readout->image parent, but the model
-    // parameters are defined wrt chip coordinates
-    pmSourceDefinePixels (source, readout, x, y, radius);
+// set this model to have the requested flux
+bool pmModelSetFlux (pmModel *model, float flux) {
 
-    return (source);
+    // set Io to be 1.0
+    model->params->data.F32[PM_PAR_I0] = 1.0;
+
+    // determine the normalized flux
+    float normFlux = model->modelFlux (model->params);
+    assert (isfinite(normFlux));
+    assert (normFlux > 0);
+
+    // set the desired normalization
+    model->params->data.F32[PM_PAR_I0] = flux / normFlux;
+
+    return true;
 }

trunk/psModules/src/objects/pmPSF.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-20 02:22:26 $
+ *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -28 +28 @@
 #include "pmMoments.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
 #include "pmSource.h"
-#include "pmGrowthCurve.h"
-#include "pmPSF.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
+#include "pmModelUtils.h"
 #include "pmSourcePhotometry.h"
 #include "pmFPAMaskWeight.h"
@@ -74 +75 @@
     psFree (psf->ApTrend);
     psFree (psf->growth);
-    psFree (psf->params_NEW);
+    psFree (psf->params);
     psFree (psf->residuals);
     return;
@@ -103 +104 @@
     psf->skyBias  = 0.0;
     psf->skySat   = 0.0;
+    psf->nPSFstars  = 0;
+    psf->nApResid   = 0;
     psf->poissonErrors = poissonErrors;
 
@@ -124 +127 @@
     psf->residuals = NULL;
 
-    Nparams = pmModelParameterCount (type);
+    Nparams = pmModelClassParameterCount (type);
     if (!Nparams) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
         return(NULL);
     }
-    psf->params_NEW = psArrayAlloc(Nparams);
+    psf->params = psArrayAlloc(Nparams);
 
     // the order of the PSF parameter (X,Y) fits is determined by the psfTrendMask polynomial
@@ -139 +142 @@
 
     if (psfTrendMask) {
-        for (int i = 0; i < psf->params_NEW->n; i++) {
+        for (int i = 0; i < psf->params->n; i++) {
             if (i == PM_PAR_SKY)
                 continue;
@@ -155 +158 @@
                 }
             }
-            psf->params_NEW->data[i] = param;
+            psf->params->data[i] = param;
         }
     }
@@ -161 +164 @@
     psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree);
     return(psf);
-}
-
-/*****************************************************************************
-pmModelFromPSF (*modelEXT, *psf):  use the model position parameters to
-construct a realization of the PSF model at the object coordinates
- *****************************************************************************/
-pmModel *pmModelFromPSF (pmModel *modelEXT, pmPSF *psf)
-{
-
-    // need to define the relationship between the modelEXT and modelPSF ?
-
-    // find function used to set the model parameters
-    pmModelFromPSFFunc modelFromPSFFunc = pmModelFromPSFFunc_GetFunction (psf->type);
-
-    // allocate a new pmModel to hold the PSF version
-    pmModel *modelPSF = pmModelAlloc (psf->type);
-
-    // set model parameters for this source based on PSF information
-    if (!modelFromPSFFunc (modelPSF, modelEXT, psf)) {
-        psError(PM_ERR_PSF, false, "Failed to set model params from PSF");
-        psFree(modelPSF);
-        return NULL;
-    }
-    // XXX note that model->residuals is just a reference
-    modelPSF->residuals = psf->residuals;
-
-    return (modelPSF);
 }
 
@@ -335 +311 @@
     va_start(ap, sxy);
 
-    pmModelType type = pmModelSetType (typeName);
+    pmModelType type = pmModelClassGetType (typeName);
     psPolynomial2D *psfTrend = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 0, 0);
     pmPSF *psf = pmPSFAlloc (type, true, psfTrend);
 
-    psVector *par = psVectorAlloc (psf->params_NEW->n, PS_TYPE_F32);
+    psVector *par = psVectorAlloc (psf->params->n, PS_TYPE_F32);
     par->data.F32[PM_PAR_SXX] = sxx;
     par->data.F32[PM_PAR_SYY] = syy;
@@ -348 +324 @@
     // set the psf shape parameters
     psPolynomial2D *poly = NULL;
-    poly = psf->params_NEW->data[PM_PAR_E0];
+    poly = psf->params->data[PM_PAR_E0];
     poly->coeff[0][0] = pol.e0;
 
-    poly = psf->params_NEW->data[PM_PAR_E1];
+    poly = psf->params->data[PM_PAR_E1];
     poly->coeff[0][0] = pol.e1;
 
-    poly = psf->params_NEW->data[PM_PAR_E2];
+    poly = psf->params->data[PM_PAR_E2];
     poly->coeff[0][0] = pol.e2;
 
-    for (int i = PM_PAR_SXY + 1; i < psf->params_NEW->n; i++) {
-        poly = psf->params_NEW->data[i];
+    for (int i = PM_PAR_SXY + 1; i < psf->params->n; i++) {
+        poly = psf->params->data[i];
         poly->coeff[0][0] = (psF32)va_arg(ap, psF64);
     }

trunk/psModules/src/objects/pmPSF.h

@@ -6 +6 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-06-20 02:22:26 $
+ * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -16 +16 @@
 /// @addtogroup Objects Object Detection / Analysis Functions
 /// @{
+
+// type of model carried by the pmModel structure
+typedef int pmModelType;
 
 typedef enum {
@@ -44 +47 @@
 typedef struct
 {
-    pmModelType type;   ///< PSF Model in use
-    psArray *params_NEW;   ///< Model parameters (psPolynomial2D)
-    // XXXXX I am changing params: we will allocate elements for the
-    // unfitted elements (So, Io, Xo, Yo) and leave them as NULL
-    // I am using a new name to catch all refs to params with gcc
+    pmModelType type;                   ///< PSF Model in use
+    psArray *params;                    ///< Model parameters (psPolynomial2D)
     psPolynomial1D *ChiTrend;         ///< Chisq vs flux fit (correction for systematic errors)
     psPolynomial4D *ApTrend;            ///< ApResid vs (x,y,rflux) (rflux = ten(0.4*mInst))
@@ -79 +79 @@
 );
 
-/**
- *
- * This function constructs a pmModel instance based on the pmPSF description
- * of the PSF. The input is a pmModel with at least the values of the centroid
- * coordinates (possibly normalization if this is needed) defined. The values of
- * the PSF-dependent parameters are specified for the specific realization based
- * on the coordinates of the object.
- *
- */
-pmModel *pmModelFromPSF(
-    pmModel *model,                     ///< Add comment
-    pmPSF *psf                          ///< Add comment
-);
-
 bool pmPSFMaskApTrend (psPolynomial4D *trend, pmPSFApTrendOptions option);
 

trunk/psModules/src/objects/pmPSF_IO.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-14 03:20:22 $
+ *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -40 +40 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmPSF_IO.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 
-bool pmFPAviewCheckDataStatusForPSFmodel (const pmFPAview *view, const pmFPAfile *file)
+bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file)
 {
     pmFPA *fpa = file->fpa;
 
     if (view->chip == -1) {
-        bool exists = pmFPACheckDataStatusForPSFmodel (fpa);
+        bool exists = pmPSFmodelCheckDataStatusForFPA (fpa);
         return exists;
     }
@@ -62 +62 @@
 
     if (view->cell == -1) {
-        bool exists = pmChipCheckDataStatusForPSFmodel (chip);
+        bool exists = pmPSFmodelCheckDataStatusForChip (chip);
         return exists;
     }
@@ -69 +69 @@
         return false;
     }
-    pmCell *cell = chip->cells->data[view->cell];
-
-    if (view->readout == -1) {
-        bool exists = pmCellCheckDataStatusForPSFmodel (cell);
-        return exists;
-    }
-
-    if (view->readout >= cell->readouts->n) {
-        psError(PS_ERR_IO, true, "Requested readout == %d >= cell->readouds->n == %ld", view->readout, cell->readouts->n);
-        return false;
-    }
-    pmReadout *readout = cell->readouts->data[view->readout];
-
-    return pmReadoutCheckDataStatusForPSFmodel (readout);
-}
-
-bool pmFPACheckDataStatusForPSFmodel (const pmFPA *fpa) {
+    psError(PS_ERR_IO, false, "PSF only valid at the chip level");
+    return false;
+}
+
+bool pmPSFmodelCheckDataStatusForFPA (const pmFPA *fpa) {
 
     for (int i = 0; i < fpa->chips->n; i++) {
         pmChip *chip = fpa->chips->data[i];
         if (!chip) continue;
-        if (pmChipCheckDataStatusForPSFmodel (chip)) return true;
+        if (pmPSFmodelCheckDataStatusForChip (chip)) return true;
     }
     return false;
 }
 
-bool pmChipCheckDataStatusForPSFmodel (const pmChip *chip) {
-
-    for (int i = 0; i < chip->cells->n; i++) {
-        pmCell *cell = chip->cells->data[i];
-        if (!cell) continue;
-        if (pmCellCheckDataStatusForPSFmodel (cell)) return true;
-    }
-    return false;
-}
-
-bool pmCellCheckDataStatusForPSFmodel (const pmCell *cell) {
-
-    for (int i = 0; i < cell->readouts->n; i++) {
-        pmReadout *readout = cell->readouts->data[i];
-        if (!readout) continue;
-        if (pmReadoutCheckDataStatusForPSFmodel (readout)) return true;
-    }
-    return false;
-}
-
-bool pmReadoutCheckDataStatusForPSFmodel (const pmReadout *readout) {
+bool pmPSFmodelCheckDataStatusForChip (const pmChip *chip) {
 
     bool status;
 
     // select the psf of interest
-    pmPSF *psf = psMetadataLookupPtr(&status, readout->analysis, "PSPHOT.PSF");
+    pmPSF *psf = psMetadataLookupPtr(&status, chip->analysis, "PSPHOT.PSF");
     return psf ? true : false;
 }
 
-bool pmFPAviewWritePSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmPSFmodelWriteForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
 
@@ -130 +98 @@
 
     if (view->chip == -1) {
-        if (!pmFPAWritePSFmodel(fpa, view, file, config)) {
+        if (!pmPSFmodelWriteFPA(fpa, view, file, config)) {
             psError(PS_ERR_IO, false, "Failed to write PSF for fpa");
             return false;
@@ -143 +111 @@
 
     if (view->cell == -1) {
-        if (!pmChipWritePSFmodel (chip, view, file, config)) {
+        if (!pmPSFmodelWriteChip (chip, view, file, config)) {
             psError(PS_ERR_IO, false, "Failed to write PSF for chip");
             return false;
@@ -150 +118 @@
     }
 
-    if (view->cell >= chip->cells->n) {
-        return false;
-    }
-    pmCell *cell = chip->cells->data[view->cell];
-
-    if (view->readout == -1) {
-        if (!pmCellWritePSFmodel (cell, view, file, config)) {
-            psError(PS_ERR_IO, false, "Failed to write PSF for cell");
-            return false;
-        }
-        return true;
-    }
-
-    if (view->readout >= cell->readouts->n) {
-        return false;
-    }
-    pmReadout *readout = cell->readouts->data[view->readout];
-
-    if (!pmReadoutWritePSFmodel (readout, view, file, config)) {
-        psError(PS_ERR_IO, false, "Failed to write PSF for readout");
-        return false;
-    }
-    return true;
+    psError(PS_ERR_IO, false, "PSF must be written at the chip level");
+    return false;
 }
 
 // read in all chip-level PSFmodel files for this FPA
-bool pmFPAWritePSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmPSFmodelWriteFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
 
     for (int i = 0; i < fpa->chips->n; i++) {
         pmChip *chip = fpa->chips->data[i];
-        if (!pmChipWritePSFmodel (chip, view, file, config)) {
+        if (!pmPSFmodelWriteChip (chip, view, file, config)) {
             psError(PS_ERR_IO, false, "Failed to write PSF for %dth chip", i);
             return false;
@@ -190 +137 @@
 
 // read in all cell-level PSFmodel files for this chip
-bool pmChipWritePSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-    for (int i = 0; i < chip->cells->n; i++) {
-        pmCell *cell = chip->cells->data[i];
-        if (!pmCellWritePSFmodel (cell, view, file, config)) {
-            psError(PS_ERR_IO, false, "Failed to write PSF for %dth cell", i);
-            return false;
-        }
+bool pmPSFmodelWriteChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+    if (!pmPSFmodelWrite (chip->analysis, view, file, config)) {
+        psError(PS_ERR_IO, false, "Failed to write PSF for chip");
+        return false;
     }
     return true;
 }
 
-// read in all readout-level PSFmodel files for this cell
-bool pmCellWritePSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-    for (int i = 0; i < cell->readouts->n; i++) {
-        pmReadout *readout = cell->readouts->data[i];
-        if (!pmReadoutWritePSFmodel (readout, view, file, config)) {
-            psError(PS_ERR_IO, false, "Failed to write PSF for %dth readout", i);
-            return false;
-        }
-    }
-    return true;
-}
-
-// for each Readout (ie, analysed image), we write out
+// for a pmPSF supplied on the analysis metadata, we write out
 // - image header        : FITS Image NAXIS = 0
 // - psf table (+header) : FITS Table
 // - psf resid (+header) : FITS Image
 // if needed, we also write out a PHU blank header
-bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
     bool status;
     pmHDU *hdu;
     char *headName, *tableName, *residName;
+
+    if (!analysis) return false;
 
     // write a PHU? (only if input image is MEF)
@@ -300 +233 @@
 
     // select the psf of interest
-    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+    pmPSF *psf = psMetadataLookupPtr (&status, analysis, "PSPHOT.PSF");
     if (!psf) {
-        psError(PS_ERR_UNKNOWN, true, "missing PSF for this readout");
+        psError(PS_ERR_UNKNOWN, true, "missing PSF for this analysis metadata");
         psFree (tableName);
         psFree (residName);
@@ -313 +246 @@
         psMetadata *header = psMetadataAlloc();
 
-        char *modelName = pmModelGetType (psf->type);
+        char *modelName = pmModelClassGetName (psf->type);
         psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName);
 
         psMetadataAddBool (header, PS_LIST_TAIL, "POISSON",  0, "Use Poisson errors in fits?", psf->poissonErrors);
 
-        int nPar = pmModelParameterCount (psf->type)    ;
+        int nPar = pmModelClassParameterCount (psf->type)    ;
         psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
 
@@ -324 +257 @@
         for (int i = 0; i < nPar; i++) {
             char name[9];
-            psPolynomial2D *poly = psf->params_NEW->data[i];
+            psPolynomial2D *poly = psf->params->data[i];
             if (poly == NULL) continue;
             snprintf (name, 9, "PAR%02d_NX", i);
@@ -344 +277 @@
         psArray *psfTable = psArrayAllocEmpty (100);
         for (int i = 0; i < nPar; i++) {
-            psPolynomial2D *poly = psf->params_NEW->data[i];
+            psPolynomial2D *poly = psf->params->data[i];
             if (poly == NULL) continue; // skip unset parameters (eg, XPOS)
             for (int ix = 0; ix <= poly->nX; ix++) {
@@ -445 +378 @@
 
 // if this file needs to have a PHU written out, write one
-bool pmPSF_WritePHU (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) {
+bool pmPSFmodelWritePHU (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) {
 
     // not needed if already written
@@ -487 +420 @@
 }
 
-bool pmFPAviewReadPSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmPSFmodelReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
 
@@ -493 +426 @@
 
     if (view->chip == -1) {
-        return pmFPAReadPSFmodel(fpa, view, file, config);
+        return pmPSFmodelReadFPA(fpa, view, file, config);
     }
 
@@ -502 +435 @@
 
     if (view->cell == -1) {
-        return pmChipReadPSFmodel(chip, view, file, config);
-    }
-
-    if (view->cell >= chip->cells->n) {
-        psAbort("Programming error: view does not apply to FPA.");
-    }
-    pmCell *cell = chip->cells->data[view->cell];
-
-    if (view->readout == -1) {
-        return pmCellReadPSFmodel(cell, view, file, config);
-    }
-
-    if (view->readout >= cell->readouts->n) {
-        psAbort("Programming error: view does not apply to FPA.");
-    }
-    pmReadout *readout = cell->readouts->data[view->readout];
-
-    return pmReadoutReadPSFmodel(readout, view, file, config);
+        return pmPSFmodelReadChip(chip, view, file, config);
+    }
+
+    psError(PS_ERR_IO, false, "PSF must be read at the chip level");
+    return false;
 }
 
 // read in all chip-level PSFmodel files for this FPA
-bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmPSFmodelReadFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
     bool success = true;                // Was everything successful?
     for (int i = 0; i < fpa->chips->n; i++) {
         pmChip *chip = fpa->chips->data[i];
-        success &= pmChipReadPSFmodel(chip, view, file, config);
+        success &= pmPSFmodelReadChip(chip, view, file, config);
     }
     return success;
@@ -534 +454 @@
 
 // read in all cell-level PSFmodel files for this chip
-bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-    bool success = true;                // Was everything successful?
-    for (int i = 0; i < chip->cells->n; i++) {
-        pmCell *cell = chip->cells->data[i];
-        success &= pmCellReadPSFmodel (cell, view, file, config);
-    }
-    return success;
-}
-
-// read in all readout-level PSFmodel files for this cell
-bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-    bool success = true;                // Was everything successful?
-    for (int i = 0; i < cell->readouts->n; i++) {
-        pmReadout *readout = cell->readouts->data[i];
-        success &= pmReadoutReadPSFmodel(readout, view, file, config);
-    }
-    return success;
+bool pmPSFmodelReadChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+    if (!pmPSFmodelRead (chip->analysis, view, file, config)) {
+        psError(PS_ERR_IO, false, "Failed to write PSF for chip");
+        return false;
+    }
+    return true;
 }
 
 // for each Readout (ie, analysed image), we write out: header + table with PSF model parameters,
 // and header + image for the PSF residual images
-bool pmReadoutReadPSFmodel(pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmPSFmodelRead (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
     bool status;
@@ -599 +507 @@
     // load the PSF model parameters from the FITS table
     char *modelName = psMetadataLookupStr (&status, header, "PSF_NAME");
-    pmModelType type = pmModelSetType (modelName);
+    pmModelType type = pmModelClassGetType (modelName);
     if (type == -1) {
         psError(PS_ERR_UNKNOWN, true, "invalid model name %s in psf file %s", modelName, file->filename);
@@ -612 +520 @@
     // check the number of expected parameters
     int nPar = psMetadataLookupS32 (&status, header, "PSF_NPAR");
-    if (nPar != pmModelParameterCount (psf->type))
+    if (nPar != pmModelClassParameterCount (psf->type))
         psAbort("mismatch model par count");
 
@@ -627 +535 @@
             return false;
         }
-        psf->params_NEW->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder);
+        psf->params->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder);
     }
 
@@ -650 +558 @@
         // XXX sanity check here
 
-        psPolynomial2D *poly = psf->params_NEW->data[iPar];
+        psPolynomial2D *poly = psf->params->data[iPar];
         if (poly == NULL) {
             psError(PS_ERR_UNKNOWN, true, "values for parameter %d, but missing NX", iPar);
@@ -698 +606 @@
     }
 
-    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
+    psMetadataAdd (analysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
     psFree (psf);
 
@@ -708 +616 @@
 }
 
-/************ old support functions, deprecate? **************/
-
+// create a psMetadata representation (human-readable) of a psf model
 psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf)
 {
@@ -717 +624 @@
     }
 
-    char *modelName = pmModelGetType (psf->type);
+    char *modelName = pmModelClassGetName (psf->type);
     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName);
 
-    int nPar = pmModelParameterCount (psf->type)    ;
+    int nPar = pmModelClassParameterCount (psf->type)    ;
     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
 
     for (int i = 0; i < nPar; i++) {
-        psPolynomial2D *poly = psf->params_NEW->data[i];
+        psPolynomial2D *poly = psf->params->data[i];
         if (poly == NULL)
             continue;
@@ -742 +649 @@
 }
 
+// parse a psMetadata representation (human-readable) of a psf model
 pmPSF *pmPSFfromMetadata (psMetadata *metadata)
 {
@@ -749 +657 @@
 
     char *modelName = psMetadataLookupPtr (&status, metadata, "PSF_MODEL_NAME");
-    pmModelType type = pmModelSetType (modelName);
+    pmModelType type = pmModelClassGetType (modelName);
 
     bool poissonErrors = psMetadataLookupPtr (&status, metadata, "PSF_POISSON_ERRORS");
@@ -759 +667 @@
 
     int nPar = psMetadataLookupS32 (&status, metadata, "PSF_MODEL_NPAR");
-    if (nPar != pmModelParameterCount (psf->type))
+    if (nPar != pmModelClassParameterCount (psf->type))
         psAbort("mismatch model par count");
 
@@ -770 +678 @@
             continue;
         psPolynomial2D *poly = psPolynomial2DfromMetadata (folder);
-        psFree (psf->params_NEW->data[i]);
-        psf->params_NEW->data[i] = poly;
+        psFree (psf->params->data[i]);
+        psf->params->data[i] = poly;
     }
     sprintf (keyword, "APTREND");
@@ -789 +697 @@
     return (psf);
 }
-
-// read in all readout-level Objects files for this cell
-bool pmReadoutWritePSFmodel_Config (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-    bool status;
-    char *filename;
-    char *realname;
-
-    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
-
-    switch (file->type) {
-      case PM_FPA_FILE_PSF:
-        filename = pmFPAfileNameFromRule (file->filerule, file, view);
-        bool create = file->mode == PM_FPA_MODE_WRITE ? true : false;
-        realname = pmConfigConvertFilename (filename, config, create);
-
-        psMetadata *psfData = pmPSFtoMetadata (NULL, psf);
-        psMetadataConfigWrite (psfData, realname);
-        psFree (psfData);
-        psFree (realname);
-        psFree (filename);
-        return true;
-
-      default:
-        fprintf (stderr, "warning: type mismatch\n");
-        break;
-    }
-    return false;
-}
-
-// read in all readout-level Objects files for this cell
-bool pmReadoutReadPSFmodel_Config (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    unsigned int Nfail;
-    char *filename;
-    char *realname;
-
-    switch (file->type) {
-      case PM_FPA_FILE_PSF:
-        filename = pmFPAfileNameFromRule (file->filerule, file, view);
-        bool create = file->mode == PM_FPA_MODE_WRITE ? true : false;
-        realname = pmConfigConvertFilename (filename, config, create);
-
-        psMetadata *psfData = psMetadataConfigRead(NULL, &Nfail, realname, FALSE);
-        pmPSF *psf = pmPSFfromMetadata (psfData);
-        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "psphot psf", psf);
-
-        psFree (psf);
-        psFree (psfData);
-        psFree (realname);
-        psFree (filename);
-
-        return true;
-
-      default:
-        fprintf (stderr, "warning: type mismatch\n");
-        break;
-    }
-    return false;
-}
-

trunk/psModules/src/objects/pmPSF_IO.h

@@ -6 +6 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-07-14 03:20:22 $
+ * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -17 +17 @@
 /// @{
 
+bool pmPSFmodelWriteForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmPSFmodelWriteFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmPSFmodelWriteChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+
+bool pmPSFmodelWritePHU (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+
+bool pmPSFmodelReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmPSFmodelReadFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmPSFmodelReadChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmPSFmodelRead (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+
+bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file);
+bool pmPSFmodelCheckDataStatusForFPA (const pmFPA *fpa);
+bool pmPSFmodelCheckDataStatusForChip (const pmChip *chip);
+
 psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf);
 pmPSF *pmPSFfromMetadata (psMetadata *metadata);
-bool pmFPAviewWritePSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmFPAWritePSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmChipWritePSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmCellWritePSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-
-bool pmPSF_WritePHU (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-
-bool pmFPAviewReadPSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmReadoutReadPSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-
-bool pmFPAviewCheckDataStatusForPSFmodel (const pmFPAview *view, const pmFPAfile *file);
-bool pmFPACheckDataStatusForPSFmodel (const pmFPA *fpa);
-bool pmChipCheckDataStatusForPSFmodel (const pmChip *chip);
-bool pmCellCheckDataStatusForPSFmodel (const pmCell *cell);
-bool pmReadoutCheckDataStatusForPSFmodel (const pmReadout *readout);
 
 /// @}

trunk/psModules/src/objects/pmPSFtry.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-20 02:22:26 $
+ *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -23 +23 @@
 #include "pmMoments.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
 #include "pmSource.h"
-#include "pmGrowthCurve.h"
-#include "pmPSF.h"
+#include "pmSourceUtils.h"
 #include "pmPSFtry.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
+#include "pmModelUtils.h"
 #include "pmSourceFitModel.h"
 #include "pmSourcePhotometry.h"
@@ -59 +61 @@
 
     // validate the requested model name
-    pmModelType type = pmModelSetType (modelName);
+    pmModelType type = pmModelClassGetType (modelName);
     if (type == -1) {
         psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName);
@@ -412 +414 @@
     // This way, the parameters masked by one of the fits will be applied to the others
     for (int i = 0; i < stats->clipIter; i++) {
-        psVectorClipFitPolynomial2D (psf->params_NEW->data[PM_PAR_E0], stats, psfTry->mask, 0xff, e0, dz, x, y);
+        psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E0], stats, psfTry->mask, 0xff, e0, dz, x, y);
         psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", stats->clippedNvalues, e0->n);
-        psVectorClipFitPolynomial2D (psf->params_NEW->data[PM_PAR_E1], stats, psfTry->mask, 0xff, e1, dz, x, y);
+        psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E1], stats, psfTry->mask, 0xff, e1, dz, x, y);
         psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", stats->clippedNvalues, e1->n);
-        psVectorClipFitPolynomial2D (psf->params_NEW->data[PM_PAR_E2], stats, psfTry->mask, 0xff, e2, dz, x, y);
+        psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E2], stats, psfTry->mask, 0xff, e2, dz, x, y);
         psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", stats->clippedNvalues, e2->n);
     }
@@ -436 +438 @@
 
     // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY)
-    for (int i = 0; i < psf->params_NEW->n; i++) {
+    for (int i = 0; i < psf->params->n; i++) {
         switch (i) {
           case PM_PAR_SKY:
@@ -461 +463 @@
         // the mask is carried from previous steps and updated with this operation
         // the weight is either the flux error or NULL, depending on 'applyWeights'
-        if (!psVectorClipFitPolynomial2D(psf->params_NEW->data[i], stats, psfTry->mask, 0xff, z, NULL, x, y)) {
+        if (!psVectorClipFitPolynomial2D(psf->params->data[i], stats, psfTry->mask, 0xff, z, NULL, x, y)) {
             psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);
             psFree(stats);
@@ -489 +491 @@
             fprintf (f, "%f %f : ", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
 
-            for (int i = 0; i < psf->params_NEW->n; i++) {
-                if (psf->params_NEW->data[i] == NULL)
+            for (int i = 0; i < psf->params->n; i++) {
+                if (psf->params->data[i] == NULL)
                     continue;
                 fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]);

trunk/psModules/src/objects/pmPeaks.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-26 01:20:29 $
+ *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -35 +35 @@
                         pmPeakType type)
 {
-    psTrace(__func__, 5, "---- begin ----\n");
+    psTrace("psModules.objects", 5, "---- begin ----\n");
 
     if (peaks == NULL) {
@@ -67 +67 @@
     psFree (peak);
 
-    psTrace(__func__, 5, "---- end ----\n");
+    psTrace("psModules.objects", 5, "---- end ----\n");
     return(peaks);
 }
@@ -118 +118 @@
 *****************************************************************************/
 static void peakFree(pmPeak *tmp)
-{} // used by pmIsPeak()
+{} // used by pmPeakTest()
 
 pmPeak *pmPeakAlloc(psS32 x,
@@ -145 +145 @@
 }
 
-bool pmIsPeak(const psPtr ptr)
+bool pmPeakTest(const psPtr ptr)
 {
     return (psMemGetDeallocator(ptr) == (psFreeFunc)peakFree);
@@ -193 +193 @@
 
 /******************************************************************************
-pmFindVectorPeaks(vector, threshold): Find all local peaks in the given vector
+pmPeaksInVector(vector, threshold): Find all local peaks in the given vector
 above the given threshold.  Returns a vector of type PS_TYPE_U32 containing
 the location (x value) of all peaks.
@@ -203 +203 @@
 Depending upon actual use, this may need to be optimized.
 *****************************************************************************/
-psVector *pmFindVectorPeaks(const psVector *vector,
-                            psF32 threshold)
+psVector *pmPeaksInVector(const psVector *vector,
+                         psF32 threshold)
 {
     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
@@ -295 +295 @@
 
 /******************************************************************************
-pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage
+pmPeaksInImage(image, threshold): Find all local peaks in the given psImage
 above the given threshold.  Returns a psArray containing location (x/y value)
 of all peaks.
@@ -309 +309 @@
 
 *****************************************************************************/
-psArray *pmFindImagePeaks(const psImage *image, psF32 threshold)
+psArray *pmPeaksInImage(const psImage *image, psF32 threshold)
 {
     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
@@ -327 +327 @@
     row = 0;
     tmpRow = getRowVectorFromImage((psImage *) image, row);
-    psVector *row1 = pmFindVectorPeaks(tmpRow, threshold);
-    // pmFindVectorPeaks returns coords in the vector, not corrected for col0
+    psVector *row1 = pmPeaksInVector(tmpRow, threshold);
+    // pmPeaksInVector returns coords in the vector, not corrected for col0
 
     for (psU32 i = 0 ; i < row1->n ; i++ ) {
@@ -382 +382 @@
     for (row = 1 ; row < (image->numRows - 1) ; row++) {
         tmpRow = getRowVectorFromImage((psImage *) image, row);
-        row1 = pmFindVectorPeaks(tmpRow, threshold);
+        row1 = pmPeaksInVector(tmpRow, threshold);
 
         // Step through all local peaks in this row.
@@ -461 +461 @@
     row = image->numRows - 1;
     tmpRow = getRowVectorFromImage((psImage *) image, row);
-    row1 = pmFindVectorPeaks(tmpRow, threshold);
+    row1 = pmPeaksInVector(tmpRow, threshold);
     for (psU32 i = 0 ; i < row1->n ; i++ ) {
         col = row1->data.U32[i];
@@ -501 +501 @@
 }
 
-
-/******************************************************************************
-psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have
-a peak value above the given maximum, or fall outside the valid region.
- 
-XXX: Should the sky value be used when comparing the maximum?
- 
-XXX: warning message if valid is NULL?
- 
-XXX: changed API to create a NEW output psArray (should change name as well)
- 
-XXX: Do we free the psList elements of those culled peaks?
- 
-XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset?
-*****************************************************************************/
-psList *pmCullPeaks(psList *peaks,
-                    psF32 maxValue,
-                    const psRegion valid)
-{
-    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(peaks, NULL);
-
-    psListElem *tmpListElem = (psListElem *) peaks->head;
-    psS32 indexNum = 0;
-
-    //    printf("pmCullPeaks(): list size is %d\n", peaks->size);
-    while (tmpListElem != NULL) {
-        pmPeak *tmpPeak = (pmPeak *) tmpListElem->data;
-        if ((tmpPeak->value > maxValue) ||
-                (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))) {
-            psListRemoveData(peaks, (psPtr) tmpPeak);
-        }
-
-        indexNum++;
-        tmpListElem = tmpListElem->next;
-    }
-
-    psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
-    return(peaks);
-}
-
-// XXX EAM: I changed this to return a new, subset array
-//          rather than alter the existing one
-// XXX: Fix the *valid pointer.
+// return a new array of peaks which are in the valid region and below threshold
+// XXX this function is unused and probably could be dropped
 psArray *pmPeaksSubset(
     psArray *peaks,
@@ -555 +513 @@
     psArray *output = psArrayAllocEmpty (200);
 
-    psTrace (".pmObjects.pmCullPeaks", 3, "list size is %ld\n", peaks->n);
+    psTrace ("psModules.objects", 3, "list size is %ld\n", peaks->n);
 
     for (int i = 0; i < peaks->n; i++) {

trunk/psModules/src/objects/pmPeaks.h

@@ -10 +10 @@
  * @author GLG, MHPCC
  *
- * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-01-24 02:54:15 $
+ * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -72 +72 @@
 );
 
-bool pmIsPeak(const psPtr ptr);
+bool pmPeakTest(const psPtr ptr);
 
-/** pmFindVectorPeaks()
+/** pmPeaksInVector()
  *
  * Find all local peaks in the given vector above the given threshold. A peak
@@ -91 +91 @@
  *
  */
-psVector *pmFindVectorPeaks(
+psVector *pmPeaksInVector(
     const psVector *vector,  ///< The input vector (float)
     float threshold   ///< Threshold above which to find a peak
@@ -97 +97 @@
 
 
-/** pmFindImagePeaks()
+/** pmPeaksInImage()
  *
  * Find all local peaks in the given image above the given threshold. This
@@ -111 +111 @@
  *
  */
-psArray *pmFindImagePeaks(
+psArray *pmPeaksInImage(
     const psImage *image,  ///< The input image where peaks will be found (float)
     float threshold   ///< Threshold above which to find a peak

trunk/psModules/src/objects/pmSource.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-08-16 18:33:00 $
+ *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27 +27 @@
 #include "pmMoments.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
 #include "pmSource.h"
@@ -113 +115 @@
 }
 
-bool pmIsSource(const psPtr ptr)
+bool pmSourceTest(const psPtr ptr)
 {
     return (psMemGetDeallocator(ptr) == (psFreeFunc)sourceFree);
@@ -337 +339 @@
             return emptyClump;
         }
-        peaks = pmFindImagePeaks (splane, stats[0].max / 2);
+        peaks = pmPeaksInImage (splane, stats[0].max / 2);
         psTrace ("psModules.objects", 2, "clump threshold is %f\n", stats[0].max/2);
 
@@ -797 +799 @@
 
 // should we call pmSourceCacheModel if it does not exist?
-bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psMaskType maskVal) {
+bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psMaskType maskVal, int dx, int dy) {
 
     bool status;
@@ -856 +858 @@
         target = source->weight;
     }
+
     if (add) {
-        status = pmModelAdd (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
+        status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
     } else {
-        status = pmModelSub (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
+        status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
     }
 
@@ -866 +869 @@
 
 bool pmSourceAdd (pmSource *source, pmModelOpMode mode, psMaskType maskVal) {
-    return pmSourceOp (source, mode, true, maskVal);
+    return pmSourceOp (source, mode, true, maskVal, 0, 0);
 }
 
 bool pmSourceSub (pmSource *source, pmModelOpMode mode, psMaskType maskVal) {
-    return pmSourceOp (source, mode, false, maskVal);
+    return pmSourceOp (source, mode, false, maskVal, 0, 0);
+}
+
+bool pmSourceAddWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy) {
+    return pmSourceOp (source, mode, true, maskVal, dx, dy);
+}
+
+bool pmSourceSubWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy) {
+    return pmSourceOp (source, mode, false, maskVal, dx, dy);
 }
 

trunk/psModules/src/objects/pmSource.h

@@ -3 +3 @@
  * @author EAM, IfA; GLG, MHPCC
  *
- * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-08-15 20:21:18 $
+ * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -57 +57 @@
  *
  */
-typedef struct
-{
+struct pmSource {
     const int id;                       ///< Unique ID for object
     pmPeak *peak;                       ///< Description of peak pixel.
@@ -81 +80 @@
     psRegion region;                    ///< area on image covered by selected pixels
     float sky, skyErr;                  ///< The sky and its error at the center of the object
-}
-pmSource;
+};
 
 /** pmPSFClump data structure
@@ -114 +112 @@
 void pmSourceFreePixels(pmSource *source);
 
-bool pmIsSource(const psPtr ptr);
+bool pmSourceTest(const psPtr ptr);
 
 /** pmSourceDefinePixels()
@@ -218 +216 @@
 bool pmSourceAdd (pmSource *source, pmModelOpMode mode, psMaskType maskVal);
 bool pmSourceSub (pmSource *source, pmModelOpMode mode, psMaskType maskVal);
-
-bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psMaskType maskVal);
+bool pmSourceAddWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy);
+bool pmSourceSubWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy);
+
+bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psMaskType maskVal, int dx, int dy);
 bool pmSourceCacheModel (pmSource *source, psMaskType maskVal);
 bool pmSourceCachePSF (pmSource *source, psMaskType maskVal);

trunk/psModules/src/objects/pmSourceContour.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -26 +26 @@
 #include "pmMoments.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
 #include "pmSource.h"

trunk/psModules/src/objects/pmSourceFitModel.c

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-21 22:58:11 $
+ *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27 +27 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceFitModel.h"
 
@@ -116 +116 @@
     psVector *dparams = model->dparams;
 
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
-    if (!modelFunc)
-        psAbort("invalid model function");
-    pmModelLimits checkLimits = pmModelLimits_GetFunction (model->type);
-    if (!checkLimits)
-        psAbort("invalid model limits function");
-
     // create the minimization constraints
     psMinConstraint *constraint = psMinConstraintAlloc();
     constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
-    constraint->checkLimits = checkLimits;
+    constraint->checkLimits = model->modelLimits;
 
     // set parameter mask based on fitting mode
@@ -156 +149 @@
     // force the floating parameters to fall within the contraint ranges
     for (int i = 0; i < params->n; i++) {
-        checkLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
-        checkLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
+        model->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
+        model->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
     }
 
@@ -174 +167 @@
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
 
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, modelFunc);
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, model->modelFunc);
     for (int i = 0; i < dparams->n; i++) {
         if (psTraceGetLevel("psModules.objects") >= 4) {
@@ -201 +194 @@
             altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1;
         }
-        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, modelFunc);
+        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, model->modelFunc);
         for (int i = 0; i < dparams->n; i++) {
             if (!constraint->paramMask->data.U8[i])
@@ -237 +230 @@
 }
 
-/********************* Source Model Set Functions ***************************/
-
-// these static variables are used by one pass of pmSourceFitSet to store
-// data for a model set.  If we are going to make psphot thread-safe, these
-// will have to go in a structure of their own or be allocated once per thread
-// sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2,
-// nPar = nSrc*(nOnePar - 1) + 1
-static pmModelFunc oneModelFunc;
-static pmModelLimits oneCheckLimits;
-static psVector *onePar;
-static psVector *oneDeriv;
-static int nOnePar;
-
-bool pmSourceFitSetInit (pmModelType type)
-{
-
-    oneModelFunc = pmModelFunc_GetFunction (type);
-    oneCheckLimits = pmModelLimits_GetFunction (type);
-    nOnePar = pmModelParameterCount (type);
-
-    onePar = psVectorAlloc (nOnePar, PS_DATA_F32);
-    oneDeriv = psVectorAlloc (nOnePar, PS_DATA_F32);
-
-    return true;
-}
-
-void pmSourceFitSetClear (void)
-{
-
-    psFree (onePar);
-    psFree (oneDeriv);
-    return;
-}
-
-bool pmSourceFitSet_CheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas)
-{
-    // convert the value of nParam into corresponding single model parameter entry
-    // convert params into corresponding single model parameter pointer
-
-    int nParamSingle = (nParam - 1) % (nOnePar - 1) + 1;
-    float *paramSingle = params + nParam - nParamSingle;
-    float *betaSingle = betas + nParam - nParamSingle;
-    bool status = oneCheckLimits (mode, nParamSingle, paramSingle, betaSingle);
-    return status;
-}
-
-psF32 pmSourceFitSet_Function(psVector *deriv,
-                              const psVector *params,
-                              const psVector *x)
-{
-
-    psF32 value;
-    psF32 model;
-
-    psF32 *PAR = onePar->data.F32;
-    psF32 *dPAR = oneDeriv->data.F32;
-
-    psF32 *pars = params->data.F32;
-    psF32 *dpars = (deriv == NULL) ? NULL : deriv->data.F32;
-
-    int nSrc = (params->n - 1) / (nOnePar - 1);
-
-    PAR[0] = model = pars[0];
-    for (int i = 0; i < nSrc; i++) {
-        int nOff = i*nOnePar - i;
-        for (int n = 1; n < nOnePar; n++) {
-            PAR[n] = pars[nOff + n];
-        }
-        if (deriv == NULL) {
-            value = oneModelFunc (NULL, onePar, x);
-        } else {
-            value = oneModelFunc (oneDeriv, onePar, x);
-            for (int n = 1; n < nOnePar; n++) {
-                dpars[nOff + n] = dPAR[n];
-            }
-        }
-        model += value;
-    }
-    if (deriv != NULL) {
-        dpars[0] = dPAR[0]*2.0;
-    }
-    return (model);
-}
-
-/*
-i:       0                   1                 2
-n:         1  2  3  4  5  6  1  2  3  4  5  6  1  2  3  4  5  6
-i*6 + n: 0 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18
-*/
-
-bool pmSourceFitSet (pmSource *source,
-                     psArray *modelSet,
-                     pmSourceFitMode mode,
-                     psMaskType maskVal)
-{
-    psTrace("psModules.objects", 3, "---- %s begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
-    PS_ASSERT_PTR_NON_NULL(source->weight, false);
-
-    psBool fitStatus = true;
-    psBool onPic     = true;
-    psBool rc        = true;
-
-    // maximum number of valid pixels
-    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
-
-    // construct the coordinate and value entries
-    psArray *x = psArrayAllocEmpty(nPix);
-    psVector *y = psVectorAllocEmpty(nPix, PS_TYPE_F32);
-    psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32);
-
-    // fill in the coordinate and value entries
-    nPix = 0;
-    for (psS32 i = 0; i < source->pixels->numRows; i++) {
-        for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            // skip masked points
-            if (source->maskObj->data.U8[i][j] & maskVal) {
-                continue;
-            }
-            // skip zero-weight points
-            if (source->weight->data.F32[i][j] == 0) {
-                continue;
-            }
-            // skip nan values in image
-            if (!isfinite(source->pixels->data.F32[i][j])) {
-                continue;
-            }
-
-            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
-
-            // Convert i/j to image space:
-            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
-            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
-            x->data[nPix] = (psPtr *) coord;
-            y->data.F32[nPix] = source->pixels->data.F32[i][j];
-
-            // psMinimizeLMChi2 takes wt = 1/dY^2.  suggestion from RHL is to use the local sky
-            // as weight to avoid the bias from systematic errors here we would just use the
-            // source sky variance
-            if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) {
-                yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
-            } else {
-                yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
-            }
-            nPix++;
-        }
-    }
-    x->n = nPix;
-    y->n = nPix;
-    yErr->n = nPix;
-
-    // base values on first model (all models must be identical)
-    pmModel *model = modelSet->data[0];
-
-    // determine number of model parameters
-    int nSrc = modelSet->n;
-    int nPar = model->params->n - 1;  // number of object parameters (excluding sky)
-
-    // define parameter vectors for source set
-    psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
-    psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
-
-    // set the static variables
-    pmSourceFitSetInit (model->type);
-
-    // create the minimization constraints
-    psMinConstraint *constraint = psMinConstraintAlloc();
-    constraint->paramMask = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8);
-    constraint->checkLimits = pmSourceFitSet_CheckLimits;
-
-    // set the parameter guesses for the multiple models
-    // first the value for the single sky parameter
-    params->data.F32[0] = model->params->data.F32[0];
-
-    // next, the values for the source parameters
-    for (int i = 0; i < nSrc; i++) {
-        model = modelSet->data[i];
-        for (int n = 1; n < nPar + 1; n++) {
-            params->data.F32[i*nPar + n] = model->params->data.F32[n];
-        }
-    }
-
-    if (psTraceGetLevel("psModules.objects") >= 5) {
-        for (int i = 0; i < params->n; i++) {
-            fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.U8[i]);
-        }
-    }
-
-    // set the parameter masks based on the fitting mode
-    int nParams = 0;
-    switch (mode) {
-    case PM_SOURCE_FIT_NORM:
-        // NORM-only model fits only source normalization (Io)
-        nParams = nSrc;
-        psVectorInit (constraint->paramMask, 1);
-        for (int i = 0; i < nSrc; i++) {
-            constraint->paramMask->data.U8[1 + i*nPar] = 0;
-        }
-        break;
-    case PM_SOURCE_FIT_PSF:
-        // PSF model only fits x,y,Io
-        nParams = 3*nSrc;
-        psVectorInit (constraint->paramMask, 1);
-        for (int i = 0; i < nSrc; i++) {
-            constraint->paramMask->data.U8[1 + i*nPar] = 0;
-            constraint->paramMask->data.U8[2 + i*nPar] = 0;
-            constraint->paramMask->data.U8[3 + i*nPar] = 0;
-        }
-        break;
-    case PM_SOURCE_FIT_EXT:
-        // EXT model fits all params (except sky)
-        nParams = nPar*nSrc;
-        psVectorInit (constraint->paramMask, 0);
-        constraint->paramMask->data.U8[0] = 1;
-        break;
-    default:
-        psAbort("invalid fitting mode");
-    }
-
-    // force the floating parameters to fall within the contraint ranges
-    for (int i = 0; i < params->n; i++) {
-        pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
-        pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
-    }
-
-    if (nPix <  nParams + 1) {
-        psTrace (__func__, 4, "insufficient valid pixels\n");
-        psTrace("psModules.objects", 3, "---- %s() end : fail pixels ----\n", __func__);
-        model->flags |= PM_MODEL_STATUS_BADARGS;
-        psFree (x);
-        psFree (y);
-        psFree (yErr);
-        psFree (params);
-        psFree (dparams);
-        psFree(constraint);
-        pmSourceFitSetClear ();
-        return(false);
-    }
-
-    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
-
-    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
-
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, pmSourceFitSet_Function);
-    if (!fitStatus) {
-        psTrace("psModules.objects", 4, "Failed to fit model (%d)\n", nSrc);
-    }
-
-    // parameter errors from the covariance matrix
-    for (int i = 0; i < dparams->n; i++) {
-        if ((constraint->paramMask != NULL) && constraint->paramMask->data.U8[i])
-            continue;
-        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
-    }
-
-    // get the Gauss-Newton distance for fixed model parameters
-    if (constraint->paramMask != NULL) {
-        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
-        psVector *altmask = psVectorAlloc (params->n, PS_TYPE_U8);
-        altmask->data.U8[0] = 1;
-        for (int i = 1; i < dparams->n; i++) {
-            altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1;
-        }
-        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSet_Function);
-        for (int i = 0; i < dparams->n; i++) {
-            if (!constraint->paramMask->data.U8[i])
-                continue;
-            // note that delta is the value *subtracted* from the parameter
-            // to get the new guess.  for dparams to represent the direction
-            // of motion, we need to take -delta
-            dparams->data.F32[i] = -delta->data.F32[i];
-        }
-        psFree (delta);
-        psFree (altmask);
-    }
-
-    // assign back the parameters to the models
-    for (int i = 0; i < nSrc; i++) {
-        model = modelSet->data[i];
-        model->params->data.F32[0] = params->data.F32[0];
-        for (int n = 1; n < nPar + 1; n++) {
-            if (psTraceGetLevel("psModules.objects") >= 4) {
-                fprintf (stderr, "%f ", params->data.F32[i*nPar + n]);
-            }
-            model->params->data.F32[n] = params->data.F32[i*nPar + n];
-            model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n];
-        }
-        psTrace ("psModules.objects", 4, " src %d", i);
-
-        // save the resulting chisq, nDOF, nIter
-        // these are not unique for any one source
-        model->chisq = myMin->value;
-        model->nIter = myMin->iter;
-        model->nDOF  = y->n - nParams;
-
-        // set the model success or failure status
-        model->flags |= PM_MODEL_STATUS_FITTED;
-        if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;
-
-        // models can go insane: reject these
-        onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->col0);
-        onPic &= (model->params->data.F32[PM_PAR_XPOS] <  source->pixels->col0 + source->pixels->numCols);
-        onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->row0);
-        onPic &= (model->params->data.F32[PM_PAR_XPOS] <  source->pixels->row0 + source->pixels->numRows);
-        if (!onPic) model->flags |= PM_MODEL_STATUS_OFFIMAGE;
-    }
-    psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
-
-    source->mode |= PM_SOURCE_MODE_FITTED;
-
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
-    psFree(myMin);
-    psFree(covar);
-    psFree(constraint);
-    psFree(params);
-    psFree(dparams);
-
-    // free static memory used by pmSourceFitSet
-    pmSourceFitSetClear ();
-
-    rc = (onPic && fitStatus);
-    psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF);
-    psTrace("psModules.objects", 5, "---- %s end (%d) ----\n", __func__, rc);
-    return(rc);
-}
-

trunk/psModules/src/objects/pmSourceIO.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-14 03:20:44 $
+ *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -34 +34 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 

trunk/psModules/src/objects/pmSourceIO_CMP.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-03 00:13:03 $
+ *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -32 +32 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 
@@ -154 +154 @@
 
     // define PSF model type
-    int modelType = pmModelSetType ("PS_MODEL_GAUSS");
+    int modelType = pmModelClassGetType ("PS_MODEL_GAUSS");
 
     char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME");
     if (PSF_NAME != NULL) {
-        modelType = pmModelSetType (PSF_NAME);
+        modelType = pmModelClassGetType (PSF_NAME);
     }
 

trunk/psModules/src/objects/pmSourceIO_OBJ.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-03 00:13:03 $
+ *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -32 +32 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 

trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-03 00:13:42 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -32 +32 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 
@@ -152 +152 @@
 
     // define PSF model type
-    int modelType = pmModelSetType ("PS_MODEL_GAUSS");
+    int modelType = pmModelClassGetType ("PS_MODEL_GAUSS");
 
     char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME");
     if (PSF_NAME != NULL) {
-        modelType = pmModelSetType (PSF_NAME);
+        modelType = pmModelClassGetType (PSF_NAME);
     }
 

trunk/psModules/src/objects/pmSourceIO_RAW.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -32 +32 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourcePhotometry.h"
 #include "pmSourceIO.h"

trunk/psModules/src/objects/pmSourceIO_SMPDATA.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-03 00:13:42 $
+ *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -30 +30 @@
 #include "pmPeaks.h"
 #include "pmMoments.h"
+#include "pmResiduals.h"
 #include "pmGrowthCurve.h"
-#include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 
@@ -132 +132 @@
 
     // define PSF model type
-    int modelType = pmModelSetType ("PS_MODEL_GAUSS");
+    int modelType = pmModelClassGetType ("PS_MODEL_GAUSS");
 
     char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME");
     if (PSF_NAME != NULL) {
-        modelType = pmModelSetType (PSF_NAME);
+        modelType = pmModelClassGetType (PSF_NAME);
     }
 

trunk/psModules/src/objects/pmSourceIO_SX.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-27 22:14:08 $
+ *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -32 +32 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 

trunk/psModules/src/objects/pmSourcePhotometry.c

@@ -3 +3 @@
  *  @author EAM, IfA; GLG, MHPCC
  *
- *  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-20 02:22:26 $
+ *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -26 +26 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourcePhotometry.h"
 
@@ -250 +250 @@
 
     // measure fitMag
-    pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
-    fitSum = modelFluxFunc (model->params);
+    fitSum = model->modelFlux (model->params);
     if (fitSum <= 0)
         return false;
@@ -324 +323 @@
 
     // the model function returns the source flux at a position
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
     psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
 
@@ -354 +352 @@
 
             // for the full model, add all points
-            value = modelFunc (NULL, params, coord) - sky;
+            value = model->modelFunc (NULL, params, coord) - sky;
             modelSum += value;
 

trunk/psModules/src/objects/pmSourcePlotApResid.c

@@ -4 +4 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-22 21:45:48 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *  Copyright 2006 IfA, University of Hawaii
  */
@@ -27 +27 @@
 #include "pmPeaks.h"
 #include "pmMoments.h"
+#include "pmResiduals.h"
 #include "pmGrowthCurve.h"
-#include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
 #include "pmSourcePlots.h"

trunk/psModules/src/objects/pmSourcePlotMoments.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-22 21:45:48 $
+ *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2006 IfA, University of Hawaii
@@ -32 +32 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
 #include "pmSourcePlots.h"

trunk/psModules/src/objects/pmSourcePlotPSFModel.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-22 21:45:48 $
+ *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2006 IfA, University of Hawaii
@@ -32 +32 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
 #include "pmSourcePlots.h"

trunk/psModules/src/objects/pmSourceSky.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-20 02:22:26 $
+ *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27 +27 @@
 #include "pmMoments.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
 #include "pmSource.h"