# include "psphotInternal.h"

bool psphotAddNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) {
  return psphotAddOrSubNoise (readout, sources, recipe, true);
}

bool psphotSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) {
  return psphotAddOrSubNoise (readout, sources, recipe, false);
}

bool psphotAddOrSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe, bool add) {

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

    PS_ASSERT (readout, false);
    PS_ASSERT (readout->parent, false);
    PS_ASSERT (readout->parent->concepts, false);

    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
    assert (maskVal);

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