/** @file  pmObjects.c
 *
 *  This file will ...
 *
 *  @author GLG, MHPCC
 *  @author EAM, IfA: significant modifications.
 *
 *  @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-12-07 20:01:05 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

#include <stdio.h>
#include <math.h>
#include <string.h>
#include "pslib.h"
#include "pmObjects.h"
#include "pmModelGroup.h"
/******************************************************************************
pmPeakAlloc(): Allocate the pmPeak data structure and set appropriate members.
*****************************************************************************/
pmPeak *pmPeakAlloc(psS32 x,
                    psS32 y,
                    psF32 counts,
                    pmPeakType class)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    pmPeak *tmp = (pmPeak *) psAlloc(sizeof(pmPeak));
    tmp->x = x;
    tmp->y = y;
    tmp->counts = counts;
    tmp->class = class;

    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    return(tmp);
}

/******************************************************************************
pmMomentsAlloc(): Allocate the pmMoments structure and initialize the members
to zero.
*****************************************************************************/
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->Sy = 0.0;
    tmp->Sxy = 0.0;
    tmp->Sum = 0.0;
    tmp->Peak = 0.0;
    tmp->Sky = 0.0;
    tmp->nPixels = 0;

    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    return(tmp);
}

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__);
}

/******************************************************************************
pmModelAlloc(): Allocate the pmModel structure, along with its parameters,
and initialize the type member.  Initialize the params to 0.0.
XXX EAM: simplifying code with pmModelParameterCount
*****************************************************************************/
pmModel *pmModelAlloc(pmModelType type)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel));

    tmp->type = type;
    tmp->chisq = 0.0;
    tmp->nIter = 0;
    tmp->radius = 0;
    tmp->status = PM_MODEL_UNTRIED;

    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++) {
        tmp->params->data.F32[i] = 0.0;
        tmp->dparams->data.F32[i] = 0.0;
    }

    psMemSetDeallocator(tmp, (psFreeFunc) modelFree);
    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    return(tmp);
}

/******************************************************************************
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->weight);
    psFree(tmp->mask);
    psFree(tmp->moments);
    psFree(tmp->modelPSF);
    psFree(tmp->modelFLT);
    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
}

/******************************************************************************
pmSourceAlloc(): Allocate the pmSource structure and initialize its members
to NULL.
*****************************************************************************/
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;
    tmp->modelPSF = NULL;
    tmp->modelFLT = NULL;
    tmp->type = 0;
    psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);

    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    return(tmp);
}

/******************************************************************************
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.
 
XXX: What types should be supported?  Only F32 is implemented.
 
XXX: We currently step through the input vector twice; once to determine the
size of the output vector, then to set the values of the output vector.
Depending upon actual use, this may need to be optimized.
*****************************************************************************/
psVector *pmFindVectorPeaks(const psVector *vector,
                            psF32 threshold)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    PS_ASSERT_VECTOR_NON_NULL(vector, NULL);
    PS_ASSERT_VECTOR_NON_EMPTY(vector, NULL);
    PS_ASSERT_VECTOR_TYPE(vector, PS_TYPE_F32, NULL);
    int count = 0;
    int n = vector->n;

    //
    // Special case: the input vector has a single element.
    //
    if (n == 1) {
        psVector *tmpVector = NULL;
        ;
        if (vector->data.F32[0] > threshold) {
            tmpVector = psVectorAlloc(1, PS_TYPE_U32);
            tmpVector->data.U32[0] = 0;
        } else {
            tmpVector = psVectorAlloc(0, PS_TYPE_U32);
        }
        psTrace(__func__, 3, "---- %s() end ----\n", __func__);
        return(tmpVector);
    }

    //
    // Determine if first pixel is a peak
    //
    if ((vector->data.F32[0] > vector->data.F32[1]) &&
            (vector->data.F32[0] > threshold)) {
        count++;
    }

    //
    // Determine if interior pixels are peaks
    //
    for (psU32 i = 1; i < n-1 ; i++) {
        if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
                (vector->data.F32[i] > vector->data.F32[i+1]) &&
                (vector->data.F32[i] > threshold)) {
            count++;
        }
    }

    //
    // Determine if last pixel is a peak
    //
    if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
            (vector->data.F32[n-1] > threshold)) {
        count++;
    }

    //
    // We know how many peaks exist, so we now allocate a psVector to store
    // those peaks.
    //
    psVector *tmpVector = psVectorAlloc(count, PS_TYPE_U32);
    count = 0;

    //
    // Determine if first pixel is a peak
    //
    if ((vector->data.F32[0] > vector->data.F32[1]) &&
            (vector->data.F32[0] > threshold)) {
        tmpVector->data.U32[count++] = 0;
    }

    //
    // Determine if interior pixels are peaks
    //
    for (psU32 i = 1; i < (n-1) ; i++) {
        if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
                (vector->data.F32[i] > vector->data.F32[i+1]) &&
                (vector->data.F32[i] > threshold)) {
            tmpVector->data.U32[count++] = i;
        }
    }

    //
    // Determine if last pixel is a peak
    //
    if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
            (vector->data.F32[n-1] > threshold)) {
        tmpVector->data.U32[count++] = n-1;
    }

    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    return(tmpVector);
}

/******************************************************************************
getRowVectorFromImage(): a private function which simply returns a
psVector containing the specified row of data from the psImage.
 
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);

    psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32);
    for (psU32 col = 0; col < image->numCols ; col++) {
        tmpVector->data.F32[col] = image->data.F32[row][col];
    }
    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    return(tmpVector);
}

/******************************************************************************
myListAddPeak(): A private function which allocates a psArray, if the list
argument is NULL, otherwise it adds the peak to that list.
XXX EAM : changed the output to psArray
XXX EAM : Switched row, col args
XXX EAM : NOTE: this was changed in the call, so the new code is consistent
*****************************************************************************/
static psArray *myListAddPeak(psArray *list,
                              psS32 row,
                              psS32 col,
                              psF32 counts,
                              pmPeakType type)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type);

    if (list == NULL) {
        list = psArrayAlloc(100);
        list->n = 0;
    }
    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);
}

/******************************************************************************
pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage
above the given threshold.  Returns a psArray containing location (x/y value)
of all peaks.
 
XXX: I'm not convinced the peak type definition in the SDRS is mutually
exclusive.  Some peaks can have multiple types.  Edges for sure.  Also, a
digonal line with the same value at each point will have a peak for every
point on that line.
 
XXX: This does not work if image has either a single row, or a single column.
 
XXX: In the output psArray elements, should we use the image row/column offsets?
     Currently, we do not.
 
XXX: Merge with CVS 1.20.  This had the proper code for images with a single
row or column.
*****************************************************************************/
psArray *pmFindImagePeaks(const psImage *image,
                          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);
    }
    psVector *tmpRow = NULL;
    psU32 col = 0;
    psU32 row = 0;
    psArray *list = NULL;

    //
    // Find peaks in row 0 only.
    //
    row = 0;
    tmpRow = getRowVectorFromImage((psImage *) image, row);
    psVector *row1 = pmFindVectorPeaks(tmpRow, threshold);

    for (psU32 i = 0 ; i < row1->n ; i++ ) {
        col = row1->data.U32[i];
        //
        // Determine if pixel (0,0) is a peak.
        //
        if (col == 0) {
            if ( (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
                    (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {

                if (image->data.F32[row][col] > threshold) {
                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
                }
            }
        } else if (col < (image->numCols - 1)) {
            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
                    (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
                    (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
                if (image->data.F32[row][col] > threshold) {
                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
                }
            }

        } else if (col == (image->numCols - 1)) {
            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
                if (image->data.F32[row][col] > threshold) {
                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
                }
            }

        } else {
            psError(PS_ERR_UNKNOWN, true, "peak specified valid column range.");
        }
    }
    psFree (tmpRow);
    psFree (row1);

    //
    // Exit if this image has a single row.
    //
    if (image->numRows == 1) {
        psTrace(__func__, 3, "---- %s() end ----\n", __func__);
        return(list);
    }

    //
    // Find peaks in interior rows only.
    //
    for (row = 1 ; row < (image->numRows - 1) ; row++) {
        tmpRow = getRowVectorFromImage((psImage *) image, row);
        row1 = pmFindVectorPeaks(tmpRow, threshold);

        // Step through all local peaks in this row.
        for (psU32 i = 0 ; i < row1->n ; i++ ) {
            pmPeakType myType = PM_PEAK_UNDEF;
            col = row1->data.U32[i];

            if (col == 0) {
                // If col==0, then we can not read col-1 pixels
                if ((image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
                        (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
                        (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
                        (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
                        (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
                    myType = PM_PEAK_EDGE;
                    list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
                }
            } else if (col < (image->numCols - 1)) {
                // This is an interior pixel
                if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
                        (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
                        (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
                        (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
                        (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
                        (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
                        (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
                        (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
                    if (image->data.F32[row][col] > threshold) {
                        if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) &&
                                (image->data.F32[row][col] > image->data.F32[row-1][col]) &&
                                (image->data.F32[row][col] > image->data.F32[row-1][col+1]) &&
                                (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
                                (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
                                (image->data.F32[row][col] > image->data.F32[row+1][col-1]) &&
                                (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
                                (image->data.F32[row][col] > image->data.F32[row+1][col+1])) {
                            myType = PM_PEAK_LONE;
                        }

                        if ((image->data.F32[row][col] == image->data.F32[row-1][col-1]) ||
                                (image->data.F32[row][col] == image->data.F32[row-1][col]) ||
                                (image->data.F32[row][col] == image->data.F32[row-1][col+1]) ||
                                (image->data.F32[row][col] == image->data.F32[row][col-1]) ||
                                (image->data.F32[row][col] == image->data.F32[row][col+1]) ||
                                (image->data.F32[row][col] == image->data.F32[row+1][col-1]) ||
                                (image->data.F32[row][col] == image->data.F32[row+1][col]) ||
                                (image->data.F32[row][col] == image->data.F32[row+1][col+1])) {
                            myType = PM_PEAK_FLAT;
                        }

                        list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
                    }
                }
            } else if (col == (image->numCols - 1)) {
                // If col==numCols - 1, then we can not read col+1 pixels
                if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
                        (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
                        (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
                        (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
                        (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
                        (image->data.F32[row][col] >= image->data.F32[row+1][col])) {
                    myType = PM_PEAK_EDGE;
                    list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
                }
            } else {
                psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
            }

        }
        psFree (tmpRow);
        psFree (row1);
    }

    //
    // Find peaks in the last row only.
    //
    row = image->numRows - 1;
    tmpRow = getRowVectorFromImage((psImage *) image, row);
    row1 = pmFindVectorPeaks(tmpRow, threshold);
    for (psU32 i = 0 ; i < row1->n ; i++ ) {
        col = row1->data.U32[i];
        if (col == 0) {
            if ( (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
                    (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
                    (image->data.F32[row][col] >  image->data.F32[row][col+1])) {
                if (image->data.F32[row][col] > threshold) {
                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
                }
            }
        } else if (col < (image->numCols - 1)) {
            if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
                    (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
                    (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
                    (image->data.F32[row][col] >  image->data.F32[row][col-1]) &&
                    (image->data.F32[row][col] >= image->data.F32[row][col+1])) {
                if (image->data.F32[row][col] > threshold) {
                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
                }
            }

        } else if (col == (image->numCols - 1)) {
            if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
                    (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
                    (image->data.F32[row][col] >  image->data.F32[row][col-1])) {
                if (image->data.F32[row][col] > threshold) {
                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
                }
            }
        } else {
            psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
        }
    }
    psFree (tmpRow);
    psFree (row1);
    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    return(list);
}

// XXX: Macro this.
static bool isItInThisRegion(const psRegion valid,
                             psS32 x,
                             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);
}


/******************************************************************************
psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have
a peak value above the given maximum, or fall outside the valid region.
 
XXX: Should the sky value be used when comparing the maximum?
 
XXX: warning message if valid is NULL?
 
XXX: changed API to create a NEW output psArray (should change name as well)
 
XXX: Do we free the psList elements of those culled peaks?
 
XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset?
*****************************************************************************/
psList *pmCullPeaks(psList *peaks,
                    psF32 maxValue,
                    const psRegion valid)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    PS_ASSERT_PTR_NON_NULL(peaks, NULL);

    psListElem *tmpListElem = (psListElem *) peaks->head;
    psS32 indexNum = 0;

    //    printf("pmCullPeaks(): list size is %d\n", peaks->size);
    while (tmpListElem != NULL) {
        pmPeak *tmpPeak = (pmPeak *) tmpListElem->data;
        if ((tmpPeak->counts > maxValue) ||
                (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))) {
            psListRemoveData(peaks, (psPtr) tmpPeak);
        }

        indexNum++;
        tmpListElem = tmpListElem->next;
    }

    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    return(peaks);
}

// XXX EAM: I changed this to return a new, subset array
//          rather than alter the existing one
// XXX: Fix the *valid pointer.
psArray *pmPeaksSubset(
    psArray *peaks,
    psF32 maxValue,
    const psRegion valid)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    PS_ASSERT_PTR_NON_NULL(peaks, NULL);

    psArray *output = psArrayAlloc (200);
    output->n = 0;

    psTrace (".pmObjects.pmCullPeaks", 3, "list size is %d\n", peaks->n);

    for (int i = 0; i < peaks->n; i++) {
        pmPeak *tmpPeak = (pmPeak *) peaks->data[i];
        if (tmpPeak->counts > maxValue)
            continue;
        if (isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))
            continue;
        psArrayAdd (output, 200, tmpPeak);
    }
    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    return(output);
}

/******************************************************************************
pmSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): this
routine creates a new pmSource data structure and sets the following members:
    ->pmPeak
    ->pmMoments->sky
 
The sky value is set from the pixels in the square annulus surrounding the
peak pixel.
 
We simply create a subSet image and mask the inner pixels, then call
psImageStats on that subImage+mask.
 
XXX: The subImage has width of 1+2*outerRadius.  Verify with IfA.
 
XXX: Use static data structures for:
     subImage
     subImageMask
     myStats
 
XXX: ensure that the inner and out radius fit in the actual image.  Should
     we generate an error, or warning?  Currently an error.
 
XXX: Sync with IfA on whether the peak x/y coords are data structure coords,
     or they use the image row/column offsets.
 
XXX: Should we simply set pmSource->peak = peak?  If so, should we increase
the reference counter?  Or, should we copy the data structure?
 
XXX: Currently the subimage always has an even number of rows/columns.  Is
     this correct?  Since there is a center pixel, maybe it should have an
     odd number of rows/columns.
 
XXX: Use psTrace() for the print statements.
 
XXX: Don't use separate structs for the subimage and mask.  Use the source->
     members.
*****************************************************************************/









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, image, mask, 0xff);
    psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED);

    psF64 tmpF64;
    p_psGetStatValue(myStats, &tmpF64);
    psFree(myStats);

    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);
}

/******************************************************************************
bool checkRadius2(): private function which simply determines if the (x, y)
point is within the radius of the specified peak.
 
XXX: macro this for performance.
XXX: this is rather inefficient - at least compute and compare against radius^2
*****************************************************************************/
static bool checkRadius2(psF32 xCenter,
                         psF32 yCenter,
                         psF32 radius,
                         psF32 x,
                         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)) {
        return(true);
    }

    psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
    return(false);
}

/******************************************************************************
pmSourceMoments(source, radius): this function takes a subImage defined in the
pmSource data structure, along with the peak location, and determines the
various moments associated with that peak.
 
Requires the following to have been created:
    pmSource
    pmSource->peak
    pmSource->pixels
    pmSource->weight
    pmSource->mask
 
XXX: The peak calculations are done in image coords, not subImage coords.
 
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);

    //
    // XXX: Verify the setting for sky if source->moments == NULL.
    //
    psF32 sky = 0.0;
    if (source->moments == NULL) {
        source->moments = pmMomentsAlloc();
    } else {
        sky = source->moments->Sky;
    }

    //
    // Sum = SUM (z - sky)
    // X1  = SUM (x - xc)*(z - sky)
    // X2  = SUM (x - xc)^2 * (z - sky)
    // XY  = SUM (x - xc)*(y - yc)*(z - sky)
    //
    psF32 peakPixel = -PS_MAX_F32;
    psS32 numPixels = 0;
    psF32 Sum = 0.0;
    psF32 X1 = 0.0;
    psF32 Y1 = 0.0;
    psF32 X2 = 0.0;
    psF32 Y2 = 0.0;
    psF32 XY = 0.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
    // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2?
    // We loop through all pixels in this subimage (source->pixels), and for each
    // pixel that is not masked, AND within the radius of the peak pixel, we
    // 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])) {
                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);

    //
    // first moment X  = X1/Sum + xc
    // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
    // Sxy             = XY / Sum
    //
    x = X1/Sum;
    y = Y1/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;
    source->moments->nPixels = numPixels;

    // XXX EAM : these values can be negative, so we need to limit the range
    source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
    source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));

    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);
}

// XXX EAM : I used
int pmComparePeakAscend (const void **a, const void **b)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    pmPeak *A = *(pmPeak **)a;
    pmPeak *B = *(pmPeak **)b;

    psF32 diff;

    diff = A->counts - B->counts;
    if (diff < FLT_EPSILON) {
        psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);
        return (-1);
    } 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);
}

int pmComparePeakDescend (const void **a, const void **b)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    pmPeak *A = *(pmPeak **)a;
    pmPeak *B = *(pmPeak **)b;

    psF32 diff;

    diff = A->counts - B->counts;
    if (diff < FLT_EPSILON) {
        psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);
        return (+1);
    } 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);
}

/******************************************************************************
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;
    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
    {
        psStats *stats  = NULL;
        psImage *splane = NULL;
        int binX, binY;
        bool status;

        psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");
        if (!status)
            SX_MAX = 10.0;
        psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");
        if (!status)
            SY_MAX = 10.0;

        // construct a sigma-plane image
        // psImageAlloc does zero the data
        splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/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?)
        for (psS32 i = 0 ; i < sources->n ; i++)
        {
            pmSource *tmpSrc = (pmSource *) sources->data[i];
            if (tmpSrc == NULL) {
                continue;
            }
            if (tmpSrc->moments == NULL) {
                continue;
            }

            // Sx,Sy are limited at 0.  a peak at 0,0 is artificial
            if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) {
                continue;
            }

            // for the moment, force splane dimensions to be 10x10 image pix
            binX = tmpSrc->moments->Sx/SCALE;
            if (binX < 0)
                continue;
            if (binX >= splane->numCols)
                continue;

            binY = tmpSrc->moments->Sy/SCALE;
            if (binY < 0)
                continue;
            if (binY >= splane->numRows)
                continue;

            splane->data.F32[binY][binX] += 1.0;
        }

        // find the peak in this image
        stats = psStatsAlloc (PS_STAT_MAX);
        stats = psImageStats (stats, splane, NULL, 0);
        peaks = pmFindImagePeaks (splane, stats[0].max / 2);
        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
    {
        pmPeak *clump;
        psF32 minSx, maxSx;
        psF32 minSy, maxSy;
        psVector *tmpSx = NULL;
        psVector *tmpSy = NULL;
        psStats *stats  = NULL;

        // XXX EAM : this lets us takes the single highest peak
        psArraySort (peaks, pmComparePeakDescend);
        clump = peaks->data[0];
        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts);

        // define section window for clump
        minSx = clump->x * SCALE - 0.2;
        maxSx = clump->x * SCALE + 0.2;
        minSy = clump->y * SCALE - 0.2;
        maxSy = clump->y * SCALE + 0.2;

        tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32);
        tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32);
        tmpSx->n = 0;
        tmpSy->n = 0;

        // XXX clip sources based on flux?
        // create vectors with Sx, Sy values in window
        for (psS32 i = 0 ; i < sources->n ; i++)
        {
            pmSource *tmpSrc = (pmSource *) sources->data[i];

            if (tmpSrc->moments->Sx < minSx)
                continue;
            if (tmpSrc->moments->Sx > maxSx)
                continue;
            if (tmpSrc->moments->Sy < minSy)
                continue;
            if (tmpSrc->moments->Sy > maxSy)
                continue;
            tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx;
            tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy;
            tmpSx->n++;
            tmpSy->n++;
            if (tmpSx->n == tmpSx->nalloc) {
                psVectorRealloc (tmpSx, tmpSx->nalloc + 100);
                psVectorRealloc (tmpSy, tmpSy->nalloc + 100);
            }
        }

        // measures stats of Sx, Sy
        stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);

        stats = psVectorStats (stats, tmpSx, NULL, NULL, 0);
        psfClump.X  = stats->clippedMean;
        psfClump.dX = stats->clippedStdev;

        stats = psVectorStats (stats, tmpSy, NULL, NULL, 0);
        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

        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;
    psRegion allArray = psRegionSet (0, 0, 0, 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?
    // make rough IDs based on clumpX,Y,DX,DY
    for (psS32 i = 0 ; i < sources->n ; i++) {

        pmSource *tmpSrc = (pmSource *) sources->data[i];

        tmpSrc->peak->class = 0;

        psF32 sigX = tmpSrc->moments->Sx;
        psF32 sigY = tmpSrc->moments->Sy;

        // 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?
        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, 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;
        }

        // 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;
        }

        // 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;
        }

        // the rest are probable stellar objects
        starsn->data.F32[starsn->n] = SN;
        starsn->n ++;
        Nstar ++;

        // 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;
        }

        // random type of star
        tmpSrc->type = PM_SOURCE_OTHER;
    }

    {
        psStats *stats  = NULL;
        stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
        stats = psVectorStats (stats, starsn, NULL, NULL, 0);
        psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);
        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: This was replaced by DefinePixels in SDRS.  Remove it.
*****************************************************************************/
bool pmSourceSetPixelsCircle(pmSource *source,
                             const psImage *image,
                             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_FLOAT_COMPARE(0.0, radius, false);

    //
    // 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).
    //
    psS32 radiusS32 = (psS32) radius;
    psS32 SubImageCenterRow = source->peak->y;
    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, SubImageCenterRow + radiusS32 + 1);
    psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - radiusS32);
    psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1);

    // XXX: Must recycle image.
    // XXX EAM: this message reflects a programming error we know about.
    //          i am setting it to a trace message which we can take out
    if (source->pixels != NULL) {
        psTrace (".psModule.pmObjects.pmSourceSetPixelsCircle", 4,
                 "WARNING: pmSourceSetPixelsCircle(): image->pixels not NULL.  Freeing and reallocating.\n");
        psFree(source->pixels);
    }
    source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,
                                   SubImageStartRow,
                                   SubImageEndCol,
                                   SubImageEndRow));

    // XXX: Must recycle image.
    if (source->mask != NULL) {
        psFree(source->mask);
    }
    source->mask = psImageAlloc(source->pixels->numCols,
                                source->pixels->numRows,
                                PS_TYPE_U8); // XXX EAM : type was F32

    //
    // Loop through the subimage mask, initialize mask to 0 or 1.
    // XXX EAM: valid pixels should have 0, not 1
    for (psS32 row = 0 ; row < source->mask->numRows; row++) {
        for (psS32 col = 0 ; col < source->mask->numCols; col++) {

            if (checkRadius2((psF32) radiusS32,
                             (psF32) radiusS32,
                             radius,
                             (psF32) col,
                             (psF32) row)) {
                source->mask->data.U8[row][col] = 0;
            } else {
                source->mask->data.U8[row][col] = 1;
            }
        }
    }
    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
    return(true);
}

/******************************************************************************
pmSourceModelGuess(source, model): This function allocates a new
pmModel structure based on the given modelType specified in the argument list.  
The corresponding pmModelGuess function is returned, and used to 
supply the values of the params array in the pmModel structure.  
 
XXX: Many parameters are based on the src->moments structure, which is in
image, not subImage coords.  Therefore, the calls to the model evaluation
functions will be in image, not subImage coords.  Remember this.
*****************************************************************************/
pmModel *pmSourceModelGuess(pmSource *source,
                            pmModelType modelType)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    PS_ASSERT_PTR_NON_NULL(source->moments, false);
    PS_ASSERT_PTR_NON_NULL(source->peak, false);

    pmModel *model = pmModelAlloc(modelType);

    pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType);
    modelGuessFunc(model, source);
    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    return(model);
}

/******************************************************************************
evalModel(source, level, row): a private function which evaluates the
source->modelPSF function at the specified coords.  The coords are subImage, not
image coords.
 
NOTE: The coords are in subImage source->pixel coords, not image coords.
 
XXX: reverse order of row,col args?
 
XXX: rename all coords in this file such that their name defines whether
the coords is in subImage or image space.
 
XXX: This should probably be a public pmModules function.
 
XXX: Use static vectors for x.
 
XXX: Figure out if it's (row, col) or (col, row) for the model functions.
 
XXX: For a while, the first psVectorAlloc() was generating a seg fault during
testing.  Try to reproduce that and debug.
*****************************************************************************/

// 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 + image->col0);
    x->data.F32[1] = (psF32) (row + image->row0);
    psF32 tmpF;
    pmModelFunc modelFunc;

    modelFunc = pmModelFunc_GetFunction (model->type);
    tmpF = modelFunc (NULL, model->params, x);
    psFree(x);
    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    return(tmpF);
}

/******************************************************************************
findValue(source, level, row, col, dir): a private function which determines
the column coordinate of the model function which has the value "level".  If
dir equals 0, then you loop leftwards from the peak pixel, otherwise,
rightwards.
 
XXX: reverse order of row,col args?
 
XXX: Input row/col are in image coords.
 
XXX: The result is returned in image coords.
*****************************************************************************/
static psF32 findValue(pmSource *source,
                       psF32 level,
                       psU32 row,
                       psU32 col,
                       psU32 dir)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    //
    // Convert coords to subImage space.
    //
    psU32 subRow = row - source->pixels->row0;
    psU32 subCol = col - source->pixels->col0;

    // Ensure that the starting column is allowable.
    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);
    }

    // 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)));
    }

    //
    // We define variables incr and lastColumn so that we can use the same loop
    // whether we are stepping leftwards, or rightwards.
    //
    psS32 incr;
    psS32 lastColumn;
    if (dir == 0) {
        incr = -1;
        lastColumn = -1;
    } else {
        incr = 1;
        lastColumn = source->pixels->numCols;
    }
    subCol+=incr;

    while (subCol != lastColumn) {
        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));
        }

        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)) );
        }

        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)) );
        }

        subCol+=incr;
    }

    psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);
    return(NAN);
}

/******************************************************************************
pmSourceContour(src, img, level, mode): For an input subImage, and model, this
routine returns a psArray of coordinates that evaluate to the specified level.
 
XXX: Probably should remove the "image" argument.
XXX: What type should the output coordinate vectors consist of?  col,row?
XXX: Why a pmArray output?
XXX: doex x,y correspond with col,row or row/col?
XXX: What is mode?
XXX: The top, bottom of the contour is not correctly determined.
XXX EAM : this function is using the model for the contour, but it should
          be using only the image counts
*****************************************************************************/
psArray *pmSourceContour(pmSource *source,
                         const psImage *image,
                         psF32 level,
                         pmContourType mode)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    PS_ASSERT_PTR_NON_NULL(source, false);
    PS_ASSERT_PTR_NON_NULL(image, 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->modelFLT, false);
    // XXX EAM : what is the purpose of modelPSF/modelFLT?

    //
    // Allocate data for x/y pairs.
    //
    psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32);
    psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32);

    //
    // Start at the row with peak pixel, then decrement.
    //
    psS32 col = source->peak->x;
    for (psS32 row = source->peak->y; row>= 0 ; row--) {
        // XXX: yVec contain no real information.  Do we really need it?
        yVec->data.F32[row] = (psF32) (source->pixels->row0 + row);
        yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row);

        // Starting at peak pixel, search leftwards for the column intercept.
        psF32 leftIntercept = findValue(source, level, row, col, 0);
        if (isnan(leftIntercept)) {
            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
            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");
        }
        xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept;

        // Starting at peak pixel, search rightwards for the column intercept.

        psF32 rightIntercept = findValue(source, level, row, col, 1);
        if (isnan(rightIntercept)) {
            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
            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");
        }
        psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept);
        xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;

        // Set starting column for next row
        col = (psS32) ((leftIntercept + rightIntercept) / 2.0);
    }
    //
    // Start at the row (+1) with peak pixel, then increment.
    //
    col = source->peak->x;
    for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) {
        // XXX: yVec contain no real information.  Do we really need it?
        yVec->data.F32[row] = (psF32) (source->pixels->row0 + row);
        yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row);

        // Starting at peak pixel, search leftwards for the column intercept.
        psF32 leftIntercept = findValue(source, level, row, col, 0);
        if (isnan(leftIntercept)) {
            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
            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");
        }
        xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept;

        // Starting at peak pixel, search rightwards for the column intercept.
        psF32 rightIntercept = findValue(source, level, row, col, 1);
        if (isnan(rightIntercept)) {
            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
            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");
        }
        xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;

        // Set starting column for next row
        col = (psS32) ((leftIntercept + rightIntercept) / 2.0);
    }

    //
    // Allocate an array for result, store coord vectors there.
    //
    psArray *tmpArray = psArrayAlloc(2);
    tmpArray->data[0] = (psPtr *) yVec;
    tmpArray->data[1] = (psPtr *) xVec;
    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    return(tmpArray);
}

// XXX EAM : these are better starting values, but should be available from metadata?
#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
/******************************************************************************
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 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->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
    // 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__);
        model->status = PM_MODEL_BADARGS;
        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];
                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++;
            }
        }
    }

    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 : 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]);
    }

    // save the resulting chisq, nDOF, nIter
    model->chisq = myMin->value;
    model->nIter = myMin->iter;
    model->nDOF  = y->n - nParams;

    // 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);
    }

    // set the model success or failure status
    if (!fitStatus) {
        model->status = PM_MODEL_NONCONVERGE;
    } else {
        model->status = PM_MODEL_SUCCESS;
    }

    // models can go insane: reject these
    onPic &= (params->data.F32[2] >= source->pixels->col0);
    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
    onPic &= (params->data.F32[3] >= source->pixels->row0);
    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
    if (!onPic) {
        model->status = PM_MODEL_OFFIMAGE;
    }

    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__);
        model->status = PM_MODEL_BADARGS;
        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 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];
        }
    }

    // set the model success or failure status
    if (!fitStatus) {
        model->status = PM_MODEL_NONCONVERGE;
    } else {
        model->status = PM_MODEL_SUCCESS;
    }

    // XXX models can go insane: reject these
    onPic &= (params->data.F32[2] >= source->pixels->col0);
    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
    onPic &= (params->data.F32[3] >= source->pixels->row0);
    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
    if (!onPic) {
        model->status = PM_MODEL_OFFIMAGE;
    }

    psFree(x);
    psFree(y);
    psFree(yErr);
    psFree(myMin);
    psFree(covar);
    psFree(paramMask);

    rc = (onPic && fitStatus);
    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    return(rc);
}

bool p_pmSourceAddOrSubModel(psImage *image,
                             psImage *mask,
                             pmModel *model,
                             bool center,
                             bool sky,
                             bool add
                                )
{
    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 = model->params;
    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    psS32 imageCol;
    psS32 imageRow;
    psF32 skyValue = params->data.F32[0];
    psF32 pixelValue;

    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;

            // XXX: Should you be adding the pixels for the entire subImage,
            // or a radius of pixels around it?

            // Convert i/j to imace coord space:
            // XXX: Make sure you have col/row order correct.
            // 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;

            // set the appropriate pixel value for this coordinate
            if (sky) {
                pixelValue = modelFunc (NULL, params, x);
            } else {
                pixelValue = modelFunc (NULL, params, x) - skyValue;
            }


            // add or subtract the value
            if (add
               ) {
                image->data.F32[i][j] += pixelValue;
            }
            else {
                image->data.F32[i][j] -= pixelValue;
            }
        }
    }
    psFree(x);
    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
    return(true);
}



/******************************************************************************
 *****************************************************************************/
bool pmSourceAddModel(psImage *image,
                      psImage *mask,
                      pmModel *model,
                      bool center,
                      bool sky)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true);
    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    return(rc);
}

/******************************************************************************
 *****************************************************************************/
bool pmSourceSubModel(psImage *image,
                      psImage *mask,
                      pmModel *model,
                      bool center,
                      bool sky)
{
    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false);
    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    return(rc);
}


