# include "psphotInternal.h"

// select objects fitted with PSF model
// re-fit all of them with the non-poisson errors
// only allow the normalization to vary
// XXX eventually, we should be able to do this with linear fitting...
// XXX UNUSED
bool psphotWeightBias (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) {

    int Nfit = 0;
    bool status;
    psF32 *PARp, *PARc;
    psF32 x, y;

    psTimerStart ("psphot");

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

    // set fitting method to use non-poisson errors (local sky error)
    float SKY_SIG = psMetadataLookupF32 (&status, recipe, "SKY_SIG");
    if (!status) {
      SKY_SIG = 1.0;
        psWarning("SKY_SIG is not set --- defaulting to %f\n", SKY_SIG);
    }
    // use poissonian errors or local-sky errors
    bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS");
    if (!status) {
        POISSON_ERRORS = true;
        psWarning("POISSON_ERRORS is not set in the recipe --- defaulting to true.\n");
    }
    pmSourceFitModelInit (15, 0.1, PS_SQR(SKY_SIG), POISSON_ERRORS);

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

    FILE *f = fopen ("bias.dat", "w");
    if (f == NULL) psAbort("can't open output file");

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

        // skip lower-quality objects
        if (source->type != PM_SOURCE_TYPE_STAR) continue;
        if (source->mode &  PM_SOURCE_MODE_POOR) continue;
        if (source->mode &  PM_SOURCE_MODE_FAIL) continue;
        if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
        if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
        if (source->mode &  PM_SOURCE_MODE_SATSTAR) continue;

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

        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;

        // replace object in image
        pmModelAdd (source->pixels, source->mask, source->modelPSF, PM_MODEL_OP_FULL);

        // make a temporary model (we don't keep the result of this analysis)
        pmModel *PSF = pmModelCopy (source->modelPSF);

        // extend source radius as needed
        psphotCheckRadiusPSF (readout, source, PSF);

        x = PSF->params->data.F32[2];
        y = PSF->params->data.F32[3];

        // fit PSF model (set/unset the pixel mask)
        psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "OR", PM_MASK_MARK);
        pmSourceFitModel (source, PSF, PM_SOURCE_FIT_NORM);
        psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));

        // re-subtract PSF for object, leave local sky
        pmModelSub (source->pixels, source->mask, source->modelPSF, PM_MODEL_OP_FULL);

        PARp = source->modelPSF->params->data.F32;
        PARc = PSF->params->data.F32;
        fprintf (f, "%7.1f %7.1f %9.2f %9.2f %10.3f\n", PARp[2], PARp[3], PARp[1], PARc[1], source->moments->dSky);
        Nfit ++;
    }
    fclose (f);
    psLogMsg ("psphot", 3, "measure PSF weighting bias for %d objects: %f sec\n", Nfit, psTimerMark ("psphot"));

    return (true);
}
