# 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, psMaskType maskVal) {

  psTimerStart ("psphot");

  // 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->x;
            modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->y;
        } 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);

    // 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?
