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

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <pslib.h>
#include <psmodules.h>
#include <psphot.h>

#include "ppSub.h"

bool ppSubBackground(pmConfig *config)
{
    psAssert(config, "Require configuration");

    bool mdok; // Status of metadata lookups

    psMetadata *ppSubRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSub
    psAssert(ppSubRecipe, "Need PPSUB recipe");
    psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE); // Recipe for psphot
    psAssert(psphotRecipe, "Need PSPHOT recipe for binning");

    psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask

    pmFPAview *view = ppSubViewReadout(); // View to readout
    pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output image
    pmReadout *modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL"); // Background model

    // Generate the background model, if required
    if (!modelRO) {
        // Create the background model
        if (!psphotModelBackgroundReadoutFileIndex(config, view, "PPSUB.OUTPUT", 0)) {
            psError(psErrorCodeLast(), false, "Unable to model background");
            psFree(view);
            return false;
        }
        // select the model readout (should now exist)
        modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL");
        if (!modelRO) {
            psError(psErrorCodeLast(), false, "Unable to find background model");
            psFree(view);
            return false;
        }
    }
    psFree(view);

    psImageBinning *binning = psMetadataLookupPtr(&mdok, modelRO->analysis,
                                                  "PSPHOT.BACKGROUND.BINNING"); // Binning for model
    psImage *modelImage = modelRO->image; // Background model
    psImage *image = outRO->image; // Image of interest
    psImage *mask = outRO->mask; // Mask of interest

    // Do the background subtraction
    int numCols = image->numCols, numRows = image->numRows; // Size of image
    for (int y = 0; y < numRows; y++) {
        for (int x = 0; x < numCols; x++) {
            if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskBad) {
                image->data.F32[y][x] = 0.0;
            } else {
                float value = psImageUnbinPixel(x + 0.5, y + 0.5, modelImage, binning); // Background value
                if (!isfinite(value)) {
                    image->data.F32[y][x] = NAN;
                    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskBad;
                } else {
                    image->data.F32[y][x] -= value;
                }
            }
        }
    }

    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");

    return true;
}
