# include "psphotInternal.h"

bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool oldWindow);
bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, psImage *smoothedPixels);


bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule)
{
    bool status = true;

    fprintf (stdout, "\n");
    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Kron Iterate ---");

    // select the appropriate recipe information
    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    psAssert (recipe, "missing recipe?");

    int num = psphotFileruleCount(config, filerule);

    // skip the chisq image (optionally?)
    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
    if (!status) chisqNum = -1;

    // loop over the available readouts
    for (int i = 0; i < num; i++) {
        if (i == chisqNum) continue; // skip chisq image

        // find the currently selected readout
        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
        psAssert (file, "missing file?");

        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
        psAssert (readout, "missing readout?");

        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
        psAssert (detections, "missing detections?");

        psArray *sources = detections->newSources ? detections->newSources : detections->allSources;
        psAssert (sources, "missing sources?");

        pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
        // psAssert (psf, "missing psf?");

        if (!psphotKronIterateReadout (config, recipe, view, filerule, readout, sources, psf, i)) {
            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
            return false;
        }
    }
    return true;
}

int psphotKapaChannel (int channel);
bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel);
bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);

bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index) {

    bool status = false;

    if (!sources->n) {
        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping masked kron");
        return true;
    }

    psTimerStart ("psphot.kron");


    // determine the number of allowed threads
    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    if (!status) {
        nThreads = 0;
    }

    float RADIUS = psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS");
    if (!status) {
        RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    }

    float MIN_KRON_RADIUS = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_MIN_KRON");
    if (!status) {
        MIN_KRON_RADIUS = 0.25*RADIUS;
    }

    int KRON_ITERATIONS = psMetadataLookupS32 (&status, recipe, "KRON_ITERATIONS");
    if (!status) {
        KRON_ITERATIONS = 1;
    }

    bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT");
    if (!status) {
        KRON_APPLY_WEIGHT = true;
    }

    bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW");
    if (!status) {
        KRON_APPLY_WINDOW = false;
    }
    bool KRON_SMOOTH = psMetadataLookupBool (&status, recipe, "KRON_SMOOTH");
    if (!status) {
        KRON_SMOOTH = false;
    }
    float KRON_SMOOTH_SIGMA = psMetadataLookupF32 (&status, recipe, "KRON_SMOOTH_SIGMA");
    if (!status) {
        KRON_SMOOTH_SIGMA = 1.7;
    }
    float KRON_SMOOTH_NSIGMA = psMetadataLookupS32 (&status, recipe, "KRON_SMOOTH_NSIGMA");
    if (!status) {
        KRON_SMOOTH_NSIGMA = 2;
    }
    /*
     *  Parameter for calculating maximum integration radius based on source's surface
     *  brightness
     *  Given minimum surface brightness SBmin and a flux the maximum radius is found from
     *
     *  SBmin = source->flux / (pi * Rmax**2)
     *  Rmax = sqrt (source->flux ) / sqrt (SBmin * pi)
     *
     *  Now what do we use for SBmin?
     *
     *  SBmin = ( some flux ) / (some area)
     *  some flux ~ flux of Ns sigma PSF source
     *  some area ~ K times area of a PSF
     *  
     * flux of Ns sigma source ~ Ns * SKY_STDEV * PSF_EFFECTIVE_AREA
     * PSF_EFFECTIVE_AREA = 4 pi sigma_PSF^2
     * (the 4 accounts for the fact that 1 sigma is not the total area, it is
     * actually a larger region).
     *
     *   SBmin = Ns * SKY_STDEV * 4 * pi * sigma_PSF^2 / (K * pi * sigma_PSF^2) 
     *         = Ns * SKY_STDEV * 4 / K
     *
     * We combine the two parameters Ns and K, and the constant 4 into a single recipe value
     * KRON_SB_MIN_FACTOR with Ns = 5 and K = 2 the corresponding value 
     * for KRON_SB_MIN_FACTOR is 5 * 4 / 2 = 10
     */
    float KRON_SB_MIN_FACTOR = psMetadataLookupF32 (&status, recipe, "KRON_SB_MIN_FACTOR");
    if (!status) {
        KRON_SB_MIN_FACTOR = 10;
    }
    float SKY_STDEV = psMetadataLookupF32 (&status, readout->analysis, "MSKY_DEV");
    float KRON_SB_MIN_DIVISOR = sqrt ( M_PI * KRON_SB_MIN_FACTOR * SKY_STDEV );

    // bit-masks to test for good/bad pixels
    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    assert (maskVal);

    // bit-mask to mark pixels not used in analysis
    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
    assert (markVal);

    // maskVal is used to test for rejected pixels, and must include markVal
    maskVal |= markVal;

    // source analysis is done in S/N order (brightest first)
    sources = psArraySort (sources, pmSourceSortByFlux);
    if (!sources->n) {
        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
        return true;
    }

    // generate the window image: multiply the flux by this to downweight neighbors
    psImage *kronWindow = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
    psImageInit (kronWindow, 1.0);

    // start with the currently known moments (Mxx, Myy, Mxy) and generate a window image
    for (int i = 0; i < sources->n; i++) {

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

	// set a window function for each source based on the moments
	// (this skips really bad sources (no peak, no moments, DEFECT)
	psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW);
    }

    // We measure the Kron Radius on a smoothed copy of the readout image
    psImage *smoothedImage = NULL;
    if (KRON_SMOOTH) {
        // Build the smoothed source image
        // Replace the subtracted sources
        psphotReplaceAllSourcesReadout(config, view, filerule, index, recipe, false);
        // Copy the image and smooth
        psTimerStart ("psphot.kron.smooth");
        smoothedImage = psImageCopy(NULL, readout->image, PS_TYPE_F32);
        psImageSmooth(smoothedImage, KRON_SMOOTH_SIGMA, KRON_SMOOTH_NSIGMA);
        psLogMsg ("psphot.kron", PS_LOG_INFO, "smoothed image %f sec\n", psTimerMark ("psphot.kron.smooth"));

        // remove the sources
        psphotRemoveAllSourcesReadout( config, view, filerule, index, recipe, false );
        // Now subtract smooth versions of the sources from the smoothed image
        psTimerStart ("psphot.kron.smooth.sources");
        for (int i=0; i< sources->n; i++) {
            pmSource *source = sources->data[i];
            // If source has been subtracted from the readout image subtract a "smoothed" version from the smoothedImage
	    if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
                // cache copy of smoothedPixels in the source
                // tmpPtr is for use by a single "module" and must be null otherwise
                psAssert(source->tmpPtr == NULL, "source->tmpPtr is not null!");

                psImage *smoothedPixels = psImageSubset(smoothedImage, source->region);
                source->tmpPtr = (psPtr) smoothedPixels;
                pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0);
            }
        }
        psLogMsg ("psphot.kron", PS_LOG_INFO, "removed %ld smoothed sources %f sec\n", sources->n, psTimerMark ("psphot.kron.smooth.sources"));

    }

    // threaded measurement of the source magnitudes
    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
    int Cx = 1, Cy = 1;
    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);

    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);

    for (int i = 0; i < cellGroups->n; i++) {

        psArray *cells = cellGroups->data[i];

        for (int j = 0; j < cells->n; j++) {

            // allocate a job -- if threads are not defined, this just runs the job
            psThreadJob *job = psThreadJobAlloc ("PSPHOT_KRON_ITERATE");

            psArrayAdd(job->args, 1, readout);
            psArrayAdd(job->args, 1, kronWindow);
            psArrayAdd(job->args, 1, cells->data[j]); // sources
            psArrayAdd(job->args, 1, smoothedImage);
            PS_ARRAY_ADD_SCALAR(job->args, markVal,            PS_TYPE_IMAGE_MASK);
            PS_ARRAY_ADD_SCALAR(job->args, maskVal,            PS_TYPE_IMAGE_MASK);
            PS_ARRAY_ADD_SCALAR(job->args, RADIUS,             PS_TYPE_F32);
            PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS,    PS_TYPE_F32);
            PS_ARRAY_ADD_SCALAR(job->args, KRON_ITERATIONS,    PS_TYPE_S32);
            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WEIGHT, PS_TYPE_S32);
            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32);
            PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA,  PS_TYPE_F32);
            PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32);

// set this to 0 to run without threading
# if (1)
            if (!psThreadJobAddPending(job)) {
                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
                return false;
            }
# else
	    if (!psphotKronIterate_Threaded(job)) {
		psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
		// psFree(AnalysisRegion);
		return false;
	    }
	    psFree(job);
# endif
        }

        // wait for the threads to finish and manage results
        if (!psThreadPoolWait (false, true)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
            return false;
        }

        // we have only supplied one type of job, so we can assume the types here
        psThreadJob *job = NULL;
        while ((job = psThreadJobGetDone()) != NULL) {
            if (job->args->n < 1) {
                fprintf (stderr, "error with job\n");
            }
            psFree(job);
        }
    }
    psFree (cellGroups);
    psFree (kronWindow);
    if (KRON_SMOOTH) {
        for (int i = 0; i < sources->n; i++) {
            pmSource *source = sources->data[i];
            psFree(source->tmpPtr);
            source->tmpPtr = NULL;
        }
    }
    psFree (smoothedImage);

    psLogMsg ("psphot.kron", PS_LOG_WARN, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
    return true;
}

bool psphotKronIterate_Threaded (psThreadJob *job) {

    pmReadout *readout              = job->args->data[0];
    psImage *kronWindow             = job->args->data[1];
    psArray *sources                = job->args->data[2];
    psImage *smoothedImage          = job->args->data[3];
    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[6],F32);
    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[7],F32);
    int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[8],S32);
    bool KRON_APPLY_WEIGHT          = PS_SCALAR_VALUE(job->args->data[9],S32);
    bool KRON_APPLY_WINDOW          = PS_SCALAR_VALUE(job->args->data[10],S32);
    float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[11],F32);
    float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[12],F32);

    for (int j = 0; j < KRON_ITERATIONS; j++) {
	for (int i = 0; i < sources->n; i++) {

	    pmSource *source = sources->data[i];
	    if (!source->peak) continue; // XXX how can we have a peak-less source?

	    // check status of this source's moments
	    if (!source->moments) continue;
	    if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
	    if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;

	    // replace object in image
	    bool reSubtract = false;
            psImage *smoothedPixels = NULL;
	    if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
                pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
                smoothedPixels = (psImage *) source->tmpPtr;
                if (smoothedPixels) {
                    pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, true, maskVal, 0, 0);
                }
		reSubtract = true;
	    }

            // On first iteration set window radius to sky radius (if valid) on second iteration
            // use a factor times the previous radial moment value up to a maximum value that
            // depends on the surface brightness of the source
	    float maxWindow;
            if (j == 0) {
                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
            } else {
                if (KRON_SB_MIN_DIVISOR) {
                    if (isfinite(source->moments->KronFlux) && (source->moments->KronFlux > 0)) {
                        // Limit window radius based on surface brightness
                        float Rmax = sqrt(source->moments->KronFlux) / KRON_SB_MIN_DIVISOR;

                        if (source->moments->Mrf > 0) {
                            maxWindow = PS_MIN(6.0*source->moments->Mrf, Rmax);
                        } else {
                            maxWindow = Rmax;
                        }
                    } else {
                        maxWindow = RADIUS;
                    }
                } else {
                    // old code
                    maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS;
                }
            }
	    float windowRadius = PS_MAX(RADIUS, maxWindow);

	    // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
	    bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
	    psAssert (source->pixels, "WTF?");
            if (extend && smoothedPixels) {
                psFree(source->tmpPtr);
                smoothedPixels = psImageSubset(smoothedImage, source->region);
                psAssert (smoothedPixels, "WTF?");
                source->tmpPtr = (psPtr) smoothedPixels ;
            }


	    // clear the window function for this source based on the moments
            // Note: this function also applies cuts on the source and returns false if it
            // does not meet the requirements for measuring the Kron Radius or Magnitude.
            // Note: that it performs the cuts even if KRON_APPLY_WINDOW is false
	    if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) {

                // this function populates moments->Mrf,KronFlux,KronFluxErr
                psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, smoothedPixels);
                psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));

                // set a window function for each source based on the moments
                psphotKronWindowSetSource (source, kronWindow, true, true, KRON_APPLY_WINDOW);
            }

	    // if we subtracted it above, re-subtract the object, leave local sky
	    if (reSubtract) {
                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
                if (smoothedPixels) {
                    pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0);
                }
	    }
	}
    }

    return true;
}

bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal,
    bool applyWeight, psImage *smoothedPixels) {

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

    psF32 R2 = PS_SQR(radius);
    float rsigma2 =  applyWeight ? 0.5 / R2 : 0;

    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    // image coordinates.  the source->pixels image has an offset relative to its parent of
    // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in
    // this subimage.  we subtract off the peak coordinates, adjusted to this subimage, to have
    // minimal round-off error in the sums.  since these values are subtracted just to minimize
    // the dynamic range and are added back below, the exact value does not matter. these are
    // (int) so they can be used in the image index below.

    // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
    // Xn  = SUM (x - xc)^n * (z - sky)


    psF32 RF = 0.0;
    psF32 RS = 0.0;

    // the peak position is less accurate but less subject to extreme deviations
    float dX = source->moments->Mx - source->peak->xf;
    float dY = source->moments->My - source->peak->yf;
    float dR = hypot(dX, dY);
    float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
    float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;

    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
    // xCM, yCM from pixel coords to pixel index here.
    psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage

    int Xwo = source->pixels->col0;
    int Ywo = source->pixels->row0;

    psF32 **vPix;
    if (smoothedPixels) {
        vPix = smoothedPixels->data.F32;
    } else {
        vPix = source->pixels->data.F32;
    }
    psF32 **vWin = kronWindow->data.F32;
    psF32 **vWgt = source->variance->data.F32;
    
    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;

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

	psF32 yDiff = row - yCM;
	if (fabs(yDiff) > radius) continue;

	// coordinate of mirror pixel
	int yFlip = yCM - yDiff;
	if (yFlip < 0) continue;
	if (yFlip >= source->pixels->numRows) continue;

	for (psS32 col = 0; col < source->pixels->numCols ; col++) {
	    // check mask and value for this pixel
	    if (vMsk && (vMsk[row][col] & maskVal)) continue;
	    if (isnan(vPix[row][col])) continue;

	    psF32 xDiff = col - xCM;
	    if (fabs(xDiff) > radius) continue;

	    // coordinate of mirror pixel
	    int xFlip = xCM - xDiff;
	    if (xFlip < 0) continue;
	    if (xFlip >= source->pixels->numCols) continue;

	    // check mask and value for mirror pixel
	    if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
	    if (isnan(vPix[yFlip][xFlip])) continue;

	    // radius is just a function of (xDiff, yDiff)
	    psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
	    if (r2 > R2) continue;

	    // flux * window
	    float z = r2 * rsigma2;
	    assert (z >= 0.0);

	    // weight by window image and wide Gaussian
	    float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
	    float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);

	    float fDiff1 = vPix[row][col]*weight1;
	    float fDiff2 = vPix[yFlip][xFlip]*weight2;

	    float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));

	    // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
	    psF32 rf = pDiff * sqrt(r2);
	    psF32 rs = 0.5 * (fDiff1 + fDiff2);

	    RF  += rf;
	    RS  += rs;
	}
    }

    // Saturate the 1st radial moment
    float Mrf = MAX(minKronRadius, RF/RS);
    if (sqrt(source->peak->detValue) < 10.0) {
	Mrf = MIN (radius, Mrf);
    }

    // Calculate the Kron magnitude (make this block optional?)
    float radKron  = 2.5*Mrf;
    float radKron2 = radKron*radKron;

    int nKronPix = 0;
    float Sum = 0.0;
    float Var = 0.0;
    float Win = 0.0;

    // set vPix to the source pixels (it may have been set to the
    // smoothed image above)
    vPix = source->pixels->data.F32;


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

	psF32 yDiff = row - yCM;
	if (fabs(yDiff) > radKron) continue;

	for (psS32 col = 0; col < source->pixels->numCols ; col++) {
	    // check mask and value for this pixel
	    if (vMsk && (vMsk[row][col] & maskVal)) continue;
	    if (isnan(vPix[row][col])) continue;

	    psF32 xDiff = col - xCM;
	    if (fabs(xDiff) > radKron) continue;

	    // radKron is just a function of (xDiff, yDiff)
	    psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
	    if (r2 > radKron2) continue;

	    float weight1  = vWin[row+Ywo][col+Xwo];
	    float fDiff1 = vPix[row][col]*weight1;

	    float pDiff = fDiff1;
	    psF32 wDiff = vWgt[row][col] * weight1;

	    Sum += pDiff;
	    Var += wDiff;
	    Win += weight1;
	    nKronPix ++;
	}
    }

    source->moments->Mrf = Mrf;
    source->moments->KronFlux    = Sum;
    source->moments->KronFluxErr = sqrt(Var);

    return true;
}

bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool applyWindow) {

    if (!source) return false;
    if (!source->peak) return false; // XXX how can we have a peak-less source?
    if (!source->moments) return false;
    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
    if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) return false;
    if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false;
    psAssert(kronWindow, "need a window");

    if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;

    if (!applyWindow) {
        return true;
    }

    // we have a source with moments Mx, My, Mxx, Mxy, Myy.  we just need to define a Gaussian that has 
    // these values (and peak of 1.0)

    int Nx = kronWindow->numCols;
    int Ny = kronWindow->numRows;

    float Xo = source->moments->Mx;
    float Yo = source->moments->My;

    psEllipseMoments moments;
    moments.x2 = source->moments->Mxx;
    moments.y2 = source->moments->Myy;
    moments.xy = source->moments->Mxy;

    psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
    if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) {
        // Shall we log a proper warning? This happens often with matched sources (forced photometry)
        // fprintf(dump, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta);
        return false;
    }

    // Why this factor of 0.5 ?
    float scale = 0.5 * source->moments->Mrf / axes.major;
    axes.major *= scale;
    axes.minor *= scale;

    psEllipseShape shape = psEllipseAxesToShape(axes);
    if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) {
        // Shall we log a proper warning? This happens often with matched sources (forced photometry)
        // fprintf(dump, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy);
        return false;
    }

    float Smajor = axes.major;

    int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1);
    int maxX = PS_MIN(PS_MAX(Xo + 5*Smajor, 0), Nx - 1);
    int minY = PS_MIN(PS_MAX(Yo - 5*Smajor, 0), Ny - 1);
    int maxY = PS_MIN(PS_MAX(Yo + 5*Smajor, 0), Ny - 1);

    float rMxx = 0.5 / PS_SQR(shape.sx);
    float rMyy = 0.5 / PS_SQR(shape.sy);
    float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
                                    // implementation. XXX: Is this correct?

    for (int iy = minY; iy < maxY; iy++) {
	for (int ix = minX; ix < maxX; ix++) {

	    float dX = (ix + 0.5 - Xo);
	    float dY = (iy + 0.5 - Yo);

	    float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;

	    float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z));

            kronWindow->data.F32[iy][ix] *= f;
	}
    }

    return true;
}
