# include "psphotInternal.h"
# include <gsl/gsl_sf_gamma.h>

typedef struct {
    psImageMaskType maskVal;
    psImageMaskType markVal;
    psImageMaskType crMask;
    float ApResid;
    float ApSysErr;
    float nSigmaApResid;
    float nSigmaMoments;
    float nSigmaCR;
    float soft;
    int grow;
} psphotSourceSizeOptions;

bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf);
bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options);
bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);

// we need to call this function after sources have been fitted to the PSF model and
// subtracted.  To determine the CR-nature, this function examines the 9 pixels in the 3x3
// square containing the peak and compares the observed flux to the model.  To determine
// the EXT-nature, this function measures the amount of positive or negative total
// deviation from the psf model at the r = FWHM/2 position

// XXX use an internal flag to mark sources which have already been measured
bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, long first)
{
    bool status;
    psphotSourceSizeOptions options;

    psTimerStart ("psphot.size");

    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    options.maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    assert (options.maskVal);

    options.markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels
    assert (options.markVal);

    // bit to mask the cosmic-ray pixels
    options.crMask  = pmConfigMaskGet("CR", config); // Mask value for cosmic rays

    options.nSigmaCR = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT");
    assert (status);

    // XXX recipe name is not great
    options.nSigmaApResid = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT");
    assert (status);

    // XXX recipe name is not great
    options.nSigmaMoments = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.MOMENTS");
    assert (status);

    options.grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs
    if (!status || options.grow < 0) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CR.GROW is not positive.");
        return false;
    }

    options.soft = psMetadataLookupF32(&status, recipe, "PSPHOT.CR.NSIGMA.SOFTEN"); // Softening parameter
    if (!status || !isfinite(options.soft) || options.soft < 0.0) {
        psWarning("PSPHOT.CR.NSIGMA.SOFTEN not set; defaulting to zero.");
        options.soft = 0.0;
    }

    // We are using the value psfMag - 2.5*log10(moment.Sum) as a measure of the extendedness
    // of and object.  We need to model this distribution for the PSF stars before we can test
    // the significance for a specific object
    // XXX move this to the code that generates the PSF?
    // XXX store the results on pmPSF?
    psphotSourceSizePSF (&options, sources, psf);

    // classify the sources based on ApResid and Moments (extended sources)
    psphotSourceClass(readout, sources, recipe, psf, &options);

    psphotSourceSizeCR (readout, sources, &options);

    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n - first, psTimerMark ("psphot.size"));

    psphotVisualPlotSourceSize (recipe, sources);
    psphotVisualShowSourceSize (readout, sources);
    psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr);

    return true;
}

bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {

    // replace the source flux
    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);

    // flag this as a CR
    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    pmPeak *peak = source->peak;
    psAssert (peak, "NULL peak");

    // grab the matching footprint
    pmFootprint *footprint = peak->footprint;
    if (!footprint) {
	// if we have not footprint, use the old code to mask by isophot
	psphotMaskCosmicRayIsophot (source, maskVal, crMask);
	return true;
    }

    if (!footprint->spans) {
	// if we have no footprint, use the old code to mask by isophot
	psphotMaskCosmicRayIsophot (source, maskVal, crMask);
	return true;
    }

    // mask all of the pixels covered by the spans of the footprint
    for (int j = 1; j < footprint->spans->n; j++) {
	pmSpan *span1 = footprint->spans->data[j];

	int iy = span1->y;
	int xs = span1->x0;
	int xe = span1->x1;

	for (int ix = xs; ix < xe; ix++) {
	    mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
	}
    }
    return true;
}

bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {

    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    pmPeak *peak = source->peak;
    psAssert (peak, "NULL peak");

    psImage *mask   = source->maskView;
    psImage *pixels = source->pixels;
    psImage *variance = source->variance;

    // XXX This should be a recipe variable
# define SN_LIMIT 5.0

    int xo = peak->x - pixels->col0;
    int yo = peak->y - pixels->row0;

    // mark the pixels in this row to the left, then the right
    for (int ix = xo; ix >= 0; ix--) {
	float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
	if (SN > SN_LIMIT) {
	    mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
	}
    }
    for (int ix = xo + 1; ix < pixels->numCols; ix++) {
	float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
	if (SN > SN_LIMIT) {
	    mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
	}
    }

    // for each of the neighboring rows, mark the high pixels if they have a marked neighbor
    // first go up:
    for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
	// mark the pixels in this row to the left, then the right
	for (int ix = 0; ix < pixels->numCols; ix++) {
	    float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
	    if (SN < SN_LIMIT) continue;

	    bool valid = false;
	    valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
	    valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
	    valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;

	    if (!valid) continue;
	    mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
	}
    }
    // next go down:
    for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
	// mark the pixels in this row to the left, then the right
	for (int ix = 0; ix < pixels->numCols; ix++) {
	    float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
	    if (SN < SN_LIMIT) continue;

	    bool valid = false;
	    valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
	    valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
	    valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;

	    if (!valid) continue;
	    mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
	}
    }
    return true;
}

// model the apmifit distribution for the psf stars:
bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) {

    // select stats from the psf stars
    psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32);
    psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
    
    psImageMaskType maskVal = options->maskVal | options->markVal;

    // XXX  why PHOT_WEIGHT??
    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;

    for (int i = 0; i < sources->n; i++) {
	pmSource *source = sources->data[i];
	if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;

        // replace object in image
        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
            pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
        }

	// clear the mask bit and set the circular mask pixels
	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
	psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);

	// XXX can we test if psfMag is set and calculate only if needed?
	pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
	
	// clear the mask bit 
	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));

        // re-subtract the object, leave local sky
        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);

	float apMag = -2.5*log10(source->moments->Sum);
	float dMag = source->psfMag - apMag;
	
	psVectorAppend (Ap, 100, dMag);
	psVectorAppend (ApErr, 100, source->errMag);
    }

    // model the distribution as a mean or median value and a systematic error from that value:
    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    psVectorStats (stats, Ap, NULL, NULL, 0);

    psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32);
    for (int i = 0; i < Ap->n; i++) {
	dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian;
    }

    options->ApResid = stats->robustMedian;
    options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05);
    psLogMsg ("psphot", PS_LOG_DETAIL, "psf - Sum: %f +/- %f\n", options->ApResid, options->ApSysErr);

    psFree (Ap);
    psFree (ApErr);
    psFree (stats);
    psFree (dAp);

    return true;
}

// classify sources based on the combination of psf-mag, Mxx, Myy
bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {

    bool status;
    pmPSFClump psfClump;
    char regionName[64];

    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nmiss", "Nskip");

    int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
    for (int i = 0; i < nRegions; i ++) {
	snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
	psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
	psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion");

	psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION");
	psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion");

	// pull FWHM_X,Y from the recipe, use to define psfClump.X,Y
	psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   psAssert (status, "missing PSF.CLUMP.X");
	psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   psAssert (status, "missing PSF.CLUMP.Y");
	psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  psAssert (status, "missing PSF.CLUMP.DX");
	psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  psAssert (status, "missing PSF.CLUMP.DY");

	if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) {
	    psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
	    continue;
	}
	
	if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, options)) {
	    psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
	    continue;
	}
    }	

    return true;
}

bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {

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

    int Nsat  = 0;
    int Next  = 0;
    int Npsf  = 0;
    int Ncr   = 0;
    int Nmiss = 0;
    int Nskip = 0;

    pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;

    psImageMaskType maskVal = options->maskVal | options->markVal;

    for (psS32 i = 0 ; i < sources->n ; i++) {

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

	// psfClumps are found for image subregions:
	// skip sources not in this region
	if (source->peak->x <  region->x0) continue;
	if (source->peak->x >= region->x1) continue;
	if (source->peak->y <  region->y0) continue;
	if (source->peak->y >= region->y1) continue;

	// skip source if it was already measured
	if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
	    psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
	    continue;
	}

	// source must have been subtracted
	if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
	    source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
	    psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
	    continue;
	}

	// we are basically classifying by moments; use the default if not found
	psAssert (source->moments, "why is this source missing moments?");
	if (source->mode & noMoments) { 
	    Nskip ++;
	    continue;
	}

	psF32 Mxx = source->moments->Mxx;
	psF32 Myy = source->moments->Myy;

        // replace object in image
        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
            pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
        }

	// clear the mask bit and set the circular mask pixels
	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
	psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);

	// XXX can we test if psfMag is set and calculate only if needed?
	pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal

	// clear the mask bit 
	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));

        // re-subtract the object, leave local sky
        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);

	float apMag = -2.5*log10(source->moments->Sum);
	float dMag = source->psfMag - apMag;
	float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);

	source->extNsigma = nSigma;
	source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;

	// Anything within this region is a probably PSF-like object. Saturated stars may land
	// in this region, but are detected elsewhere on the basis of their peak value.
	bool isPSF = (fabs(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY);
	if (isPSF) {
	    Npsf ++;
	    continue;
	}

	// Defects may not always match CRs from peak curvature analysis
	// Defects may also be marked as SATSTAR -- XXX deactivate this flag?
	// XXX this rule is not great
	if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
	    source->mode |= PM_SOURCE_MODE_DEFECT;
	    Ncr ++;
	    continue;
	}

	// saturated star (determined in PSF fit).  These may also be saturated galaxies, or
	// just large saturated regions.
	if (source->mode & PM_SOURCE_MODE_SATSTAR) { 
	    Nsat ++;
	    continue;
	}

	// XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
	bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y));
	if (isEXT) {
	    source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
	    Next ++;
	    continue;
	}

	psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma);
	Nmiss ++;
    }

    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip);

    return true;
}

// given the PSF ellipse parameters, navigate around the 1sigma contour, return the total
// deviation in sigmas.  This is measured on the residual image - should we ignore negative
// deviations?  NOTE: This function was an early attempt to classify extended objects, and is
// no longer used by psphot.
float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
			 psImageMaskType maskVal, const pmModel *model, float Ro)
{
    psF32 *PAR = model->params->data.F32; // Model parameters
    float sxx = PAR[PM_PAR_SXX], sxy = PAR[PM_PAR_SXY], syy = PAR[PM_PAR_SYY]; // Ellipse parameters

    // We treat the contour as an ellipse:
    // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY
    // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0;
    // This is a quadratic, Ay^2 + By + C with A = 1/SYY^2, B = x*SXY, C = (x / SXX)^2 - Ro
    // The solution is y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A], so:
    // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2]

    // min/max value of x is where B^2 - 4AC = 0; solve this for x
    float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0);
    if (Q < 0.0) {
	// ellipse is imaginary
	return NAN;
    }

    int radius = sqrtf(Q) + 0.5;        // Radius of ellipse
    int nPts = 0;                       // Number of points in ellipse
    float nSigma = 0.0;                 //

    for (int x = -radius; x <= radius; x++) {
	// Polynomial coefficients
	// XXX Should we be using the centre of the pixel as x or x+0.5?
	float A = PS_SQR (1.0 / syy);
	float B = x * sxy;
	float C = PS_SQR (x / sxx) - Ro;
	float T = PS_SQR(B) - 4*A*C;
	if (T < 0.0) {
	    continue;
	}

	// y position in source frame
	float yP = (-B + sqrt (T)) / (2.0 * A);
	float yM = (-B - sqrt (T)) / (2.0 * A);

	// Get the closest pixel positions (image frame)
	int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
	int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
	int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;

	if (xPix < 0 || xPix >= image->numCols) {
	    continue;
	}

	if (yPixM >= 0 && yPixM < image->numRows &&
	    !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
	    float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
	    nSigma += dSigma;
	    nPts++;
	}

	if (yPixM == yPixP) {
	    continue;
	}

	if (yPixP >= 0 && yPixP < image->numRows &&
	    !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
	    float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
	    nSigma += dSigma;
	    nPts++;
	}
    }
    nSigma /= nPts;
    return nSigma;
}

bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {

    // classify the sources based on the CR test (place this in a function?)
    // XXX use an internal flag to mark sources which have already been measured
    for (int i = 0; i < sources->n; i++) {
	pmSource *source = sources->data[i];

	// skip source if it was already measured
	if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
	    psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
	    continue;
	}

	// source must have been subtracted
	if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
	    source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
	    psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
	    continue;
	}

	psF32 **resid  = source->pixels->data.F32;
	psF32 **variance = source->variance->data.F32;
	psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;

	// Integer position of peak
	int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
	int yPeak = source->peak->yf - source->pixels->row0 + 0.5;

	// Skip sources which are too close to a boundary.  These are mostly caught as DEFECT
	if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
	    yPeak < 1 || yPeak > source->pixels->numRows - 2) {
	    psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
	    continue;
	}

	// Skip sources with masked pixels.  These are mostly caught as DEFECT
	bool keep = true;
	for (int iy = -1; (iy <= +1) && keep; iy++) {
	    for (int ix = -1; (ix <= +1) && keep; ix++) {
		if (mask[yPeak+iy][xPeak+ix] & options->maskVal) {
		    keep = false;
		}
	    }
	}
	if (!keep) {
	    psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
	    continue;
	}

	// Compare the central pixel with those on either side, for the four possible lines through it.

	// Soften variances (add systematic error)
	float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances

	// Across the middle: y = 0
	float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
	float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
	float nX = cX / sqrtf(dcX + softening);

	// Up the centre: x = 0
	float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
	float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
	float nY = cY / sqrtf(dcY + softening);

	// Diagonal: x = y
	float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
	float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
	float nL = cL / sqrtf(dcL + softening);

	// Diagonal: x = - y
	float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
	float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
	float nR = cR / sqrtf(dcR + softening);

	// P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
	// Ndof = 4 ? (four measurements, no free parameters)
	// XXX this value is going to be biased low because of systematic errors.
	// we need to calibrate it somehow
	// source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);

	// not strictly accurate: overcounts the chisq contribution from the center pixel (by
	// factor of 4); also biases a bit low if any pixels are masked
	// XXX I am not sure I want to keep this value...
	source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR);

	float fCR = 0.0;
	int nCR = 0;
	if (nX > 0.0) {
	    fCR += nX;
	    nCR ++;
	}
	if (nY > 0.0) {
	    fCR += nY;
	    nCR ++;
	}
	if (nL > 0.0) {
	    fCR += nL;
	    nCR ++;
	}
	if (nR > 0.0) {
	    fCR += nR;
	    nCR ++;
	}
	source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
	source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;

	if (!isfinite(source->crNsigma)) {
	    continue;
	}

	// this source is thought to be a cosmic ray.  flag the detection and mask the pixels
	if (source->crNsigma > options->nSigmaCR) {
	    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
	    // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
	    // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask);
	}
    }

    // now that we have masked pixels associated with CRs, we can grow the mask
    if (options->grow > 0) {
	bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask
	psImage *newMask = psImageConvolveMask(NULL, readout->mask, options->crMask, options->crMask, -options->grow, options->grow, -options->grow, options->grow);
	psImageConvolveSetThreads(oldThreads);
	if (!newMask) {
	    psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask");
	    return false;
	}
	psFree(readout->mask);
	readout->mask = newMask;
    }
    return true;
}
