# include "psphotInternal.h"

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

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

    psTimerStart ("psphot.fit.nonlinear");

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

    // 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->mode &  PM_SOURCE_MODE_CR_LIMIT) 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
        if (source->peak->SN < FIT_SN_LIM) continue;

        // XXX this should use peak?
        if (source->peak->xf < AnalysisRegion.x0) continue;
        if (source->peak->yf < AnalysisRegion.y0) continue;
        if (source->peak->xf > AnalysisRegion.x1) continue;
        if (source->peak->yf > 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 ++;

	// XXX TEST
	// if ((fabs(source->peak->x - 1202) < 2) && (fabs(source->peak->y - 1065) < 2)) {
	// psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6);
	// }

        // try fitting PSFs, then try extended sources
        // these functions subtract the resulting fitted source (XXX and update the modelFlux?)
	// XXX re-consider conditions under which the source has EXT fit:
	// I should try EXT if the source size measurement says it is large
	if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
	    if (psphotFitBlob (readout, source, sources, psf, maskVal, markVal)) {
		source->type = PM_SOURCE_TYPE_EXTENDED;
		psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
		Next ++;
		continue;
	    }
	} else {
	    if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
		source->type = PM_SOURCE_TYPE_STAR;
		psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
		Npsf ++;
		continue;
	    }
	}

        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);
        Nfail ++;

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

    if (psTraceGetLevel("psphot") >= 6) {
      psphotSaveImage (NULL, readout->image,  "image.v2.fits");
    }

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

    psphotVisualShowResidualImage (readout);
    psphotVisualShowFlags (sources);

    return true;
}
