Changeset 5844

Timestamp:

Dec 23, 2005, 3:24:38 PM (20 years ago)

Author:

desonia

Message:

changes from eam_rel9_b1

Location:

trunk/psModules

Files:

• src/objects/Makefile.am (modified) (1 diff)
• src/objects/pmModelGroup.c (modified) (10 diffs)
• src/objects/pmModelGroup.h (modified) (2 diffs)
• src/objects/pmObjects.c (modified) (26 diffs)
• src/objects/pmObjects.h (modified) (56 diffs)
• src/objects/pmPSF.c (modified) (1 diff)
• src/objects/pmPSFtry.c (modified) (4 diffs)
• src/objects/pmPSFtry.h (modified) (8 diffs)
• test/imsubtract/tst_pmSubtractBias.c (modified) (5 diffs)
• test/objects/verified/tst_pmObjects01.stdout (modified) (1 diff)

trunk/psModules/src/objects/Makefile.am

@@ -15 +15 @@
         models/pmModel_QGAUSS.c \
         models/pmModel_SGAUSS.c 
-    
+
 psmoduleincludedir = $(includedir)
 psmoduleinclude_HEADERS = \

trunk/psModules/src/objects/pmModelGroup.c

@@ -10 +10 @@
 #include "models/pmModel_SGAUSS.c"
 
-static pmModelGroup models[] = {
-                                   {"PS_MODEL_GAUSS",        7, pmModelFunc_GAUSS,   pmModelFlux_GAUSS,   pmModelRadius_GAUSS,   pmModelLimits_GAUSS,   pmModelGuess_GAUSS,  pmModelFromPSF_GAUSS,  pmModelFitStatus_GAUSS},
-                                   {"PS_MODEL_PGAUSS",       7, pmModelFunc_PGAUSS,  pmModelFlux_PGAUSS,  pmModelRadius_PGAUSS,  pmModelLimits_PGAUSS,  pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelFitStatus_PGAUSS},
-                                   {"PS_MODEL_QGAUSS",       8, pmModelFunc_QGAUSS,  pmModelFlux_QGAUSS,  pmModelRadius_QGAUSS,  pmModelLimits_QGAUSS,  pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelFitStatus_QGAUSS},
-                                   {"PS_MODEL_SGAUSS",       9, pmModelFunc_SGAUSS,  pmModelFlux_SGAUSS,  pmModelRadius_SGAUSS,  pmModelLimits_SGAUSS,  pmModelGuess_SGAUSS, pmModelFromPSF_SGAUSS, pmModelFitStatus_SGAUSS},
-                               };
+static pmModelGroup defaultModels[] = {
+                                          {"PS_MODEL_GAUSS",        7, pmModelFunc_GAUSS,   pmModelFlux_GAUSS,   pmModelRadius_GAUSS,   pmModelLimits_GAUSS,   pmModelGuess_GAUSS,  pmModelFromPSF_GAUSS,  pmModelFitStatus_GAUSS},
+                                          {"PS_MODEL_PGAUSS",       7, pmModelFunc_PGAUSS,  pmModelFlux_PGAUSS,  pmModelRadius_PGAUSS,  pmModelLimits_PGAUSS,  pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelFitStatus_PGAUSS},
+                                          {"PS_MODEL_QGAUSS",       8, pmModelFunc_QGAUSS,  pmModelFlux_QGAUSS,  pmModelRadius_QGAUSS,  pmModelLimits_QGAUSS,  pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelFitStatus_QGAUSS},
+                                          {"PS_MODEL_SGAUSS",       9, pmModelFunc_SGAUSS,  pmModelFlux_SGAUSS,  pmModelRadius_SGAUSS,  pmModelLimits_SGAUSS,  pmModelGuess_SGAUSS, pmModelFromPSF_SGAUSS, pmModelFitStatus_SGAUSS},
+                                      };
+
+static pmModelGroup *models = NULL;
+static int Nmodels = 0;
+
+static void ModelGroupFree (pmModelGroup *modelGroup)
+{
+
+    if (modelGroup == NULL)
+        return;
+    psFree (modelGroup);
+    return;
+}
+
+pmModelGroup *pmModelGroupAlloc (int nModels)
+{
+
+    pmModelGroup *modelGroup = (pmModelGroup *) psAlloc (nModels * sizeof(pmModelGroup));
+    psMemSetDeallocator(modelGroup, (psFreeFunc) ModelGroupFree);
+    return (modelGroup);
+}
+
+void pmModelGroupAdd (pmModelGroup *model)
+{
+
+    if (models == NULL) {
+        pmModelGroupInit ();
+    }
+
+    Nmodels ++;
+    models = (pmModelGroup *) psRealloc (models, Nmodels*sizeof(pmModelGroup));
+    models[Nmodels-1] = model[0];
+    return;
+}
+
+void pmModelGroupInit (void)
+{
+
+    int Nnew = sizeof (defaultModels) / sizeof (pmModelGroup);
+
+    models = pmModelGroupAlloc (Nnew);
+    for (int i = 0; i < Nnew; i++) {
+        models[i] = defaultModels[i];
+    }
+    Nmodels = Nnew;
+    return;
+}
 
 pmModelFunc pmModelFunc_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -29 +74 @@
 pmModelFlux pmModelFlux_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -39 +83 @@
 pmModelRadius pmModelRadius_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -49 +92 @@
 pmModelLimits pmModelLimits_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -59 +101 @@
 pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -69 +110 @@
 pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -79 +119 @@
 pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -89 +128 @@
 psS32 pmModelParameterCount (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -99 +137 @@
 psS32 pmModelSetType (char *name)
 {
-
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     for (int i = 0; i < Nmodels; i++) {
         if (!strcmp(models[i].name, name)) {
@@ -111 +147 @@
 char *pmModelGetType (pmModelType type)
 {
-
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");

trunk/psModules/src/objects/pmModelGroup.h

@@ -9 +9 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-10-10 19:53:40 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -209 +209 @@
 pmModelGroup;
 
+// allocate a pmModelGroup to hold nModels entries
+pmModelGroup *pmModelGroupAlloc (int nModels);
+
+// initialize the internal (static) model group with the default models
+void pmModelGroupInit (void);
+
+// add a new model to the internal (static) model group
+void pmModelGroupAdd (pmModelGroup *model);
+
 # endif

trunk/psModules/src/objects/pmObjects.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -85 +85 @@
 psVector containing the specified row of data from the psImage.
  
-XXX: Is there a better way to do this?  
+XXX: Is there a better way to do this?
 XXX EAM: does this really need to alloc a new vector???
 *****************************************************************************/
@@ -271 +271 @@
     tmp->chisq = 0.0;
     tmp->nIter = 0;
+    tmp->radius = 0;
+    tmp->status = PM_MODEL_UNTRIED;
+
     psS32 Nparams = pmModelParameterCount(type);
     if (Nparams == 0) {
@@ -850 +853 @@
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
         for (psS32 col = 0; col < source->pixels->numCols ; col++) {
-            if ((source->mask != NULL) && (source->mask->data.U8[row][col]))
+            if ((source->mask != NULL) && (source->mask->data.U8[row][col])) {
                 continue;
+            }
 
             psF32 xDiff = col + source->pixels->col0 - xPeak;
@@ -858 +862 @@
             // XXX EAM : calculate xDiff, yDiff up front;
             //           radius is just a function of (xDiff, yDiff)
-            if (!VALID_RADIUS(xDiff, yDiff, R2))
+            if (!VALID_RADIUS(xDiff, yDiff, R2)) {
                 continue;
+            }
 
             psF32 pDiff = source->pixels->data.F32[row][col] - sky;
@@ -865 +870 @@
             // XXX EAM : check for valid S/N in pixel
             // XXX EAM : should this limit be user-defined?
-            if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1)
+            if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) {
                 continue;
+            }
 
             Sum += pDiff;
@@ -970 +976 @@
 
 /******************************************************************************
-pmSourcePSFClump(source, metadata): Find the likely PSF clump in the 
-sigma-x, sigma-y plane. return 0,0 clump in case of error. 
+pmSourcePSFClump(source, metadata): Find the likely PSF clump in the
+sigma-x, sigma-y plane. return 0,0 clump in case of error.
 *****************************************************************************/
 
@@ -994 +1000 @@
         psImage *splane = NULL;
         int binX, binY;
+        bool status;
+
+        psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");
+        if (!status)
+            SX_MAX = 10.0;
+        psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");
+        if (!status)
+            SY_MAX = 10.0;
 
         // construct a sigma-plane image
         // psImageAlloc does zero the data
-        splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);
+        splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);
         for (int i = 0; i < splane->numRows; i++)
         {
@@ -1131 +1145 @@
 XXX: How can this function ever return FALSE?
  
-XXX EAM : add the saturated mask value to metadata 
+XXX EAM : add the saturated mask value to metadata
 *****************************************************************************/
 
@@ -1146 +1160 @@
     int Ncr      = 0;
     int Nsatstar = 0;
+    psRegion allArray = psRegionSet (0, 0, 0, 0);
 
     psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
@@ -1177 +1192 @@
 
         // XXX EAM : can we use the value of SATURATE if mask is NULL?
-        //
-        // XXX: Must verify this region (the region argument was added to psImageCountPixelMask()
-        // after EAM wrote this code.
-        //
-        psRegion tmpRegion = psRegionSet(0, tmpSrc->mask->numCols, 0, tmpSrc->mask->numRows);
-        int Nsatpix = psImageCountPixelMask(tmpSrc->mask, tmpRegion, PSPHOT_MASK_SATURATED);
+        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, PSPHOT_MASK_SATURATED);
 
         // saturated star (size consistent with PSF or larger)
@@ -1254 +1264 @@
 
 /** pmSourceDefinePixels()
- * 
+ *
  * Define psImage subarrays for the source located at coordinates x,y on the
  * image set defined by readout. The pixels defined by this operation consist of
@@ -1266 +1276 @@
  * example). This function should be used to define a region of interest around a
  * source, including both source and sky pixels.
- * 
+ *
  * XXX: must code this.
- * 
+ *
  */
 bool pmSourceDefinePixels(
@@ -1360 +1370 @@
 /******************************************************************************
 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.  
+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
@@ -1598 +1608 @@
         psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
         psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        model->status = PM_MODEL_BADARGS;
         return(false);
     }
@@ -1649 +1660 @@
     }
 
-    // XXX EAM: we need to do something (give an error?) if rc is false
-    // XXX EAM: psMinimizeLMChi2 does not check convergence
-
-    // XXX models can go insane: reject these
-    onPic &= (params->data.F32[2] >= source->pixels->col0);
-    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
-    onPic &= (params->data.F32[3] >= source->pixels->row0);
-    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
-
-    // XXX EAM: save the resulting chisq, nDOF, nIter
+    // save the resulting chisq, nDOF, nIter
     model->chisq = myMin->value;
     model->nIter = myMin->iter;
     model->nDOF  = y->n - nParams;
 
-    // XXX EAM get the Gauss-Newton distance for fixed model parameters
+    // get the Gauss-Newton distance for fixed model parameters
     if (paramMask != NULL) {
         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
@@ -1673 +1675 @@
         }
         psFree (delta);
+    }
+
+    // set the model success or failure status
+    if (!fitStatus) {
+        model->status = PM_MODEL_NONCONVERGE;
+    } else {
+        model->status = PM_MODEL_SUCCESS;
+    }
+
+    // models can go insane: reject these
+    onPic &= (params->data.F32[2] >= source->pixels->col0);
+    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
+    onPic &= (params->data.F32[3] >= source->pixels->row0);
+    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
+    if (!onPic) {
+        model->status = PM_MODEL_OFFIMAGE;
     }
 
@@ -1732 +1750 @@
         psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
         psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        model->status = PM_MODEL_BADARGS;
         return(false);
     }
@@ -1810 +1829 @@
     // XXX EAM: psMinimizeLMChi2 does not check convergence
 
-    // XXX models can go insane: reject these
-    onPic &= (params->data.F32[2] >= source->pixels->col0);
-    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
-    onPic &= (params->data.F32[3] >= source->pixels->row0);
-    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
-
     // XXX EAM: save the resulting chisq, nDOF, nIter
     model->chisq = myMin->value;
@@ -1832 +1845 @@
     }
 
-    psFree(paramMask);
+    // set the model success or failure status
+    if (!fitStatus) {
+        model->status = PM_MODEL_NONCONVERGE;
+    } else {
+        model->status = PM_MODEL_SUCCESS;
+    }
+
+    // XXX models can go insane: reject these
+    onPic &= (params->data.F32[2] >= source->pixels->col0);
+    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
+    onPic &= (params->data.F32[3] >= source->pixels->row0);
+    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
+    if (!onPic) {
+        model->status = PM_MODEL_OFFIMAGE;
+    }
+
     psFree(x);
     psFree(y);
+    psFree(yErr);
     psFree(myMin);
+    psFree(covar);
+    psFree(paramMask);
 
     rc = (onPic && fitStatus);
@@ -1846 +1877 @@
                              pmModel *model,
                              bool center,
-                             psS32 flag)
+                             bool sky,
+                             bool add
+                                )
 {
     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
@@ -1859 +1892 @@
     psS32 imageCol;
     psS32 imageRow;
+    psF32 skyValue = params->data.F32[0];
+    psF32 pixelValue;
 
     for (psS32 i = 0; i < image->numRows; i++) {
@@ -1864 +1899 @@
             if ((mask != NULL) && mask->data.U8[i][j])
                 continue;
-            psF32 pixelValue;
+
             // XXX: Should you be adding the pixels for the entire subImage,
             // or a radius of pixels around it?
@@ -1882 +1917 @@
             x->data.F32[0] = (float) imageCol;
             x->data.F32[1] = (float) imageRow;
-            pixelValue = modelFunc (NULL, params, x);
-            // fprintf (stderr, "%f %f  %d %d  %f\n", x->data.F32[0], x->data.F32[1], i, j, pixelValue);
-
-            if (flag == 1) {
-                pixelValue = -pixelValue;
-            }
-
-            // XXX: Must figure out how to calculate the image coordinates and
-            // how to use the boolean "center" flag.
-
-            image->data.F32[i][j]+= pixelValue;
+
+            // set the appropriate pixel value for this coordinate
+            if (sky) {
+                pixelValue = modelFunc (NULL, params, x);
+            } else {
+                pixelValue = modelFunc (NULL, params, x) - skyValue;
+            }
+
+
+            // add or subtract the value
+            if (add
+               ) {
+                image->data.F32[i][j] += pixelValue;
+            }
+            else {
+                image->data.F32[i][j] -= pixelValue;
+            }
         }
     }
@@ -1907 +1948 @@
                       psImage *mask,
                       pmModel *model,
-                      bool center)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 0);
+                      bool center,
+                      bool sky)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true);
     psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);
@@ -1920 +1962 @@
                       psImage *mask,
                       pmModel *model,
-                      bool center)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 1);
+                      bool center,
+                      bool sky)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false);
     psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);

trunk/psModules/src/objects/pmObjects.h

@@ -4 +4 @@
  * images is one of the critical tasks of the IPP or any astronomical software
  * system. This file will define structures and functions related to the task
- * of source detection and measurement. The elements defined in this section 
+ * of source detection and measurement. The elements defined in this section
  * are generally low-level components which can be connected together to
  * construct a complete object measurement suite.
@@ -10 +10 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -45 +45 @@
 
 /** pmPeakType
- * 
+ *
  *  A peak pixel may have several features which may be determined when the
  *  peak is found or measured. These are specified by the pmPeakType enum.
@@ -52 +52 @@
  *  edge. The PM_PEAK_FLAT represents a peak pixel which has more than a specific
  *  number of neighbors at the same value, within some tolarence:
- * 
+ *
  */
 typedef enum {
@@ -63 +63 @@
 
 /** pmPeak data structure
- *  
+ *
  *  A source has the capacity for several types of measurements. The
  *  simplest measurement of a source is the location and flux of the peak pixel
  *  associated with the source:
- *  
+ *
  */
 typedef struct
@@ -80 +80 @@
 
 /** pmMoments data structure
- *  
+ *
  * One of the simplest measurements which can be made quickly for an object
  * are the object moments. We specify a structure to carry the moment information
  * for a specific source:
- *  
+ *
  */
 typedef struct
@@ -103 +103 @@
 
 /** pmPSFClump data structure
- * 
+ *
  * A collection of object moment measurements can be used to determine
  * approximate object classes. The key to this analysis is the location and
  * statistics (in the second-moment plane,
- *  
+ *
  */
 typedef struct
@@ -118 +118 @@
 pmPSFClump;
 
+// type of model carried by the pmModel structure
 typedef int pmModelType;
-#define PS_MODEL_GAUSS 0
-#define PS_MODEL_PGAUSS 1
-#define PS_MODEL_QGAUSS 2
-#define PS_MODEL_SGAUSS 3
-
+
+typedef enum {
+    PM_MODEL_UNTRIED,               ///< model fit not yet attempted
+    PM_MODEL_SUCCESS,               ///< model fit succeeded
+    PM_MODEL_NONCONVERGE,           ///< model fit did not converge
+    PM_MODEL_OFFIMAGE,              ///< model fit drove out of range
+    PM_MODEL_BADARGS                ///< model fit called with invalid args
+} pmModelStatus;
 
 /** pmModel data structure
- * 
+ *
  * Every source may have two types of models: a PSF model and a FLT (floating)
  * model. The PSF model represents the best fit of the image PSF to the specific
@@ -132 +136 @@
  * object by the PSF, not by the fit. The FLT model represents the best fit of
  * the given model to the object, with all parameters floating in the fit.
- *  
+ *
  */
 typedef struct
@@ -142 +146 @@
     int nDOF;    ///< number of degrees of freedom
     int nIter;    ///< number of iterations to reach min
+    int status;         ///< fit status
     float radius;   ///< fit radius actually used
 }
@@ -147 +152 @@
 
 /** pmSourceType enumeration
- * 
+ *
  * A given source may be identified as most-likely to be one of several source
  * types. The pmSource entry pmSourceType defines the current best-guess for this
  * source.
- * 
+ *
  * XXX: The values given below are currently illustrative and will require
  * some modification as the source classification code is developed. (TBD)
- * 
+ *
  */
 typedef enum {
@@ -178 +183 @@
 
 /** pmSource data structure
- *  
+ *
  *  This source has the capacity for several types of measurements. The
  *  simplest measurement of a source is the location and flux of the peak pixel
  *  associated with the source:
- *  
+ *
  */
 typedef struct
@@ -194 +199 @@
     pmModel *modelFLT;   ///< FLT (floating) Model fit (parameters and type).
     pmSourceType type;   ///< Best identification of object.
+    float apMag;
+    float fitMag;
 }
 pmSource;
@@ -211 +218 @@
 
 /** pmMomentsAlloc()
- * 
+ *
  */
 pmMoments *pmMomentsAlloc();
@@ -217 +224 @@
 
 /** pmModelAlloc()
- * 
+ *
  */
 pmModel *pmModelAlloc(pmModelType type);
@@ -223 +230 @@
 
 /** pmSourceAlloc()
- * 
+ *
  */
 pmSource  *pmSourceAlloc();
@@ -229 +236 @@
 
 /** pmFindVectorPeaks()
- * 
+ *
  * Find all local peaks in the given vector above the given threshold. A peak
  * is defined as any element with a value greater than its two neighbors and with
@@ -243 +250 @@
  * a vector containing the coordinates (element number) of the detected peaks
  * (type psU32).
- * 
+ *
  */
 psVector *pmFindVectorPeaks(
@@ -252 +259 @@
 
 /** pmFindImagePeaks()
- * 
+ *
  * Find all local peaks in the given image above the given threshold. This
  * function should find all row peaks using pmFindVectorPeaks, then test each row
@@ -263 +270 @@
  * peaks at the (+x,+y) corner. When selecting the peaks, their type must also be
  * set. The result of this function is an array of pmPeak entries.
- * 
+ *
  */
 psArray *pmFindImagePeaks(
@@ -272 +279 @@
 
 /** pmCullPeaks()
- * 
+ *
  * Eliminate peaks from the psList that have a peak value above the given
  * maximum, or fall outside the valid region.
- * 
+ *
  */
 psList *pmCullPeaks(
@@ -285 +292 @@
 
 /** pmPeaksSubset()
- * 
+ *
  * Create a new peaks array, removing certain types of peaks from the input
  * array of peaks based on the given criteria. Peaks should be eliminated if they
@@ -291 +298 @@
  * the valid region.  The result of the function is a new array with a reduced
  * number of peaks.
- * 
+ *
  */
 psArray *pmPeaksSubset(
@@ -301 +308 @@
 
 /** pmSourceDefinePixels()
- * 
+ *
  * Define psImage subarrays for the source located at coordinates x,y on the
  * image set defined by readout. The pixels defined by this operation consist of
@@ -313 +320 @@
  * example). This function should be used to define a region of interest around a
  * source, including both source and sky pixels.
- * 
+ *
  * XXX: must code this.
- * 
+ *
  */
 // XXX: Uncommenting the pmReadout causes compile errors.
@@ -328 +335 @@
 
 /** pmSourceLocalSky()
- * 
+ *
  * Measure the local sky in the vicinity of the given source. The Radius
  * defines the square aperture in which the moments will be measured. This
@@ -338 +345 @@
  * This function allocates the pmMoments structure. The resulting sky is used to
  * set the value of the pmMoments.sky element of the provided pmSource structure.
- * 
+ *
  */
 bool pmSourceLocalSky(
@@ -348 +355 @@
 
 /** pmSourceMoments()
- * 
+ *
  * Measure source moments for the given source, using the value of
  * source.moments.sky provided as the local background value and the peak
@@ -355 +362 @@
  * are measured within the given circular radius of the source.peak coordinates.
  * The return value indicates the success (TRUE) of the operation.
- * 
+ *
  */
 bool pmSourceMoments(
@@ -364 +371 @@
 
 /** pmSourcePSFClump()
- * 
+ *
  * We use the source moments to make an initial, approximate source
  * classification, and as part of the information needed to build a PSF model for
@@ -373 +380 @@
  * similar. The function returns a pmPSFClump structure, representing the
  * centroid and size of the clump in the sigma_x, sigma_y second-moment plane.
- * 
+ *
  * The goal is to identify and characterize the stellar clump within the
  * sigma_x, sigma_y second-moment plane.  To do this, an image is constructed to
@@ -384 +391 @@
  *  * used to calculate the median and standard deviation of the sigma_x, sigma_y
  * values. These resulting values are returned via the pmPSFClump structure.
- * 
+ *
  * The return value indicates the success (TRUE) of the operation.
- * 
+ *
  * XXX: Limit the S/N of the candidate sources (part of Metadata)? (TBD).
  * XXX: Save the clump parameters on the Metadata (TBD)
- * 
+ *
  */
 pmPSFClump pmSourcePSFClump(
@@ -398 +405 @@
 
 /** pmSourceRoughClass()
- * 
+ *
  * Based on the specified data values, make a guess at the source
  * classification. The sources are provides as a psArray of pmSource entries.
@@ -405 +412 @@
  * can be extracted from the metadata using the given keywords. Except as noted,
  * the data type for these parameters are psF32.
- * 
+ *
  */
 bool pmSourceRoughClass(
@@ -415 +422 @@
 
 /** pmSourceModelGuess()
- * 
+ *
  * Convert available data to an initial guess for the given model. This
  * function allocates a pmModel entry for the pmSource structure based on the
@@ -421 +428 @@
  * are specified for each model below. The guess values are placed in the model
  * parameters. The function returns TRUE on success or FALSE on failure.
- * 
+ *
  */
 pmModel *pmSourceModelGuess(
@@ -430 +437 @@
 
 /** pmContourType
- * 
+ *
  * Only one type is defined at present.
- * 
+ *
  */
 typedef enum {
@@ -442 +449 @@
 
 /** pmSourceContour()
- * 
+ *
  * Find points in a contour for the given source at the given level. If type
  * is PM_CONTOUR_CRUDE, the contour is found by starting at the source peak,
@@ -452 +459 @@
  * This function may be used as part of the model guess inputs.  Other contour
  * types may be specified in the future for more refined contours (TBD)
- * 
+ *
  */
 psArray *pmSourceContour(
@@ -463 +470 @@
 
 /** pmSourceFitModel()
- * 
+ *
  * Fit the requested model to the specified source. The starting guess for the
  * model is given by the input source.model parameter values. The pixels of
@@ -469 +476 @@
  * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
  * on success or FALSE on failure.
- * 
+ *
  */
 bool pmSourceFitModel(
@@ -479 +486 @@
 
 /** pmModelFitStatus()
- * 
+ *
  * This function wraps the call to the model-specific function returned by
  * pmModelFitStatusFunc_GetFunction.  The model-specific function examines the
  * model parameters, parameter errors, Chisq, S/N, and other parameters available
  * from model to decide if the particular fit was successful or not.
- * 
+ *
  * XXX: Must code this.
- * 
+ *
  */
 bool pmModelFitStatus(
@@ -494 +501 @@
 
 /** pmSourceAddModel()
- * 
+ *
  * Add the given source model flux to/from the provided image. The boolean
  * option center selects if the source is re-centered to the image center or if
@@ -501 +508 @@
  * is at most the pixel range specified by the source.pixels image. The success
  * status is returned.
- * 
+ *
  */
 bool pmSourceAddModel(
@@ -507 +514 @@
     psImage *mask,   ///< The image pixel mask (valid == 0)
     pmModel *model,   ///< The input pmModel
-    bool center    ///< A boolean flag that determines whether pixels are centered
+    bool center,    ///< A boolean flag that determines whether pixels are centered
+    bool sky        ///< A boolean flag that determines if the sky is subtracted
 );
 
 
 /** pmSourceSubModel()
- * 
+ *
  * Subtract the given source model flux to/from the provided image. The
  * boolean option center selects if the source is re-centered to the image center
@@ -519 +527 @@
  * image is at most the pixel range specified by the source.pixels image. The
  * success status is returned.
- * 
+ *
  */
 bool pmSourceSubModel(
@@ -525 +533 @@
     psImage *mask,   ///< The image pixel mask (valid == 0)
     pmModel *model,   ///< The input pmModel
-    bool center    ///< A boolean flag that determines whether pixels are centered
+    bool center,    ///< A boolean flag that determines whether pixels are centered
+    bool sky        ///< A boolean flag that determines if the sky is subtracted
 );
 
 
 /**
- * 
+ *
  * The function returns both the magnitude of the fit, defined as -2.5log(flux),
  * where the flux is integrated under the model, theoretically from a radius of 0
@@ -537 +546 @@
  * not excluded by the aperture mask. The model flux is calculated by calling the
  * model-specific function provided by pmModelFlux_GetFunction.
- * 
+ *
  * XXX: must code this.
- * 
+ *
  */
 bool pmSourcePhotometry(
@@ -551 +560 @@
 
 /**
- * 
+ *
  * This function converts the source classification into the closest available
  * approximation to the Dophot classification scheme:
- * 
+ *
  * PM_SOURCE_DEFECT: 8
  * PM_SOURCE_SATURATED: 8
@@ -569 +578 @@
  * PM_SOURCE_POOR_FIT_GAL: 2
  * PM_SOURCE_OTHER: ?
- * 
+ *
  */
 int pmSourceDophotType(
@@ -577 +586 @@
 
 /** pmSourceSextractType()
- * 
+ *
  * This function converts the source classification into the closest available
  * approximation to the Sextractor classification scheme. the correspondence is
  * not yet defined (TBD) .
- * 
+ *
  * XXX: Must code this.
- * 
+ *
  */
 int pmSourceSextractType(
@@ -590 +599 @@
 
 /** pmSourceFitModel_v5()
- * 
+ *
  * Fit the requested model to the specified source. The starting guess for the
  * model is given by the input source.model parameter values. The pixels of
@@ -596 +605 @@
  * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
  * on success or FALSE on failure.
- * 
+ *
  */
 bool pmSourceFitModel_v5(
@@ -606 +615 @@
 
 /** pmSourceFitModel_v7()
- * 
+ *
  * Fit the requested model to the specified source. The starting guess for the
  * model is given by the input source.model parameter values. The pixels of
@@ -612 +621 @@
  * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
  * on success or FALSE on failure.
- * 
+ *
  */
 bool pmSourceFitModel_v7(
@@ -622 +631 @@
 
 /** pmSourcePhotometry()
- * 
+ *
  * XXX: Need descriptions
- * 
+ *
  */
 bool pmSourcePhotometry(

trunk/psModules/src/objects/pmPSF.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/objects/pmPSFtry.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -176 +176 @@
 
     // XXX this function wants aperture radius for pmSourcePhotometry
-    pmPSFtryMetric (psfTry, RADIUS);
+    pmPSFtryMetric_Alt (psfTry, RADIUS);
     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
               modelName,
@@ -288 +288 @@
 
     // linear fit to rfBin, daBin
-    psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1);
+    psPolynomial1D *poly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD);
+    psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
 
     // XXX EAM : this is the intended API (cycle 7? cycle 8?)
-    poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
+    // poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
 
     // XXX EAM : replace this when the above version is implemented
@@ -314 +316 @@
     psFree (poly);
     psFree (stats);
+    psFree (fitstat);
 
     return true;
 }
+
+bool pmPSFtryMetric_Alt (pmPSFtry *try
+                         , float RADIUS)
+{
+
+    // the measured (aperture - fit) magnitudes (dA == try->metric)
+    //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
+    //     dA = dAo + dsky/flux
+    //   where flux is the flux of the star
+    // we fit this trend to find the infinite flux aperture correction (dAo),
+    //   the nominal sky bias (dsky), and the error on dAo
+    // the values of dA are contaminated by stars with close neighbors in the aperture
+    //   we use an outlier rejection to avoid this bias
+
+    // rflux = ten(0.4*fitMag);
+    psVector *rflux = psVectorAlloc (try
+                                     ->sources->n, PS_TYPE_F64);
+    for (int i = 0; i < try
+                ->sources->n; i++) {
+            if (try
+                    ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
+            rflux->data.F64[i] = pow(10.0, 0.4*try
+                                     ->fitMag->data.F64[i]);
+        }
+
+    // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    // XXX EAM
+    stats->min = 1.0;
+    stats->max = 3.0;
+    stats->clipIter = 3;
+
+    // linear clipped fit to rfBin, daBin
+    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+    poly = psVectorClipFitPolynomial1D (poly, stats, try
+                                            ->mask, PSFTRY_MASK_ALL, try
+                                                ->metric, NULL, rflux);
+    fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+
+    try
+        ->psf->ApResid = poly->coeff[0];
+    try
+        ->psf->dApResid = stats->sampleStdev;
+    try
+        ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+
+    fprintf (stderr, "*******************************************************************************\n");
+
+    FILE *f;
+    f = fopen ("apresid.dat", "w");
+    if (f == NULL)
+        psAbort ("pmPSFtry", "can't open output file");
+
+    for (int i = 0; i < try
+                ->sources->n; i++) {
+            fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try
+                         ->fitMag->data.F64[i], rflux->data.F64[i], try
+                             ->metric->data.F64[i], try
+                                 ->mask->data.U8[i]);
+        }
+    fclose (f);
+
+    psFree (rflux);
+    psFree (poly);
+    psFree (stats);
+
+    // psFree (daFit);
+    // psFree (daResid);
+
+    return true;
+}

trunk/psModules/src/objects/pmPSFtry.h

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 20:32:44 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -18 +18 @@
 
 /**
- * 
+ *
  * This structure contains a pointer to the collection of sources which will
  * be used to test the PSF model form. It lists the pmModelType type of model
@@ -38 +38 @@
  * ultimate metric to intercompare multiple types of PSF models is the value of
  * the aperture correction scatter.
- * 
+ *
  * XXX: There are many more members in the SDRS then in the prototype code.
  * I stuck with the prototype code.
- * 
- * 
+ *
+ *
  */
 typedef struct
@@ -58 +58 @@
 
 /** pmPSFtryMaskValues
- * 
+ *
  * The following datatype defines the masks used by the pmPSFtry analysis to
  * identify sources which should or should not be included in the analysis.
- * 
+ *
  */
 enum {
@@ -74 +74 @@
 
 /** pmPSFtryAlloc()
- * 
+ *
  * Allocate a pmPSFtry data structure.
- * 
+ *
  */
 pmPSFtry *pmPSFtryAlloc(
@@ -85 +85 @@
 
 /** pmPSFtryModel()
- * 
+ *
  * This function takes the input collection of sources and performs a complete
  * analysis to determine a PSF model of the given type (specified by model name).
  * The result is a pmPSFtry with the results of the analysis.
- * 
+ *
  */
 pmPSFtry *pmPSFtryModel(
@@ -99 +99 @@
 
 /** pmPSFtryMetric()
- * 
+ *
  * This function is used to measure the PSF model metric for the set of
  * results contained in the pmPSFtry structure.
- * 
+ *
  */
 bool pmPSFtryMetric(
@@ -109 +109 @@
 );
 
+/** pmPSFtryMetric_Alt()
+ *
+ * This function is used to measure the PSF model metric for the set of
+ * results contained in the pmPSFtry structure (alternative implementation).
+ *
+ */
+bool pmPSFtryMetric_Alt(
+    pmPSFtry *try
+    ,                      ///< Add comment.
+    float RADIUS                        ///< Add comment.
+);
+
 # endif

trunk/psModules/test/imsubtract/tst_pmSubtractBias.c

@@ -19 +19 @@
  *  test03: Calculate a row overscan vector and subtract it from each
  *  row in the input image.
- * test05: 
+ * test05:
  *
  *  @author GLG, MHPCC
@@ -25 +25 @@
  *  XXX: Memory leaks are not being detected.
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-11-23 23:54:30 $
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:36 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -326 +326 @@
         rc = pmSubtractBias(myReadout, NULL, NULL, PM_OVERSCAN_NONE, stat,
                             0, PM_FIT_NONE, myBias);
-     
+ 
         for (i=0;i<numRows;i++) {
             for (j=0;j<numCols;j++) {
@@ -335 +335 @@
                     testStatus = 1;
                 }
-     
+ 
                 // Restore myReadout for next test.
                 myReadout->image->data.F32[i][j] = (float) (i + j);
@@ -627 +627 @@
     // Set overscan axis in the metadata.
     //
-    psBool rc;
+    psBool rc = false;
     if (overscanaxis == PM_OVERSCAN_ROWS) {
         rc = psMetadataAddS32(myReadout->parent->concepts, PS_LIST_HEAD, "CELL.READDIR", 0, NULL, 1);

trunk/psModules/test/objects/verified/tst_pmObjects01.stdout

@@ -1 +1 @@
 Testing pmPeakAlloc()...
 Testing pmMomentsAlloc()...
-Testing pmModelAlloc(PS_MODEL_GAUSS)...
-Testing pmModelAlloc(PS_MODEL_GAUSS)...
-Testing pmModelAlloc(PS_MODEL_GAUSS)...
-Testing pmModelAlloc(PS_MODEL_GAUSS)...
 ----------------------------------------------------------------------------------
 Calling pmFindVectorPeaks with NULL psVector.  Should generate error and return NULL.