# include "psphotInternal.h"
static int npass = 0;

// generate the median in NxN boxes, clipping heavily
// linear interpolation to generate full-scale model
bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filename) 
{
    bool status = true;
    pmReadout *background = NULL;
    pmReadout *backSub = NULL;

    psTimerStart ("psphot.background");

    // find the currently selected readout
    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, filename);
    pmFPA *inFPA = file->fpa;
    pmReadout *readout = pmFPAviewThisReadout (view, inFPA);
    psImage *image = readout->image;
    psImage *mask  = readout->mask;

    // find the currently selected readout
    pmReadout *model = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKMDL");
    assert (model);

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

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

    psImageBinning *binning = psMetadataLookupPtr(&status, recipe, "PSPHOT.BACKGROUND.BINNING");
    assert (binning);

    // select background pixels, from output background file, or create
    file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND");
    if (file) {
        // we are using PSPHOT.BACKGND as an I/O file: select readout or create
        if (file->mode == PM_FPA_MODE_INTERNAL) {
            background = file->readout;
        } else {
            background = pmFPAviewThisReadout (view, file->fpa);
        }
        if (background == NULL) {
            // readout does not yet exist: create from input
            pmFPAfileCopyStructureView (file->fpa, inFPA, 1, 1, view);
            background = pmFPAviewThisReadout (view, file->fpa);
            if ((image->numCols != background->image->numCols) || (image->numRows != background->image->numRows)) {
                psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for background dimensions");
                return false;
            }
        }
    } else {
        background = pmFPAfileDefineInternal (config->files, "PSPHOT.BACKGND", image->numCols, image->numRows, PS_TYPE_F32);
    }
    psF32 **backData = background->image->data.F32;

    // linear interpolation to full-scale
    if (!psImageUnbin (background->image, model->image, binning)) {
        psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
        return false;
    }

    psLogMsg ("psphot", PS_LOG_MINUTIA, "build resampled image: %f sec\n", psTimerMark ("psphot.background"));

    // back-sub image pixels, from output background file (don't create if not requested)
    file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKSUB");
    if (file) {
        // we are using PSPHOT.BACKSUB as an I/O file: select readout or create
        backSub = pmFPAviewThisReadout (view, file->fpa);
        if (backSub == NULL) {
            // readout does not yet exist: create from input
            pmFPAfileCopyStructureView (file->fpa, inFPA, 1, 1, view);
            backSub = pmFPAviewThisReadout (view, file->fpa);
        }
    }

    if (psTraceGetLevel("psphot") > 5) {
        char name[256];
        sprintf (name, "image.%02d.fits", npass);
        psphotSaveImage (NULL, image, name);
        sprintf (name, "back.%02d.fits", npass);
        psphotSaveImage (NULL, background->image, name);
        sprintf (name, "mask.%02d.fits", npass);
        psphotSaveImage (NULL, mask, name);
        sprintf (name, "backmdl.%02d.fits", npass);
        psphotSaveImage (NULL, model->image, name);
    }

    // subtract the background model (save in backSub, if requested)
    // XXX if needed, multithread this (fairly trivial)
    for (int j = 0; j < image->numRows; j++) {
        for (int i = 0; i < image->numCols; i++) {
            image->data.F32[j][i] -= backData[j][i];
            if (backSub) {
                backSub->image->data.F32[j][i] = image->data.F32[j][i];
            }
        }
    }

    if (psTraceGetLevel("psphot") > 5) {
        char name[256];
        sprintf (name, "backsub.%02d.fits", npass);
        psphotSaveImage (NULL, image, name);
    }
    psLogMsg ("psphot", PS_LOG_INFO, "subtracted background model: %f sec\n", psTimerMark ("psphot.background"));

    // the pmReadout selected in this function are all view on entries in config->files

    npass ++;
    return true;
}
