/** @file  pmSource.c
 *
 *  Functions to define and manipulate sources on images
 *
 *  @author GLG, MHPCC
 *  @author EAM, IfA: significant modifications.
 *
 *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2007-04-26 01:20:29 $
 *
 *  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 <pslib.h>
#include "pmHDU.h"
#include "pmFPA.h"
#include "pmFPAMaskWeight.h"
#include "pmPeaks.h"
#include "pmMoments.h"
#include "pmResiduals.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->moments);
    psFree(tmp->modelPSF);
    psFree(tmp->modelEXT);
    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);

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

/******************************************************************************
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->peak = NULL;
    source->pixels = NULL;
    source->weight = NULL;
    source->maskObj = NULL;
    source->maskView = NULL;
    source->modelFlux = NULL;
    source->moments = NULL;
    source->blends = NULL;
    source->modelPSF = NULL;
    source->modelEXT = NULL;
    source->type = PM_SOURCE_TYPE_UNKNOWN;
    source->mode = PM_SOURCE_MODE_DEFAULT;
    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;

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

bool pmIsSource(const psPtr ptr)
{
    return (psMemGetDeallocator(ptr) == (psFreeFunc)sourceFree);
}

/******************************************************************************
pmSourceCopy(): copy the pmSource structure and contents
XXX : are we OK with incrementing the ID?
*****************************************************************************/
pmSource *pmSourceCopy(pmSource *in)
{
    // 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?
    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
    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 = psImageCopyView (NULL, in->maskView);

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

    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)
{
    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   = psMemIncrRefCounter(psImageSubset(readout->image, srcRegion));
    mySource->weight   = psMemIncrRefCounter(psImageSubset(readout->weight, srcRegion));
    mySource->maskView = psMemIncrRefCounter(psImageSubset(readout->mask,  srcRegion));
    mySource->region   = srcRegion;

    // the object mask is a copy, and used to define the source pixels
    mySource->maskObj = psImageCopy (NULL, mySource->maskView, PS_TYPE_MASK);

    return true;
}

bool pmSourceRedefinePixels(pmSource *mySource,
                            const pmReadout *readout,
                            psF32 x,
                            psF32 y,
                            psF32 Radius)
{
    bool extend;
    psRegion newRegion;

    if (Radius == 0)
        return false;

    // 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) {
        // Grab a new subimage
        psFree (mySource->pixels);
        psFree (mySource->weight);
        psFree (mySource->maskView);

        mySource->pixels   = psMemIncrRefCounter(psImageSubset(readout->image,  newRegion));
        mySource->weight   = psMemIncrRefCounter(psImageSubset(readout->weight, newRegion));
        mySource->maskView = psMemIncrRefCounter(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;
    }
    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 = pmFindImagePeaks (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)
{
    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;
            source->mode |= PM_SOURCE_MODE_DEFAULT;
            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, PM_MASK_SAT);

        // 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_DEFAULT;
            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_DEFAULT;
            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;
            source->mode |= PM_SOURCE_MODE_DEFAULT;
            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;
        source->mode |= PM_SOURCE_MODE_DEFAULT;
    }

    psStats *stats = NULL;
    stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    psVectorStats (stats, starsn_moments, NULL, NULL, 0);
    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);
    psVectorStats (stats, starsn_peaks, NULL, NULL, 0);
    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_PTR_NON_NULL(source->maskObj, 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 != NULL) && *vMsk) {
                vMsk++;
                continue;
            }

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

            // radius is just a function of (xDiff, yDiff)
            if (!VALID_RADIUS(xDiff, yDiff, R2)) {
                if (vMsk != NULL)
                    vMsk++;
                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 (vMsk != NULL)
                    vMsk++;
                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 (vMsk != NULL)
                vMsk++;
        }
    }

    // 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, "no valid pixels for source\n");
        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
    source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
    source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));

    psTrace ("psModules.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) {

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

// should we call pmSourceCacheModel if it does not exist?
bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add) {

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

	// 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]) continue;
		float value = Io*source->modelFlux->data.F32[iy][ix];
		if (add) {
		    source->pixels->data.F32[oy][ox] += value;
		} else {
		    source->pixels->data.F32[oy][ox] -= value;
		}
	    }
	}
	return true;
    } 

    if (add) {
	status = pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
    } else {
	status = pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
    }

    return true;
}

bool pmSourceAdd (pmSource *source, pmModelOpMode mode) {
    bool status = pmSourceOp (source, mode, true);
    return status;
}

bool pmSourceSub (pmSource *source, pmModelOpMode mode) {
    bool status = pmSourceOp (source, mode, false);
    return status;
}

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

    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;

    case PM_SOURCE_TYPE_EXTENDED:
        model = source->modelEXT;
        if (!model && source->modelPSF) {
            if (isPSF) {
                *isPSF = true;
	    }
	    return source->modelPSF;
        }
        return (model);
        break;

    default:
        return NULL;
    }
    return NULL;
}
