/** @file ppSubVarianceFactors.c
 *
 *  @brief
 *
 *  @ingroup ppSub
 *
 *  @author IfA
 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-02-18 00:31:20 $
 *  Copyright 2009 Institute for Astronomy, University of Hawaii
 */

#include "ppSub.h"

/**
 * Calculate the variance factor for the output image based on the input images
 */
bool ppSubVarianceFactors (pmConfig *config, psMetadata *stats, const pmFPAview *view) {

    // Look up recipe values
    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
    psAssert(recipe, "We checked this earlier, so it should be here.");

    // convolved input images
    pmCell *inCell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Input cell
    pmCell *refCell = pmFPAfileThisCell(config->files, view, "PPSUB.REF.CONV"); // Reference cell

    pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input cell
    pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference cell

    pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output readout

    // Update variance factors
    // It's not possible to do this perfectly, since we're combining different images:
    // vf_sum sigma_sum^2 = vf_1 * sigma_1^2 + vf_2 * sigma_2^2
    // It's not possible to write vf_sum as a function of vf_1 and vf_2 with no reference to the sigmas.
    // Instead, we're going to cheat.

    bool mdok;                      // Status of MD lookup
    pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, outRO->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
    if (!mdok) {
        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find subtraction kernels.");
        return false;
    }

    float vfIn = psMetadataLookupF32(NULL, inCell->concepts, "CELL.VARFACTOR"); // Variance factor for input
    if (!isfinite(vfIn)) {
        vfIn = 1.0;
    }
    float vfRef = psMetadataLookupF32(NULL, refCell->concepts, "CELL.VARFACTOR"); // Variance factor for ref
    if (!isfinite(vfRef)) {
        vfRef = 1.0;
    }

    if (kernels->mode == PM_SUBTRACTION_MODE_1) {
        vfIn *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
    }
    if (kernels->mode == PM_SUBTRACTION_MODE_2) {
        vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
    }
    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
        vfIn *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
        vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, true);
    }

    psStats *vfStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics
    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator

    // mask these values when examining the background
    psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask going in to pmSubtractionMatch
    psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels

    // measure the mean of the convolved input variance image
    if (!psImageBackground(vfStats, NULL, inConv->variance, inConv->mask, maskVal | maskBad, rng)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved input");
        psFree(vfStats);
        psFree(rng);
        return false;
    }
    float inMeanVar = vfStats->robustMedian; // Mean variance of input

    // measure the mean of the convolved reference variance image
    if (!psImageBackground(vfStats, NULL, refConv->variance, refConv->mask, maskVal | maskBad, rng)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved reference");
        psFree(vfStats);
        psFree(rng);
        return false;
    }
    float refMeanVar = vfStats->robustMedian; // Mean variance of reference

    psFree(rng);
    psFree(vfStats);

    float vf = (vfIn * inMeanVar + vfRef * refMeanVar) / (inMeanVar + refMeanVar); // Resulting var factor
    psMetadataItem *vfItem = psMetadataLookup(outRO->parent->concepts, "CELL.VARFACTOR");
    vfItem->data.F32 = vf;

    return true;
}
