// Copyright (C) 2006  Paul A. Price (price@ifa.hawaii.edu)
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

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

#include <stdio.h>
#include <assert.h>
#include "pslib.h"
#include "pois.h"

#define MAXCHAR 80

psVector *poisCalculateDeviations(psVector *deviations, // Output array of deviations, or NULL
                                  psArray *stamps, // Array of stamps
                                  psImage *refImage, // Reference image
                                  psImage *inImage, // Input image
                                  psImage *mask, // Mask image
                                  psArray *kernelParams, // Array of kernel parameters
                                  psVector *solution, // Solution vector
                                  poisConfig *config // Configuration
    )
{
    // Check inputs
    assert(!deviations || deviations->type.type == PS_TYPE_F32);
    assert(!deviations || deviations->n == stamps->n);
    assert(stamps);
    assert(refImage->type.type == PS_TYPE_F32);
    assert(inImage->type.type == PS_TYPE_F32);
    assert(refImage->numRows == inImage->numRows);
    assert(refImage->numCols == inImage->numCols);
    assert(mask->type.type == PS_TYPE_U8);
    assert(refImage->numRows == mask->numRows);
    assert(refImage->numCols == mask->numCols);
    assert(solution->n == kernelParams->n + 1);
    assert(solution->type.type == PS_TYPE_F64);
    assert(config);

    if (!deviations) {
        deviations = psVectorAlloc(stamps->n, PS_TYPE_F32);
    }

    int footprint = config->footprint;  // Size of stamp footprint
    int xSize = footprint + config->xKernel; // (Half) size of subimage in x
    int ySize = footprint + config->yKernel; // (Half) size of subimage in y
    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics

    // Pre-allocate for optimisation:
    psImage *subStamp = psImageAlloc(2 * xSize, 2 * ySize, PS_TYPE_F32); // Subtraction of stamp

    for (int s = 0; s < stamps->n; s++) {
        poisStamp *stamp = stamps->data[s]; // The coordinates of the stamp of interest
        int x = stamp->x;               // Stamp x coord
        int y = stamp->y;               // Stamp y coord
        if (stamp->status == POIS_STAMP_USED) {
            psRegion stampRegion = {x - xSize, x + xSize, y - ySize, y + ySize};
            psImage *refStamp = psImageSubset(refImage, stampRegion);
            psImage *inStamp = psImageSubset(inImage, stampRegion);
            psImage *maskStamp = psImageSubset(mask, stampRegion);
            psImage *convRefStamp = poisConvolveImage(refStamp, maskStamp, solution, kernelParams, config);
            // Calculate chi^2
            (void)psBinaryOp(subStamp, inStamp, "-", convRefStamp);
            (void)psBinaryOp(subStamp, subStamp, "*", subStamp);
            (void)psBinaryOp(subStamp, subStamp, "/", inStamp);

            // Trim regions.  Need separate regions for the subtraction and the mask, since the subtraction
            // isn't a child image (it's the subtraction of two child images, but placed in a pre-allocated
            // image which isn't a child image), but the mask *is* a child image.  Therefore the regions are
            // different, since they refer to the region on the parent image (if it exists).
            psRegion stampTrim = { config->xKernel, config->xKernel + 2 * footprint, config->yKernel,
                                   config->yKernel + 2 * footprint };
#if 1
            psRegion maskTrim = { x - xSize + config->xKernel, x + xSize - config->xKernel,
                                  y - ySize + config->yKernel, y + ySize - config->yKernel };
            psImage *maskStampTrim = psImageSubset(maskStamp, maskTrim);
#else
            psImage *maskStampTrim = psImageSubset(maskStamp, stampTrim);
#endif
            psImage *subStampTrim = psImageSubset(subStamp, stampTrim);
            (void)psImageStats(stats, subStampTrim, maskStampTrim, POIS_MASK_BAD | POIS_MASK_NEAR_BAD);

            deviations->data.F32[s] = sqrtf(stats->sampleMean / 2.0);
            psTrace("pois.calculateDeviations", 5, "Deviation of stamp %d (%d,%d) is %f\n", s, x, y,
                    deviations->data.F32[s]);

#ifdef TESTING
            char stampName[MAXCHAR];            // File name for stamp
            snprintf(stampName, MAXCHAR, "stamp%d.fits", s);
            psFits *stampFile = psFitsOpen(stampName, "w");
            if (!psFitsWriteImage(stampFile, NULL, subStampTrim, 0, NULL)) {
                psErrorStackPrint(stderr, "Unable to write stamp: %s\n", stampName);
            }
            psFitsClose(stampFile);
#endif

#if 0
            psFree(convRefStamp);
            psFree(maskStampTrim);
            psFree(subStampTrim);
            psFree(maskStamp);
            psFree(refStamp);
            psFree(inStamp);
#endif
        }
    }

    psFree(stats);
    psFree(subStamp);

    return deviations;
}
