# include "psphotInternal.h"

// given a pmSource which has been fitted using modelPSF, evaluate the
// resulting fit: did the fit succeed? is this object PSF-like? is this object 
// extended?  return status is TRUE for a valid PSF, FALSE otherwise.

// identify objects consistent with PSF shape/magnitude distribution
// we expect dparams[PM_PAR_SXX],dparams[PM_PAR_SXX] to have a scatter of:
// sigma_x / (S/N)

// any objects which is consistent with the PSF should have 
// abs(dparams[PM_PAR_SXX]) < N * dsxLine(mag) & abs(dparams[PM_PAR_SYY]) < N * dsyLine(mag)
// this includes a minimum buffer (DS) for the brighter objects

// EAM : 2006.10.20 : I have re-defined params[PM_PAR_SXX] : it is now sqrt(2)*sigma_x, not 1/sigma_x

// saturated stars should fall outside (larger), but have peaks above SATURATION
// extended sources should be larger, cosmic rays smaller
// we also reject objects with S/N too low or ChiSquare to high

// floor for DS value 
// XXX EAM : add to configuration?

static float SATURATION;
static float PSF_MIN_SN;
static float PSF_MIN_DS;
static float PSF_MAX_CHI;
static float PSF_SHAPE_NSIGMA;
static float PSF_DSX_MEAN;
static float PSF_DSY_MEAN;
static float PSF_DSX_STDEV;
static float PSF_DSY_STDEV;

bool psphotInitLimitsPSF (psMetadata *recipe, pmReadout *readout) {

    bool status;

    // XXX do we need to set this differently from the value used to mark saturated pixels?
    pmCell *cell     = readout->parent;
    SATURATION       = psMetadataLookupF32 (&status, cell->concepts, "CELL.SATURATION");
    PSF_MIN_SN       = psMetadataLookupF32 (&status, recipe, "PSF_MIN_SN");
    PSF_MAX_CHI      = psMetadataLookupF32 (&status, recipe, "PSF_MAX_CHI");
    PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_SHAPE_NSIGMA");
    PSF_MIN_DS       = psMetadataLookupF32 (&status, recipe, "PSF_MIN_DS");
    if (!status) PSF_MIN_DS = 0.01;

    PSF_DSX_MEAN     = psMetadataLookupF32 (&status, recipe, "DSX_MEAN");
    PSF_DSY_MEAN     = psMetadataLookupF32 (&status, recipe, "DSY_MEAN");
    PSF_DSX_STDEV    = psMetadataLookupF32 (&status, recipe, "DSX_STDV");
    PSF_DSY_STDEV    = psMetadataLookupF32 (&status, recipe, "DSY_STDV");

    return true;
}

// examine the model->status, fit parameters, etc and decide if the model succeeded
// set the source->type and source->mode appropriately
bool psphotEvalPSF (pmSource *source, pmModel *model) { 

    int keep;
    float dSX, dSY, SX, SY, SN;
    float nSx, nSy, Chi;

    // do we actually have a valid PSF model?
    if (model == NULL) {
	source->mode &= ~PM_SOURCE_MODE_FITTED;
	return false;
    }

    // did the model fit fail for one or another reason?
    switch (model->status) {
      case PM_MODEL_SUCCESS:
	break;
      case PM_MODEL_UNTRIED:
	source->mode &= ~PM_SOURCE_MODE_FITTED; 
	if (source->mode & PM_SOURCE_MODE_SATSTAR) {
	    psTrace ("psphot", 5, "satstar was not fitted");
	}
	return false;
      case PM_MODEL_BADARGS:
      case PM_MODEL_NONCONVERGE:
      case PM_MODEL_OFFIMAGE:
      default:
	source->mode |= PM_SOURCE_MODE_FAIL;
	if (source->mode & PM_SOURCE_MODE_SATSTAR) {
	    psTrace ("psphot", 5, "satstar failed fit");
	}
	return false;
    }

    // unless we prove otherwise, this object is a star.
    source->type = PM_SOURCE_TYPE_STAR;

    // the following source->mode information pertains to modelPSF:
    source->mode |= PM_SOURCE_MODE_PSFMODEL;

    // if the object has fitted peak above saturation, label as SATSTAR
    // this is a valid PSF object, but ignore the other quality tests
    // remember: fit does not use saturated pixels (masked)
    // XXX no extended object can saturate and stay extended...
    if (model->params->data.F32[PM_PAR_I0] >= SATURATION) {
	if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
	    psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n");
	}
	source->mode |=  PM_SOURCE_MODE_SATSTAR;
	return true;
    } 

    // if the object has a fitted peak below 0, the fit did not converge cleanly
    if (model->params->data.F32[PM_PAR_I0] <= 0) {
	source->mode |= PM_SOURCE_MODE_FAIL;
	if (source->mode & PM_SOURCE_MODE_SATSTAR) {
	    psTrace ("psphot", 5, "satstar failed fit (peak below 0)");
	}
	return false;
    } 

    // if the source was predicted to be a SATSTAR, but it fitted below saturation, 
    // make a note to the user
    if (source->mode & PM_SOURCE_MODE_SATSTAR) {
	psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n");
	source->mode &= ~PM_SOURCE_MODE_SATSTAR;
    }

    SN  = model->params->data.F32[PM_PAR_I0]/model->dparams->data.F32[PM_PAR_I0];
    SX  = model->params->data.F32[PM_PAR_SXX];
    SY  = model->params->data.F32[PM_PAR_SYY];
    dSX = model->dparams->data.F32[PM_PAR_SXX];
    dSY = model->dparams->data.F32[PM_PAR_SYY];
    Chi = model->chisqNorm / model->nDOF;

    // swing of sigma_x,y in sigmas
    nSx = (dSX - PSF_DSX_MEAN) / (PSF_DSX_STDEV * PS_MAX (PSF_MIN_DS, (SX / SN)));
    nSy = (dSY - PSF_DSY_MEAN) / (PSF_DSY_STDEV * PS_MAX (PSF_MIN_DS, (SY / SN)));

    // assign PM_SOURCE_MODE_GOODSTAR to bright objects within PSF region of dparams[]
    keep = TRUE;
    keep &= (fabs(nSx) < PSF_SHAPE_NSIGMA);
    keep &= (fabs(nSy) < PSF_SHAPE_NSIGMA);
    keep &= (SN > PSF_MIN_SN);
    keep &= (Chi < PSF_MAX_CHI);

    if (source->mode & PM_SOURCE_MODE_SATSTAR) {
	psTrace ("psphot", 5, "satstar fit results: %f, %f  %d :  %f %f : %f %f : %f %f\n", 
		     model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], keep,
		 dSX, nSx, dSY, nSy, SN, Chi);
    }

    if (keep) {
	if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
	    psTrace ("psphot", 7, "PSFSTAR kept (%f, %f  :  %f %f : %f %f : %f %f)\n", 
		     model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], dSX, nSx, dSY, nSy, SN, Chi);
	}
	return true;
    }

    // this source is not a star, warn if it was a PSFSTAR
    if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
	psTrace ("psphot", 5, "PSFSTAR demoted based on fit quality   (%f, %f  :  %f %f : %f %f : %f %f)\n", 
		  model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], dSX, nSx, dSY, nSy, SN, Chi);
    } else {
	psTrace ("psphot", 5, "fails PSF fit (%f, %f  :  %f %f : %f %f : %f %f)\n", 
		  model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], dSX, nSx, dSY, nSy, SN, Chi);
    }

    // object appears to be small, suspected defect.  this is a weak indication for a defect: a
    // large positive or negative swing can also come from the bright sources for which we are 
    // more sensitive to the wings than for the faint sources.  
    // XXX allow -NSIGMA and +NSIGMA to have different values?
    if ((nSx <= -PSF_SHAPE_NSIGMA) || (nSy <= -PSF_SHAPE_NSIGMA)) {
	source->type = PM_SOURCE_TYPE_DEFECT;
	return false;
    }

    // object appears to be large, suspected extended source
    if ((nSx >= PSF_SHAPE_NSIGMA) || (nSy >= PSF_SHAPE_NSIGMA)) {
	source->type = PM_SOURCE_TYPE_EXTENDED;
	return false;
    }

    // poor-quality fit; only keep if nothing else works...
    source->mode |= PM_SOURCE_MODE_POOR;
    return false;
}	

// examine the model->status, fit parameters, etc and decide if the model succeeded
// set the source->type and source->mode appropriately
bool psphotEvalDBL (pmSource *source, pmModel *model) { 

    // do we actually have a valid PSF model?
    if (model == NULL) {
	source->mode &= ~PM_SOURCE_MODE_FITTED;
	return false;
    }

    // did the model fit fail for one or another reason?
    switch (model->status) {
      case PM_MODEL_SUCCESS:
	break;
      case PM_MODEL_UNTRIED:
	source->mode &= ~PM_SOURCE_MODE_FITTED; 
	return false;
      case PM_MODEL_BADARGS:
      case PM_MODEL_NONCONVERGE:
      case PM_MODEL_OFFIMAGE:
      default:
	source->mode |= PM_SOURCE_MODE_FAIL;
	return false;
    }

    // unless we prove otherwise, this object is a star.
    source->type = PM_SOURCE_TYPE_STAR;

    // the following source->mode information pertains to modelPSF:
    source->mode |= PM_SOURCE_MODE_PSFMODEL;

    // if the object has fitted peak above saturation, label as SATSTAR
    // this is a valid PSF object, but ignore the other quality tests
    // remember: fit does not use saturated pixels (masked)
    // XXX no extended object can saturate and stay extended...
    if (model->params->data.F32[PM_PAR_I0] >= SATURATION) {
	if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
	    psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n");
	}
	source->mode |=  PM_SOURCE_MODE_SATSTAR;
	return true;
    } 

    // if the object has a fitted peak below 0, the fit did not converge cleanly
    if (model->params->data.F32[PM_PAR_I0] <= 0) {
	source->mode |= PM_SOURCE_MODE_FAIL;
	return false;
    } 

    // if the source was predicted to be a SATSTAR, but it fitted below saturation, 
    // make a note to the user
    if (source->mode & PM_SOURCE_MODE_SATSTAR) {
	psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n");
	source->mode &= ~PM_SOURCE_MODE_SATSTAR;
    }
    return true;
}	
