/** @file  pmSourcePhotometry.c
 *
 *  @author EAM, IfA; GLG, MHPCC
 *
 *  @version $Revision: 1.50 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-02-16 22:28:54 $
 *
 *  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 "pmErrorCodes.h"
#include "pmHDU.h"
#include "pmFPA.h"
#include "pmFPAMaskWeight.h"
#include "pmSpan.h"
#include "pmFootprint.h"
#include "pmPeaks.h"
#include "pmMoments.h"
#include "pmGrowthCurve.h"
#include "pmResiduals.h"
#include "pmTrend2D.h"
#include "pmPSF.h"
#include "pmModel.h"
#include "pmSource.h"
#include "pmModelClass.h"
#include "pmSourcePhotometry.h"

# define DO_SKY 0

static float AP_MIN_SN = 0.0;

bool pmSourceMagnitudesInit (psMetadata *config)
{
    PS_ASSERT_PTR_NON_NULL(config, false);
    bool status;

    float limit = psMetadataLookupF32 (&status, config, "AP_MIN_SN");
    if (status) {
        AP_MIN_SN = limit;
    }
    return true;
}

/**
    this function is used to calculate the three defined source magnitudes:
    - apMag  : only if S/N > AP_MIN_SN
             : is optionally corrected for curve-of-growth if:
        - the source is a STAR (PSF)
        - the option is selected (mode & PM_SOURCE_PHOT_GROWTH)
    - psfMag : all sources with non-NULL modelPSF
             : is optionally corrected for aperture residual if:
        - the source is a STAR (PSF)
        - the option is selected (mode & PM_SOURCE_PHOT_APCORR)

    - extMag : all sources with non-NULL modelEXT
**/

// XXX masked region should be (optionally) elliptical
bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal)
{
    PS_ASSERT_PTR_NON_NULL(source, false);
    PS_ASSERT_PTR_NON_NULL(psf, false);

    int status = false;
    bool isPSF;
    float x, y;
    float rflux;
    float SN;
    pmModel *model;

    source->psfMag = NAN;
    source->extMag = NAN;
    source->errMag = NAN;
    source->apMag  = NAN;

    // we must have a valid model
    // XXX allow aperture magnitudes for sources without a model
    model = pmSourceGetModel (&isPSF, source);
    if (model == NULL) {
        psTrace ("psModules.objects", 3, "fail mag : no valid model");
        return false;
    }

    // XXX handle negative flux, low-significance
    if (model->dparams->data.F32[PM_PAR_I0] > 0) {
        SN = model->params->data.F32[PM_PAR_I0] / model->dparams->data.F32[PM_PAR_I0];
        source->errMag = 1.0 / SN;
    } else {
        SN = NAN;
        source->errMag = NAN;
    }
    x = model->params->data.F32[PM_PAR_XPOS];
    y = model->params->data.F32[PM_PAR_YPOS];

    // measure PSF model photometry
    // XXX TEST: do not use flux scale
    if (0 && psf->FluxScale) {
        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
        double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
	psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
	psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
	source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
	source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
	source->psfMag = -2.5*log10(source->psfFlux);
    } else {
        status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF);
	source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
    }

    // if we have a collection of model fits, check if one of them is a pointer to modelEXT
    if (source->modelFits) {
        bool foundEXT = false;
	for (int i = 0; i < source->modelFits->n; i++) {
	    pmModel *model = source->modelFits->data[i];
	    status = pmSourcePhotometryModel (&model->mag, NULL, model);
	    if (model == source->modelEXT) foundEXT = true;
	}
	if (foundEXT) {
	    source->extMag = source->modelEXT->mag;
	} else {
	    status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
	}
    } else {
	if (source->modelEXT) {
	    status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
	}
    }

    // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
    if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf && psf->ApTrend) {
        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
        double apTrend = pmTrend2DEval (psf->ApTrend, (float)source->peak->x, (float)source->peak->y);
        source->psfMag += apTrend;
    }

    // measure the contribution of included pixels
    if (mode & PM_SOURCE_PHOT_WEIGHT) {
        pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
    }

    // measure the contribution of included pixels
    if (mode & PM_SOURCE_PHOT_DIFFSTATS) {
        pmSourceMeasureDiffStats (source, maskVal);
    }

    // measure the aperture magnitude, if (SN > AP_MIN_SN)
    if (!isfinite(SN)) {
        psTrace ("psModules.objects", 3, "fail mag : bad SN: %f (limit: %f)", SN, AP_MIN_SN);
        return false;
    }

    // measure the aperture magnitude, if (SN > AP_MIN_SN)
    if (SN < AP_MIN_SN) {
        psTrace ("psModules.objects", 3, "skip ap mag : SN < limit : %f vs %f)", SN, AP_MIN_SN);
        return true;
    }

    // if we measure aperture magnitudes, the source must not currently be subtracted!
    psAssert (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED), "cannot measure ap mags if source is subtracted!");

    // if we are measuring aperture photometry and applying the growth correction,
    // we need to shift the flux in the selected pixels (but not the mask)
    psImage *flux = NULL, *mask = NULL; // Star flux and mask images, to photometer
    if (mode & PM_SOURCE_PHOT_INTERP) {
        float dx = 0.5 - x + (int)x;
        float dy = 0.5 - y + (int)y;
        x += dx;
        y += dy;

        if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy,
                              NAN, 0xff, PS_INTERPOLATE_BIQUADRATIC)) {
            // Not much we can do about it
            psErrorClear();
            psTrace ("psModules.objects", 3, "fail shift");
            return false;
        }

    } else {
        flux = source->pixels;
        mask = source->maskObj;
    }

    // measure object aperture photometry
    status = pmSourcePhotometryAper  (&source->apMag, model, flux, mask, maskVal);
    if (!status) {
        psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag");
    }

    // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
    // if the aper mag is NAN, the flux < 0.  this can happen for sources near the
    // detection limits (esp near bright neighbors)
    if (isfinite (source->apMag) && isPSF && psf) {
        if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
            source->apMag += pmGrowthCurveCorrect (psf->growth, source->apRadius);
        }
        if (mode & PM_SOURCE_PHOT_APCORR) {
	    // XXX this should be removed -- we no longer fit for the 'sky bias'
            rflux   = pow (10.0, 0.4*source->psfMag);
            source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux;
        }
    }
    if (mode & PM_SOURCE_PHOT_INTERP) {
        psFree(flux);
        psFree(mask);
    }

    return status;
}

/*
aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y)
(aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y)
(aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y)

*/

// return source model magnitude
bool pmSourcePhotometryModel (float *fitMag, float *fitFlux, pmModel *model)
{
    psAssert (fitMag || fitFlux, "at least one of magnitude or flux must be requested (not NULL)");
    if (model == NULL) return false;

    float mag  = NAN;
    float flux = NAN;

    // measure fitMag
    flux = model->modelFlux (model->params);
    if (flux > 0) {
	mag = -2.5*log10(flux);
    }
    if (fitMag) {
	*fitMag = mag;
    }
    if (fitFlux) {
	*fitFlux = flux;
    }

    if (flux <= 0) return false;
    if (!isfinite(flux)) return false;

    return (true);
}

// return source aperture magnitude
bool pmSourcePhotometryAper (float *apMag, pmModel *model, psImage *image, psImage *mask, psImageMaskType maskVal)
{
    PS_ASSERT_PTR_NON_NULL(apMag, false);
    PS_ASSERT_PTR_NON_NULL(image, false);
    PS_ASSERT_PTR_NON_NULL(mask, false);
    PS_ASSERT_PTR_NON_NULL(model, false);

    float apSum = 0;
    float sky = 0;
    *apMag = NAN;

    if (DO_SKY) {
        sky = model->params->data.F32[PM_PAR_SKY];
    } else {
        sky = 0;
    }

    psF32 **imData = image->data.F32;
    psImageMaskType **mkData = mask->data.PS_TYPE_IMAGE_MASK_DATA;

    // measure apMag
    for (int ix = 0; ix < image->numCols; ix++) {
        for (int iy = 0; iy < image->numRows; iy++) {
            if (mkData[iy][ix] & maskVal)
                continue;
            apSum += imData[iy][ix] - sky;
        }
    }
    if (apSum <= 0) {
        *apMag = NAN;
        return true;
    }

    *apMag = -2.5*log10(apSum);
    return true;
}

// return source aperture magnitude
bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *mask, psImageMaskType maskVal)
{
    PS_ASSERT_PTR_NON_NULL(pixWeight, false);
    PS_ASSERT_PTR_NON_NULL(mask, false);
    PS_ASSERT_PTR_NON_NULL(model, false);

    float modelSum = 0;
    float validSum = 0;
    float sky = 0;
    float value;

    int Xo, Yo, dP;
    int dX, DX, NX;
    int dY, DY, NY;

    *pixWeight = 0.0;

    // we only care about the value of the object model, not the local sky
    if (DO_SKY) {
        sky = model->params->data.F32[PM_PAR_SKY];
    } else {
        sky = 0;
    }

    // the model function returns the source flux at a position
    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);

    psVector *params = model->params;

    Xo = params->data.F32[PM_PAR_XPOS];
    Yo = params->data.F32[PM_PAR_YPOS];

    dX = Xo - mask->col0;
    dP = mask->numCols - dX;
    DX = PS_MAX(dX, dP);
    NX = mask->numCols;

    dY = Yo - mask->row0;
    dP = mask->numRows - dY;
    DY = PS_MAX(dY, dP);
    NY = mask->numRows;

    // measure modelSum and validSum.  this function is applied to a sources' subimage.  the
    // value of DX is chosen (see above) to cover the full possible size of the subimage if it
    // were not by an edge; ie, if the source is cut in half by an image edge, we correctly
    // count the virtual pixels off the edge in normalizing the value of the pixWeight
    for (int ix = -DX; ix < DX + 1; ix++) {
        int mx = ix + dX;
        for (int iy = -DY; iy < DY + 1; iy++) {
            int my = iy + dY;

            coord->data.F32[0] = (psF32) (ix + Xo);
            coord->data.F32[1] = (psF32) (iy + Yo);

            // for the full model, add all points
            value = model->modelFunc (NULL, params, coord) - sky;
            modelSum += value;

            // include count only the unmasked pixels within the image area
            if (mx < 0)
                continue;
            if (my < 0)
                continue;
            if (mx >= NX)
                continue;
            if (my >= NY)
                continue;
            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskVal)
                continue;

            validSum += value;
        }
    }
    psFree (coord);

    if (validSum <= 0)
        return false;

    *pixWeight = validSum / modelSum;
    return (true);
}

# define FLUX_LIMIT 3.0

// return source aperture magnitude
bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal)
{
    PS_ASSERT_PTR_NON_NULL(source, false);

    if (source->diffStats == NULL) {
	source->diffStats = pmSourceDiffStatsAlloc();
    }

    float fGood = 0.0;
    float fBad  = 0.0;
    int   nGood = 0;
    int   nMask = 0;
    int   nBad  = 0;
    
    psImage *flux     = source->pixels;
    psImage *variance = source->variance;
    psImage *mask     = source->maskObj;

    for (int iy = 0; iy < flux->numRows; iy++) {
	for (int ix = 0; ix < flux->numCols; ix++) {
            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) {
		nMask ++;
                continue;
	    }

	    float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);

	    if (SN > +FLUX_LIMIT) { 
		nGood ++;
		fGood += fabs(flux->data.F32[iy][ix]);
	    }

	    if (SN < -FLUX_LIMIT) { 
		nBad ++;
		fBad += fabs(flux->data.F32[iy][ix]);
	    }
	}
    }

    source->diffStats->nGood      = nGood;
    source->diffStats->fRatio     = (fGood + fBad         == 0.0) ? NAN : fGood / (fGood + fBad);	   
    source->diffStats->nRatioBad  = (nGood + nBad         == 0)   ? NAN : nGood / (float) (nGood + nBad);	   
    source->diffStats->nRatioMask = (nGood + nMask        == 0)   ? NAN : nGood / (float) (nGood + nMask);	   
    source->diffStats->nRatioAll  = (nGood + nMask + nBad == 0)   ? NAN : nGood / (float) (nGood + nMask + nBad);

    return (true);
}

# if (0)
double pmSourceCrossProduct (const pmSource *Mi,
                             const pmSource *Mj,
                             const bool unweighted_sum) // should the cross product be weighted?
{
    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    PS_ASSERT_PTR_NON_NULL(Mj, NAN);

    int Xs, Xe, Ys, Ye;
    int xi, xj, yi, yj;
    int xIs, xJs, yIs, yJs;
    int xIe, yIe;
    double flux, wt;

    const psImage *Pi = Mi->pixels;
    assert (Pi != NULL);
    const psImage *Pj = Mj->pixels;
    assert (Pj != NULL);

    const psImage *Wi = Mi->variance;
    if (!unweighted_sum) {
        assert (Wi != NULL);
    }

    const psImage *Ti = Mi->maskObj;
    assert (Ti != NULL);
    const psImage *Tj = Mj->maskObj;
    assert (Tj != NULL);

    Xs = PS_MAX (Pi->col0, Pj->col0);
    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);

    Ys = PS_MAX (Pi->row0, Pj->row0);
    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);

    xIs = Xs - Pi->col0;
    xJs = Xs - Pj->col0;
    yIs = Ys - Pi->row0;
    yJs = Ys - Pj->row0;

    xIe = Xe - Pi->col0;
    yIe = Ye - Pi->row0;

    // note that this is addressing the same image pixels,
    // though only if both are source not model images
    flux = 0;
    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
                continue;
            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj])
                continue;

            if (unweighted_sum) {
                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
            } else {
                wt = Wi->data.F32[yi][xi];
                if (wt > 0) {
                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
                }
            }
        }
    }
    return flux;
}

double pmSourceCrossWeight(const pmSource *Mi,
                           const pmSource *Mj,
                           const bool unweighted_sum) // should the cross product be weighted?
{
    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    PS_ASSERT_PTR_NON_NULL(Mj, NAN);

    int Xs, Xe, Ys, Ye;
    int xi, xj, yi, yj;
    int xIs, xJs, yIs, yJs;
    int xIe, yIe;
    double flux, wt;

    const psImage *Pi = Mi->pixels;
    assert (Pi != NULL);
    const psImage *Pj = Mj->pixels;
    assert (Pj != NULL);

    const psImage *Wi = Mi->variance;
    if (!unweighted_sum) {
        assert (Wi != NULL);
    }

    const psImage *Ti = Mi->maskObj;
    assert (Ti != NULL);
    const psImage *Tj = Mj->maskObj;
    assert (Tj != NULL);

    Xs = PS_MAX (Pi->col0, Pj->col0);
    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);

    Ys = PS_MAX (Pi->row0, Pj->row0);
    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);

    xIs = Xs - Pi->col0;
    xJs = Xs - Pj->col0;
    yIs = Ys - Pi->row0;
    yJs = Ys - Pj->row0;

    xIe = Xe - Pi->col0;
    yIe = Ye - Pi->row0;

    // note that this is addressing the same image pixels,
    // though only if both are source not model images
    flux = 0;
    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
                continue;
            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj])
                continue;

            if (unweighted_sum) {
                flux++;
            } else {
                wt = Wi->data.F32[yi][xi];
                if (wt > 0) {
                    flux += 1.0 / wt;
                }
            }
        }
    }
    return flux;
}

double pmSourceWeight(const pmSource *Mi,
                      int term,
                      const bool unweighted_sum) // should the cross product be weighted?
{
    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    double flux = 0, wt = 0, factor = 0;

    const psImage *Pi = Mi->pixels;
    assert (Pi != NULL);
    const psImage *Wi = Mi->variance;
    if (!unweighted_sum) {
        assert (Wi != NULL);
    }
    const psImage *Ti = Mi->maskObj;
    assert (Ti != NULL);

    // note that this is addressing the same image pixels,
    // though only if both are source not model images
    for (int yi = 0; yi < Pi->numRows; yi++) {
        for (int xi = 0; xi < Pi->numCols; xi++) {
            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
                continue;
            if (!unweighted_sum) {
                wt = Wi->data.F32[yi][xi];
                if (wt == 0)
                    continue;
            }

            switch (term) {
            case 0:
                factor = 1;
                break;
            case 1:
                factor = xi + Pi->col0;
                break;
            case 2:
                factor = yi + Pi->row0;
                break;
            default:
                psAbort("invalid term for pmSourceWeight");
            }

            if (unweighted_sum) {
                flux += (factor * Pi->data.F32[yi][xi]);
            } else {
                flux += (factor * Pi->data.F32[yi][xi]) / wt;
            }
            // fprintf (stderr, "Pi: %f, flux: %f\n", Pi->data.F32[yi][xi], flux);
        }
    }
    return flux;
}
# endif

// determine chisq, etc for linear normalization-only fit
bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal, const float covarFactor)
{
    PS_ASSERT_PTR_NON_NULL(model, false);
    PS_ASSERT_PTR_NON_NULL(image, false);
    PS_ASSERT_PTR_NON_NULL(mask, false);
    PS_ASSERT_PTR_NON_NULL(variance, false);

    double dC = 0.0;
    int Npix = 0;
    for (int j = 0; j < image->numRows; j++) {
        for (int i = 0; i < image->numCols; i++) {
            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] & maskVal)
                continue;
            if (variance->data.F32[j][i] <= 0)
                continue;
            dC += PS_SQR (image->data.F32[j][i]) / (covarFactor * variance->data.F32[j][i]);
            Npix ++;
        }
    }
    model->nPix = Npix; 
    model->nDOF = Npix - 1;
    model->chisq = dC;

    return (true);
}


double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor)
{
    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    double flux = 0, wt = 0, factor = 0;

    const psImage *Pi = Mi->modelFlux;
    assert (Pi != NULL);
    const psImage *Wi = Mi->variance;
    if (!unweighted_sum) {
        assert (Wi != NULL);
    }
    const psImage *Ti = Mi->maskObj;
    assert (Ti != NULL);

    for (int yi = 0; yi < Pi->numRows; yi++) {
        for (int xi = 0; xi < Pi->numCols; xi++) {
            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
                continue;
            if (!unweighted_sum) {
                wt = covarFactor * Wi->data.F32[yi][xi];
                if (wt == 0)
                    continue;
            }

            switch (term) {
            case 0:
                factor = 1;
                break;
            case 1:
                factor = xi + Pi->col0;
                break;
            case 2:
                factor = yi + Pi->row0;
                break;
            default:
                psAbort("invalid term for pmSourceWeight");
            }

            if (unweighted_sum) {
                flux += (factor * Pi->data.F32[yi][xi]);
            } else {
                flux += (factor * Pi->data.F32[yi][xi]) / wt;
            }
        }
    }
    return flux;
}

double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
{
    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    PS_ASSERT_PTR_NON_NULL(Mj, NAN);
    int Xs, Xe, Ys, Ye;
    int xi, xj, yi, yj;
    int xIs, xJs, yIs, yJs;
    int xIe, yIe;
    double flux, wt;

    const psImage *Pi = Mi->modelFlux;
    assert (Pi != NULL);
    const psImage *Pj = Mj->modelFlux;
    assert (Pj != NULL);

    const psImage *Wi = Mi->variance;
    if (!unweighted_sum) {
        assert (Wi != NULL);
    }

    const psImage *Ti = Mi->maskObj;
    assert (Ti != NULL);
    const psImage *Tj = Mj->maskObj;
    assert (Tj != NULL);

    Xs = PS_MAX (Pi->col0, Pj->col0);
    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);

    Ys = PS_MAX (Pi->row0, Pj->row0);
    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);

    xIs = Xs - Pi->col0;
    xJs = Xs - Pj->col0;
    yIs = Ys - Pi->row0;
    yJs = Ys - Pj->row0;

    xIe = Xe - Pi->col0;
    yIe = Ye - Pi->row0;

    // note that weight is addressing the same image pixels
    flux = 0;
    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
                continue;
            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj])
                continue;

            // XXX skip the nonsense weight pixels?
            if (unweighted_sum) {
                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
            } else {
                wt = covarFactor * Wi->data.F32[yi][xi];
                if (wt > 0) {
                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
                }
            }
        }
    }
    return flux;
}

double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
{
    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    PS_ASSERT_PTR_NON_NULL(Mj, NAN);
    int Xs, Xe, Ys, Ye;
    int xi, xj, yi, yj;
    int xIs, xJs, yIs, yJs;
    int xIe, yIe;
    double flux, wt;

    const psImage *Pi = Mi->pixels;
    assert (Pi != NULL);
    const psImage *Pj = Mj->modelFlux;
    assert (Pj != NULL);

    const psImage *Wi = Mi->variance;
    if (!unweighted_sum) {
        assert (Wi != NULL);
    }

    const psImage *Ti = Mi->maskObj;
    assert (Ti != NULL);
    const psImage *Tj = Mj->maskObj;
    assert (Tj != NULL);

    Xs = PS_MAX (Pi->col0, Pj->col0);
    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);

    Ys = PS_MAX (Pi->row0, Pj->row0);
    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);

    xIs = Xs - Pi->col0;
    xJs = Xs - Pj->col0;
    yIs = Ys - Pi->row0;
    yJs = Ys - Pj->row0;

    xIe = Xe - Pi->col0;
    yIe = Ye - Pi->row0;

    // note that weight is addressing the same image pixels,
    flux = 0;
    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
                continue;
            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj])
                continue;

            // XXX skip the nonsense weight pixels?
            if (unweighted_sum) {
                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
            } else {
                wt = covarFactor * Wi->data.F32[yi][xi];
                if (wt > 0) {
                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
                }
            }
        }
    }
    return flux;
}
