# include "psphotInternal.h"

bool psphotAddNoise (pmConfig *config, const pmFPAview *view) {
    return psphotAddOrSubNoise (config, view, true);
}

bool psphotSubNoise (pmConfig *config, const pmFPAview *view) {
    return psphotAddOrSubNoise (config, view, false);
}

// for now, let's store the detections on the readout->analysis for each readout
bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, bool add)
{
    bool status = true;

    // select the appropriate recipe information
    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    psAssert (recipe, "missing recipe?");

    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");

    // loop over the available readouts
    for (int i = 0; i < num; i++) {
	if (!psphotAddOrSubNoiseReadout (config, view, "PSPHOT.INPUT", i, recipe, add)) {
            psError (PSPHOT_ERR_CONFIG, false, "failed on to modify noise for PSPHOT.INPUT entry %d", i);
	    return false;
	}
    }
    return true;
}

bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool add) {

    bool status = false;
    psEllipseShape oldshape;
    psEllipseShape newshape;
    psEllipseAxes axes;

    // find the currently selected readout
    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
    psAssert (file, "missing file?");

    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    psAssert (readout, "missing readout?");
    psAssert (readout->parent, "missing cell?");
    psAssert (readout->parent->concepts, "missing concepts?");

    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    psAssert (detections, "missing detections?");

    psArray *sources = detections->allSources;
    psAssert (sources, "missing sources?");

    psTimerStart ("psphot.noise");

    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    psAssert (maskVal, "missing mask value?");

    // increase variance by factor*(object noise):
    // weight = flux/gain + rn^2/g^2
    // we are adding factor*flux/gain

    float FACTOR = psMetadataLookupF32 (&status, recipe, "NOISE.FACTOR");
    PS_ASSERT (status, false);
    float SIZE = psMetadataLookupF32 (&status, recipe, "NOISE.SIZE");
    PS_ASSERT (status, false);

    if (SIZE <= 0) {
       return true;
    }

    float GAIN = psMetadataLookupF32(&status, readout->parent->concepts, "CELL.GAIN"); // Cell gain
    PS_ASSERT (status, false);
    if (isfinite(GAIN)) {
        FACTOR /= GAIN;
    }

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

        // skip sources which were not subtracted
	// NOTE: this bit is not modified when pmSourceOp applies to noise
        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;

        // select appropriate model
        pmModel *model = pmSourceGetModel (NULL, source);
        if (model == NULL) continue;  // model must be defined

        if (add) {
            psTrace ("psphot", 4, "adding noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
        } else {
            psTrace ("psphot", 4, "remove noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
        }

        psF32 *PAR = model->params->data.F32;

        // save original values
        float oldI0  = PAR[PM_PAR_I0];
        oldshape.sx  = PAR[PM_PAR_SXX];
        oldshape.sy  = PAR[PM_PAR_SYY];
        oldshape.sxy = PAR[PM_PAR_SXY];

        // XXX can this be done more intelligently?
        if (oldI0 == 0.0) continue;
        if (!isfinite(oldI0)) continue;

        // increase size and height of source
        axes = psEllipseShapeToAxes (oldshape, 20.0);
        axes.major *= SIZE;
        axes.minor *= SIZE;
        newshape = psEllipseAxesToShape (axes);
        PAR[PM_PAR_I0]  = FACTOR*oldI0;
        PAR[PM_PAR_SXX] = newshape.sx;
        PAR[PM_PAR_SYY] = newshape.sy;
        PAR[PM_PAR_SXY] = newshape.sxy;

        // XXX if we use pmSourceOp, the size (and possibly Io) will not be respected
        pmSourceOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, add, maskVal, 0, 0);

        // restore original values
        PAR[PM_PAR_I0]  = oldI0;
        PAR[PM_PAR_SXX] = oldshape.sx;
        PAR[PM_PAR_SYY] = oldshape.sy;
        PAR[PM_PAR_SXY] = oldshape.sxy;
    }
    if (add) {
        psLogMsg ("psphot.noise", PS_LOG_INFO, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
    } else {
        psLogMsg ("psphot.noise", PS_LOG_INFO, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
    }
    return true;
}
