// 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 SQUARE(x) ((x)*(x))

bool poisRejectStamps(psArray *stamps,  // Array of stamps
                      psImage *mask,    // Mask image
                      const psVector *deviations, // Vector of deviations for the stamps
                      const poisConfig *config // Configuration
    )
{
    assert(stamps);
    assert(mask);
    assert(deviations);
    assert(config);
    assert(stamps->n == deviations->n);
    assert(mask->numCols == config->xImage && mask->numRows == config->yImage);
    assert(mask->type.type == PS_TYPE_U8);
    assert(deviations->type.type == PS_TYPE_F32);

    bool masked = false;                // Have we masked any stamps?

    // Don't want the standard deviation, but rather the deviation from zero
    double sum = 0.0;
    int numDev = 0;
    for (int i = 0; i < deviations->n; i++) {
        poisStamp *stamp = stamps->data[i];     // The stamp
        // Only interested in the stamps that we're actually using
        if (stamp->status == POIS_STAMP_USED) {
            sum += SQUARE(deviations->data.F32[i]);
            numDev++;
        }
    }
    float meanDev = sqrtf(sum / (float)numDev);
    float limit = meanDev * config->sigmaRej;
    psTrace("pois.rejectStamps", 2, "Mean RMS deviation: %f\n", meanDev);
    psTrace("pois.rejectStamps", 2, "Rejection limit: %f\n", limit);

    // Reject stamps
    for (int s = 0; s < deviations->n; s++) {
        poisStamp *stamp = stamps->data[s];
        if (stamp->status == POIS_STAMP_USED &&
            (fabsf(deviations->data.F32[s]) > limit || ! isfinite(deviations->data.F32[s]))) {
            masked = true;
            psTrace("pois.rejectStamps", 1, "Rejecting stamp %d (%d,%d): %f\n", s, stamp->x, stamp->y,
                    deviations->data.F32[s]);

            // Mask out the stamp in the image so you don't find it again
            for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) {
                for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) {
                    mask->data.U8[y][x] |= POIS_MASK_STAMP;
                }
            }

            // Set stamp for replacement
            stamp->x = 0;
            stamp->y = 0;
            stamp->status = POIS_STAMP_RESET;

        } // Bad stamps
    } // Iterating over stamps

    return masked;
}
