/** @file  pmObjects.c
 *
 *  This file will ...
 *
 *  @author GLG, MHPCC
 *
 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
 *  @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.
 *****************************************************************************/
psModel *pmModelAlloc(psModelType type)
{
    psModel *tmp = (psModel *) psAlloc(sizeof(psModel));

    tmp->type = type;
    tmp->chisq = 0.0;
    switch (type) {
    case PS_MODEL_GAUSS:
        tmp->Nparams = 7;
        tmp->params = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32));
        tmp->dparams = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32));
        break;
    case PS_MODEL_PGAUSS:
        tmp->Nparams = 7;
        tmp->params = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32));
        tmp->dparams = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32));
        break;
    case PS_MODEL_TWIST_GAUSS:
        tmp->Nparams = 11;
        tmp->params = (psF32 *) psAlloc(11 * sizeof(PS_TYPE_F32));
        tmp->dparams = (psF32 *) psAlloc(11 * sizeof(PS_TYPE_F32));
        break;
    case PS_MODEL_WAUSS:
        tmp->Nparams = 9;
        tmp->params = (psF32 *) psAlloc(9 * sizeof(PS_TYPE_F32));
        tmp->dparams = (psF32 *) psAlloc(9 * sizeof(PS_TYPE_F32));
        break;
    case PS_MODEL_SERSIC:
        tmp->Nparams = 8;
        tmp->params = (psF32 *) psAlloc(8 * sizeof(PS_TYPE_F32));
        tmp->dparams = (psF32 *) psAlloc(8 * sizeof(PS_TYPE_F32));
        break;
    case PS_MODEL_SERSIC_CORE:
        tmp->Nparams = 12;
        tmp->params = (psF32 *) psAlloc(12 * sizeof(PS_TYPE_F32));
        tmp->dparams = (psF32 *) psAlloc(12 * sizeof(PS_TYPE_F32));
        break;
    default:
        psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
        return(NULL);
    }

    for (psS32 i = 0 ; i < tmp->Nparams ; i++) {
        tmp->params[i] = 0.0;
        tmp->dparams[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;
    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 psList, if the list
argument is NULL, otherwise it adds the peak to that list.
 
XXX: Switch row, col args?
 *****************************************************************************/
psList *MyListAddPeak(psList *list,
                      psS32 row,
                      psS32 col,
                      psF32 counts,
                      psPeakType type)
{
    psPeak *tmpPeak = pmPeakAlloc(row, col, counts, type);

    if (list == NULL) {
        list = psListAlloc(tmpPeak);
    } else {
        psListAdd(list, PS_LIST_HEAD, tmpPeak);
    }

    return(list);
}

/******************************************************************************
pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage
above the given threshold.  Returns a psList 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 psList elements, should we use the image row/column offsets?
     Currently, we do not.
 *****************************************************************************/
psList *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;
    psList *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, 0);
        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 psList 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?
 *****************************************************************************/
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);
}

/******************************************************************************
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).
    //
    psS32 SubImageCenterRow = peak->y;
    psS32 SubImageCenterCol = peak->x;
    psS32 SubImageStartRow = SubImageCenterRow - outerRadiusS32;
    psS32 SubImageEndRow = SubImageCenterRow + outerRadiusS32;
    psS32 SubImageStartCol = SubImageCenterCol - outerRadiusS32;
    psS32 SubImageEndCol = 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);

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

    //
    // 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.
    //
    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.
  *****************************************************************************/
bool CheckRadius2(psF32 xCenter,
                  psF32 yCenter,
                  psF32 radius,
                  psF32 x,
                  psF32 y)
{
    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;
    //
    // 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.
    //
    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;
                    X2+= PS_SQR(xDiff) * pDiff;
                    Y2+= PS_SQR(yDiff) * pDiff;
                    XY+= xDiff * 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
    //
    source->moments->x = X1/Sum + ((psF32) source->peak->x);
    source->moments->y = Y1/Sum + ((psF32) source->peak->y);
    source->moments->Sx = sqrt(X2/Sum - PS_SQR(X1/Sum));
    source->moments->Sy = sqrt(Y2/Sum - PS_SQR(Y1/Sum));
    source->moments->Sxy = XY/Sum;
    source->moments->Peak = peakPixel;
    source->moments->nPixels = numPixels;

    return(source);
}

/******************************************************************************
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 SATURATE 0.0
#define FAINT_SN_LIM 0.0
#define PSF_SN_LIM 0.0
#define SATURATE 0.0
#define SATURATE 0.0

bool pmSourceRoughClass(psArray *source,
                        psMetadata *metadata)
{
    PS_PTR_CHECK_NULL(source, false);
    PS_PTR_CHECK_NULL(metadata, false);
    psBool rc = true;

    for (psS32 i = 0 ; i < source->n ; i++) {
        psSource *tmpSrc = (psSource *) source->data[i];
        PS_PTR_CHECK_NULL(tmpSrc->moments, false);
        tmpSrc->peak->class = 0;

        psF32 sigX = 0.0;
        psF32 sigY = 0.0;
        // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE.
        psF32 clumpX = 0.0;
        psF32 clumpDX = 0.0;
        psF32 clumpY = 0.0;
        psF32 clumpDY = 0.0;

        if (tmpSrc->moments->Peak > SATURATE) {
            tmpSrc->peak->class|= PS_SOURCE_SATURATED;
        } else {
            // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE.
            psF32 gain = 0.0;
            psF32 readNoise = 0.0;
            psF32 S = tmpSrc->moments->Sum;
            psF32 A = PS_PI * tmpSrc->moments->Sx * tmpSrc->moments->Sy;
            psF32 B = tmpSrc->moments->Sky;
            psF32 SN = (PS_SQRT_F32(gain) * S) /
                       PS_SQRT_F32(S + (A * B) + ((A * readNoise * readNoise) / PS_SQRT_F32(gain)));
            if (SN < FAINT_SN_LIM) {
                tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR;
            }
            if (SN < PSF_SN_LIM) {
                tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR;
            }
            // XXX: The SDRS is not real clear on how to calculate sigX, sigY.
            if ((fabs(sigX - clumpX) < clumpDX) &&
                    (fabs(sigY - clumpY) < clumpDY)) {
                tmpSrc->peak->class|= PS_SOURCE_PSFSTAR;
            }

            if ((sigX < (clumpX - clumpDX)) &&
                    (sigY < (clumpY - clumpDY)))
                tmpSrc->peak->class|= PS_SOURCE_DEFECT;
        }

        if ((sigX > (clumpX + clumpDX)) &&
                (sigY > (clumpY + clumpDY))) {
            tmpSrc->peak->class|= PS_SOURCE_GALAXY;
        }

        if (tmpSrc->peak->class == 0) {
            tmpSrc->peak->class|= PS_SOURCE_OTHER;
        }
    }

    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;
    psS32 SubImageStartRow = SubImageCenterRow - radiusS32;
    psS32 SubImageEndRow = SubImageCenterRow + radiusS32;
    psS32 SubImageStartCol = SubImageCenterCol - radiusS32;
    psS32 SubImageEndCol = SubImageCenterCol + radiusS32;

    if (SubImageStartRow < 0) {
        psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
                SubImageStartRow);
        return(false);
    }
    if (SubImageEndRow+1 >= 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+1 >= image->numCols) {
        psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
                SubImageEndCol);
        return(false);
    }

    // XXX: Must recycle image.
    if (source->pixels != NULL) {
        psLogMsg(__func__, PS_LOG_WARN,
                 "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL.  Freeing and reallocating.\n");
        psFree(source->pixels);
    }
    source->pixels = psImageSubset((psImage *) image,
                                   SubImageStartCol,
                                   SubImageStartRow,
                                   SubImageEndCol+1,
                                   SubImageEndRow+1);

    // XXX: Must recycle image.
    if (source->mask != NULL) {
        psFree(source->mask);
    }
    source->mask = psImageAlloc(1 + 2 * radiusS32,
                                1 + 2 * radiusS32,
                                PS_TYPE_F32);

    //
    // Loop through the subimage mask, initialize mask to 0 or 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] = 1;
            } else {
                source->mask->data.U8[row][col] = 0;
            }
        }
    }
    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);

    switch (model) {
    case PS_MODEL_GAUSS:
        source->models->params[0] = source->moments->Sky;
        source->models->params[1] = source->peak->counts - source->moments->Sky;
        source->models->params[2] = source->moments->x;
        source->models->params[3] = source->moments->y;
        source->models->params[4] = sqrt(2.0) / source->moments->Sx;
        source->models->params[5] = sqrt(2.0) / source->moments->Sy;
        source->models->params[6] = source->moments->Sxy;
        return(true);

    case PS_MODEL_PGAUSS:
        source->models->params[0] = source->moments->Sky;
        source->models->params[1] = source->peak->counts - source->moments->Sky;
        source->models->params[2] = source->moments->x;
        source->models->params[3] = source->moments->y;
        source->models->params[4] = sqrt(2.0) / source->moments->Sx;
        source->models->params[5] = sqrt(2.0) / source->moments->Sy;
        source->models->params[6] = source->moments->Sxy;
        return(true);

    case PS_MODEL_TWIST_GAUSS:
        source->models->params[0] = source->moments->Sky;
        source->models->params[1] = source->peak->counts - source->moments->Sky;
        source->models->params[2] = source->moments->x;
        source->models->params[3] = source->moments->y;
        // XXX: What are these?
        // source->models->params[4] = SxInner;
        // source->models->params[5] = SyInner;
        // source->models->params[6] = SxyInner;
        // source->models->params[7] = SxOuter;
        // source->models->params[8] = SyOuter;
        // source->models->params[9] = SxyOuter;
        // source->models->params[10] = N;
        return(true);

    case PS_MODEL_WAUSS:
        source->models->params[0] = source->moments->Sky;
        source->models->params[1] = source->peak->counts - source->moments->Sky;
        source->models->params[2] = source->moments->x;
        source->models->params[3] = source->moments->y;
        source->models->params[4] = sqrt(2.0) / source->moments->Sx;
        source->models->params[5] = sqrt(2.0) / source->moments->Sy;
        source->models->params[6] = source->moments->Sxy;
        // XXX: What are these?
        // source->models->params[7] = B2;
        // source->models->params[8] = B3;
        return(true);

    case PS_MODEL_SERSIC:
        source->models->params[0] = source->moments->Sky;
        source->models->params[1] = source->peak->counts - source->moments->Sky;
        source->models->params[2] = source->moments->x;
        source->models->params[3] = source->moments->y;
        source->models->params[4] = sqrt(2.0) / source->moments->Sx;
        source->models->params[5] = sqrt(2.0) / source->moments->Sy;
        source->models->params[6] = source->moments->Sxy;
        // XXX: What are these?
        //source->models->params[7] = Nexp;
        return(true);

    case PS_MODEL_SERSIC_CORE:
        source->models->params[0] = source->moments->Sky;
        source->models->params[1] = source->peak->counts - source->moments->Sky;
        source->models->params[2] = source->moments->x;
        source->models->params[3] = source->moments->y;
        // XXX: What are these?
        //source->models->params[4] SxInner;
        //source->models->params[5] SyInner;
        //source->models->params[6] SxyInner;
        //source->models->params[7] Zd;
        //source->models->params[8] SxOuter;
        //source->models->params[9] SyOuter;
        //source->models->params[10] = SxyOuter;
        //source->models->params[11] = Nexp;
        return(true);

    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.
 *****************************************************************************/
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.
    psVector *params = psVectorAlloc(7, PS_TYPE_F32);
    for (psS32 i = 0 ; i < src->models->Nparams ; i++) {
        params->data.F32[i] = src->models->params[i];
    }

    //
    // 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(params);
    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.
 *****************************************************************************/
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: What should these values be?
#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 100
#define PM_SOURCE_FIT_MODEL_TOLERANCE 1.0
/******************************************************************************
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;
    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++;
            }
        }
    }
    psArray *x = psArrayAlloc(count);
    psVector *y = 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:
                coord->data.F32[0] = (psF32) (i + source->pixels->row0);
                coord->data.F32[1] = (psF32) (j + source->pixels->col0);
                x->data[tmpCnt] = (psPtr *) coord;
                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
                tmpCnt++;
            }
        }
    }

    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
                            PM_SOURCE_FIT_MODEL_TOLERANCE);

    psVector *params = psVectorAlloc(source->models->Nparams, PS_TYPE_F32);

    switch (source->models->type) {
    case PS_MODEL_GAUSS:

        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec);
        break;
    case PS_MODEL_PGAUSS:
        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_PsuedoGauss2D_Vec);
        break;
    case PS_MODEL_TWIST_GAUSS:
        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_Wauss2D_Vec);
        break;
    case PS_MODEL_WAUSS:
        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_TwistGauss2D_Vec);
        break;
    case PS_MODEL_SERSIC:
        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_Sersic_Vec);
        break;
    case PS_MODEL_SERSIC_CORE:
        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_SersicCore_Vec);
        break;
    default:
        psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
        rc = false;
    }

    psFree(x);
    psFree(y);
    psFree(myMin);
    psFree(params);
    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 *params = psVectorAlloc(src->models->Nparams, PS_TYPE_F32);
    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
    for (psS32 i = 0 ; i < src->models->Nparams ; i++) {
        params->data.F32[i] = src->models->params[i];
    }

    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);
                psFree(params);
                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);
    psFree(params);
    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); \
}


/******************************************************************************
pmMinLM_Gauss2D(*deriv, *params, *x): the argument "x" contains a single "x,
y" coordinate pair.  This function computes the gaussian, specified by the
parameters in "params" at that x,y point and returns the value.  The
derivatives are also caculated and returned in the "deriv" argument.
 
    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;
 
XXX: Consider getting rid of the parameter checks since this might consume
a significant fraction of this function CPU time.
 
XXX: I added the following.  Must conform with IfA.  If deriv==NULL, then
we simply don't perform the derivative calculations.
 
XXX: It is not clear whether the subImage coords, or the image coords should
be used when calling these functions.  I don't think it really matters, as
long as we are consistent.
 *****************************************************************************/
psF32 pmMinLM_Gauss2D(psVector *deriv,
                      psVector *params,
                      psVector *x)
{
    PS_VECTOR_CHECK_NULL(params, NAN);
    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    PS_VECTOR_CHECK_SIZE(params, 7, NAN);
    PS_VECTOR_CHECK_NULL(x, NAN);
    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    PS_VECTOR_CHECK_SIZE(x, 2, NAN);

    psF32 X = x->data.F32[0] - params->data.F32[2];
    psF32 Y = x->data.F32[1] - params->data.F32[3];
    psF32 px = params->data.F32[4]*X;
    psF32 py = params->data.F32[5]*Y;
    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
    psF32 r = exp(-z);
    psF32 f = params->data.F32[1]*r + params->data.F32[0];

    if (deriv != NULL) {
        PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
        PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);

        psF32 q = params->data.F32[1]*r;
        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);
}

/******************************************************************************
p_pmMinLM_Gauss2D_Vec(*deriv, *params, *x): this function wraps the above
function in a form that is usable in the LM minimization routines.
 *****************************************************************************/
psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv,
                                psVector *params,
                                psArray *x)
{
    PS_IMAGE_CHECK_NULL(deriv, NULL);
    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    PS_VECTOR_CHECK_NULL(params, NULL);
    PS_VECTOR_CHECK_EMPTY(params, NULL);
    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    PS_PTR_CHECK_NULL(x, NULL);
    if (deriv->numRows != x->n) {
        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    }
    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    // XXX: use static memory here.
    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);

    for (psS32 i = 0 ; i < x->n ; i++) {
        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
            tmpRow->data.F32[j] = deriv->data.F32[i][j];
        }

        psVector *tmpVec2 = (psVector *) x->data[i];
        tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow,
                                              params,
                                              tmpVec2);

        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
            deriv->data.F32[i][j] = tmpRow->data.F32[j];
        }
    }

    psFree(tmpRow);
    return(tmpVec);
}


/******************************************************************************
    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;
 *****************************************************************************/
psF32 pmMinLM_PsuedoGauss2D(psVector *deriv,
                            psVector *params,
                            psVector *x)
{
    PS_VECTOR_CHECK_NULL(params, NAN);
    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    PS_VECTOR_CHECK_SIZE(params, 7, NAN);
    PS_VECTOR_CHECK_NULL(x, NAN);
    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    PS_VECTOR_CHECK_SIZE(x, 2, NAN);

    psF32 X = x->data.F32[0] - params->data.F32[2];
    psF32 Y = x->data.F32[1] - params->data.F32[3];
    psF32 px = params->data.F32[4]*X;
    psF32 py = params->data.F32[5]*Y;
    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
    psF32 t = 1 + z + 0.5*z*z;
    psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */
    psF32 f = params->data.F32[1]*r + params->data.F32[0];

    if (deriv != NULL) {
        PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
        PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);

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

/******************************************************************************
p_pmMinLM_PsuedoGauss2D_Vec(*deriv, *params, *x): this function wraps the
above function in a form that is usable in the LM minimization routines.
 *****************************************************************************/
psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv,
                                      psVector *params,
                                      psArray *x)
{
    PS_IMAGE_CHECK_NULL(deriv, NULL);
    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    PS_VECTOR_CHECK_NULL(params, NULL);
    PS_VECTOR_CHECK_EMPTY(params, NULL);
    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    PS_PTR_CHECK_NULL(x, NULL);
    if (deriv->numRows != x->n) {
        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    }
    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    // XXX: use static memory here.
    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);

    for (psS32 i = 0 ; i < x->n ; i++) {
        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
            tmpRow->data.F32[j] = deriv->data.F32[i][j];
        }

        tmpVec->data.F32[i] = pmMinLM_PsuedoGauss2D(tmpRow,
                              params,
                              (psVector *) x->data[i]);

        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
            deriv->data.F32[i][j] = tmpRow->data.F32[j];
        }
    }

    psFree(tmpRow);
    return(tmpVec);
}




/******************************************************************************
    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;
 *****************************************************************************/
psF32 pmMinLM_Wauss2D(psVector *deriv,
                      psVector *params,
                      psVector *x)
{
    PS_VECTOR_CHECK_NULL(params, NAN);
    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    PS_VECTOR_CHECK_SIZE(params, 9, NAN);
    PS_VECTOR_CHECK_NULL(x, NAN);
    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    PS_VECTOR_CHECK_SIZE(x, 2, NAN);

    psF32 X = x->data.F32[0] - params->data.F32[2];
    psF32 Y = x->data.F32[1] - params->data.F32[2];
    psF32 px = params->data.F32[4]*X;
    psF32 py = params->data.F32[5]*Y;
    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
    psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0);
    psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */
    psF32 f = params->data.F32[1]*r + params->data.F32[0];

    if (deriv != NULL) {
        PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
        PS_VECTOR_CHECK_SIZE(deriv, 9, NAN);
        //
        // 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);
}

/******************************************************************************
p_pmMinLM_Wauss2D_Vec(*deriv, *params, *x): this function wraps the above
function in a form that is usable in the LM minimization routines.
 *****************************************************************************/
psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv,
                                psVector *params,
                                psArray *x)
{
    PS_IMAGE_CHECK_NULL(deriv, NULL);
    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    PS_VECTOR_CHECK_NULL(params, NULL);
    PS_VECTOR_CHECK_EMPTY(params, NULL);
    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    PS_PTR_CHECK_NULL(x, NULL);
    if (deriv->numRows != x->n) {
        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    }
    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    // XXX: use static memory here.
    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);

    for (psS32 i = 0 ; i < x->n ; i++) {
        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
            tmpRow->data.F32[j] = deriv->data.F32[i][j];
        }

        tmpVec->data.F32[i] = pmMinLM_Wauss2D(tmpRow,
                                              params,
                                              (psVector *) x->data[i]);

        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
            deriv->data.F32[i][j] = tmpRow->data.F32[j];
        }
    }

    psFree(tmpRow);
    return(tmpVec);
}






// 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;
 *****************************************************************************/
psF32 pmMinLM_TwistGauss2D(psVector *deriv,
                           psVector *params,
                           psVector *x)
{
    PS_VECTOR_CHECK_NULL(params, NAN);
    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    PS_VECTOR_CHECK_SIZE(params, 11, NAN);
    PS_VECTOR_CHECK_NULL(x, NAN);
    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    PS_VECTOR_CHECK_SIZE(x, 2, NAN);

    psF32 X = x->data.F32[0] - params->data.F32[2];
    psF32 Y = x->data.F32[1] - params->data.F32[3];
    psF32 px1 = params->data.F32[4]*X;
    psF32 py1 = params->data.F32[5]*Y;
    psF32 px2 = params->data.F32[7]*X;
    psF32 py2 = params->data.F32[8]*Y;
    psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y;
    psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y;
    psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10]));


    psF32 f = params->data.F32[5]*r + params->data.F32[6];

    if (deriv != NULL) {
        PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
        PS_VECTOR_CHECK_SIZE(deriv, 11, NAN);

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

/******************************************************************************
p_pmMinLM_TwistGauss2D_Vec(*deriv, *params, *x): this function wraps the above
function in a form that is usable in the LM minimization routines.
 *****************************************************************************/
psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv,
                                     psVector *params,
                                     psArray *x)
{
    PS_IMAGE_CHECK_NULL(deriv, NULL);
    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    PS_VECTOR_CHECK_NULL(params, NULL);
    PS_VECTOR_CHECK_EMPTY(params, NULL);
    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    PS_PTR_CHECK_NULL(x, NULL);
    if (deriv->numRows != x->n) {
        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    }
    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    // XXX: use static memory here.
    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);

    for (psS32 i = 0 ; i < x->n ; i++) {
        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
            tmpRow->data.F32[j] = deriv->data.F32[i][j];
        }

        tmpVec->data.F32[i] = pmMinLM_TwistGauss2D(tmpRow,
                              params,
                              (psVector *) x->data[i]);

        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
            deriv->data.F32[i][j] = tmpRow->data.F32[j];
        }
    }

    psFree(tmpRow);
    return(tmpVec);
}



/******************************************************************************
    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;
 *****************************************************************************/
psF32 pmMinLM_Sersic(psVector *deriv,
                     psVector *params,
                     psVector *x)
{
    PS_VECTOR_CHECK_NULL(params, NAN);
    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    PS_VECTOR_CHECK_SIZE(params, 8, NAN);
    PS_VECTOR_CHECK_NULL(x, NAN);
    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    PS_VECTOR_CHECK_SIZE(x, 2, NAN);

    if (deriv != NULL) {
        PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
        PS_VECTOR_CHECK_SIZE(deriv, 8, NAN);
    }

    psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
    return(0.0);
}
/******************************************************************************
p_pmMinLM_Sersic_Vec(*deriv, *params, *x): this function wraps the above
function in a form that is usable in the LM minimization routines.
 *****************************************************************************/
psVector *p_pmMinLM_Sersic_Vec(psImage *deriv,
                               psVector *params,
                               psArray *x)
{
    PS_IMAGE_CHECK_NULL(deriv, NULL);
    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    PS_VECTOR_CHECK_NULL(params, NULL);
    PS_VECTOR_CHECK_EMPTY(params, NULL);
    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    PS_PTR_CHECK_NULL(x, NULL);
    if (deriv->numRows != x->n) {
        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    }
    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    // XXX: use static memory here.
    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);

    for (psS32 i = 0 ; i < x->n ; i++) {
        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
            tmpRow->data.F32[j] = deriv->data.F32[i][j];
        }

        tmpVec->data.F32[i] = pmMinLM_Sersic(tmpRow,
                                             params,
                                             (psVector *) x->data[i]);

        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
            deriv->data.F32[i][j] = tmpRow->data.F32[j];
        }
    }

    psFree(tmpRow);
    return(tmpVec);
}

/******************************************************************************
    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;
 *****************************************************************************/
psF32 pmMinLM_SersicCore(psVector *deriv,
                         psVector *params,
                         psVector *x)
{
    PS_VECTOR_CHECK_NULL(params, NAN);
    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    PS_VECTOR_CHECK_SIZE(params, 12, NAN);
    PS_VECTOR_CHECK_NULL(x, NAN);
    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    PS_VECTOR_CHECK_SIZE(x, 2, NAN);

    if (deriv != NULL) {
        PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
        PS_VECTOR_CHECK_SIZE(deriv, 12, NAN);
    }

    psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
    return(0.0);
}
/******************************************************************************
p_pmMinLM_SersicCore_Vec(*deriv, *params, *x): this function wraps the above
function in a form that is usable in the LM minimization routines.
 *****************************************************************************/
psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv,
                                   psVector *params,
                                   psArray *x)
{
    PS_IMAGE_CHECK_NULL(deriv, NULL);
    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    PS_VECTOR_CHECK_NULL(params, NULL);
    PS_VECTOR_CHECK_EMPTY(params, NULL);
    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    PS_PTR_CHECK_NULL(x, NULL);
    if (deriv->numRows != x->n) {
        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    }
    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    // XXX: use static memory here.
    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);

    for (psS32 i = 0 ; i < x->n ; i++) {
        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
            tmpRow->data.F32[j] = deriv->data.F32[i][j];
        }

        tmpVec->data.F32[i] = pmMinLM_SersicCore(tmpRow,
                              params,
                              (psVector *) x->data[i]);

        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
            deriv->data.F32[i][j] = tmpRow->data.F32[j];
        }
    }

    psFree(tmpRow);
    return(tmpVec);
}



/******************************************************************************
 *****************************************************************************/
psF32 pmMinLM_PsuedoSersic(psVector *deriv,
                           psVector *params,
                           psVector *x)
{
    return(0.0);
}
