/** @file  pmSourcePhotometry.c
 *
 *  @author EAM, IfA; GLG, MHPCC
 *
 *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2008-02-03 21:58:56 $
 *
 *  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 "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,
                         psMaskType maskVal, psMaskType mark)
{
    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
    model = pmSourceGetModel (&isPSF, source);
    if (model == NULL) {
	psTrace ("psModules.objects", 3, "fail mag : no valid model");
        return false;
    }

    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
    if (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);
	if (!isfinite(fluxScale)) {
	    // XXX objects on the edge can be slightly outside -- if we get an 
	    // error, use the full fit.
	    psErrorClear();
	    status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
	} else {
	    if (isfinite(fluxScale) && (fluxScale > 0.0)) {
		source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]);
	    } else {
		source->psfMag = NAN;
	    }
	}
    } else {
	status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
    }

    // measure EXT model photometry (do both modelEXT and modelConv or just the one?)
    status = pmSourcePhotometryModel (&source->extMag, 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;
    }

    if (!isfinite(SN) || (SN < AP_MIN_SN)) {
	psTrace ("psModules.objects", 3, "fail mag : bad SN: %f (limit: %f)", SN, AP_MIN_SN);
        return false;
    }

    // replace source flux
    // XXX need to be certain which model and size of mask for prior subtraction
    // XXX full model or just analytical?
    // XXX use pmSourceAdd instead?
    if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
        pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
    }

    // 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)
    // psImageShift ();
    psImage *flux;
    if (mode & PM_SOURCE_PHOT_INTERP) {
        float dx = 0.5 - x + (int)x;
        float dy = 0.5 - y + (int)y;
        x += dx;
        y += dy;
        flux = psImageShift (NULL, source->pixels, dx, dy, 0.0, PS_INTERPOLATE_BICUBE);
        // flux = psImageShift (NULL, source->pixels, dx, dy, 0.0, PS_INTERPOLATE_BILINEAR);

        // XXX this is test code to verify the shift is doing the right thing (seems to be)
        # if (0)
            // measure centroid of unshifted gaussian (should be 16.0,16.0)
        {
            psImage *image = source->pixels;
            float xo = 0.0;
            float yo = 0.0;
            float xo2 = 0.0;
            float yo2 = 0.0;
            float no = 0.0;
            for (int j = 0; j < image->numRows; j++)
            {
                for (int i = 0; i < image->numCols; i++) {
                    xo += i*image->data.F32[j][i];
                    yo += j*image->data.F32[j][i];
                    xo2 += i*i*image->data.F32[j][i];
                    yo2 += j*j*image->data.F32[j][i];
                    no += image->data.F32[j][i];
                }
            }
            xo /= no;
            yo /= no;
            xo2 = sqrt (xo2/no - xo*xo);
            yo2 = sqrt (yo2/no - yo*yo);
            fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
        }

        // measure centroid of unshifted gaussian (should be 16.0,16.0)
        {
            psImage *image = flux;
            float xo = 0.0;
            float yo = 0.0;
            float xo2 = 0.0;
            float yo2 = 0.0;
            float no = 0.0;
            for (int j = 0; j < image->numRows; j++)
            {
                for (int i = 0; i < image->numCols; i++) {
                    xo += i*image->data.F32[j][i];
                    yo += j*image->data.F32[j][i];
                    xo2 += i*i*image->data.F32[j][i];
                    yo2 += j*j*image->data.F32[j][i];
                    no += image->data.F32[j][i];
                }
            }
            xo /= no;
            yo /= no;
            xo2 = sqrt (xo2/no - xo*xo);
            yo2 = sqrt (yo2/no - yo*yo);
            fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
        }
        # endif

    } else {
        flux = source->pixels;
    }
    // set aperture mask circle to model radius
    // XXX use a different radius for apertures and fits...
    // XXX can I remove this?  the source should have the mask defined when it is constructed or
    // when the fit / aperture radius is changed...
    psImageKeepCircle (source->maskObj, x, y, model->radiusFit, "OR", mark);

    // measure the weight of included pixels
    // XXX is this supposed to use the weight or the flux?
    if (mode & PM_SOURCE_PHOT_WEIGHT) {
        pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->maskObj, maskVal);
    }

    // measure object aperture photometry
    status = pmSourcePhotometryAper  (&source->apMag, model, flux, source->maskObj, maskVal);
    if (!status) {
	psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag");
        psErrorCode last = psErrorCodeLast();
	// XXX flag the object?
        if (last == PM_ERR_PHOTOM) {
            // the aper mag was undefined (flux < 0).  bad object, but
            // don't keep the error on the stack
            psErrorClear();
        }
    }

    // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
    if (status && isPSF && psf) {
        if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
            source->apMag += pmGrowthCurveCorrect (psf->growth, model->radiusFit);
        }
        if (mode & PM_SOURCE_PHOT_APCORR) {
            rflux   = pow (10.0, 0.4*source->psfMag);
            source->apMag  -= PS_SQR(model->radiusFit)*rflux * psf->skyBias + psf->skySat / rflux;
        }
    }
    if (mode & PM_SOURCE_PHOT_INTERP) {
        psFree (flux);
    }

    // unmask aperture
    // XXX can I remove this?  this will probably break things downstream...
    psImageKeepCircle (source->maskObj, x, y, model->radiusFit, "AND", PS_NOT_U8(mark));

    // if source was originally subtracted, re-subtract object, leave local sky
    // XXX replace with pmSourceSub...
    if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
        pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
    }

    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, pmModel *model)
{
    PS_ASSERT_PTR_NON_NULL(fitMag, false);
    if (model == NULL) {
        return false;
    }

    float fitSum = 0;
    *fitMag = NAN;

    // measure fitMag
    fitSum = model->modelFlux (model->params);
    if (fitSum <= 0)
        return false;
    if (!isfinite(fitSum))
        return false;
    *fitMag = -2.5*log10(fitSum);

    return (true);
}

// return source aperture magnitude
bool pmSourcePhotometryAper (float *apMag, pmModel *model, psImage *image, psImage *mask,
                             psMaskType 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;
    psU8 **mkData = mask->data.U8;

    // 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) {
        psError(PM_ERR_PHOTOM, true, "apSum is negative");
        return false;
    }

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

// return source aperture magnitude
bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *image, psImage *mask,
                          psMaskType maskVal)
{
    PS_ASSERT_PTR_NON_NULL(pixWeight, false);
    PS_ASSERT_PTR_NON_NULL(image, 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 - image->col0;
    dP = image->numCols - dX;
    DX = PS_MAX(dX, dP);
    NX = image->numCols;

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

    // measure modelSum and validSum
    // XXX this does not work for sources near the edge: we need to measure for
    // a full square region
    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.U8[my][mx] & maskVal)
                continue;

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

    if (validSum <= 0)
        return false;

    *pixWeight = validSum / modelSum;
    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->weight;
    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.U8[yi][xi])
                continue;
            if (Tj->data.U8[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->weight;
    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.U8[yi][xi])
                continue;
            if (Tj->data.U8[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->weight;
    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.U8[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

bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight,
                    psMaskType maskVal)
{
    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(weight, 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.U8[j][i] & maskVal)
                continue;
            if (weight->data.F32[j][i] <= 0)
                continue;
            dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i];
            Npix ++;
        }
    }
    model->nDOF = Npix - 1;
    model->chisq = dC;

    return (true);
}


double pmSourceModelWeight(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->modelFlux;
    assert (Pi != NULL);
    const psImage *Wi = Mi->weight;
    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.U8[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;
            }
        }
    }
    return flux;
}

double pmSourceModelDotModel (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->modelFlux;
    assert (Pi != NULL);
    const psImage *Pj = Mj->modelFlux;
    assert (Pj != NULL);

    const psImage *Wi = Mi->weight;
    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.U8[yi][xi])
                continue;
            if (Tj->data.U8[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 = 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) // 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->modelFlux;
    assert (Pj != NULL);

    const psImage *Wi = Mi->weight;
    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.U8[yi][xi])
                continue;
            if (Tj->data.U8[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 = Wi->data.F32[yi][xi];
                if (wt > 0) {
                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
                }
            }
        }
    }
    return flux;
}

