# include "psphotInternal.h"

// XXX I don't like this name
bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal) {

    int Nfit = 0;
    int Npsf = 0;
    int Next = 0;
    int Nfail = 0;
    bool status;

    psTimerStart ("psphot");

    // source analysis is done in S/N order (brightest first)
    sources = psArraySort (sources, pmSourceSortBySN);

    // S/N limit to perform full non-linear fits
    float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM");

    psphotInitLimitsPSF (recipe, readout);
    psphotInitLimitsEXT (recipe);
    psphotInitRadiusPSF (recipe, psf->type);

    psphotFitInit ();

    // option to limit analysis to a specific region
    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");

    for (int i = 0; i < sources->n; i++) {
        // if (i%100 == 0) psphotFitSummary ();

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

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

        // skip DBL second sources (ie, added by psphotFitBlob)
        if (source->mode &  PM_SOURCE_MODE_PAIR) continue;

        // limit selection to some SN limit
        // XXX this should use peak?
        if (source->moments == NULL) continue;
        if (source->moments->SN < FIT_SN_LIM) continue;

        // XXX this should use peak?
        if (source->moments->x < AnalysisRegion.x0) continue;
        if (source->moments->y < AnalysisRegion.y0) continue;
        if (source->moments->x > AnalysisRegion.x1) continue;
        if (source->moments->y > AnalysisRegion.y1) continue;

        // if model is NULL, we don't have a starting guess
        if (source->modelPSF == NULL) continue;

        // skip sources which are insignificant flux?
        if (source->modelPSF->params->data.F32[1] < 0.1) {
            psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
                     source->modelPSF->params->data.F32[1],
                     source->modelPSF->params->data.F32[2],
                     source->modelPSF->params->data.F32[3]);
            continue;
        }

        // replace object in image
        if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
        }
        Nfit ++;

        // try fitting PSFs, then try extended sources
        // these functions subtract the resulting fitted source (XXX and update the modelFlux?)
        if (psphotFitBlend (readout, source, psf, maskVal)) {
            psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y);
            Npsf ++;
            continue;
        }
        if (psphotFitBlob (readout, source, sources, psf, maskVal)) {
            psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y);
            Next ++;
            continue;
        }

        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->moments->x, source->moments->y);
        Nfail ++;

        // re-subtract the object, leave local sky
        pmSourceCacheModel (source, maskVal);
        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
        source->mode |= PM_SOURCE_MODE_SUBTRACTED;
        source->mode |= PM_SOURCE_MODE_TEMPSUB;
    }

    psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot"), Nfit, Npsf, Next, Nfail, sources->n - Nfit);
    return true;
}
