Ticket #386: pmObjects.c

/** @file  pmObjects.c
 *
 *  This file will ...
 *
 *  @author GLG, MHPCC
 *
 *  @version $Revision: $ $Name: $
 *  @date $Date: 2005/04/01 20:47:40 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

#include<stdio.h>
#include<math.h>
#include "pslib.h"
#include "psConstants.h"
#include "pmObjects.h"

/******************************************************************************
pmPeakAlloc(): Allocate the psPeak data structure and set appropriate members.
*****************************************************************************/
psPeak *pmPeakAlloc(psS32 x,
                    psS32 y,
                    psF32 counts,
                    psPeakType class)
{
    psPeak *tmp = (psPeak *) psAlloc(sizeof(psPeak));
    tmp->x = x;
    tmp->y = y;
    tmp->counts = counts;
    tmp->class = class;

    return(tmp);
}

/******************************************************************************
pmMomentsAlloc(): Allocate the psMoments structure and initialize the members
to zero.
*****************************************************************************/
psMoments *pmMomentsAlloc()
{
    psMoments *tmp = (psMoments *) psAlloc(sizeof(psMoments));
    tmp->x = 0.0;
    tmp->y = 0.0;
    tmp->Sx = 0.0;
    tmp->Sx = 0.0;
    tmp->Sxy = 0.0;
    tmp->Sum = 0.0;
    tmp->Peak = 0.0;
    tmp->Sky = 0.0;
    tmp->nPixels = 0;

    return(tmp);
}

static void p_psModelFree(psModel *tmp)
{
    psFree(tmp->params);
    psFree(tmp->dparams);
}

/******************************************************************************
pmModelAlloc(): Allocate the psModel structure, along with its parameters,
and initialize the type member.  Initialize the params to 0.0.
XXX EAM: changing params and dparams to psVector 
*****************************************************************************/
psModel *pmModelAlloc(psModelType type)
{
    psModel *tmp = (psModel *) psAlloc(sizeof(psModel));

    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:
        psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
        return(NULL);
    }

    for (psS32 i = 0; i < tmp->params->n; i++) {
        tmp->params->data.F32[i] = 0.0;
        tmp->dparams->data.F32[i] = 0.0;
    }

    p_psMemSetDeallocator(tmp, (psFreeFcn) p_psModelFree);
    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.
*****************************************************************************/
static void p_psSourceFree(psSource *tmp)
{
    psFree(tmp->peak);
    //    psFree(tmp->pixels);
    //    psFree(tmp->mask);
    psFree(tmp->moments);
    psFree(tmp->models);
}

/******************************************************************************
pmSourceAlloc(): Allocate the psSource structure and initialize its members
to NULL.
*****************************************************************************/
psSource *pmSourceAlloc()
{
    psSource *tmp = (psSource *) psAlloc(sizeof(psSource));
    tmp->peak = NULL;
    tmp->pixels = NULL;
    tmp->mask = NULL;
    tmp->moments = NULL;
    tmp->models = NULL;
    tmp->type = 0;
    p_psMemSetDeallocator(tmp, (psFreeFcn) p_psSourceFree);

    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)
{
    PS_VECTOR_CHECK_NULL(vector, NULL);
    PS_VECTOR_CHECK_EMPTY(vector, NULL);
    PS_VECTOR_CHECK_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);
        }
        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;
    }

    return(tmpVector);
}

/******************************************************************************
p_psGetRowVectorFromImage(): 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?
*****************************************************************************/
psVector *p_psGetRowVectorFromImage(psImage *image,
                                    psU32 row)
{
    PS_IMAGE_CHECK_NULL(image, NULL);
    PS_IMAGE_CHECK_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];
    }
    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
*****************************************************************************/
psArray *MyListAddPeak(psArray *list,
                       psS32 row,
                       psS32 col,
                       psF32 counts,
                       psPeakType type)
{
    psPeak *tmpPeak = pmPeakAlloc(col, row, counts, type);

    if (list == NULL) {
        list = psArrayAlloc(100);
        list->n = 0;
    } 
    psArrayAdd(list, 100, tmpPeak);

    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.
*****************************************************************************/
psArray *pmFindImagePeaks(const psImage *image,
                          psF32 threshold)
{
    PS_IMAGE_CHECK_NULL(image, NULL);
    PS_IMAGE_CHECK_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.");
    }
    psVector *tmpRow = NULL;
    psU32 col = 0;
    psU32 row = 0;
    psArray *list = NULL;

    //
    // Find peaks in row 0 only.
    //
    row = 0;
    tmpRow = p_psGetRowVectorFromImage((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.");
        }
    }

    //
    // Exit if this image has a single row.
    //
    if (image->numRows == 1) {
        return(list);
    }

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

        // Step through all local peaks in this row.
        for (psU32 i = 0 ; i < row1->n ; i++ ) {
            psPeakType 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 valid column range.");
            }

        }
    }

    //
    // Find peaks in the last row only.
    //
    row = image->numRows - 1;
    tmpRow = p_psGetRowVectorFromImage((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 valid colum range.");
        }
    }
    return(list);
}

// XXX: Macro this.
bool IsItInThisRegion(const psRegion *valid,
                      psS32 x,
                      psS32 y)
{

    if ((x >= valid->x0) &&
        (x <= valid->x1) &&
        (y >= valid->y0) &&
        (y <= valid->y1)) {
        return(true);
    }

    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)
*****************************************************************************/
psList *pmCullPeaks(psList *peaks,
                    psF32 maxValue,
                    const psRegion *valid)
{
    PS_PTR_CHECK_NULL(peaks, NULL);
    //    PS_PTR_CHECK_NULL(valid, NULL);

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

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

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

    return(peaks);
}

// 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)
{
    PS_PTR_CHECK_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++) {
        psPeak *tmpPeak = (psPeak *) peaks->data[i];
        if (tmpPeak->counts > maxValue) continue;
        if (valid != NULL) {
            if (IsItInThisRegion(valid, tmpPeak->x, tmpPeak->y)) continue;
        }
        psArrayAdd (output, 200, tmpPeak);
    }
    return(output);
}

/******************************************************************************
psSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): this
routine creates a new psSource data structure and sets the following members:
    ->psPeak
    ->psMoments->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 psSource->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.
*****************************************************************************/
psSource *pmSourceLocalSky(const psImage *image,
                           const psPeak *peak,
                           psStatsOptions statsOptions,
                           psF32 innerRadius,
                           psF32 outerRadius)
{
    PS_IMAGE_CHECK_NULL(image, NULL);
    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL);
    PS_PTR_CHECK_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).
    //
    psImage *subImage = psImageSubset((psImage *) image,
                                      SubImageStartCol,
                                      SubImageStartRow,
                                      SubImageEndCol,
                                      SubImageEndRow);
    //    printf("pmSourceLocalSky: subimage width/length is (%d, %d)\n", subImage->numCols, subImage->numRows);
    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.
    //
    psStats *myStats = psStatsAlloc(statsOptions);
    myStats = psImageStats(myStats, subImage, subImageMask, 1);

    //
    // Create the output mySource, and set appropriate members.
    //
    psSource *mySource = pmSourceAlloc();
    mySource->peak = (psPeak *) peak;
    mySource->moments = pmMomentsAlloc();
    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.
*****************************************************************************/
bool CheckRadius(psPeak *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);
}

/******************************************************************************
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
*****************************************************************************/
bool CheckRadius2(psF32 xCenter,
                  psF32 yCenter,
                  psF32 radius,
                  psF32 x,
                  psF32 y)
{
    /// XXX EAM should compare with hypot (x,y) for speed
    if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) {
        return(true);
    }

    return(false);
}

/******************************************************************************
pmSourceMoments(source, radius): this function takes a subImage defined in the
psSource data structure, along with the peak location, and determines the
various moments associated with that peak.
 
Requires the following to have been created:
    psSource
    psSource->peak
    psSource->pixels
 
XXX: The peak calculations are done in image coords, not subImage coords.
 
XXX: mask values?
*****************************************************************************/
psSource *pmSourceMoments(psSource *source,
                          psF32 radius)
{
    PS_PTR_CHECK_NULL(source, NULL);
    PS_PTR_CHECK_NULL(source->peak, NULL);
    PS_PTR_CHECK_NULL(source->pixels, NULL);
    PS_FLOAT_COMPARE(0.0, radius, NULL);

    //
    // 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 Sum = 0.0;
    psF32 peakPixel = -PS_MAX_F32;
    psS32 numPixels = 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;
    //
    // 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.
    //
    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++;
                }
            }
        }
    }

    //
    // first moment X  = X1/Sum + xc
    // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
    // Sxy             = XY / Sum
    //
    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;
    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));
    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
}

// XXX EAM : I used 
int pmComparePeakAscend (const void **a, const void **b) 
{
    psPeak *A = *(psPeak **)a;
    psPeak *B = *(psPeak **)b;
    
    psF32 diff;

    diff = A->counts - B->counts;
    if (diff < FLT_EPSILON) return (-1);
    if (diff > FLT_EPSILON) return (+1);
    return (0);
}

int pmComparePeakDescend (const void **a, const void **b) 
{
    psPeak *A = *(psPeak **)a;
    psPeak *B = *(psPeak **)b;

    psF32 diff;

    diff = A->counts - B->counts;
    if (diff < FLT_EPSILON) return (+1);
    if (diff > FLT_EPSILON) return (-1);
    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_PTR_CHECK_NULL(sources, false);
    PS_PTR_CHECK_NULL(metadata, false);
    psBool rc = true;
    psArray *peaks  = NULL;
    psF32 clumpX = 0.0;
    psF32 clumpDX = 0.0;
    psF32 clumpY = 0.0;
    psF32 clumpDY = 0.0;
        
    // find the sigmaX, sigmaY clump
    { 
        psStats *stats  = NULL;
        psImage *splane = NULL;
        int binX, binY;

        // construct a sigma-plane image
        splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);

        // place the sources in the sigma-plane image (ignore 0,0 values?)
        for (psS32 i = 0 ; i < sources->n ; i++) {
            psSource *tmpSrc = (psSource *) sources->data[i];
            PS_PTR_CHECK_NULL(tmpSrc, false); // just skip this one?
            PS_PTR_CHECK_NULL(tmpSrc->moments, false); // just skip this one?

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

    }

    // measure statistics on Sx, Sy if Sx, Sy within range of clump
    { 
        psPeak *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\n", clump->x, clump->y);

        // 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++) {
            psSource *tmpSrc = (psSource *) 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);
        clumpX  = stats->clippedMean;
        clumpDX = 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);
        // 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;
    psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
    starsn->n = 0;

    // 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++) {

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

        tmpSrc->peak->class = 0;

        psF32 sigX = tmpSrc->moments->Sx;
        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;
            Nsat ++;
            continue;
        } 
        
        // too small to be stellar
        if ((sigX < (clumpX - clumpDX)) || (sigY < (clumpY - clumpDY))) {
            tmpSrc->type |= PS_SOURCE_DEFECT;
            Ncr ++;
            continue;
        }

        // possible galaxy
        if ((sigX > (clumpX + clumpDX)) || (sigY > (clumpY + clumpDY))) {
            tmpSrc->type |= PS_SOURCE_GALAXY;
            Ngal ++;
            continue;
        }

        // the rest are probable stellar objects
        psF32 S  = tmpSrc->moments->Sum;
        psF32 A  = PS_PI * sigX * sigY;
        psF32 B  = tmpSrc->moments->Sky;
        psF32 RT = PS_SQRT_F32(S + (A * B) + (A * PS_SQR(RDNOISE) / PS_SQRT_F32(GAIN)));
        psF32 SN = (S * PS_SQRT_F32(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;
            Npsf ++;
            continue;
        }

        // random type of star
        tmpSrc->type |= PS_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);
    }

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

    return(rc);
}

/******************************************************************************
pmSourceSetPixelCircle(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.
*****************************************************************************/
bool pmSourceSetPixelCircle(psSource *source,
                            const psImage *image,
                            psF32 radius)
{
    PS_IMAGE_CHECK_NULL(image, false);
    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
    PS_PTR_CHECK_NULL(source, false);
    PS_PTR_CHECK_NULL(source->moments, false);
    // PS_PTR_CHECK_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
    psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - radiusS32);
    psS32 SubImageEndRow    = PS_MIN (image->numRows - 1, SubImageCenterRow + radiusS32);
    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

    // 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.pmSourceSetPixelCircle", 4,
                 "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL.  Freeing and reallocating.\n");
        psFree(source->pixels);
    }
    source->pixels = psImageSubset((psImage *) image,
                                   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_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;
            }
        }
    }
    return(true);
}

/******************************************************************************
pmSourceModelGuess(source, image, model): This function allocates a new
psModel structure and stores it in the psSource data structure specified in
the argument list.  The model type is specified in the argument list.  The
params array in that psModel 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.
 
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.
*****************************************************************************/
bool pmSourceModelGuess(psSource *source,
                        const psImage *image,
                        psModelType model)
{
    PS_PTR_CHECK_NULL(source, false);
    PS_PTR_CHECK_NULL(source->moments, false);
    PS_PTR_CHECK_NULL(source->peak, false);
    PS_IMAGE_CHECK_NULL(image, false);
    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
    if (source->models != NULL) {
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->models was non-NULL; calling psFree(source->models).\n");
        psFree(source->models);
    }
    source->models = pmModelAlloc(model);

    psVector *params = source->models->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->models->params[7] = B2;
        // source->models->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);
# endif
      default:
        psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
        return(false);
    }
}

/******************************************************************************
evalModel(source, level, row): a private function which evaluates the
source->model 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 psModules 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.
*****************************************************************************/
static psF32 evalModel(psSource *src,
                psU32 row,
                psU32 col)
{
    PS_PTR_CHECK_NULL(src, false);
    PS_PTR_CHECK_NULL(src->models, false);
    PS_PTR_CHECK_NULL(src->models->params, false);

    // XXX: The following step will not be necessary if the models->params
    // member is a psVector.  Suggest to IfA.  

    // XXX EAM: done: models->params is now a vector
    psVector *params = src->models->params;

    //
    // 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);
    psF32 tmpF;

    switch (src->models->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 psModelType");
        return(NAN);
    }
    psFree(x);
    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(psSource *source,
                psF32 level,
                psU32 row,
                psU32 col,
                psU32 dir)
{
    //
    // 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");
        return(NAN);
    }
    if (!((0 <= subRow) && (subRow < source->pixels->numRows))) {
        psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range");
        return(NAN);
    }

    psF32 oldValue = evalModel(source, subRow, subCol);
    if (oldValue == level) {
        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 = evalModel(source, subRow, subCol);
        if (oldValue == level) {
            return((psF32) (subCol + source->pixels->col0));
        }

        if ((newValue <= level) && (level <= oldValue)) {
            // This is simple linear interpolation.
            return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) );
        }

        if ((oldValue <= level) && (level <= newValue)) {
            // This is simple linear interpolation.
            return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) );
        }

        subCol+=incr;
    }

    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 momde?
XXX: The top, bottom of the contour is not correctly determined.
*****************************************************************************/
psArray *pmSourceContour(psSource *source,
                         const psImage *image,
                         psF32 level,
                         pmContourType mode)
{
    PS_PTR_CHECK_NULL(source, false);
    PS_PTR_CHECK_NULL(image, false);
    PS_PTR_CHECK_NULL(source->moments, false);
    PS_PTR_CHECK_NULL(source->peak, false);
    PS_PTR_CHECK_NULL(source->pixels, false);
    PS_PTR_CHECK_NULL(source->models, false);

    //
    // 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);
            return(NULL);
        }
        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);
            return(NULL);
        }
        //printf("HERE 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);
            return(NULL);
        }
        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);
            return(NULL);
        }
        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;
    return(tmpArray);
}

psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x);
psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x);

// 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_TOLERANCE 0.1
/******************************************************************************
pmSourceFitModel(source, image): 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(psSource *source,
                      const psImage *image)
{
    PS_PTR_CHECK_NULL(source, false);
    PS_PTR_CHECK_NULL(source->moments, false);
    PS_PTR_CHECK_NULL(source->peak, false);
    PS_PTR_CHECK_NULL(source->pixels, false);
    PS_PTR_CHECK_NULL(source->models, false);
    PS_IMAGE_CHECK_NULL(image, false);
    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
    psBool rc;

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

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

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

    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
                                                PM_SOURCE_FIT_MODEL_TOLERANCE);

    psVector *params = source->models->params;

    switch (source->models->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 psModelType");
        rc = false;
    }
    // XXX EAM: we need to do something (give an error?) if rc is false
    // XXX EAM: save the resulting chisq, nDOF, nIter
    source->models->chisq = myMin->value;
    source->models->nDOF  = y->n - params->n;
    source->models->nIter = myMin->iter;

    psFree(x);
    psFree(y);
    psFree(myMin);
    return(rc);
}

bool p_pmSourceAddOrSubModel(psImage *image,
                             psSource *src,
                             bool center,
                             psS32 flag)
{
    PS_PTR_CHECK_NULL(src, false);
    PS_PTR_CHECK_NULL(src->moments, false);
    PS_PTR_CHECK_NULL(src->peak, false);
    PS_PTR_CHECK_NULL(src->pixels, false);
    PS_PTR_CHECK_NULL(src->models, false);
    PS_IMAGE_CHECK_NULL(image, false);
    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);

    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
    psVector *params = src->models->params;

    for (psS32 i = 0; i < src->pixels->numRows; i++) {
        for (psS32 j = 0; j < src->pixels->numCols; j++) {
            psF32 pixelValue;
            // 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.
            psS32 imageRow = i + src->pixels->row0;
            psS32 imageCol = j + src->pixels->col0;

            x->data.F32[0] = (float) imageCol;
            x->data.F32[1] = (float) imageRow;
            switch (src->models->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 psModelType");
                psFree(x);
                return(false);
            }
            if (flag == 1) {
                pixelValue = -pixelValue;
            }

            // XXX: Must figure out how to calculate the image coordinates and
            // how to use the boolean "center" flag.

            image->data.F32[imageRow][imageCol]+= pixelValue;
        }
    }
    psFree(x);
    return(true);
}



/******************************************************************************
 *****************************************************************************/
bool pmSourceAddModel(psImage *image,
                      psSource *src,
                      bool center)
{
    return(p_pmSourceAddOrSubModel(image, src, center, 0));
}

/******************************************************************************
 *****************************************************************************/
bool pmSourceSubModel(psImage *image,
                      psSource *src,
                      bool center)
{
    return(p_pmSourceAddOrSubModel(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;
*****************************************************************************/
psF64 pmMinLM_Gauss2D(psVector *deriv,
                      psVector *params,
                      psVector *x)
{
    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;
*****************************************************************************/
psF64 pmMinLM_PsuedoGauss2D(psVector *deriv,
                            psVector *params,
                            psVector *x)
{
    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;
*****************************************************************************/
psF64 pmMinLM_Wauss2D(psVector *deriv,
                      psVector *params,
                      psVector *x)
{
    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;
*****************************************************************************/
psF64 pmMinLM_TwistGauss2D(psVector *deriv,
                           psVector *params,
                           psVector *x)
{
    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;
*****************************************************************************/
psF64 pmMinLM_Sersic(psVector *deriv,
                     psVector *params,
                     psVector *x)
{
    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;
*****************************************************************************/
psF64 pmMinLM_SersicCore(psVector *deriv,
                         psVector *params,
                         psVector *x)
{
    psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
    return(0.0);
}