/* @file pmSubtraction.h
 *
 * PSF-matched image subtraction, based on the Alard & Lupton (1998) and Alard (2000) methods.
 *
 * @author Paul Price, IfA
 * @author GLG, MHPCC
 *
 * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
 * @date $Date: 2007-08-30 21:45:26 $
 *
 * Copyright 2004-207 Institute for Astronomy, University of Hawaii
 */

#ifndef PM_SUBTRACTION_H
#define PM_SUBTRACTION_H

#include <pslib.h>
#include "pmSubtractionKernels.h"

/// @addtogroup imcombine Image Combinations
/// @{

typedef enum {
    PM_SUBTRACTION_MASK_CLEAR     = 0x00, // No masking
    PM_SUBTRACTION_MASK_REF       = 0x01, // Reference image is bad
    PM_SUBTRACTION_MASK_INPUT     = 0x02, // Input image is bad
    PM_SUBTRACTION_MASK_CONVOLVE  = 0x04, // If convolved, would be bad
    PM_SUBTRACTION_MASK_FOOTPRINT = 0x08, // Bad pixel within the stamp footprint
    PM_SUBTRACTION_MASK_BORDER    = 0x10, // Image border
    PM_SUBTRACTION_MASK_REJ       = 0x20, // Previously tried as a stamp, and rejected
} pmSubtractionMasks;

/// Generate a mask for use in the subtraction process
psImage *pmSubtractionMask(const psImage *refMask, ///< Mask for the reference image (will be convolved)
                           const psImage *inMask, ///< Mask for the input image, or NULL
                           psMaskType maskVal, ///< Value to mask out
                           int size, ///< Half-size of the kernel (pmSubtractionKernels.size)
                           int footprint ///< Half-size of the kernel footprint
    );

/// Calculate the least-squares equation to match the image quality
bool pmSubtractionCalculateEquation(psArray *stamps, ///< The stamps for which to calculate the equation,
                                    const pmSubtractionKernels *kernels, ///< Kernel parameters
                                    int footprint ///< Half-size of region over which to calculate equation
                                    );

/// Solve the least-squares equation to match the image quality
psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL
                                     const psArray *stamps ///< Array of stamps
                                     );

/// Reject stamps
int pmSubtractionRejectStamps(psArray *stamps, ///< Array of stamps to check for rejection
                              psImage *subMask, ///< Subtraction mask
                              const psVector *solution, ///< Solution vector
                              int footprint, ///< Region to mask if stamp is bad
                              float sigmaRej, ///< Number of RMS deviations above zero at which to reject
                              const pmSubtractionKernels *kernels ///< Kernel parameters
    );

/// Generate an image of the convolution kernel
psImage *pmSubtractionKernelImage(const psVector *solution, ///< Solution vector
                                  const pmSubtractionKernels *kernels, ///< Kernel parameters
                                  float x, float y ///< Normalised position [-1,1] for which to generate image
                                  );

/// Generate images of the convolution kernel elements
psArray *pmSubtractionKernelSolutions(const psVector *solution, ///< Solution vector
                                      const pmSubtractionKernels *kernels, ///< Kernel parameters
                                      float x, float y ///< Normalised position [-1,1] for images
    );

/// Convolve image in preparation for subtraction
bool pmSubtractionConvolve(psImage **outImage, ///< Output image
                           psImage **outWeight, ///< Output weight map (or NULL)
                           psImage **outMask, ///< Output mask (or NULL)
                           const psImage *inImage, ///< Input image
                           const psImage *inWeight, ///< Input weight map (or NULL)
                           const psImage *subMask, ///< Subtraction mask (or NULL)
                           psMaskType blank, ///< Mask value for blank regions
                           const psRegion *region, ///< Region to convolve (or NULL)
                           const psVector *solution, ///< The solution vector
                           const pmSubtractionKernels *kernels, ///< Kernel parameters
                           bool useFFT  ///< Use Fast Fourier Transform for the convolution?
    );

/// Mark the non-convolved part of the image as blank
bool pmSubtractionBorder(psImage *image,///< Image
                         psImage *weight, ///< Weight map (or NULL)
                         psImage *mask, ///< Mask (or NULL)
                         const pmSubtractionKernels *kernels, ///< Kernel parameters (for kernel size)
                         psMaskType blank ///< Mask value for blank regions
    );

/// Generate the convolution of an image, given a precalculated kernel
///
/// The 'image' is a kernel for convenience --- intended to be a stamp
psKernel *p_pmSubtractionConvolveStampPrecalc(const psKernel *image, ///< Image to convolve
                                              const psKernel *kernel ///< Kernel by which to convolve
    );

/// @}
#endif
