# include "psphotInternal.h"

// A guess for when the moments aren't available
static pmModel *wildGuess(pmSource *source, // Source for which to guess
                          pmPSF *psf    // The point-spread function
    )
{
    pmModel *model = pmModelAlloc(psf->type);
    psF32 *PAR = model->params->data.F32;
    PAR[PM_PAR_SKY]  = 0;
    // XXX get this from the image pixels
    PAR[PM_PAR_I0]   = source->peak->flux;
    PAR[PM_PAR_XPOS] = source->peak->xf;
    PAR[PM_PAR_YPOS] = source->peak->yf;
    return model;
}

// construct an initial PSF model for each object
bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) {

    bool status;

    psTimerStart ("psphot");

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

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

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

    // setup the PSF fit radius details
    psphotInitRadiusPSF (recipe, psf->type);

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

	// skip non-astronomical objects (very likely defects)
	if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
	if (source->type == PM_SOURCE_TYPE_SATURATED) continue;

	// XXX if a source is faint, it will not have moments measured.
	// it must be modelled as a PSF.  In this case, we need to use
	// the peak centroid to get the coordinates and get the peak flux
	// from the image?
	pmModel *modelEXT;
	if (!source->moments) {
	    modelEXT = wildGuess(source, psf);
	} else {
	    // use the source moments, etc to guess basic model parameters
	    modelEXT = pmSourceModelGuess (source, psf->type);
	    if (!modelEXT) {
		modelEXT = wildGuess(source, psf);
	    }
	    // these valuse are set in pmSourceModelGuess, should this rule be in there as well?
	    if (source->mode &  PM_SOURCE_MODE_SATSTAR) {
		modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->Mx;
		modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->My;
	    } else {
		modelEXT->params->data.F32[PM_PAR_XPOS] = source->peak->xf;
		modelEXT->params->data.F32[PM_PAR_YPOS] = source->peak->yf;
	    }
	}

	// set PSF parameters for this model (apply 2D shape model)
	pmModel *modelPSF = pmModelFromPSF (modelEXT, psf);
	if (modelPSF == NULL) {
	    psError(PSPHOT_ERR_PSF, false,
		    "Failed to determine PSF model at r,c = (%d,%d); trying centre of image",
		    source->peak->y, source->peak->x);
	    //
	    // Try the centre of the image
	    //
	    modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols;
	    modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows;
	    modelPSF = pmModelFromPSF (modelEXT, psf);
	    if (modelPSF == NULL) {
		psError(PSPHOT_ERR_PSF, false,
			"Failed to determine PSF model at centre of image");
		psFree(modelEXT);
		return false;
	    }

	    source->mode |= PM_SOURCE_MODE_BADPSF;
	}
	psFree (modelEXT);

	// XXX need to define the guess flux?
	// set the fit radius based on the object flux limit and the model
	psphotCheckRadiusPSF (readout, source, modelPSF, markVal);

	// set the source PSF model
	source->modelPSF = modelPSF;
	source->modelPSF->residuals = psf->residuals;

	pmSourceCacheModel (source, maskVal);
    }
    psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot"));
    return true;
}

// XXX do we always know which model is supposed to be used?
