/** @file  pmSource.c
 *
 *  Functions to define and manipulate sources on images
 *
 *  @author GLG, MHPCC
 *  @author EAM, IfA: significant modifications.
 *
 *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2008-06-10 02:26:17 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <strings.h>
#include <pslib.h>
#include "pmHDU.h"
#include "pmFPA.h"
#include "pmFPAMaskWeight.h"
#include "pmPeaks.h"
#include "pmMoments.h"
#include "pmResiduals.h"
#include "pmGrowthCurve.h"
#include "pmTrend2D.h"
#include "pmPSF.h"
#include "pmModel.h"
#include "pmSource.h"

static void sourceFree(pmSource *tmp)
{
    if (!tmp)
        return;

    psTrace("psModules.objects", 5, "---- begin ----\n");
    psFree(tmp->peak);
    psFree(tmp->pixels);
    psFree(tmp->weight);
    psFree(tmp->maskObj);
    psFree(tmp->maskView);
    psFree(tmp->modelFlux);
    psFree(tmp->psfFlux);
    psFree(tmp->moments);
    psFree(tmp->modelPSF);
    psFree(tmp->modelEXT);
    psFree(tmp->modelFits);
    psFree(tmp->extpars);
    psFree(tmp->blends);
    psTrace("psModules.objects", 5, "---- end ----\n");
}

// free only the pixel data associated with this source
void pmSourceFreePixels(pmSource *source)
{

    if (!source)
        return;

    psFree (source->pixels);
    psFree (source->weight);
    psFree (source->maskObj);
    psFree (source->maskView);
    psFree (source->modelFlux);
    psFree (source->psfFlux);

    source->pixels = NULL;
    source->weight = NULL;
    source->maskObj = NULL;
    source->maskView = NULL;
    source->modelFlux = NULL;
    source->psfFlux = NULL;
    return;
}

bool psMemCheckSource(psPtr ptr)
{
    PS_ASSERT_PTR(ptr, false);
    return ( psMemGetDeallocator(ptr) == (psFreeFunc) sourceFree);
}

/******************************************************************************
pmSourceAlloc(): Allocate the pmSource structure and initialize its members
to NULL.
*****************************************************************************/

pmSource *pmSourceAlloc()
{
    psTrace("psModules.objects", 5, "---- begin ----\n");
    static int id = 1;
    pmSource *source = (pmSource *) psAlloc(sizeof(pmSource));
    *(int *)&source->id = id++;
    source->seq = -1;
    source->peak = NULL;
    source->pixels = NULL;
    source->weight = NULL;
    source->maskObj = NULL;
    source->maskView = NULL;
    source->modelFlux = NULL;
    source->psfFlux = NULL;
    source->moments = NULL;
    source->blends = NULL;
    source->modelPSF = NULL;
    source->modelEXT = NULL;
    source->modelFits = NULL;
    source->type = PM_SOURCE_TYPE_UNKNOWN;
    source->mode = PM_SOURCE_MODE_DEFAULT;
    source->extpars = NULL;
    source->region = psRegionSet(NAN, NAN, NAN, NAN);
    psMemSetDeallocator(source, (psFreeFunc) sourceFree);

    // default values are NAN
    source->psfMag = NAN;
    source->extMag = NAN;
    source->errMag = NAN;
    source->apMag  = NAN;
    source->sky    = NAN;
    source->skyErr = NAN;
    source->pixWeight = NAN;

    source->psfChisq = NAN;
    source->crNsigma = NAN;
    source->extNsigma = NAN;

    psTrace("psModules.objects", 5, "---- end ----\n");
    return(source);
}

/******************************************************************************
pmSourceCopy(): copy the pmSource structure and contents
XXX : are we OK with incrementing the ID?
*****************************************************************************/
pmSource *pmSourceCopy(pmSource *in)
{
    if (in == NULL) {
        return(NULL);
    }
    // this copy is used to allow multiple fit attempts on the
    // same object.  the pixels, weight, and mask arrays all point to
    // the same original subarrays.  the peak and moments point at
    // the original values.
    pmSource *source = pmSourceAlloc ();

    // this is actually the same peak as the original, is this the correct way to handle this?
    if (in->peak != NULL) {
        source->peak = pmPeakAlloc (in->peak->x, in->peak->y, in->peak->value, in->peak->type);
        source->peak->xf = in->peak->xf;
        source->peak->yf = in->peak->yf;
        source->peak->flux = in->peak->flux;
        source->peak->SN = in->peak->SN;
    }

    // copy the values in the moments structure
    if (in->moments != NULL) {
        source->moments  =  pmMomentsAlloc();
        *source->moments = *in->moments;
    }

    // These images are all views to the parent.
    // We want a new view, but pointing at the same pixels.
    source->pixels   = psImageCopyView(NULL, in->pixels);
    source->weight   = psImageCopyView(NULL, in->weight);
    source->maskView = in->maskView ? psImageCopyView(NULL, in->maskView) : NULL;

    // the maskObj is a unique mask array; create a new mask image
    source->maskObj = in->maskObj ? psImageCopy (NULL, in->maskObj, PS_TYPE_MASK) : NULL;

    source->type = in->type;
    source->mode = in->mode;

    return(source);
}

// x,y are defined in the parent image coords of readout->image
bool pmSourceDefinePixels(pmSource *mySource,
                          const pmReadout *readout,
                          psF32 x,
                          psF32 y,
                          psF32 Radius)
{
    PS_ASSERT_PTR_NON_NULL(mySource, false);
    PS_ASSERT_PTR_NON_NULL(readout, false);
    PS_ASSERT_PTR_NON_NULL(readout->image, false);
    PS_ASSERT_INT_POSITIVE(Radius, false);

    psRegion srcRegion;

    // Grab a subimage of the original image of size (2 * outerRadius).
    srcRegion = psRegionForSquare (x, y, Radius);
    srcRegion = psRegionForImage (readout->image, srcRegion);

    // these images are subset images of the equivalent parents
    mySource->pixels = psImageSubset(readout->image, srcRegion);
    if (readout->weight) {
        mySource->weight = psImageSubset(readout->weight, srcRegion);
    }
    if (readout->mask) {
        mySource->maskView = psImageSubset(readout->mask,  srcRegion);
        // the object mask is a copy, and used to define the source pixels
        mySource->maskObj = psImageCopy(NULL, mySource->maskView, PS_TYPE_MASK);
    }
    mySource->region   = srcRegion;

    return true;
}

bool pmSourceRedefinePixels(pmSource *mySource,
                            const pmReadout *readout,
                            psF32 x,
                            psF32 y,
                            psF32 Radius)
{
    PS_ASSERT_PTR_NON_NULL(mySource, false);
    PS_ASSERT_PTR_NON_NULL(readout, false);
    PS_ASSERT_PTR_NON_NULL(readout->image, false);
    PS_ASSERT_INT_POSITIVE(Radius, false);

    bool extend;
    psRegion newRegion;

    // check to see if new region is completely contained within old region
    newRegion = psRegionForSquare (x, y, Radius);
    newRegion = psRegionForImage (readout->image, newRegion);

    extend = false;
    extend |= (int)(newRegion.x0) < (int)(mySource->region.x0);
    extend |= (int)(newRegion.x1) > (int)(mySource->region.x1);
    extend |= (int)(newRegion.y0) < (int)(mySource->region.y0);
    extend |= (int)(newRegion.y1) > (int)(mySource->region.y1);

    extend |= (mySource->pixels == NULL);
    extend |= (mySource->weight == NULL);
    extend |= (mySource->maskObj == NULL);
    extend |= (mySource->maskView == NULL);

    // extend = true;
    if (extend) {
        // re-create the subimage
        psFree (mySource->pixels);
        psFree (mySource->weight);
        psFree (mySource->maskView);

        mySource->pixels   = psImageSubset(readout->image,  newRegion);
        mySource->weight   = psImageSubset(readout->weight, newRegion);
        mySource->maskView = psImageSubset(readout->mask,   newRegion);
        mySource->region   = newRegion;

        // re-copy the main mask pixels.  NOTE: the user will need to reset the object mask
        // pixels (eg, with psImageKeepCircle)
        mySource->maskObj = psImageCopy (mySource->maskObj, mySource->maskView, PS_TYPE_MASK);

        // drop the old modelFlux pixels and force the user to re-create
        psFree (mySource->modelFlux);
        mySource->modelFlux = NULL;

        // drop the old psfFlux pixels and force the user to re-create
        psFree (mySource->psfFlux);
        mySource->psfFlux = NULL;
    }
    return extend;
}

/******************************************************************************
    pmSourcePSFClump(source, recipe): Find the likely PSF clump in the
    sigma-x, sigma-y plane. return 0,0 clump in case of error.
*****************************************************************************/

// XXX EAM include a S/N cutoff in selecting the sources?
// XXX this function should probably accept the values, not a recipe. wrap with a
// psphot-specific function which applies the recipe values
pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *recipe)
{
    psTrace("psModules.objects", 5, "---- begin ----\n");

    # define NPIX 10
    # define SCALE 0.1

    psArray *peaks  = NULL;
    pmPSFClump errorClump = {-1.0, -1.0, 0.0, 0.0};
    pmPSFClump emptyClump = {+0.0, +0.0, 0.0, 0.0};
    pmPSFClump psfClump;

    PS_ASSERT_PTR_NON_NULL(sources, errorClump);
    PS_ASSERT_PTR_NON_NULL(recipe, errorClump);

    bool status = true;
    float PSF_CLUMP_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_SN_LIM");
    if (!status) {
        PSF_CLUMP_SN_LIM = 0;
    }

    // find the sigmaX, sigmaY clump
    {
        psStats *stats  = NULL;
        psImage *splane = NULL;
        int binX, binY;
        bool status;

        psF32 SX_MAX = psMetadataLookupF32 (&status, recipe, "MOMENTS_SX_MAX");
        if (!status)
            SX_MAX = 10.0;
        psF32 SY_MAX = psMetadataLookupF32 (&status, recipe, "MOMENTS_SY_MAX");
        if (!status)
            SY_MAX = 10.0;
        psF32 AR_MAX = psMetadataLookupF32 (&status, recipe, "MOMENTS_AR_MAX");
        if (!status)
            AR_MAX =  3.0;
        psF32 AR_MIN = 1.0 / AR_MAX;

        // construct a sigma-plane image
        splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);
        psImageInit(splane, 0);  // psImageAlloc doesn't zero the data

        // place the sources in the sigma-plane image (ignore 0,0 values?)
        int nValid = 0;
        for (psS32 i = 0 ; i < sources->n ; i++)
        {
            pmSource *tmpSrc = (pmSource *) sources->data[i];
            if (tmpSrc == NULL)
                continue;
            if (tmpSrc->moments == NULL)
                continue;
            if (tmpSrc->moments->SN < PSF_CLUMP_SN_LIM)
                continue;

            // Sx,Sy are limited at 0.  a peak at 0,0 is artificial
            if (fabs(tmpSrc->moments->Sx) < 0.05)
                continue;
            if (fabs(tmpSrc->moments->Sy) < 0.05)
                continue;
            if ((tmpSrc->moments->Sx / tmpSrc->moments->Sy) > AR_MAX)
                continue;
            if ((tmpSrc->moments->Sx / tmpSrc->moments->Sy) < AR_MIN)
                continue;
            if (tmpSrc->mode & PM_SOURCE_MODE_BLEND)
                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;
            nValid ++;
        }

        // find the peak in this image
        stats = psStatsAlloc (PS_STAT_MAX);
        if (!psImageStats (stats, splane, NULL, 0)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
            psFree(stats);
            psFree(splane);
            return emptyClump;
        }
        peaks = pmPeaksInImage (splane, stats[0].max / 2);
        psTrace ("psModules.objects", 2, "clump threshold is %f\n", stats[0].max/2);

        const bool keep_psf_clump = psMetadataLookupBool(NULL, recipe, "KEEP_PSF_CLUMP");
        if (keep_psf_clump)
        {
            psMetadataAdd(recipe, PS_LIST_TAIL,
                          "PSF_CLUMP", PS_DATA_IMAGE, "Image of PSF coefficients", splane);
        }
        psFree (splane);
        psFree (stats);

        // if we failed to find a valid peak, return the empty clump (failure signal)
        if (peaks == NULL)
        {
            psLogMsg ("psphot", 3, "failed to find a peak in the PSF clump image\n");
            if (nValid == 0) {
                psLogMsg ("psphot", 3, "no valid sources kept for PSF search\n");
            } else {
                psLogMsg ("psphot", 3, "no significant peak\n");
            }
            return (emptyClump);
        }
    }

    // measure statistics on Sx, Sy if Sx, Sy within range of clump
    {
        pmPeak *clump;
        psF32 minSx, maxSx;
        psF32 minSy, maxSy;
        psVector *tmpSx = NULL;
        psVector *tmpSy = NULL;
        psStats *stats  = NULL;

        // select the single highest peak
        psArraySort (peaks, pmPeaksCompareDescend);
        clump = peaks->data[0];
        psTrace ("psModules.objects", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->value);

        // 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 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
        tmpSy = psVectorAllocEmpty (sources->n, PS_TYPE_F32);

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

            if (tmpSrc == NULL)
                continue;
            if (tmpSrc->moments == NULL)
                continue;
            if (tmpSrc->moments->SN < PSF_CLUMP_SN_LIM)
                continue;

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

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

        psVectorStats (stats, tmpSy, NULL, NULL, 0);
        psfClump.Y  = stats->clippedMean;
        psfClump.dY = stats->clippedStdev;

        psTrace ("psModules.objects", 2, "clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
        psTrace ("psModules.objects", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);

        psFree (stats);
        psFree (peaks);
        psFree (tmpSx);
        psFree (tmpSy);
    }

    psTrace("psModules.objects", 5, "---- end ----\n");
    return (psfClump);
}

/******************************************************************************
    pmSourceRoughClass(source, recipe): make a guess at the source
    classification.
    XXX: How can this function ever return FALSE?
*****************************************************************************/

bool pmSourceRoughClass(psArray *sources, psMetadata *recipe, pmPSFClump clump, psMaskType maskSat)
{
    psTrace("psModules.objects", 5, "---- begin ----");

    PS_ASSERT_PTR_NON_NULL(sources, false);
    PS_ASSERT_PTR_NON_NULL(recipe, false);

    int Nsat     = 0;
    int Next     = 0;
    int Nstar    = 0;
    int Npsf     = 0;
    int Ncr      = 0;
    int Nsatstar = 0;
    psRegion inner;

    // report stats on S/N values for star-like objects
    psVector *starsn_peaks = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    psVector *starsn_moments = psVectorAllocEmpty (sources->n, PS_TYPE_F32);

    // get basic parameters, or set defaults
    bool status;
    float PSF_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_SN_LIM");
    if (!status)
        PSF_SN_LIM = 20.0;
    float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
    if (!status)
        PSF_CLUMP_NSIGMA = 1.5;
    float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");

    // XXX allow clump size to be scaled relative to sigmas?
    // make rough IDs based on clumpX,Y,DX,DY
    for (psS32 i = 0 ; i < sources->n ; i++) {

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

        source->peak->type = 0;

        // we are basically classifying by moments; use the default if not found
        if (!source->moments) {
            source->type = PM_SOURCE_TYPE_STAR;
            Nstar++;
            continue;
        }

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

        // XXX EAM : can we use the value of SATURATE if mask is NULL?
        inner = psRegionForSquare (source->peak->x, source->peak->y, 2);
        inner = psRegionForImage (source->maskView, inner);
        int Nsatpix = psImageCountPixelMask (source->maskView, inner, maskSat);

        // saturated star (size consistent with PSF or larger)
        // Nsigma should be user-configured parameter
        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
        big = true;
        if ((Nsatpix > 1) && big) {
            source->type = PM_SOURCE_TYPE_STAR;
            source->mode |= PM_SOURCE_MODE_SATSTAR;
            // recalculate moments here with larger box?
            pmSourceMoments (source, INNER_RADIUS);
            Nsatstar ++;
            continue;
        }

        // saturated object (not a star, eg bleed trails, hot pixels)
        if (Nsatpix > 1) {
            source->type = PM_SOURCE_TYPE_SATURATED;
            source->mode |= PM_SOURCE_MODE_SATURATED;
            Nsat ++;
            continue;
        }

        // likely defect (too small to be stellar) (push out to 3 sigma)
        // low S/N objects which are small are probably stellar
        // only set candidate defects if
        // XXX these limits are quite arbitrary
        if ((sigX < 0.05) || (sigY < 0.05)) {
            source->type = PM_SOURCE_TYPE_DEFECT;
            source->mode |= PM_SOURCE_MODE_DEFECT;
            Ncr ++;
            continue;
        }

        // likely unsaturated extended source (too large to be stellar)
        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
            source->type = PM_SOURCE_TYPE_EXTENDED;
            Next ++;
            continue;
        }

        // the rest are probable stellar objects
        starsn_moments->data.F32[starsn_moments->n] = source->moments->SN;
        starsn_moments->n ++;
        starsn_peaks->data.F32[starsn_peaks->n] = source->peak->SN;
        starsn_peaks->n ++;
        Nstar ++;

        // PSF star (within 1.5 sigma of clump center, S/N > limit)
        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
        if ((source->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) {
            source->type = PM_SOURCE_TYPE_STAR;
            source->mode |= PM_SOURCE_MODE_PSFSTAR;
            Npsf ++;
            continue;
        }

        // random type of star
        source->type = PM_SOURCE_TYPE_STAR;
    }

    psStats *stats = NULL;
    stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    if (!psVectorStats (stats, starsn_moments, NULL, NULL, 0)) {
        // Don't care about this error
        psErrorClear();
    }
    psLogMsg ("pmObjects", 3, "SN range (moments): %f - %f\n", stats[0].min, stats[0].max);
    psFree (starsn_moments);
    psFree (stats);

    stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    if (!psVectorStats (stats, starsn_peaks, NULL, NULL, 0)) {
        // Don't care about this error
        psErrorClear();
    }
    psLogMsg ("pmObjects", 3, "SN range (peaks)  : %f - %f\n", stats[0].min, stats[0].max);
    psFree (starsn_peaks);
    psFree (stats);

    psTrace ("psModules.objects", 2, "Nstar:    %3d\n", Nstar);
    psTrace ("psModules.objects", 2, "Npsf:     %3d\n", Npsf);
    psTrace ("psModules.objects", 2, "Next:     %3d\n", Next);
    psTrace ("psModules.objects", 2, "Nsatstar: %3d\n", Nsatstar);
    psTrace ("psModules.objects", 2, "Nsat:     %3d\n", Nsat);
    psTrace ("psModules.objects", 2, "Ncr:      %3d\n", Ncr);

    psTrace("psModules.objects", 5, "---- end ----\n");
    return true;
}

/******************************************************************************
pmSourceMoments(source, radius): this function takes a subImage defined in the
pmSource data structure, along with the peak location, and determines the
various moments associated with that peak.

Requires the following to have been created:
    pmSource
    pmSource->peak
    pmSource->pixels
    pmSource->weight
    pmSource->mask

XXX: The peak calculations are done in image coords, not subImage coords.

XXX EAM : this version clips input pixels on S/N
XXX EAM : this version returns false for several reasons
*****************************************************************************/
# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)

bool pmSourceMoments(pmSource *source,
                     psF32 radius)
{
    psTrace("psModules.objects", 5, "---- begin ----\n");
    PS_ASSERT_PTR_NON_NULL(source, false);
    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);

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

    //
    // Sum = SUM (z - sky)
    // X1  = SUM (x - xc)*(z - sky)
    // X2  = SUM (x - xc)^2 * (z - sky)
    // XY  = SUM (x - xc)*(y - yc)*(z - sky)
    //
    psF32 peakPixel = -PS_MAX_F32;
    psS32 numPixels = 0;
    psF32 Sum = 0.0;
    psF32 Var = 0.0;
    psF32 X1 = 0.0;
    psF32 Y1 = 0.0;
    psF32 X2 = 0.0;
    psF32 Y2 = 0.0;
    psF32 XY = 0.0;
    psF32 x  = 0;
    psF32 y  = 0;
    psF32 R2 = PS_SQR(radius);

    psF32 xPeak = source->peak->x;
    psF32 yPeak = source->peak->y;
    psF32 xOff = source->pixels->col0 - source->peak->x;
    psF32 yOff = source->pixels->row0 - source->peak->y;

    // XXX why do I get different results for these two methods of finding Sx?
    // XXX Sx, Sy would be better measured if we clip pixels close to sky
    // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2?
    // We loop through all pixels in this subimage (source->pixels), and for each
    // pixel that is not masked, AND within the radius of the peak pixel, we
    // proceed with the moments calculation.  need to do two loops for a
    // numerically stable result.  first loop: get the sums.
    // XXX EAM : mask == 0 is valid

    for (psS32 row = 0; row < source->pixels->numRows ; row++) {

        psF32 *vPix = source->pixels->data.F32[row];
        psF32 *vWgt = source->weight->data.F32[row];
        psU8  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.U8[row];

        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
            if (vMsk) {
                if (*vMsk) {
                    vMsk++;
                    psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to mask: %d\n",
                            col, row, (int)*vMsk);
                    continue;
                }
                vMsk++;
            }
	    if (isnan(*vPix)) continue;

            psF32 xDiff = col + xOff;
            psF32 yDiff = row + yOff;

            // radius is just a function of (xDiff, yDiff)
            if (!VALID_RADIUS(xDiff, yDiff, R2)) {
#if 1
                psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to position: %f %f\n",
                        col, row, xDiff, yDiff);
#endif
                continue;
            }

            psF32 pDiff = *vPix - sky;
            psF32 wDiff = *vWgt;

            // XXX EAM : check for valid S/N in pixel
            // XXX EAM : should this limit be user-defined?
            if (PS_SQR(pDiff) < wDiff) {
#if 1
                psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to insignificance: %f, %f\n",
                        col, row, pDiff, wDiff);
#endif
                continue;
            }

            Var += wDiff;
            Sum += pDiff;

            psF32 xWght = xDiff * pDiff;
            psF32 yWght = yDiff * pDiff;

            X1  += xWght;
            Y1  += yWght;

            XY  += xDiff * yWght;
            X2  += xDiff * xWght;
            Y2  += yDiff * yWght;

            peakPixel = PS_MAX (*vPix, peakPixel);
            numPixels++;
        }
    }

    // if we have less than (1/4) of the possible pixels, force a retry
    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
        psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n",
                 numPixels, (int)(0.75*R2), Sum);
        psTrace("psModules.objects", 5, "---- end (false) ----\n");
        return (false);
    }

    psTrace ("psModules.objects", 4, "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
             sky, Sum, X1, Y1, X2, Y2, XY, numPixels);

    //
    // first moment X  = X1/Sum + xc
    // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
    // Sxy             = XY / Sum
    //
    x = X1/Sum;
    y = Y1/Sum;
    if ((fabs(x) > radius) || (fabs(y) > radius)) {
        psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n",
                 source->peak->x, source->peak->y);
        psTrace("psModules.objects", 5, "---- end(false)  ----\n");
        return (false);
    }

    source->moments->x = x + xPeak;
    source->moments->y = y + yPeak;

    // XXX EAM : Sxy needs to have x*y subtracted
    source->moments->Sxy = XY/Sum - x*y;
    source->moments->Sum = Sum;
    source->moments->SN  = Sum / sqrt(Var);
    source->moments->Peak = peakPixel;
    source->moments->nPixels = numPixels;

    // XXX EAM : these values can be negative, so we need to limit the range
    // XXX EAM : make the use of this consistent: should this be the second moment or sqrt?
    // source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
    // source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
    source->moments->Sx = PS_MAX(X2/Sum - PS_SQR(x), 0);
    source->moments->Sy = PS_MAX(Y2/Sum - PS_SQR(y), 0);

    psTrace ("psModules.objects", 4,
             "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n",
             sky, Sum, source->moments->x, source->moments->y,
             source->moments->Sx, source->moments->Sy, source->moments->Sxy);

    psTrace("psModules.objects", 5, "---- end ----\n");
    return(true);
}

// construct a realization of the source model
bool pmSourceCacheModel (pmSource *source, psMaskType maskVal) {
    PS_ASSERT_PTR_NON_NULL(source, false);
    // select appropriate model
    pmModel *model = pmSourceGetModel (NULL, source);
    if (model == NULL) return false;  // model must be defined

    // if we already have a cached image, re-use that memory
    source->modelFlux = psImageCopy (source->modelFlux, source->pixels, PS_TYPE_F32);
    psImageInit (source->modelFlux, 0.0);

    // in some places (psphotEnsemble), we need a normalized version
    // in others, we just want the model.  which is more commonly used?
    // modelFlux always has unity normalization (I0 = 1.0)
    pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
    return true;
}

// construct a realization of the source model
// XXX this function should optionally save an existing psf image from modelFlux
bool pmSourceCachePSF (pmSource *source, psMaskType maskVal) {
    PS_ASSERT_PTR_NON_NULL(source, false);

    // select appropriate model
    if (source->modelPSF == NULL) return false;  // model must be defined

    // if we already have a cached image, re-use that memory
    source->psfFlux = psImageCopy (source->psfFlux, source->pixels, PS_TYPE_F32);
    psImageInit (source->psfFlux, 0.0);

    // in some places (psphotEnsemble), we need a normalized version
    // in others, we just want the model.  which is more commonly used?
    // psfFlux always has unity normalization (I0 = 1.0)
    pmModelAdd (source->psfFlux, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
    return true;
}

// should we call pmSourceCacheModel if it does not exist?
bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add,
                 psMaskType maskVal, int dx, int dy)
{
    PS_ASSERT_PTR_NON_NULL(source, false);
    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    bool status;

    if (add) {
        psTrace ("psphot", 3, "replacing object at %f,%f\n", source->peak->xf, source->peak->yf);
    } else {
        psTrace ("psphot", 3, "removing object at %f,%f\n", source->peak->xf, source->peak->yf);
    }

    pmModel *model = pmSourceGetModel (NULL, source);
    if (model == NULL) return false;  // model must be defined

    if (source->modelFlux) {
        // add in the pixels from the modelFlux image
        int dX = source->modelFlux->col0 - source->pixels->col0;
        int dY = source->modelFlux->row0 - source->pixels->row0;
        assert (dX >= 0);
        assert (dY >= 0);
        assert (dX + source->modelFlux->numCols <= source->pixels->numCols);
        assert (dY + source->modelFlux->numRows <= source->pixels->numRows);

        // modelFlux has unity normalization
        float Io = model->params->data.F32[PM_PAR_I0];
        if (mode & PM_MODEL_OP_NORM) {
            Io = 1.0;
        }

        psU8 **mask = NULL;
        if (source->maskObj) {
            mask = source->maskObj->data.U8;
        }

        psF32 **target = source->pixels->data.F32;
        if (mode & PM_MODEL_OP_NOISE) {
	    // XXX need to scale by the gain...
            target = source->weight->data.F32;
        }

        // XXX need to respect the source and model masks?
        for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
            int oy = iy + dY;
            for (int ix = 0; ix < source->modelFlux->numCols; ix++) {
                int ox = ix + dX;
                if (mask && (mask[iy][ix] & maskVal)) continue;
                float value = Io*source->modelFlux->data.F32[iy][ix];
                if (add) {
                    target[oy][ox] += value;
                } else {
                    target[oy][ox] -= value;
                }
            }
        }
        return true;
    }

    psImage *target = source->pixels;
    if (mode & PM_MODEL_OP_NOISE) {
        target = source->weight;
    }

    if (add) {
        status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
    } else {
        status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
    }

    return true;
}

bool pmSourceAdd (pmSource *source, pmModelOpMode mode, psMaskType maskVal) {
    return pmSourceOp (source, mode, true, maskVal, 0, 0);
}

bool pmSourceSub (pmSource *source, pmModelOpMode mode, psMaskType maskVal) {
    return pmSourceOp (source, mode, false, maskVal, 0, 0);
}

bool pmSourceAddWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy) {
    return pmSourceOp (source, mode, true, maskVal, dx, dy);
}

bool pmSourceSubWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy) {
    return pmSourceOp (source, mode, false, maskVal, dx, dy);
}

// given a source, which model is currently appropriate?
// choose PSF or EXT based on source->type, but fall back on PSF
// if the EXT model is NULL
pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source)
{
    PS_ASSERT_PTR_NON_NULL(source, NULL);

    pmModel *model;

    if (isPSF) {
        *isPSF = false;
    }

    switch (source->type) {
      case PM_SOURCE_TYPE_STAR:
        model = source->modelPSF;
        if (model == NULL)
            return NULL;
        if (isPSF) {
            *isPSF = true;
        }
        return model;

	// the 'best' extended model is saved in source->modelEXT (may be a pointer to one of
	// the elements of source->modelFits)
      case PM_SOURCE_TYPE_EXTENDED:
	model = source->modelEXT;
        if (!model && source->modelPSF) {
	    // XXX raise an error or warning here?
            if (isPSF) {
                *isPSF = true;
            }
            return source->modelPSF;
        }
        return (model);
        break;

      default:
        return NULL;
    }
    return NULL;
}

// sort by SN (descending)
int pmSourceSortBySN (const void **a, const void **b)
{
    pmSource *A = *(pmSource **)a;
    pmSource *B = *(pmSource **)b;

    psF32 fA = (A->peak == NULL) ? 0 : A->peak->SN;
    psF32 fB = (B->peak == NULL) ? 0 : B->peak->SN;
    if (isnan (fA)) fA = 0;
    if (isnan (fB)) fB = 0;

    psF32 diff = fA - fB;
    if (diff > FLT_EPSILON) return (-1);
    if (diff < FLT_EPSILON) return (+1);
    return (0);
}

// sort by Y (ascending)
int pmSourceSortByY (const void **a, const void **b)
{
    pmSource *A = *(pmSource **)a;
    pmSource *B = *(pmSource **)b;

    psF32 fA = (A->peak == NULL) ? 0 : A->peak->y;
    psF32 fB = (B->peak == NULL) ? 0 : B->peak->y;

    psF32 diff = fA - fB;
    if (diff > FLT_EPSILON) return (+1);
    if (diff < FLT_EPSILON) return (-1);
    return (0);
}

pmSourceMode pmSourceModeFromString (const char *name) {
  if (!strcasecmp (name, "DEFAULT"   )) return PM_SOURCE_MODE_DEFAULT;
  if (!strcasecmp (name, "PSFMODEL"  )) return PM_SOURCE_MODE_PSFMODEL;
  if (!strcasecmp (name, "EXTMODEL"  )) return PM_SOURCE_MODE_EXTMODEL;
  if (!strcasecmp (name, "FITTED"    )) return PM_SOURCE_MODE_FITTED;
  if (!strcasecmp (name, "FAIL"      )) return PM_SOURCE_MODE_FAIL;
  if (!strcasecmp (name, "POOR"      )) return PM_SOURCE_MODE_POOR;
  if (!strcasecmp (name, "PAIR"      )) return PM_SOURCE_MODE_PAIR;
  if (!strcasecmp (name, "PSFSTAR"   )) return PM_SOURCE_MODE_PSFSTAR;
  if (!strcasecmp (name, "SATSTAR"   )) return PM_SOURCE_MODE_SATSTAR;
  if (!strcasecmp (name, "BLEND"     )) return PM_SOURCE_MODE_BLEND;
  if (!strcasecmp (name, "EXTERNAL"  )) return PM_SOURCE_MODE_EXTERNAL;
  if (!strcasecmp (name, "BADPSF"    )) return PM_SOURCE_MODE_BADPSF;
  if (!strcasecmp (name, "DEFECT"    )) return PM_SOURCE_MODE_DEFECT;
  if (!strcasecmp (name, "SATURATED" )) return PM_SOURCE_MODE_SATURATED;
  if (!strcasecmp (name, "CRLIMIT"   )) return PM_SOURCE_MODE_CR_LIMIT;
  if (!strcasecmp (name, "EXTLIMIT"  )) return PM_SOURCE_MODE_EXT_LIMIT;
  if (!strcasecmp (name, "SUBTRACTED")) return PM_SOURCE_MODE_SUBTRACTED;
  return PM_SOURCE_MODE_DEFAULT;
}

char *pmSourceModeToString (const pmSourceMode mode) {
  switch (mode) {
    case PM_SOURCE_MODE_DEFAULT    : return psStringCopy ("DEFAULT"   );
    case PM_SOURCE_MODE_PSFMODEL   : return psStringCopy ("PSFMODEL"  );
    case PM_SOURCE_MODE_EXTMODEL   : return psStringCopy ("EXTMODEL"  );
    case PM_SOURCE_MODE_FITTED     : return psStringCopy ("FITTED"    );
    case PM_SOURCE_MODE_FAIL       : return psStringCopy ("FAIL"      );
    case PM_SOURCE_MODE_POOR       : return psStringCopy ("POOR"      );
    case PM_SOURCE_MODE_PAIR       : return psStringCopy ("PAIR"      );
    case PM_SOURCE_MODE_PSFSTAR    : return psStringCopy ("PSFSTAR"   );
    case PM_SOURCE_MODE_SATSTAR    : return psStringCopy ("SATSTAR"   );
    case PM_SOURCE_MODE_BLEND      : return psStringCopy ("BLEND"     );
    case PM_SOURCE_MODE_EXTERNAL   : return psStringCopy ("EXTERNAL"  );
    case PM_SOURCE_MODE_BADPSF     : return psStringCopy ("BADPSF"    );
    case PM_SOURCE_MODE_DEFECT     : return psStringCopy ("DEFECT"    );
    case PM_SOURCE_MODE_SATURATED  : return psStringCopy ("SATURATED" );
    case PM_SOURCE_MODE_CR_LIMIT   : return psStringCopy ("CRLIMIT"   );
    case PM_SOURCE_MODE_EXT_LIMIT  : return psStringCopy ("EXTLIMIT"  );
    case PM_SOURCE_MODE_SUBTRACTED : return psStringCopy ("SUBTRACTED");
    default:
      return NULL;
  }
  return NULL;
}
