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

// we need to call this function after sources have been fitted to the PSF model and
// subtracted.  this function examines the 9 pixels in the 3x3 square containing the peak and
// compares the observed flux to the model.

bool psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe) {

    // loop over all source
    for (int i = 0; i < sources->n; i++) {
	pmSource *source = sources->data[i];

	// source must have been subtracted
	if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;

	psF32 **resid = source->pixels->data.F32;
	psF32 **weight = source->weight->data.F32;
	psU8 **mask = source->maskObj->data.U8;

	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
	if (xPeak < 1) continue;
	if (xPeak > source->pixels->numCols - 2) continue;
	if (yPeak < 1) continue;
	if (yPeak > source->pixels->numRows - 2) continue;

	// XXX for now, just skip any sources with masked pixels
	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]) { keep &= false; }
	    }
	}
	if (!keep) continue;

	// XXX need to deal with edge peaks... and mask
	float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1];
	float cY = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak+0] - resid[yPeak+1][xPeak+0];
	float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1];
	float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1];

	float dcX = 4*weight[yPeak][xPeak] + weight[yPeak+0][xPeak-1] + weight[yPeak+0][xPeak+1];
	float dcY = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak+0] + weight[yPeak+1][xPeak+0];
	float dcL = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak-1] + weight[yPeak+1][xPeak+1];
	float dcR = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak-1] + weight[yPeak-1][xPeak+1];

	float nX = cX / sqrt(dcX);
	float nY = cY / sqrt(dcY);
	float nL = cL / sqrt(dcL);
	float nR = cR / sqrt(dcR);

	// 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)
	source->psfChisq = PS_SQR (nX) + PS_SQR (nY) + PS_SQR (nX) + PS_SQR (nR);

	float fCR = 0.0;
	float fEXT = 0.0;
	int nCR = 0;
	int nEXT = 0;
	if (nX > 0.0) {
	    fCR += nX;
	    nCR ++;
	} else {
	    fEXT += nX;
	    nEXT ++;
	}
	if (nY > 0.0) {
	    fCR += nY;
	    nCR ++;
	} else {
	    fEXT += nY;
	    nEXT ++;
	}
	if (nL > 0.0) {
	    fCR += nL;
	    nCR ++;
	} else {
	    fEXT += nL;
	    nEXT ++;
	}
	if (nR > 0.0) {
	    fCR += nR;
	    nCR ++;
	} else {
	    fEXT += nR;
	    nEXT ++;
	}
	source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
	source->extNsigma = (nEXT > 0) ? fabs(fEXT) / nEXT : 0.0;
	// NOTE: abs needed to make the Nsigma value positive 
    }
    return true;
}


// some possible classes:
// all 4 curvatures are highly negative : CR
// all 4 curvatures are highly positive : EXT

// at least 2 are significantly negative, none are significantly positive : CR
// at least 2 are significantly positive, none are significantly negative : EXT

// any are significantly negative, some may be significantly positive : CR
// any are significantly positive, none may be significantly positive : EXT
	
/* Nn  Np  No
   4   0   0   CR_1
   3   1   0   CR_1
   3   0   1   CR_1
   2   2   0   CR_2
   2   1   1   CR_2
   2   0   2   CR_2
   1   3   0   CR_3
   1   2   1   CR_3
   1   1   2   CR_3
   1   0   3   CR_3
   0   4   0   EXT
   0   3   1   EXT
   0   2   2   EXT
   0   1   3   PSF
   0   0   4   PSF
*/

/* Alternatively, write a f(CR) = Sum(nX,etc if >0) */
	  

/* I can write the formal probability that the 4 measurements are consistent with a PSF
 * based on how large the error values are (Nsigma -> erf -> P(Nsigma)).
 * I should examine this value as a function of flux and also examine the distribution of
 * the probabilities for the PSF sources.
 */

