Changeset 5255 for trunk/psModules/src/objects

Timestamp:

Oct 10, 2005, 9:53:54 AM (21 years ago)

Author:

gusciora

Message:

A fairly large check-in. This incorporates must of Eugene's mods to
the object detection routines.

Location:

trunk/psModules/src/objects

Files:

• Makefile.am (modified) (1 diff)
• pmModelGroup.c (added)
• pmModelGroup.h (added)
• pmObjects.c (modified) (78 diffs)
• pmObjects.h (modified) (6 diffs)
• pmPSF.c (added)
• pmPSF.h (added)
• pmPSFtry.c (added)
• pmPSFtry.h (added)
• psEllipse.c (added)
• psEllipse.h (added)
• psLibUtils.c (added)
• psLibUtils.h (added)
• psModulesUtils.c (added)
• psModulesUtils.h (added)

trunk/psModules/src/objects/Makefile.am

@@ -3 +3 @@
 libpsmoduleobjects_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS)
 libpsmoduleobjects_la_LDFLAGS  = -release $(PACKAGE_VERSION)
-libpsmoduleobjects_la_SOURCES  = pmObjects.c
+libpsmoduleobjects_la_SOURCES  = \
+    pmObjects.c \
+    pmPSF.c \
+    pmPSFtry.c \
+    pmModelGroup.c \
+    psModulesUtils.c \
+    psLibUtils.c \
+    psEllipse.c
 
 psmoduleincludedir = $(includedir)
 psmoduleinclude_HEADERS = \
-  pmObjects.h
+    pmObjects.h \
+    pmPSF.h \
+    pmPSFtry.h \
+    pmModelGroup.h \
+    psModulesUtils.h \
+    psLibUtils.h \
+    psEllipse.h

trunk/psModules/src/objects/pmObjects.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-28 20:43:52 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-10-10 19:53:40 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -13 +13 @@
  */
 
-#include<stdio.h>
-#include<math.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
 #include "pslib.h"
-#include "psConstants.h"
 #include "pmObjects.h"
-
+#include "pmModelGroup.h"
 /******************************************************************************
 pmPeakAlloc(): Allocate the pmPeak data structure and set appropriate members.
@@ -27 +27 @@
                     pmPeakType class)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     pmPeak *tmp = (pmPeak *) psAlloc(sizeof(pmPeak));
     tmp->x = x;
@@ -33 +34 @@
     tmp->class = class;
 
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmp);
 }
@@ -42 +44 @@
 pmMoments *pmMomentsAlloc()
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     pmMoments *tmp = (pmMoments *) psAlloc(sizeof(pmMoments));
     tmp->x = 0.0;
     tmp->y = 0.0;
     tmp->Sx = 0.0;
-    tmp->Sx = 0.0;
+    tmp->Sy = 0.0;
     tmp->Sxy = 0.0;
     tmp->Sum = 0.0;
@@ -53 +56 @@
     tmp->nPixels = 0;
 
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmp);
 }
@@ -58 +62 @@
 static void modelFree(pmModel *tmp)
 {
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     psFree(tmp->params);
     psFree(tmp->dparams);
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
 }
 
@@ -65 +71 @@
 pmModelAlloc(): Allocate the pmModel structure, along with its parameters,
 and initialize the type member.  Initialize the params to 0.0.
-XXX EAM: changing params and dparams to psVector
+XXX EAM: simplifying code with pmModelParameterCount
 *****************************************************************************/
 pmModel *pmModelAlloc(pmModelType type)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel));
 
     tmp->type = type;
     tmp->chisq = 0.0;
-    switch (type) {
-    case PS_MODEL_GAUSS:
-        tmp->params  = psVectorAlloc(7, PS_TYPE_F32);
-        tmp->dparams = psVectorAlloc(7, PS_TYPE_F32);
-        break;
-    case PS_MODEL_PGAUSS:
-        tmp->params  = psVectorAlloc(7, PS_TYPE_F32);
-        tmp->dparams = psVectorAlloc(7, PS_TYPE_F32);
-        break;
-    case PS_MODEL_TWIST_GAUSS:
-        tmp->params  = psVectorAlloc(11, PS_TYPE_F32);
-        tmp->dparams = psVectorAlloc(11, PS_TYPE_F32);
-        break;
-    case PS_MODEL_WAUSS:
-        tmp->params  = psVectorAlloc(9, PS_TYPE_F32);
-        tmp->dparams = psVectorAlloc(9, PS_TYPE_F32);
-        break;
-    case PS_MODEL_SERSIC:
-        tmp->params  = psVectorAlloc(8, PS_TYPE_F32);
-        tmp->dparams = psVectorAlloc(8, PS_TYPE_F32);
-        break;
-    case PS_MODEL_SERSIC_CORE:
-        tmp->params  = psVectorAlloc(12, PS_TYPE_F32);
-        tmp->dparams = psVectorAlloc(12, PS_TYPE_F32);
-        break;
-    default:
+    tmp->nIter = 0;
+    psS32 Nparams = pmModelParameterCount(type);
+    if (Nparams == 0) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
         return(NULL);
     }
+
+    tmp->params  = psVectorAlloc(Nparams, PS_TYPE_F32);
+    tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32);
 
     for (psS32 i = 0; i < tmp->params->n; i++) {
@@ -109 +96 @@
 
     psMemSetDeallocator(tmp, (psFreeFunc) modelFree);
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmp);
 }
 
 /******************************************************************************
-XXX: We don't free pixels and mask since that caused a memory error.
-We might need to increase the reference counter and decrease it here.
+XXX EAM : we can now free these pixels - memory ref is incremented now
 *****************************************************************************/
 static void sourceFree(pmSource *tmp)
 {
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     psFree(tmp->peak);
-    //    psFree(tmp->pixels);
-    //    psFree(tmp->mask);
+    psFree(tmp->pixels);
+    psFree(tmp->weight);
+    psFree(tmp->mask);
     psFree(tmp->moments);
     psFree(tmp->modelPSF);
     psFree(tmp->modelFLT);
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
 }
 
@@ -132 +122 @@
 pmSource *pmSourceAlloc()
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     pmSource *tmp = (pmSource *) psAlloc(sizeof(pmSource));
     tmp->peak = NULL;
     tmp->pixels = NULL;
+    tmp->weight = NULL;
     tmp->mask = NULL;
     tmp->moments = NULL;
@@ -142 +134 @@
     psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);
 
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmp);
 }
@@ -159 +152 @@
                             psF32 threshold)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_VECTOR_NON_NULL(vector, NULL);
     PS_ASSERT_VECTOR_NON_EMPTY(vector, NULL);
@@ -177 +171 @@
             tmpVector = psVectorAlloc(0, PS_TYPE_U32);
         }
+        psTrace(__func__, 3, "---- %s() end ----\n", __func__);
         return(tmpVector);
     }
@@ -241 +236 @@
     }
 
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmpVector);
 }
@@ -248 +244 @@
 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???
 *****************************************************************************/
 static psVector *getRowVectorFromImage(psImage *image,
                                        psU32 row)
 {
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
     PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
@@ -260 +258 @@
         tmpVector->data.F32[col] = image->data.F32[row][col];
     }
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
     return(tmpVector);
 }
@@ -276 +275 @@
                               pmPeakType type)
 {
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type);
 
@@ -283 +283 @@
     }
     psArrayAdd(list, 100, tmpPeak);
-
+    psFree (tmpPeak);
+    // XXX EAM : is this free appropriate?  (does psArrayAdd increment memory counter?)
+
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
     return(list);
 }
@@ -308 +311 @@
                           psF32 threshold)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
     PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
     if ((image->numRows == 1) || (image->numCols == 1)) {
         psError(PS_ERR_UNKNOWN, true, "Currently, input image must have at least 2 rows and 2 columns.");
+        psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
         return(NULL);
     }
@@ -364 +369 @@
         }
     }
+    psFree (tmpRow);
+    psFree (row1);
 
     //
@@ -369 +376 @@
     //
     if (image->numRows == 1) {
+        psTrace(__func__, 3, "---- %s() end ----\n", __func__);
         return(list);
     }
@@ -442 +450 @@
                 }
             } else {
-                psError(PS_ERR_UNKNOWN, true, "peak specified valid column range.");
+                psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
             }
 
         }
+        psFree (tmpRow);
+        psFree (row1);
     }
 
@@ -484 +494 @@
             }
         } else {
-            psError(PS_ERR_UNKNOWN, true, "peak specified valid colum range.");
-        }
-    }
+            psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
+        }
+    }
+    psFree (tmpRow);
+    psFree (row1);
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(list);
 }
@@ -495 +508 @@
                              psS32 y)
 {
-
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     if ((x >= valid.x0) &&
             (x <= valid.x1) &&
             (y >= valid.y0) &&
             (y <= valid.y1)) {
+        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
         return(true);
     }
-
+    psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
     return(false);
 }
@@ -516 +530 @@
  
 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,
@@ -521 +539 @@
                     const psRegion valid)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(peaks, NULL);
-    //    PS_ASSERT_PTR_NON_NULL(valid, NULL);
 
     psListElem *tmpListElem = (psListElem *) peaks->head;
@@ -539 +557 @@
     }
 
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(peaks);
 }
@@ -544 +563 @@
 // XXX EAM: I changed this to return a new, subset array
 //          rather than alter the existing one
-psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion valid)
-{
+// XXX: Fix the *valid pointer.
+psArray *pmPeaksSubset(
+    psArray *peaks,
+    psF32 maxValue,
+    const psRegion valid)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(peaks, NULL);
 
@@ -561 +585 @@
         psArrayAdd (output, 200, tmpPeak);
     }
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(output);
 }
@@ -601 +626 @@
      members.
 *****************************************************************************/
-pmSource *pmSourceLocalSky(const psImage *image,
-                           const pmPeak *peak,
-                           psStatsOptions statsOptions,
-                           psF32 innerRadius,
-                           psF32 outerRadius)
-{
-    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
-    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
-    PS_ASSERT_PTR_NON_NULL(peak, NULL);
-    PS_FLOAT_COMPARE(0.0, innerRadius, NULL);
-    PS_FLOAT_COMPARE(innerRadius, outerRadius, NULL);
-    psS32 innerRadiusS32 = (psS32) innerRadius;
-    psS32 outerRadiusS32 = (psS32) outerRadius;
-
-    //
-    // We define variables for code readability.
-    //
-    // XXX: Since the peak->xy coords are in image, not subImage coords,
-    // these variables should be renamed for clarity (imageCenterRow, etc).
-    //
-    // peak->x,y is guaranteed to be on image
-    psS32 SubImageCenterRow = peak->y;
-    psS32 SubImageCenterCol = peak->x;
-
-    // XXX EAM : I added this code to stay on the image. So did George
-    psS32 SubImageStartRow  = PS_MAX(0, SubImageCenterRow - outerRadiusS32);
-    psS32 SubImageEndRow    = PS_MIN(image->numRows - 1, SubImageCenterRow + outerRadiusS32);
-    psS32 SubImageStartCol  = PS_MAX(0, SubImageCenterCol - outerRadiusS32);
-    psS32 SubImageEndCol    = PS_MIN(image->numCols - 1, SubImageCenterCol + outerRadiusS32);
-    // AnulusWidth == number of pixels width in the annulus.  We add one since
-    // the pixels at the inner AND outher radius are included.
-    psS32 AnulusWidth = 1 + (outerRadiusS32 - innerRadiusS32);
-    // Example: assume an outer/inner radius of 20/10.  Then the subimage
-    // should have width/length of 40.  An 18-by-18 interior region will
-    // be masked.
-    //    printf("pmSourceLocalSky(): innerRadiusS32 is %d\n", innerRadiusS32);
-    //    printf("pmSourceLocalSky(): outerRadiusS32 is %d\n", outerRadiusS32);
-    //    printf("pmSourceLocalSky(): AnulusWidth is %d\n", AnulusWidth);
-
-    // XXX EAM : these tests should not be needed: we can never hit this error because of above
-    # if (1)
-
-        if (SubImageStartRow < 0) {
-            psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
-                    SubImageStartRow);
-            return(NULL);
-        }
-    if (SubImageEndRow >= image->numRows) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
-                SubImageEndRow);
-        return(NULL);
-    }
-    if (SubImageStartCol < 0) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n",
-                SubImageStartCol);
-        return(NULL);
-    }
-    if (SubImageEndCol >= image->numCols) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
-                SubImageEndCol);
-        return(NULL);
-    }
-    # endif
-
-    //
-    // Grab a subimage of the original image of size (2 * outerRadius).
-    //
-    // XXX: Must fix for new psImageSubset
-    //    psImage *subImage = psImageSubset((psImage *) image,
-    //                                      SubImageStartCol,
-    //                                      SubImageStartRow,
-    //                                      SubImageEndCol,
-    //                                      SubImageEndRow);
-    //    printf("pmSourceLocalSky: subimage width/length is (%d, %d)\n", subImage->numCols, subImage->numRows);
-    psRegion tmpRegion = psRegionSet(SubImageStartCol,
-                                     SubImageEndCol,
-                                     SubImageStartRow,
-                                     SubImageEndRow);
-    psImage *subImage = psImageSubset((psImage *) image, tmpRegion);
-
-    psImage *subImageMask = psImageAlloc(subImage->numCols,
-                                         subImage->numRows,
-                                         PS_TYPE_U8);
-
-    //
-    // Loop through the subimage mask, initialize mask to 0.
-    //
-    for (psS32 row = 0 ; row < subImageMask->numRows; row++) {
-        for (psS32 col = 0 ; col < subImageMask->numCols; col++) {
-            subImageMask->data.U8[row][col] = 0;
-        }
-    }
-
-    //
-    // Loop through the subimage, mask off pixels in the inner square.
-    // XXX this uses a static mask value of 1
-    //
-    for (psS32 row = AnulusWidth; row <= (subImageMask->numRows - AnulusWidth) - 1; row++) {
-        for (psS32 col = AnulusWidth; col <= (subImageMask->numCols - AnulusWidth) - 1; col++) {
-            subImageMask->data.U8[row][col] = 1;
-        }
-    }
-
-
-    //    for (psS32 row = 0 ; row < subImage->numRows; row++) {
-    //        for (psS32 col = 0 ; col < subImage->numCols; col++) {
-    //            printf("(%d) ", subImageMask->data.U8[row][col]);
-    //        }
-    //        printf("\n");
-    //    }
-
-    //
-    // Allocate the myStats structure, then call psImageStats(), which will
-    // calculate the specified statistic.
-    //
+
+
+
+
+
+
+
+
+
+bool pmSourceLocalSky(
+    pmSource *source,
+    psStatsOptions statsOptions,
+    psF32 Radius)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_IMAGE_NON_NULL(source->pixels, false);
+    PS_ASSERT_IMAGE_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_INT_POSITIVE(Radius, false);
+    PS_ASSERT_INT_NONNEGATIVE(Radius, false);
+
+    psImage *image = source->pixels;
+    psImage *mask  = source->mask;
+    pmPeak *peak  = source->peak;
+    psRegion srcRegion;
+
+    srcRegion = psRegionForSquare(peak->x, peak->y, Radius);
+    srcRegion = psRegionForImage(mask, srcRegion);
+
+    psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED);
     psStats *myStats = psStatsAlloc(statsOptions);
-    myStats = psImageStats(myStats, subImage, subImageMask, 1);
-
-    //
-    // Create the output mySource, and set appropriate members.
-    //
-    pmSource *mySource = pmSourceAlloc();
-    mySource->peak = (pmPeak *) peak;
-    mySource->moments = pmMomentsAlloc();
+    myStats = psImageStats(myStats, image, mask, 0xff);
+    psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED);
+
     psF64 tmpF64;
     p_psGetStatValue(myStats, &tmpF64);
-    mySource->moments->Sky = (psF32) tmpF64;
-    mySource->pixels = subImage;
-    mySource->mask = subImageMask;
-
-    //
-    // Free things.  XXX: This should be static memory.
-    //
     psFree(myStats);
 
-    return(mySource);
-}
-
-/******************************************************************************
-bool checkRadius(*peak, radius, x, y): private function which simply
-determines if the (x, y) point is within the radius of the specified peak.
- 
-XXX: macro this for performance.
-*****************************************************************************/
-static bool checkRadius(pmPeak *peak,
-                        psF32 radius,
-                        psS32 x,
-                        psS32 y)
-{
-    if (PS_SQR(radius) >= (psF32) (PS_SQR(x - peak->x) + PS_SQR(y - peak->y))) {
-        return(true);
-    }
-
-    return(false);
+    if (isnan(tmpF64)) {
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return(false);
+    }
+    source->moments = pmMomentsAlloc();
+    source->moments->Sky = (psF32) tmpF64;
+    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
+    return (true);
 }
 
@@ -770 +688 @@
                          psF32 y)
 {
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     /// XXX EAM should compare with hypot (x,y) for speed
     if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) {
@@ -775 +694 @@
     }
 
+    psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
     return(false);
 }
@@ -787 +707 @@
     pmSource->peak
     pmSource->pixels
+    pmSource->weight
+    pmSource->mask
  
 XXX: The peak calculations are done in image coords, not subImage coords.
  
-XXX: mask values?
-*****************************************************************************/
-pmSource *pmSourceMoments(pmSource *source,
-                          psF32 radius)
-{
-    PS_ASSERT_PTR_NON_NULL(source, NULL);
-    PS_ASSERT_PTR_NON_NULL(source->peak, NULL);
-    PS_ASSERT_PTR_NON_NULL(source->pixels, NULL);
-    PS_FLOAT_COMPARE(0.0, radius, NULL);
+XXX EAM : this version clips input pixels on S/N
+XXX EAM : this version returns false for several reasons
+*****************************************************************************/
+# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
+
+bool pmSourceMoments(pmSource *source,
+                     psF32 radius)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
+    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
 
     //
@@ -816 +743 @@
     // XY  = SUM (x - xc)*(y - yc)*(z - sky)
     //
-    psF32 Sum = 0.0;
     psF32 peakPixel = -PS_MAX_F32;
     psS32 numPixels = 0;
+    psF32 Sum = 0.0;
     psF32 X1 = 0.0;
     psF32 Y1 = 0.0;
@@ -824 +751 @@
     psF32 Y2 = 0.0;
     psF32 XY = 0.0;
-    psF32 x = 0;
-    psF32 y = 0;
-    //
+    psF32 x  = 0;
+    psF32 y  = 0;
+    psF32 R2 = PS_SQR(radius);
+
+    psF32 xPeak = source->peak->x;
+    psF32 yPeak = source->peak->y;
+
     // XXX why do I get different results for these two methods of finding Sx?
     // XXX Sx, Sy would be better measured if we clip pixels close to sky
@@ -834 +765 @@
     // proceed with the moments calculation.  need to do two loops for a
     // numerically stable result.  first loop: get the sums.
-    //
+    // XXX EAM : mask == 0 is valid
+
     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] != 0)) {
-                psS32 imgColCoord = col + source->pixels->col0;
-                psS32 imgRowCoord = row + source->pixels->row0;
-                if (checkRadius(source->peak,
-                                radius,
-                                imgColCoord,
-                                imgRowCoord)) {
-                    psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
-                    psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
-                    psF32 pDiff = source->pixels->data.F32[row][col] - sky;
-
-                    Sum+= pDiff;
-                    X1+= xDiff * pDiff;
-                    Y1+= yDiff * pDiff;
-                    XY+= xDiff * yDiff * pDiff;
-
-                    X2+= PS_SQR(xDiff) * pDiff;
-                    Y2+= PS_SQR(yDiff) * pDiff;
-
-                    if (source->pixels->data.F32[row][col] > peakPixel) {
-                        peakPixel = source->pixels->data.F32[row][col];
-                    }
-                    numPixels++;
-                }
-            }
-        }
-    }
+            if ((source->mask != NULL) && (source->mask->data.U8[row][col]))
+                continue;
+
+            psF32 xDiff = col + source->pixels->col0 - xPeak;
+            psF32 yDiff = row + source->pixels->row0 - yPeak;
+
+            // XXX EAM : calculate xDiff, yDiff up front;
+            //           radius is just a function of (xDiff, yDiff)
+            if (!VALID_RADIUS(xDiff, yDiff, R2))
+                continue;
+
+            psF32 pDiff = source->pixels->data.F32[row][col] - sky;
+
+            // 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)
+                continue;
+
+            Sum += pDiff;
+            X1  += xDiff * pDiff;
+            Y1  += yDiff * pDiff;
+            XY  += xDiff * yDiff * pDiff;
+
+            X2  += PS_SQR(xDiff) * pDiff;
+            Y2  += PS_SQR(yDiff) * pDiff;
+
+            peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel);
+            numPixels++;
+        }
+    }
+    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
+    if ((numPixels < 3) || (Sum <= 0)) {
+        psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return (false);
+    }
+
+    psTrace (".psModules.pmSourceMoments", 5,
+             "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
+             sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
 
     //
@@ -872 +817 @@
     x = X1/Sum;
     y = Y1/Sum;
-    source->moments->x = x + ((psF32) source->peak->x);
-    source->moments->y = y + ((psF32) source->peak->y);
-
-    source->moments->Sxy = XY/Sum;
+    if ((fabs(x) > radius) || (fabs(y) > radius)) {
+        psTrace (".psModules.pmSourceMoments", 5,
+                 "large centroid swing; invalid peak %d, %d\n",
+                 source->peak->x, source->peak->y);
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return (false);
+    }
+
+    source->moments->x = x + xPeak;
+    source->moments->y = y + yPeak;
+
+    // XXX EAM : Sxy needs to have x*y subtracted
+    source->moments->Sxy = XY/Sum - x*y;
     source->moments->Sum = Sum;
     source->moments->Peak = peakPixel;
@@ -883 +837 @@
     source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
     source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
-    return(source);
-
-    // XXX EAM : the following code should be the same as above, but it is not very stable: ignore it
-    # if (0)
-        //
-        // second loop: get the difference sums
-        //
-        X2 = Y2 = 0;
-    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] != 0)) {
-                psS32 imgColCoord = col + source->pixels->col0;
-                psS32 imgRowCoord = row + source->pixels->row0;
-                if (checkRadius(source->peak,
-                                radius,
-                                imgColCoord,
-                                imgRowCoord)) {
-                    psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
-                    psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
-                    psF32 pDiff = source->pixels->data.F32[row][col] - sky;
-
-                    Sum+= pDiff;
-                    X2+= PS_SQR(xDiff - x) * pDiff;
-                    Y2+= PS_SQR(yDiff - y) * pDiff;
-                }
-            }
-        }
-    }
-
-    //
-    // second moment X = sqrt (X2/Sum)
-    //
-    source->moments->Sx = (X2/Sum);
-    source->moments->Sy = (Y2/Sum);
-    return(source);
-    # endif
+
+    psTrace (".psModules.pmSourceMoments", 4,
+             "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n",
+             sky, Sum, source->moments->x, source->moments->y,
+             source->moments->Sx, source->moments->Sy, source->moments->Sxy);
+
+    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
+    return(true);
 }
 
@@ -924 +850 @@
 int pmComparePeakAscend (const void **a, const void **b)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     pmPeak *A = *(pmPeak **)a;
     pmPeak *B = *(pmPeak **)b;
@@ -930 +857 @@
 
     diff = A->counts - B->counts;
-    if (diff < FLT_EPSILON)
+    if (diff < FLT_EPSILON) {
+        psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);
         return (-1);
-    if (diff > FLT_EPSILON)
+    } else if (diff > FLT_EPSILON) {
+        psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);
         return (+1);
+    }
+    psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);
     return (0);
 }
@@ -939 +870 @@
 int pmComparePeakDescend (const void **a, const void **b)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     pmPeak *A = *(pmPeak **)a;
     pmPeak *B = *(pmPeak **)b;
@@ -945 +877 @@
 
     diff = A->counts - B->counts;
-    if (diff < FLT_EPSILON)
+    if (diff < FLT_EPSILON) {
+        psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);
         return (+1);
-    if (diff > FLT_EPSILON)
+    } else if (diff > FLT_EPSILON) {
+        psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);
         return (-1);
+    }
+    psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);
     return (0);
 }
 
 /******************************************************************************
-pmSourceRoughClass(source, metadata): make a guess at the source
-classification.
- 
-XXX: This is not useable code, as of the release date.  There remains a fair
-bit of coding to be completed.
- 
-XXX: The sigX and sigY stuff in the SDRS is unclear.
- 
-XXX: How can this function ever return FALSE?
-*****************************************************************************/
-
-# define NPIX 10
-# define SCALE 0.1
-
-// XXX I am ignore memory freeing issues (EAM)
-bool pmSourceRoughClass(psArray *sources, psMetadata *metadata)
-{
-    PS_ASSERT_PTR_NON_NULL(sources, false);
-    PS_ASSERT_PTR_NON_NULL(metadata, false);
-    psBool rc = true;
+pmSourcePSFClump(source, metadata): Find the likely PSF clump in the 
+sigma-x, sigma-y plane. return 0,0 clump in case of error. 
+*****************************************************************************/
+
+// XXX EAM include a S/N cutoff in selecting the sources?
+pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+
+    # define NPIX 10
+    # define SCALE 0.1
+
     psArray *peaks  = NULL;
-    psF32 clumpX = 0.0;
-    psF32 clumpDX = 0.0;
-    psF32 clumpY = 0.0;
-    psF32 clumpDY = 0.0;
+    pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0};
+    pmPSFClump psfClump = emptyClump;
+
+    PS_ASSERT_PTR_NON_NULL(sources, emptyClump);
+    PS_ASSERT_PTR_NON_NULL(metadata, emptyClump);
 
     // find the sigmaX, sigmaY clump
@@ -986 +915 @@
 
         // construct a sigma-plane image
+        // psImageAlloc does zero the data
         splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);
+        for (int i = 0; i < splane->numRows; i++)
+        {
+            memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32));
+        }
 
         // place the sources in the sigma-plane image (ignore 0,0 values?)
@@ -992 +926 @@
         {
             pmSource *tmpSrc = (pmSource *) sources->data[i];
-            PS_ASSERT_PTR_NON_NULL(tmpSrc, false); // just skip this one?
-            PS_ASSERT_PTR_NON_NULL(tmpSrc->moments, false); // just skip this one?
+            if (tmpSrc == NULL) {
+                continue;
+            }
+            if (tmpSrc->moments == NULL) {
+                continue;
+            }
 
             // Sx,Sy are limited at 0.  a peak at 0,0 is artificial
@@ -1022 +960 @@
         psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2);
 
-    }
+        psFree (splane);
+        psFree (stats);
+
+    }
+    // XXX EAM : possible errors:
+    //           1) no peak in splane
+    //           2) no significant peak in splane
 
     // measure statistics on Sx, Sy if Sx, Sy within range of clump
@@ -1036 +980 @@
         psArraySort (peaks, pmComparePeakDescend);
         clump = peaks->data[0];
-        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d\n", clump->x, clump->y);
+        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts);
 
         // define section window for clump
@@ -1077 +1021 @@
 
         stats = psVectorStats (stats, tmpSx, NULL, NULL, 0);
-        clumpX  = stats->clippedMean;
-        clumpDX = stats->clippedStdev;
+        psfClump.X  = stats->clippedMean;
+        psfClump.dX = stats->clippedStdev;
 
         stats = psVectorStats (stats, tmpSy, NULL, NULL, 0);
-        clumpY  = stats->clippedMean;
-        clumpDY = stats->clippedStdev;
-
-        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump  X,  Y: %f, %f\n", clumpX, clumpY);
-        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", clumpDX, clumpDY);
+        psfClump.Y  = stats->clippedMean;
+        psfClump.dY = stats->clippedStdev;
+
+        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
+        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);
         // these values should be pushed on the metadata somewhere
-    }
-
-    int Nsat   = 0;
-    int Ngal   = 0;
-    int Nfaint = 0;
-    int Nstar  = 0;
-    int Npsf   = 0;
-    int Ncr    = 0;
+
+        psFree (stats);
+        psFree (peaks);
+        psFree (tmpSx);
+        psFree (tmpSy);
+    }
+
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
+    return (psfClump);
+}
+
+/******************************************************************************
+pmSourceRoughClass(source, metadata): make a guess at the source
+classification.
+ 
+XXX: push the clump info into the metadata?
+ 
+XXX: How can this function ever return FALSE?
+ 
+XXX EAM : add the saturated mask value to metadata 
+*****************************************************************************/
+
+bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+
+    psBool rc = true;
+
+    int Nsat     = 0;
+    int Ngal     = 0;
+    int Nstar    = 0;
+    int Npsf     = 0;
+    int Ncr      = 0;
+    int Nsatstar = 0;
+
     psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
     starsn->n = 0;
+
+    // check return status value (do these exist?)
+    bool status;
+    psF32 RDNOISE    = psMetadataLookupF32 (&status, metadata, "RDNOISE");
+    psF32 GAIN       = psMetadataLookupF32 (&status, metadata, "GAIN");
+    psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
+    // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
 
     // XXX allow clump size to be scaled relative to sigmas?
@@ -1109 +1087 @@
         psF32 sigY = tmpSrc->moments->Sy;
 
-        // check return status value (do these exist?)
-        bool status;
-        psF32 RDNOISE  = psMetadataLookupF32 (&status, metadata, "RDNOISE");
-        psF32 GAIN     = psMetadataLookupF32 (&status, metadata, "GAIN");
-        psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
-
-        psF32 PSF_SN_LIM   = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
-        psF32 FAINT_SN_LIM = psMetadataLookupF32 (&status, metadata, "FAINT_SN_LIM");
-
-        // saturated object (star or single pixel not distinguished)
-        if (tmpSrc->moments->Peak > SATURATE) {
-            tmpSrc->type |= PS_SOURCE_SATURATED;
+        // calculate and save signal-to-noise estimates
+        psF32 S  = tmpSrc->moments->Sum;
+        psF32 A  = 4 * M_PI * sigX * sigY;
+        psF32 B  = tmpSrc->moments->Sky;
+        psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN)));
+        psF32 SN = (S * sqrt(GAIN) / RT);
+        tmpSrc->moments->SN = SN;
+
+        // 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);
+
+        // saturated star (size consistent with PSF or larger)
+        // Nsigma should be user-configured parameter
+        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
+        if ((Nsatpix > 1) && big) {
+            tmpSrc->type |= PM_SOURCE_SATSTAR;
+            Nsatstar ++;
+            continue;
+        }
+
+        // saturated object (not a star, eg bleed trails, hot pixels)
+        if (Nsatpix > 1) {
+            tmpSrc->type |= PM_SOURCE_SATURATED;
             Nsat ++;
             continue;
         }
 
-        // too small to be stellar
-        if ((sigX < (clumpX - clumpDX)) || (sigY < (clumpY - clumpDY))) {
-            tmpSrc->type |= PS_SOURCE_DEFECT;
+        // likely defect (too small to be stellar) (push out to 3 sigma)
+        // low S/N objects which are small are probably stellar
+        // only set candidate defects if
+        if ((sigX < 0.05) || (sigY < 0.05)) {
+            tmpSrc->type |= PM_SOURCE_DEFECT;
             Ncr ++;
             continue;
         }
 
-        // possible galaxy
-        if ((sigX > (clumpX + clumpDX)) || (sigY > (clumpY + clumpDY))) {
-            tmpSrc->type |= PS_SOURCE_GALAXY;
+        // likely unsaturated galaxy (too large to be stellar)
+        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
+            tmpSrc->type |= PM_SOURCE_GALAXY;
             Ngal ++;
             continue;
@@ -1140 +1136 @@
 
         // the rest are probable stellar objects
-        psF32 S  = tmpSrc->moments->Sum;
-        psF32 A  = M_PI * sigX * sigY;
-        psF32 B  = tmpSrc->moments->Sky;
-        psF32 RT = sqrtf(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrtf(GAIN)));
-        psF32 SN = (S * sqrtf(GAIN) / RT);
-
         starsn->data.F32[starsn->n] = SN;
         starsn->n ++;
         Nstar ++;
 
-        // faint star
-        if (SN < FAINT_SN_LIM) {
-            tmpSrc->type |= PS_SOURCE_FAINTSTAR;
-            Nfaint ++;
-            continue;
-        }
-
-        // PSF star
-        if (SN > PSF_SN_LIM) {
-            tmpSrc->type |= PS_SOURCE_PSFSTAR;
+        // PSF star (within 1.5 sigma of clump center
+        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
+        if ((SN > PSF_SN_LIM) && (radius < 1.5)) {
+            tmpSrc->type |= PM_SOURCE_PSFSTAR;
             Npsf ++;
             continue;
@@ -1165 +1149 @@
 
         // random type of star
-        tmpSrc->type |= PS_SOURCE_OTHER;
+        tmpSrc->type |= PM_SOURCE_OTHER;
     }
 
@@ -1173 +1157 @@
         stats = psVectorStats (stats, starsn, NULL, NULL, 0);
         psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);
-    }
-
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:  %3d\n", Nstar);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:   %3d\n", Npsf);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nfaint: %3d\n", Nfaint);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal:   %3d\n", Ngal);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:   %3d\n", Nsat);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr:    %3d\n", Ncr);
-
+        psFree (stats);
+        psFree (starsn);
+    }
+
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:    %3d\n", Nstar);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:     %3d\n", Npsf);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal:     %3d\n", Ngal);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:     %3d\n", Nsat);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr:      %3d\n", Ncr);
+
+    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);
 }
 
+/** 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
+ * a square window (of full width 2Radius+1) centered on the pixel which contains
+ * the given coordinate, in the frame of the readout. The window is defined to
+ * have limits which are valid within the boundary of the readout image, thus if
+ * the radius would fall outside the image pixels, the subimage is truncated to
+ * only consist of valid pixels. If readout->mask or readout->weight are not
+ * NULL, matching subimages are defined for those images as well. This function
+ * fails if no valid pixels can be defined (x or y less than Radius, for
+ * 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(
+    pmSource *mySource,                 ///< Add comment.
+    pmReadout *readout,                 ///< Add comment.
+    psF32 x,                            ///< Add comment.
+    psF32 y,                            ///< Add comment.
+    psF32 Radius)                       ///< Add comment.
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented.  Returning FALSE.\n");
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
+    return(false);
+}
+
 /******************************************************************************
 pmSourceSetPixelsCircle(source, image, radius)
  
-XXX: Why boolean output?
- 
-XXX: Why are we checking source->moments for NULL?  Should the circle be
-     centered on the centroid or the peak?
- 
-XXX: The circle will have a diameter of (1+radius).  This is different from
-     the pmSourceSetLocal() function.
+XXX: This was replaced by DefinePixels in SDRS.  Remove it.
 *****************************************************************************/
 bool pmSourceSetPixelsCircle(pmSource *source,
@@ -1200 +1211 @@
                              psF32 radius)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_IMAGE_NON_NULL(image, false);
     PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_PTR_NON_NULL(source->moments, false);
-    // PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_FLOAT_COMPARE(0.0, radius, false);
 
@@ -1217 +1229 @@
     psS32 SubImageCenterCol = source->peak->x;
     // XXX EAM : for the circle to stay on the image
+    // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
     psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - radiusS32);
-    psS32 SubImageEndRow    = PS_MIN (image->numRows - 1, SubImageCenterRow + radiusS32);
+    psS32 SubImageEndRow    = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1);
     psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - radiusS32);
-    psS32 SubImageEndCol    = PS_MIN (image->numCols - 1, SubImageCenterCol + radiusS32);
-
-    // XXX EAM : this should not be needed: we can never hit this error
-    # if (1)
-
-        if (SubImageStartRow < 0) {
-            psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
-                    SubImageStartRow);
-            return(false);
-        }
-    if (SubImageEndRow >= image->numRows) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
-                SubImageEndRow);
-        return(false);
-    }
-    if (SubImageStartCol < 0) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n",
-                SubImageStartCol);
-        return(false);
-    }
-    if (SubImageEndCol >= image->numCols) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
-                SubImageEndCol);
-        return(false);
-    }
-    # endif
+    psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1);
 
     // XXX: Must recycle image.
@@ -1255 +1243 @@
         psFree(source->pixels);
     }
-    // XXX: Must fix this.  psImageSubset() has different parameters in latest CVS.
-    //    source->pixels = psImageSubset((psImage *) image,
-    //                                   SubImageStartCol,
-    //                                   SubImageStartRow,
-    //                                   SubImageEndCol,
-    //                                   SubImageEndRow);
+    source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,
+                                   SubImageStartRow,
+                                   SubImageEndCol,
+                                   SubImageEndRow));
 
     // XXX: Must recycle image.
@@ -1268 +1254 @@
     source->mask = psImageAlloc(source->pixels->numCols,
                                 source->pixels->numRows,
-                                PS_TYPE_F32);
+                                PS_TYPE_U8); // XXX EAM : type was F32
 
     //
@@ -1287 +1273 @@
         }
     }
+    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
     return(true);
 }
 
 /******************************************************************************
-pmSourceModelGuess(source, image, model): This function allocates a new
-pmModel structure and stores it in the pmSource data structure specified in
-the argument list.  The model type is specified in the argument list.  The
-params array in that pmModel structure are allocated, and then set to the
-appropriate values.  This function returns true if everything was successful.
- 
-XXX: Many of the initial parameters are set to 0.0 since I don't know what
-the appropiate initial guesses are.
- 
-XXX: The image argument is redundant.
+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.
- 
-XXX: The source->models member used to be allocated here.  Now I assume
-->modelPSF should be allocated
-*****************************************************************************/
-bool pmSourceModelGuess(pmSource *source,
-                        const psImage *image,
-                        pmModelType model)
-{
-    PS_ASSERT_PTR_NON_NULL(source, false);
+*****************************************************************************/
+pmModel *pmSourceModelGuess(pmSource *source,
+                            pmModelType modelType)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(source->moments, false);
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
-    PS_ASSERT_IMAGE_NON_NULL(image, false);
-    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
-    if (source->modelPSF != NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->modelPSF was non-NULL; calling psFree(source->modelPSF).\n");
-        psFree(source->modelPSF);
-    }
-    if (!((model == PS_MODEL_GAUSS) ||
-            (model == PS_MODEL_PGAUSS) ||
-            (model == PS_MODEL_WAUSS) ||
-            (model == PS_MODEL_TWIST_GAUSS) ||
-            (model == PS_MODEL_SERSIC) ||
-            (model == PS_MODEL_SERSIC_CORE))) {
-        psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
-        return(false);
-    }
-
-    source->modelPSF = pmModelAlloc(model);
-
-    psVector *params = source->modelPSF->params;
-
-    switch (model) {
-    case PS_MODEL_GAUSS:
-        params->data.F32[0] = source->moments->Sky;
-        params->data.F32[1] = source->peak->counts - source->moments->Sky;
-        params->data.F32[2] = source->moments->x;
-        params->data.F32[3] = source->moments->y;
-        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
-        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
-        params->data.F32[6] = source->moments->Sxy;
-        return(true);
-
-    case PS_MODEL_PGAUSS:
-        params->data.F32[0] = source->moments->Sky;
-        params->data.F32[1] = source->peak->counts - source->moments->Sky;
-        params->data.F32[2] = source->moments->x;
-        params->data.F32[3] = source->moments->y;
-        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
-        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
-        params->data.F32[6] = source->moments->Sxy;
-        return(true);
-
-    case PS_MODEL_WAUSS:
-        params->data.F32[0] = source->moments->Sky;
-        params->data.F32[1] = source->peak->counts - source->moments->Sky;
-        params->data.F32[2] = source->moments->x;
-        params->data.F32[3] = source->moments->y;
-        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
-        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
-        params->data.F32[6] = source->moments->Sxy;
-        // XXX: What are these?
-        // source->modelPSF->params[7] = B2;
-        // source->modelPSF->params[8] = B3;
-        return(true);
-
-        // XXX EAM : I might drop this model (or rather, replace it)
-    case PS_MODEL_TWIST_GAUSS:
-        params->data.F32[0] = source->moments->Sky;
-        params->data.F32[1] = source->peak->counts - source->moments->Sky;
-        params->data.F32[2] = source->moments->x;
-        params->data.F32[3] = source->moments->y;
-        // XXX: What are these?
-        // params->data.F32[4] = SxInner;
-        // params->data.F32[5] = SyInner;
-        // params->data.F32[6] = SxyInner;
-        // params->data.F32[7] = SxOuter;
-        // params->data.F32[8] = SyOuter;
-        // params->data.F32[9] = SxyOuter;
-        // params->data.F32[10] = N;
-        return(true);
-
-    case PS_MODEL_SERSIC:
-        params->data.F32[0] = source->moments->Sky;
-        params->data.F32[1] = source->peak->counts - source->moments->Sky;
-        params->data.F32[2] = source->moments->x;
-        params->data.F32[3] = source->moments->y;
-        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
-        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
-        params->data.F32[6] = source->moments->Sxy;
-        // XXX: What are these?
-        //params->data.F32[7] = Nexp;
-        return(true);
-
-    case PS_MODEL_SERSIC_CORE:
-        params->data.F32[0] = source->moments->Sky;
-        params->data.F32[1] = source->peak->counts - source->moments->Sky;
-        params->data.F32[2] = source->moments->x;
-        params->data.F32[3] = source->moments->y;
-        // XXX: What are these?
-        // params->data.F32[4] SxInner;
-        // params->data.F32[5] SyInner;
-        // params->data.F32[6] SxyInner;
-        // params->data.F32[7] Zd;
-        // params->data.F32[8] SxOuter;
-        // params->data.F32[9] SyOuter;
-        // params->data.F32[10] = SxyOuter;
-        // params->data.F32[11] = Nexp;
-        return(true);
-
-    default:
-        psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
-        return(false);
-    }
+
+    pmModel *model = pmModelAlloc(modelType);
+
+    pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType);
+    modelGuessFunc(model, source);
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
+    return(model);
 }
 
@@ -1441 +1323 @@
 testing.  Try to reproduce that and debug.
 *****************************************************************************/
-static psF32 evalModel(pmSource *src,
-                       psU32 row,
-                       psU32 col)
-{
-    PS_ASSERT_PTR_NON_NULL(src, false);
-    PS_ASSERT_PTR_NON_NULL(src->modelPSF, false);
-    PS_ASSERT_PTR_NON_NULL(src->modelPSF->params, false);
-
-    // XXX: The following step will not be necessary if the modelPSF->params
-    // member is a psVector.  Suggest to IfA.
-
-    // XXX EAM: done: modelPSF->params is now a vector
-    psVector *params = src->modelPSF->params;
-
-    //
+
+// XXX EAM : I have made this a public function
+// XXX EAM : this now uses a pmModel as the input
+// XXX EAM : it was using src->type to find the model, not model->type
+psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row)
+{
+    psTrace(__func__, 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 + src->pixels->col0);
-    x->data.F32[1] = (psF32) (row + src->pixels->row0);
+    x->data.F32[0] = (psF32) (col + image->col0);
+    x->data.F32[1] = (psF32) (row + image->row0);
     psF32 tmpF;
-
-    switch (src->modelPSF->type) {
-    case PS_MODEL_GAUSS:
-        tmpF = pmMinLM_Gauss2D(NULL, params, x);
-        break;
-    case PS_MODEL_PGAUSS:
-        tmpF = pmMinLM_PsuedoGauss2D(NULL, params, x);
-        break;
-    case PS_MODEL_TWIST_GAUSS:
-        tmpF = pmMinLM_TwistGauss2D(NULL, params, x);
-        break;
-    case PS_MODEL_WAUSS:
-        tmpF = pmMinLM_Wauss2D(NULL, params, x);
-        break;
-    case PS_MODEL_SERSIC:
-        tmpF = pmMinLM_Sersic(NULL, params, x);
-        break;
-    case PS_MODEL_SERSIC_CORE:
-        tmpF = pmMinLM_SersicCore(NULL, params, x);
-        break;
-    default:
-        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
-        return(NAN);
-    }
+    pmModelFunc modelFunc;
+
+    modelFunc = pmModelFunc_GetFunction (model->type);
+    tmpF = modelFunc (NULL, model->params, x);
     psFree(x);
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmpF);
 }
@@ -1508 +1367 @@
                        psU32 dir)
 {
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     //
     // Convert coords to subImage space.
@@ -1517 +1377 @@
     if (!((0 <= subCol) && (subCol < source->pixels->numCols))) {
         psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range");
+        psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);
         return(NAN);
     }
     if (!((0 <= subRow) && (subRow < source->pixels->numRows))) {
+        psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);
         psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range");
         return(NAN);
     }
 
-    psF32 oldValue = evalModel(source, subRow, subCol);
+    // XXX EAM : i changed this to match pmModelEval above, but see
+    // XXX EAM   the note below in pmSourceContour
+    psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
     if (oldValue == level) {
+        psTrace(__func__, 4, "---- %s() end ----\n", __func__);
         return(((psF32) (subCol + source->pixels->col0)));
     }
@@ -1545 +1410 @@
 
     while (subCol != lastColumn) {
-        psF32 newValue = evalModel(source, subRow, subCol);
+        psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
         if (oldValue == level) {
+            psTrace(__func__, 4, "---- %s() end ----\n", __func__);
             return((psF32) (subCol + source->pixels->col0));
         }
@@ -1552 +1418 @@
         if ((newValue <= level) && (level <= oldValue)) {
             // This is simple linear interpolation.
+            psTrace(__func__, 4, "---- %s() end ----\n", __func__);
             return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) );
         }
@@ -1557 +1424 @@
         if ((oldValue <= level) && (level <= newValue)) {
             // This is simple linear interpolation.
+            psTrace(__func__, 4, "---- %s() end ----\n", __func__);
             return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) );
         }
@@ -1563 +1431 @@
     }
 
+    psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);
     return(NAN);
 }
@@ -1576 +1445 @@
 XXX: What is mode?
 XXX: The top, bottom of the contour is not correctly determined.
+XXX EAM : this functions is using the model for the contour, but it should
+          be using only the image counts
 *****************************************************************************/
 psArray *pmSourceContour(pmSource *source,
@@ -1582 +1453 @@
                          pmContourType mode)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_PTR_NON_NULL(image, false);
@@ -1587 +1459 @@
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->modelPSF, false);
+    PS_ASSERT_PTR_NON_NULL(source->modelFLT, false);
+    // XXX EAM : what is the purpose of modelPSF/modelFLT?
 
     //
@@ -1610 +1483 @@
             psFree(xVec);
             psFree(yVec);
+            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
             return(NULL);
             //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
@@ -1622 +1496 @@
             psFree(xVec);
             psFree(yVec);
+            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
             return(NULL);
             //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
         }
-        //printf("The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept);
+        psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept);
         xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;
 
@@ -1646 +1521 @@
             psFree(xVec);
             psFree(yVec);
+            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
             return(NULL);
             //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
@@ -1657 +1533 @@
             psFree(xVec);
             psFree(yVec);
+            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
             return(NULL);
             //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
@@ -1672 +1549 @@
     tmpArray->data[0] = (psPtr *) yVec;
     tmpArray->data[1] = (psPtr *) xVec;
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmpArray);
 }
 
-#if 0
-static psVector *minLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
-static psVector *minLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
-static psVector *minLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
-static psVector *minLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
-static psVector *minLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x);
-static psVector *minLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x);
-#endif
-
 // XXX EAM : these are better starting values, but should be available from metadata?
-#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 20
+#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
 #define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
 /******************************************************************************
-pmSourceFitModel(source, image): must create the appropiate arguments to the
+pmSourceFitModel(source, model): must create the appropiate arguments to the
 LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions.
  
 XXX: should there be a mask value?
-XXX: Probably should remove the "image" argument.
-*****************************************************************************/
-bool pmSourceFitModel(pmSource *source,
-                      const psImage *image)
-{
+XXX EAM : fit the specified model (not necessarily the one in source)
+*****************************************************************************/
+bool pmSourceFitModel_v5(pmSource *source,
+                         pmModel *model,
+                         const bool PSF)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_PTR_NON_NULL(source->moments, false);
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->modelPSF, false);
-    PS_ASSERT_IMAGE_NON_NULL(image, false);
-    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
-    psBool rc;
+    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->weight, false);
+    psBool fitStatus = true;
+    psBool onPic     = true;
+    psBool rc        = true;
+
+    // XXX EAM : is it necessary for the mask & weight to exist?  the
+    //           tests below could be conditions (!NULL)
+
+    psVector *params = model->params;
+    psVector *dparams = model->dparams;
+    psVector *paramMask = NULL;
+
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+
+    int nParams = PSF ? params->n - 4 : params->n;
 
     // find the number of valid pixels
@@ -1716 +1599 @@
             }
         }
+    }
+    if (count <  nParams + 1) {
+        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return(false);
     }
 
@@ -1733 +1621 @@
                 x->data[tmpCnt] = (psPtr *) coord;
                 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
-
-                // XXX EAM : this is approximate: need to apply the gain and rdnoise
-                yErr->data.F32[tmpCnt] = sqrt(PS_MAX(1, source->pixels->data.F32[i][j]));
+                yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]);
+                // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then
+                //           the minimization function calculates sq()
                 tmpCnt++;
             }
@@ -1744 +1632 @@
                             PM_SOURCE_FIT_MODEL_TOLERANCE);
 
-    psVector *params = source->modelPSF->params;
-
-    switch (source->modelPSF->type) {
-    case PS_MODEL_GAUSS:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Gauss2D);
-        break;
-    case PS_MODEL_PGAUSS:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_PsuedoGauss2D);
-        break;
-    case PS_MODEL_TWIST_GAUSS:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Wauss2D);
-        break;
-    case PS_MODEL_WAUSS:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_TwistGauss2D);
-        break;
-    case PS_MODEL_SERSIC:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Sersic);
-        break;
-    case PS_MODEL_SERSIC_CORE:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_SersicCore);
-        break;
-    default:
-        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
-        rc = false;
-    }
+    // PSF model only fits first 4 parameters, FLT model fits all
+    if (PSF) {
+        paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
+        for (int i = 0; i < 4; i++) {
+            paramMask->data.U8[i] = 0;
+        }
+        for (int i = 4; i < paramMask->n; i++) {
+            paramMask->data.U8[i] = 1;
+        }
+    }
+
+    // XXX EAM : covar must be F64?
+    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+
+    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
+    for (int i = 0; i < dparams->n; i++) {
+        if ((paramMask != NULL) && paramMask->data.U8[i])
+            continue;
+        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
+    }
+
     // 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
-    source->modelPSF->chisq = myMin->value;
-    source->modelPSF->nDOF  = y->n - params->n;
-    source->modelPSF->nIter = myMin->iter;
-
+    model->chisq = myMin->value;
+    model->nIter = myMin->iter;
+    model->nDOF  = y->n - nParams;
+
+    // XXX EAM get the Gauss-Newton distance for fixed model parameters
+    if (paramMask != NULL) {
+        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
+        psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
+        for (int i = 0; i < dparams->n; i++) {
+            if (!paramMask->data.U8[i])
+                continue;
+            dparams->data.F32[i] = delta->data.F64[i];
+        }
+        psFree (delta);
+    }
+
+    psFree(x);
+    psFree(y);
+    psFree(yErr);
+    psFree(myMin);
+    psFree(covar);
+    psFree(paramMask);
+
+    rc = (onPic && fitStatus);
+    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
+    return(rc);
+}
+
+// XXX EAM : new version with parameter range limits and weight enhancement
+bool pmSourceFitModel (pmSource *source,
+                       pmModel *model,
+                       const bool PSF)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->moments, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
+    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->weight, false);
+
+    // XXX EAM : is it necessary for the mask & weight to exist?  the
+    //           tests below could be conditions (!NULL)
+
+    psBool fitStatus = true;
+    psBool onPic     = true;
+    psBool rc        = true;
+    psF32  Ro, ymodel;
+
+    psVector *params = model->params;
+    psVector *dparams = model->dparams;
+    psVector *paramMask = NULL;
+
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+
+    // XXX EAM : I need to use the sky value to constrain the weight model
+    int nParams = PSF ? params->n - 4 : params->n;
+    psF32 So = params->data.F32[0];
+
+    // find the number of valid pixels
+    // XXX EAM : this loop and the loop below could just be one pass
+    //           using the psArrayAdd and psVectorExtend functions
+    psS32 count = 0;
+    for (psS32 i = 0; i < source->pixels->numRows; i++) {
+        for (psS32 j = 0; j < source->pixels->numCols; j++) {
+            if (source->mask->data.U8[i][j] == 0) {
+                count++;
+            }
+        }
+    }
+    if (count <  nParams + 1) {
+        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return(false);
+    }
+
+    // construct the coordinate and value entries
+    psArray *x = psArrayAlloc(count);
+    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
+    psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
+    psS32 tmpCnt = 0;
+    for (psS32 i = 0; i < source->pixels->numRows; i++) {
+        for (psS32 j = 0; j < source->pixels->numCols; j++) {
+            if (source->mask->data.U8[i][j] == 0) {
+                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+                // XXX: Convert i/j to image space:
+                // XXX EAM: coord order is (x,y) == (col,row)
+                coord->data.F32[0] = (psF32) (j + source->pixels->col0);
+                coord->data.F32[1] = (psF32) (i + source->pixels->row0);
+                x->data[tmpCnt] = (psPtr *) coord;
+                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
+
+                // compare observed flux to model flux to adjust weight
+                ymodel = modelFunc (NULL, model->params, coord);
+
+                // this test enhances the weight based on deviation from the model flux
+                Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
+
+                // psMinimizeLMChi2_EAM takes wt = 1/dY^2
+                if (source->weight->data.F32[i][j] == 0) {
+                    yErr->data.F32[tmpCnt] = 0.0;
+                } else {
+                    yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro);
+                }
+                tmpCnt++;
+            }
+        }
+    }
+
+    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
+                            PM_SOURCE_FIT_MODEL_TOLERANCE);
+
+    // PSF model only fits first 4 parameters, FLT model fits all
+    if (PSF) {
+        paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
+        for (int i = 0; i < 4; i++) {
+            paramMask->data.U8[i] = 0;
+        }
+        for (int i = 4; i < paramMask->n; i++) {
+            paramMask->data.U8[i] = 1;
+        }
+    }
+
+    // XXX EAM : I've added three types of parameter range checks
+    // XXX EAM : this requires my new psMinimization functions
+    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
+    psVector *beta_lim = NULL;
+    psVector *params_min = NULL;
+    psVector *params_max = NULL;
+
+    // XXX EAM : in this implementation, I pass in the limits with the covar matrix.
+    //           in the SDRS, I define a new psMinimization which will take these in
+    psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
+    modelLimits (&beta_lim, &params_min, &params_max);
+    for (int i = 0; i < params->n; i++) {
+        covar->data.F64[0][i] = beta_lim->data.F32[i];
+        covar->data.F64[1][i] = params_min->data.F32[i];
+        covar->data.F64[2][i] = params_max->data.F32[i];
+    }
+
+    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
+    for (int i = 0; i < dparams->n; i++) {
+        if ((paramMask != NULL) && paramMask->data.U8[i])
+            continue;
+        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
+    }
+
+    // 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
+    model->chisq = myMin->value;
+    model->nIter = myMin->iter;
+    model->nDOF  = y->n - nParams;
+
+    // XXX EAM get the Gauss-Newton distance for fixed model parameters
+    if (paramMask != NULL) {
+        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
+        psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc);
+        for (int i = 0; i < dparams->n; i++) {
+            if (!paramMask->data.U8[i])
+                continue;
+            dparams->data.F32[i] = delta->data.F64[i];
+        }
+    }
+
+    psFree(paramMask);
     psFree(x);
     psFree(y);
     psFree(myMin);
+
+    rc = (onPic && fitStatus);
+    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);
 }
 
-static bool sourceAddOrSubModel(psImage *image,
-                                pmSource *src,
-                                bool center,
-                                psS32 flag)
-{
-    PS_ASSERT_PTR_NON_NULL(src, false);
-    PS_ASSERT_PTR_NON_NULL(src->moments, false);
-    PS_ASSERT_PTR_NON_NULL(src->peak, false);
-    PS_ASSERT_PTR_NON_NULL(src->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(src->modelPSF, false);
+bool p_pmSourceAddOrSubModel(psImage *image,
+                             psImage *mask,
+                             pmModel *model,
+                             bool center,
+                             psS32 flag)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+
+    PS_ASSERT_PTR_NON_NULL(model, false);
     PS_ASSERT_IMAGE_NON_NULL(image, false);
     PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
 
     psVector *x = psVectorAlloc(2, PS_TYPE_F32);
-    psVector *params = src->modelPSF->params;
-
-    for (psS32 i = 0; i < src->pixels->numRows; i++) {
-        for (psS32 j = 0; j < src->pixels->numCols; j++) {
+    psVector *params = model->params;
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+    psS32 imageCol;
+    psS32 imageRow;
+
+    for (psS32 i = 0; i < image->numRows; i++) {
+        for (psS32 j = 0; j < image->numCols; j++) {
+            if ((mask != NULL) && mask->data.U8[i][j])
+                continue;
             psF32 pixelValue;
             // XXX: Should you be adding the pixels for the entire subImage,
@@ -1805 +1875 @@
             // Convert i/j to imace coord space:
             // XXX: Make sure you have col/row order correct.
-            psS32 imageRow = i + src->pixels->row0;
-            psS32 imageCol = j + src->pixels->col0;
+            // XXX EAM : 'center' option changes this
+            // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2]
+            if (center) {
+                imageCol = j - 0.5*image->numCols + model->params->data.F32[2];
+                imageRow = i - 0.5*image->numRows + model->params->data.F32[3];
+            } else {
+                imageCol = j + image->col0;
+                imageRow = i + image->row0;
+            }
 
             x->data.F32[0] = (float) imageCol;
             x->data.F32[1] = (float) imageRow;
-            switch (src->modelPSF->type) {
-            case PS_MODEL_GAUSS:
-                pixelValue = pmMinLM_Gauss2D(NULL, params, x);
-                break;
-            case PS_MODEL_PGAUSS:
-                pixelValue = pmMinLM_PsuedoGauss2D(NULL, params, x);
-                break;
-            case PS_MODEL_TWIST_GAUSS:
-                pixelValue = pmMinLM_TwistGauss2D(NULL, params, x);
-                break;
-            case PS_MODEL_WAUSS:
-                pixelValue = pmMinLM_Wauss2D(NULL, params, x);
-                break;
-            case PS_MODEL_SERSIC:
-                pixelValue = pmMinLM_Sersic(NULL, params, x);
-                break;
-            case PS_MODEL_SERSIC_CORE:
-                pixelValue = pmMinLM_SersicCore(NULL, params, x);
-                break;
-            default:
-                psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
-                psFree(x);
-                return(false);
-            }
+            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;
@@ -1841 +1897 @@
             // how to use the boolean "center" flag.
 
-            image->data.F32[imageRow][imageCol]+= pixelValue;
+            image->data.F32[i][j]+= pixelValue;
         }
     }
     psFree(x);
+    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
     return(true);
 }
@@ -1853 +1910 @@
  *****************************************************************************/
 bool pmSourceAddModel(psImage *image,
-                      pmSource *src,
+                      psImage *mask,
+                      pmModel *model,
                       bool center)
 {
-    return(sourceAddOrSubModel(image, src, center, 0));
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 0);
+    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
+    return(rc);
 }
 
@@ -1862 +1923 @@
  *****************************************************************************/
 bool pmSourceSubModel(psImage *image,
-                      pmSource *src,
+                      psImage *mask,
+                      pmModel *model,
                       bool center)
 {
-    return(sourceAddOrSubModel(image, src, center, 1));
-}
-
-
-// XXX: Put this is psConstants.h
-#define PS_VECTOR_CHECK_SIZE(VEC1, N, RVAL) \
-if (VEC1->n != N) { \
-    psError(PS_ERR_BAD_PARAMETER_SIZE, true, \
-            "psVector %s has size %d, should be %d.", \
-            #VEC1, VEC1->n, N); \
-    return(RVAL); \
-}
-
-
-/**
-   all of these object representation functions have the same form : func(*deriv, *params, *x)
- 
-   the argument "x" contains a single "x,y" coordinate pair.  The function computes the object
-   model, based on the parameters in "params" at the x,y point specified by *x, and returns the value.
-   The derivatives are also caculated and returned in the "deriv" argument.  parameter error checking is
-   skipped because speed is most important.
-**/
-
-/******************************************************************************
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = sqrt(2.0) / SigmaX;
-    params->data.F32[5] = sqrt(2.0) / SigmaY;
-    params->data.F32[6] = Sxy;
-*****************************************************************************/
-float pmMinLM_Gauss2D(
-    psVector *deriv,
-    const psVector *params,
-    const psVector *x)
-{
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
-    psF32 X  = x->data.F32[0] - params->data.F32[2];
-    psF32 Y  = x->data.F32[1] - params->data.F32[3];
-    psF32 px = params->data.F32[4]*X;
-    psF32 py = params->data.F32[5]*Y;
-    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
-    psF32 r  = exp(-z);
-    psF32 q  = params->data.F32[1]*r;
-    psF32 f  = q + params->data.F32[0];
-
-    if (deriv != NULL) {
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y);
-        deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-    }
-    return(f);
-}
-
-/******************************************************************************
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = sqrt(2) / SigmaX;
-    params->data.F32[5] = sqrt(2) / SigmaY;
-    params->data.F32[6] = Sxy;
-*****************************************************************************/
-float pmMinLM_PsuedoGauss2D(
-    psVector *deriv,
-    const psVector *params,
-    const psVector *x)
-{
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
-    psF32 X  = x->data.F32[0] - params->data.F32[2];
-    psF32 Y  = x->data.F32[1] - params->data.F32[3];
-    psF32 px = params->data.F32[4]*X;
-    psF32 py = params->data.F32[5]*Y;
-    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
-    psF32 t  = 1 + z + 0.5*z*z;
-    psF32 r  = 1.0 / (t*(1 + z/3)); /* exp (-Z) */
-    psF32 f  = params->data.F32[1]*r + params->data.F32[0];
-
-    if (deriv != NULL) {
-        // note difference from a pure gaussian: q = params->data.F32[1]*r
-        psF32 q = params->data.F32[1]*r*r*t;
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
-        deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-    }
-    return(f);
-}
-
-/******************************************************************************
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = Sx;
-    params->data.F32[5] = Sy;
-    params->data.F32[6] = Sxy;
-    params->data.F32[7] = B2;
-    params->data.F32[8] = B3;
-*****************************************************************************/
-float pmMinLM_Wauss2D(
-    psVector *deriv,
-    const psVector *params,
-    const psVector *x)
-{
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
-    psF32 X = x->data.F32[0] - params->data.F32[2];
-    psF32 Y = x->data.F32[1] - params->data.F32[2];
-    psF32 px = params->data.F32[4]*X;
-    psF32 py = params->data.F32[5]*Y;
-    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
-    psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0);
-    psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */
-    psF32 f = params->data.F32[1]*r + params->data.F32[0];
-
-    if (deriv != NULL) {
-        // note difference from gaussian: q = params->data.F32[1]*r
-        psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0));
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
-        deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-        deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t;
-        deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0;
-        // The values of 100 dampen the swing of params->data.F32[7,8] */
-    }
-    return(f);
-}
-
-// XXX: What should these be?
-#define FFACTOR 1.0
-#define FSCALE 1.0
-/******************************************************************************
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = SxInner;
-    params->data.F32[5] = SyInner;
-    params->data.F32[6] = SxyInner;
-    params->data.F32[7] = SxOuter;
-    params->data.F32[8] = SyOuter;
-    params->data.F32[9] = SxyOuter;
-    params->data.F32[10] = N;
-*****************************************************************************/
-float pmMinLM_TwistGauss2D(
-    psVector *deriv,
-    const psVector *params,
-    const psVector *x)
-{
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
-    psF32 X = x->data.F32[0] - params->data.F32[2];
-    psF32 Y = x->data.F32[1] - params->data.F32[3];
-    psF32 px1 = params->data.F32[4]*X;
-    psF32 py1 = params->data.F32[5]*Y;
-    psF32 px2 = params->data.F32[7]*X;
-    psF32 py2 = params->data.F32[8]*Y;
-    psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y;
-    psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y;
-    psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10]));
-
-    psF32 f = params->data.F32[5]*r + params->data.F32[6];
-
-    if (deriv != NULL) {
-        psF32 q1 = params->data.F32[5]*PS_SQR(r);
-        psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0));
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q1*(2.0*px1*params->data.F32[4] + params->data.F32[6]*Y) + q2*(2*px2*params->data.F32[7] + params->data.F32[9]*Y);
-        deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X);
-
-        // These fudge factors impede the growth of params->data.F32[4] beyond
-        // params->data.F32[7].
-        psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]);
-        psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
-        deriv->data.F32[4] = -2.0*q1*px1*X*f2;
-
-        // These fudge factors impede the growth of params->data.F32[5] beyond
-        // params->data.F32[8].
-        f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]);
-        f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
-        deriv->data.F32[5] = -2.0*q1*py1*Y*f2;
-        deriv->data.F32[6] = -q1*X*Y;
-        deriv->data.F32[7] = -2.0*q2*px2*X;
-        deriv->data.F32[8] = -2.0*q2*py2*Y;
-        deriv->data.F32[9] = -q2*X*Y;
-        deriv->data.F32[10] = -q1*log(z2);
-    }
-
-    return(f);
-}
-
-/******************************************************************************
-    float Sersic()
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = Sx;
-    params->data.F32[5] = Sy;
-    params->data.F32[6] = Sxy;
-    params->data.F32[7] = Nexp;
-*****************************************************************************/
-float pmMinLM_Sersic(
-    psVector *deriv,
-    const psVector *params,
-    const psVector *x)
-{
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
-    psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
-    return(0.0);
-}
-
-/******************************************************************************
-    float SersicBulge()
-    params->data.F32[0] So;
-    params->data.F32[1] Zo;
-    params->data.F32[2] Xo;
-    params->data.F32[3] Yo;
-    params->data.F32[4] SxInner;
-    params->data.F32[5] SyInner;
-    params->data.F32[6] SxyInner;
-    params->data.F32[7] Zd;
-    params->data.F32[8] SxOuter;
-    params->data.F32[9] SyOuter;
-    params->data.F32[10] = SxyOuter;
-    params->data.F32[11] = Nexp;
-*****************************************************************************/
-float pmMinLM_SersicCore(
-    psVector *deriv,
-    const psVector *params,
-    const psVector *x)
-{
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
-    psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
-    return(0.0);
-}
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 1);
+    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
+    return(rc);
+}
+
+

trunk/psModules/src/objects/pmObjects.h

@@ -1 +1 @@
 /** @file  pmObjects.h
  *
- *  This file will ...
+ * The process of finding, measuring, and classifying astronomical sources on
+ * 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 
+ * are generally low-level components which can be connected together to
+ * construct a complete object measurement suite.
  *
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-28 20:43:52 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-10-10 19:53:40 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -19 +24 @@
 #endif
 
-#include<stdio.h>
-#include<math.h>
+#include <stdio.h>
+#include <math.h>
 #include "pslib.h"
+#include "pmAstrometry.h"
+/**
+ * In the object analysis process, we will use specific mask values to mark the
+ * image pixels. The following structure defines the relevant mask values.
+ *
+ * XXX: This is probably a bad solution: we will want to set mask values
+ * outside of the PSPHOT code.  Perhaps we can set up a registered set of mask
+ * values with specific meanings that other functions can add to or define?
+ */
+enum {
+    PSPHOT_MASK_CLEAR     = 0x00,
+    PSPHOT_MASK_INVALID   = 0x01,
+    PSPHOT_MASK_SATURATED = 0x02,
+    PSPHOT_MASK_MARKED    = 0x08,
+} psphotMaskValues;
+
 
 /** pmPeakType
@@ -40 +61 @@
 } pmPeakType;
 
+
 /** 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:
  *  
  */
@@ -52 +78 @@
 pmPeak;
 
+
 /** 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
 {
-    float x;                            ///< X-coord of centroid.
-    float y;                            ///< Y-coord of centroid.
-    float Sx;                           ///< x-second moment.
-    float Sy;                           ///< y-second moment.
-    float Sxy;                          ///< xy cross moment.
-    float Sum;                          ///< Pixel sum above sky (background).
-    float Peak;                         ///< Peak counts above sky.
-    float Sky;                          ///< Sky level (background).
-    int nPixels;                        ///< Number of pixels used.
+    float x;    ///< X-coord of centroid.
+    float y;    ///< Y-coord of centroid.
+    float Sx;    ///< x-second moment.
+    float Sy;    ///< y-second moment.
+    float Sxy;    ///< xy cross moment.
+    float Sum;    ///< Pixel sum above sky (background).
+    float Peak;    ///< Peak counts above sky.
+    float Sky;    ///< Sky level (background).
+    float SN;    ///< approx signal-to-noise
+    int nPixels;   ///< Number of pixels used.
 }
 pmMoments;
 
-/** pmModelType enumeration
- *  
- */
-typedef enum {
-    PS_MODEL_GAUSS,                     ///< Regular 2-D Gaussian
-    PS_MODEL_PGAUSS,                    ///< Psuedo 2-D Gaussian
-    PS_MODEL_TWIST_GAUSS,               ///< 2-D Twisted Gaussian
-    PS_MODEL_WAUSS,                     ///< 2-D Waussian
-    PS_MODEL_SERSIC,                    ///< Sersic
-    PS_MODEL_SERSIC_CORE,               ///< Sersic Core
-    PS_MODEL_UNDEFINED                  ///< Undefined
-} pmModelType;
-
-/** pmModel data structure
- *  
- */
-// XXX: The SDRS has the "type" member of type psS32.
+
+/** 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
 {
-    pmModelType type;                   ///< Model to be used.
-    psVector *params;                   ///< Paramater values.
-    psVector *dparams;                  ///< Parameter errors.
-    float chisq;                        ///< Fit chi-squared.
-    int nDOF;                           ///< number of degrees of freedom
-    int nIter;                          ///< number of iterations to reach min
+    float X;
+    float dX;
+    float Y;
+    float dY;
+}
+pmPSFClump;
+
+typedef int pmModelType;
+#define PS_MODEL_GAUSS 0
+#define PS_MODEL_PGAUSS 1
+#define PS_MODEL_QGAUSS 2
+#define PS_MODEL_SGAUSS 3
+
+
+/** 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
+ * object. In this case, the PSF-dependent parameters are specified for the
+ * 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
+{
+    pmModelType type;   ///< Model to be used.
+    psVector *params;   ///< Paramater values.
+    psVector *dparams;   ///< Parameter errors.
+    float chisq;   ///< Fit chi-squared.
+    int nDOF;    ///< number of degrees of freedom
+    int nIter;    ///< number of iterations to reach min
+    float radius;   ///< fit radius actually used
 }
 pmModel;
 
 /** 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 {
-    PS_SOURCE_PSFSTAR,
-    PS_SOURCE_GALAXY,
-    PS_SOURCE_DEFECT,
-    PS_SOURCE_SATURATED,
-    PS_SOURCE_SATSTAR,
-    PS_SOURCE_FAINTSTAR,
-    PS_SOURCE_BRIGHTSTAR,
-    PS_SOURCE_OTHER
+    PM_SOURCE_DEFECT,                   ///< a cosmic-ray
+    PM_SOURCE_SATURATED,                ///< random saturated pixels
+
+    PM_SOURCE_SATSTAR,                  ///< a saturated star
+    PM_SOURCE_PSFSTAR,                  ///< a PSF star
+    PM_SOURCE_GOODSTAR,                 ///< a good-quality star
+
+    PM_SOURCE_POOR_FIT_PSF,             ///< poor quality PSF fit
+    PM_SOURCE_FAIL_FIT_PSF,             ///< failed to get a good PSF fit
+    PM_SOURCE_FAINTSTAR,                ///< below S/N cutoff
+
+    PM_SOURCE_GALAXY,                   ///< an extended object (galaxy)
+    PM_SOURCE_FAINT_GALAXY,             ///< a galaxy below S/N cutoff
+    PM_SOURCE_DROP_GALAXY,              ///< ?
+    PM_SOURCE_FAIL_FIT_GAL,             ///< failed on the galaxy fit
+    PM_SOURCE_POOR_FIT_GAL,             ///< poor quality galaxy fit
+
+    PM_SOURCE_OTHER,                    ///< unidentified
 } pmSourceType;
 
@@ -122 +186 @@
 typedef struct
 {
-    pmPeak *peak;                       ///< Description of peak pixel.
-    psImage *pixels;                    ///< Rectangular region including object pixels.
-    psImage *mask;                      ///< Mask which marks pixels associated with objects.
-    pmMoments *moments;                 ///< Basic moments measure for the object.
-    pmModel *modelPSF;                  ///< PSF model parameters and type
-    pmModel *modelFLT;                  ///< FLT model parameters and type
-    pmSourceType type;                  ///< Best identification of object.
+    pmPeak *peak;   ///< Description of peak pixel.
+    psImage *pixels;   ///< Rectangular region including object pixels.
+    psImage *weight;   ///< Image variance.
+    psImage *mask;   ///< Mask which marks pixels associated with objects.
+    pmMoments *moments;   ///< Basic moments measure for the object.
+    pmModel *modelPSF;   ///< PSF Model fit (parameters and type)
+    pmModel *modelFLT;   ///< FLT (floating) Model fit (parameters and type).
+    pmSourceType type;   ///< Best identification of object.
 }
 pmSource;
 
-/** pmPeak data structure
- *  
- *  
- *  
- */
-typedef struct
-{
-    psS32 type;                         ///< PSF Model in use
-    psArray *params;                    ///< Model parameters (psPolynomial2D)
-    psF32 chisq;                        ///< PSF goodness statistic
-    psS32 nPSFstars;                    ///< number of stars used to measure PSF
-}
-pmPSF;
-
-
-
+
+/** pmPeakAlloc()
+ *
+ *  @return pmPeak*    newly allocated pmPeak with all internal pointers set to NULL
+ */
 pmPeak *pmPeakAlloc(
-    int x,                              ///< Row-coordinate in image space
-    int y,                              ///< Col-coordinate in image space
-    float counts,                       ///< The value of the peak pixel
-    pmPeakType class                    ///< The type of peak pixel
-);
-
+    int x,    ///< Row-coordinate in image space
+    int y,    ///< Col-coordinate in image space
+    float counts,   ///< The value of the peak pixel
+    pmPeakType class   ///< The type of peak pixel
+);
+
+
+/** pmMomentsAlloc()
+ * 
+ */
 pmMoments *pmMomentsAlloc();
+
+
+/** pmModelAlloc()
+ * 
+ */
 pmModel *pmModelAlloc(pmModelType type);
-pmSource *pmSourceAlloc();
-
-/******************************************************************************
-pmFindVectorPeaks(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.
- *****************************************************************************/
+
+
+/** pmSourceAlloc()
+ * 
+ */
+pmSource  *pmSourceAlloc();
+
+
+/** 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
+ * a value above the threshold. Two types of special cases must be addressed.
+ * Equal value elements: If an element has the same value as the following
+ * element, it is not considered a peak. If an element has the same value as the
+ * preceding element (but not the following), then it is considered a peak. Note
+ * that this rule (arbitrarily) identifies flat regions by their trailing edge.
+ * Edge cases: At start of the vector, the element must be higher than its
+ * neighbor. At the end of the vector, the element must be higher or equal to its
+ * neighbor. These two rules again places the peak associated with a flat region
+ * which touches the image edge at the image edge. The result of this function is
+ * a vector containing the coordinates (element number) of the detected peaks
+ * (type psU32).
+ * 
+ */
 psVector *pmFindVectorPeaks(
-    const psVector *vector,             ///< The input vector (float)
-    float threshold                     ///< Threshold above which to find a peak
-);
-
-/******************************************************************************
-pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage
-above the given threshold.  Returns a psList containing the location (x/y
-value) of all peaks.
- *****************************************************************************/
+    const psVector *vector,  ///< The input vector (float)
+    float threshold   ///< Threshold above which to find a peak
+);
+
+
+/** 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
+ * peak and exclude peaks which are not local peaks. A peak is a local peak if it
+ * has a higher value than all 8 neighbors. If the peak has the same value as its
+ * +y neighbor or +x neighbor, it is NOT a local peak. If any other neighbors
+ * have an equal value, the peak is considered a valid peak. Note two points:
+ * first, the +x neighbor condition is already enforced by pmFindVectorPeaks.
+ * Second, these rules have the effect of making flat-topped regions have single
+ * 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(
-    const psImage *image,               ///< The input image where peaks will be found (float)
-    float threshold                     ///< Threshold above which to find a peak
-);
-
-/******************************************************************************
-psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psList that have
-a peak value above the given maximum, or fall outside the valid region.
- *****************************************************************************/
+    const psImage *image,  ///< The input image where peaks will be found (float)
+    float threshold   ///< Threshold above which to find a peak
+);
+
+
+/** pmCullPeaks()
+ * 
+ * Eliminate peaks from the psList that have a peak value above the given
+ * maximum, or fall outside the valid region.
+ * 
+ */
 psList *pmCullPeaks(
-    psList *peaks,                      ///< The psList of peaks to be culled
-    float maxValue,                     ///< Cull peaks above this value
-    const psRegion valid               ///< Cull peaks otside this psRegion
-);
-
-/******************************************************************************
-pmSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius):
- 
- *****************************************************************************/
-pmSource *pmSourceLocalSky(
-    const psImage *image,               ///< The input image (float)
-    const pmPeak *peak,                 ///< The peak for which the psSource struct is created.
-    psStatsOptions statsOptions,        ///< The statistic used in calculating the background sky
-    float innerRadius,                  ///< The inner radius of the suqare annulus for calculating sky
-    float outerRadius                   ///< The outer radius of the suqare annulus for calculating sky
-);
-
-/******************************************************************************
- *****************************************************************************/
-pmSource *pmSourceMoments(
-    pmSource *source,                   ///< The input pmSource for which moments will be computed
-    float radius                        ///< Use a circle of pixels around the peak
-);
-
-/******************************************************************************
-pmSourceRoughClass(pmArray *source, psMetaDeta *metadata): make a guess at the
-source classification.
- *****************************************************************************/
+    psList *peaks,   ///< The psList of peaks to be culled
+    float maxValue,   ///< Cull peaks above this value
+    const psRegion valid                ///< Cull peaks otside this psRegion
+);
+
+
+/** 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
+ * have a peak value above the given maximum value limit or if the fall outside
+ * the valid region.  The result of the function is a new array with a reduced
+ * number of peaks.
+ * 
+ */
+psArray *pmPeaksSubset(
+    psArray *peaks,                     ///< Add comment.
+    float maxvalue,                     ///< Add comment.
+    const psRegion valid                ///< Add comment.
+);
+
+
+/** 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
+ * a square window (of full width 2Radius+1) centered on the pixel which contains
+ * the given coordinate, in the frame of the readout. The window is defined to
+ * have limits which are valid within the boundary of the readout image, thus if
+ * the radius would fall outside the image pixels, the subimage is truncated to
+ * only consist of valid pixels. If readout->mask or readout->weight are not
+ * NULL, matching subimages are defined for those images as well. This function
+ * fails if no valid pixels can be defined (x or y less than Radius, for
+ * 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.
+bool pmSourceDefinePixels(
+    pmSource *mySource,                 ///< Add comment.
+    pmReadout *readout,                 ///< Add comment.
+    psF32 x,                            ///< Add comment.
+    psF32 y,                            ///< Add comment.
+    psF32 Radius                        ///< Add comment.
+);
+
+
+/** 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
+ * function assumes the source pixels have been defined, and that the value of
+ * Radius here is smaller than the value of Radius used to define the pixels. The
+ * annular region not contained within the radius defined here is used to measure
+ * the local background in the vicinity of the source. The local background
+ * measurement uses the specified statistic passed in via the statsOptions entry.
+ * 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(
+    pmSource *source,   ///< The input image (float)
+    psStatsOptions statsOptions, ///< The statistic used in calculating the background sky
+    float Radius   ///< The inner radius of the square annulus to exclude
+);
+
+
+/** pmSourceMoments()
+ * 
+ * Measure source moments for the given source, using the value of
+ * source.moments.sky provided as the local background value and the peak
+ * coordinates as the initial source location. The resulting moment values are
+ * applied to the source.moments entry, and the source is returned. The moments
+ * are measured within the given circular radius of the source.peak coordinates.
+ * The return value indicates the success (TRUE) of the operation.
+ * 
+ */
+bool pmSourceMoments(
+    pmSource *source,   ///< The input pmSource for which moments will be computed
+    float radius   ///< Use a circle of pixels around the peak
+);
+
+
+/** 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
+ * the image. As long as the PSF shape does not vary excessively across the
+ * image, the sources which are represented by a PSF (the start) will have very
+ * similar second moments. The function pmSourcePSFClump searches a collection of
+ * sources with measured moments for a group with moments which are all very
+ * 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
+ * represent this plane.  The units of sigma_x and sigma_y are in image pixels. A
+ * pixel in this analysis image represents 0.1 pixels in the input image. The
+ * dimensions of the image need only be 10 pixels. The peak pixel in this image
+ * (above a threshold of half of the image maximum) is found. The coordinates of
+ * this peak pixel represent the 2D mode of the sigma_x, sigma_y distribution.
+ * The sources with sigma_x, sigma_y within 0.2 pixels of this value are then
+ *  * 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(
+    psArray *source,   ///< The input pmSource
+    psMetadata *metadata  ///< Contains classification parameters
+);
+
+
+/** pmSourceRoughClass()
+ * 
+ * Based on the specified data values, make a guess at the source
+ * classification. The sources are provides as a psArray of pmSource entries.
+ * Definable parameters needed to make the classification are provided to the
+ * routine with the psMetadata structure. The rules (in SDRS) refer to values which
+ * can be extracted from the metadata using the given keywords. Except as noted,
+ * the data type for these parameters are psF32.
+ * 
+ */
 bool pmSourceRoughClass(
-    psArray *source,                    ///< The input pmSource
-    psMetadata *metadata                ///< Contains classification parameters
-);
-/******************************************************************************
-pmSourceSetPixelCircle(source, image, radius)
- *****************************************************************************/
-bool pmSourceSetPixelsCircle(
-    pmSource *source,                   ///< The input pmSource
-    const psImage *image,               ///< The input image (float)
-    float radius                        ///< The radius of the circle
-);
-
-/******************************************************************************
- *****************************************************************************/
-bool pmSourceModelGuess(
-    pmSource *source,                   ///< The input pmSource
-    const psImage *image,               ///< The input image (float)
-    pmModelType model                   ///< The type of model to be created.
-);
-
-/******************************************************************************
- *****************************************************************************/
+    psArray *source,   ///< The input pmSource
+    psMetadata *metadata,  ///< Contains classification parameters
+    pmPSFClump clump   ///< Statistics about the PSF clump
+);
+
+
+/** 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
+ * provided model selection. The method of defining the model parameter guesses
+ * 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(
+    pmSource *source,   ///< The input pmSource
+    pmModelType model   ///< The type of model to be created.
+);
+
+
+/** pmContourType
+ * 
+ * Only one type is defined at present.
+ * 
+ */
 typedef enum {
     PS_CONTOUR_CRUDE,
+    PS_CONTOUR_UNKNOWN01,
+    PS_CONTOUR_UNKNOWN02
 } pmContourType;
 
+
+/** 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,
+ * running along each pixel row until the level is crossed, then interpolating to
+ * the level coordinate for that row. This is done for each row, with the
+ * starting point determined by the midpoint of the previous row, until the
+ * starting point has a value below the contour level. The returned contour
+ * consists of two vectors giving the x and y coordinates of the contour levels.
+ * 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(
-    pmSource *source,                   ///< The input pmSource
-    const psImage *image,               ///< The input image (float) (this arg should be removed)
-    float level,                        ///< The level of the contour
-    pmContourType mode                  ///< Currently this must be PS_CONTOUR_CRUDE
-);
-
-/******************************************************************************
- *****************************************************************************/
+    pmSource *source,   ///< The input pmSource
+    const psImage *image,  ///< The input image (float) (this arg should be removed)
+    float level,   ///< The level of the contour
+    pmContourType mode   ///< Currently this must be PS_CONTOUR_CRUDE
+);
+
+
+/** 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
+ * interest are specified by the source.pixelsand source.maskentries. This
+ * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
+ * on success or FALSE on failure.
+ * 
+ */
 bool pmSourceFitModel(
-    pmSource *source,                   ///< The input pmSource
-    const psImage *image                ///< The input image (float)
-);
-
-/******************************************************************************
- *****************************************************************************/
+    pmSource *source,   ///< The input pmSource
+    pmModel *model,   ///< model to be fitted
+    const bool PSF   ///< Treat model as PSF or FLT?
+);
+
+
+/** 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(
+    pmModel *model                      ///< Add comment.
+);
+
+
+/** 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
+ * it is placed at its centroid location. The boolean option sky selects if the
+ * background sky is applied (TRUE) or not. The pixel range in the target image
+ * is at most the pixel range specified by the source.pixels image. The success
+ * status is returned.
+ * 
+ */
 bool pmSourceAddModel(
-    psImage *image,                     ///< The opuut image (float)
-    pmSource *source,                   ///< The input pmSource
-    bool center                         ///< A boolean flag that determines whether pixels are centered
-);
-
-/******************************************************************************
- *****************************************************************************/
+    psImage *image,   ///< The output image (float)
+    psImage *mask,   ///< The image pixel mask (valid == 0)
+    pmModel *model,   ///< The input pmModel
+    bool center    ///< A boolean flag that determines whether pixels are centered
+);
+
+
+/** 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
+ * or if it is placed at its centroid location. The boolean option sky selects if
+ * the background sky is applied (TRUE) or not. The pixel range in the target
+ * image is at most the pixel range specified by the source.pixels image. The
+ * success status is returned.
+ * 
+ */
 bool pmSourceSubModel(
-    psImage *image,                     ///< The output image (float)
-    pmSource *source,                   ///< The input pmSource
-    bool center                         ///< A boolean flag that determines whether pixels are centered
-);
-
-/******************************************************************************
-XXX: Why only *x argument?
-XXX EAM: psMinimizeLMChi2Func returns psF64, not float
- *****************************************************************************/
-float pmMinLM_Gauss2D(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
-);
-
-/******************************************************************************
- *****************************************************************************/
-float pmMinLM_PsuedoGauss2D(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
-);
-
-/******************************************************************************
- *****************************************************************************/
-float pmMinLM_Wauss2D(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
-);
-
-/******************************************************************************
- *****************************************************************************/
-float pmMinLM_TwistGauss2D(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
-);
-
-/******************************************************************************
- *****************************************************************************/
-float pmMinLM_Sersic(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
-);
-
-/******************************************************************************
- *****************************************************************************/
-float pmMinLM_SersicCore(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
-);
-
-/******************************************************************************
- *****************************************************************************/
-float pmMinLM_PsuedoSersic(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
+    psImage *image,   ///< The output image (float)
+    psImage *mask,   ///< The image pixel mask (valid == 0)
+    pmModel *model,   ///< The input pmModel
+    bool center    ///< A boolean flag that determines whether pixels are centered
 );
 
@@ -330 +531 @@
 /**
  * 
- *  The object model functions 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 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
+ * to infinity. In practice, we integrate the model beyond 50sigma.  The aperture magnitude is
+ * defined as -2.5log(flux) , where the flux is summed for all pixels which are
+ * 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(
+    float *fitMag,                      ///< integrated fit magnitude
+    float *obsMag,   ///< aperture flux magnitude
+    pmModel *model,                     ///< model used for photometry
+    psImage *image,                     ///< image pixels to be used
+    psImage *mask                       ///< mask of pixels to ignore
+);
+
 
 /**
  * 
- *  This function is the model chi-square minimization function for this model.
- * 
- */
-typedef psMinimizeLMChi2Func pmModelFunc;
-
-
-/**
- * 
- * This function returns the integrated flux for the given model parameters.
- */
-typedef psF64 (*pmModelFlux)(const psVector *params);
-
-
-/**
- * 
- *  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 FLT model for the object.
- * 
- */
-typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf);
-
-
-/**
- * 
- *  This function returns the radius at which the given model and parameters
- *  achieves the given flux.
- * 
- */
-typedef psF64 (*pmModelRadius)(const psVector *params, double flux);
-
-
-/**
- * 
- *  Each of the function types above has a corresponding function which returns
- *  the function given the model type:
- * 
- */
-pmModelFunc pmModelFunc_GetFunction (pmModelType type);
-pmModelFlux pmModelFlux_GetFunction (pmModelType type);
-pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type);
-pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type);
-pmModelRadius pmModelRadius_GetFunction (pmModelType type);
-psS32 pmModelParameterCount(pmModelType type);
-psS32 pmModelSetType(char *name);
-char *pmModelGetType(pmModelType type);
+ * This function converts the source classification into the closest available
+ * approximation to the Dophot classification scheme:
+ * 
+ * PM_SOURCE_DEFECT: 8
+ * PM_SOURCE_SATURATED: 8
+ * PM_SOURCE_SATSTAR: 10
+ * PM_SOURCE_PSFSTAR: 1
+ * PM_SOURCE_GOODSTAR: 1
+ * PM_SOURCE_POOR_FIT_PSF: 7
+ * PM_SOURCE_FAIL_FIT_PSF: 4
+ * PM_SOURCE_FAINTSTAR: 4
+ * PM_SOURCE_GALAXY: 2
+ * PM_SOURCE_FAINT_GALAXY: 2
+ * PM_SOURCE_DROP_GALAXY: 2
+ * PM_SOURCE_FAIL_FIT_GAL: 2
+ * PM_SOURCE_POOR_FIT_GAL: 2
+ * PM_SOURCE_OTHER: ?
+ * 
+ */
+int pmSourceDophotType(
+    pmSource *source                    ///< Add comment.
+);
+
+
+/** 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(
+    pmSource *source                    ///< Add comment.
+);
+
+/** 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
+ * interest are specified by the source.pixelsand source.maskentries. This
+ * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
+ * on success or FALSE on failure.
+ * 
+ */
+bool pmSourceFitModel_v5(
+    pmSource *source,   ///< The input pmSource
+    pmModel *model,   ///< model to be fitted
+    const bool PSF   ///< Treat model as PSF or FLT?
+);
+
+
+/** 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
+ * interest are specified by the source.pixelsand source.maskentries. This
+ * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
+ * on success or FALSE on failure.
+ * 
+ */
+bool pmSourceFitModel_v7(
+    pmSource *source,   ///< The input pmSource
+    pmModel *model,   ///< model to be fitted
+    const bool PSF   ///< Treat model as PSF or FLT?
+);
 
 #endif