#ifndef PM_SUBTRACTION_STAMPS_H
#define PM_SUBTRACTION_STAMPS_H

#include <pslib.h>

#include "pmSubtractionKernels.h"

/// Status of stamp
typedef enum {
    PM_SUBTRACTION_STAMP_INIT,          ///< Initial state
    PM_SUBTRACTION_STAMP_FOUND,         ///< Found a suitable source for this stamp
    PM_SUBTRACTION_STAMP_CALCULATE,     ///< Calculate matrix and vector values for this stamp
    PM_SUBTRACTION_STAMP_USED,          ///< Use this stamp
    PM_SUBTRACTION_STAMP_REJECTED,      ///< This stamp has been rejected
    PM_SUBTRACTION_STAMP_NONE           ///< No stamp in this region
} pmSubtractionStampStatus;

/// A list of stamps
typedef struct {
    long num;                           ///< Number of stamps
    psArray *stamps;                    ///< The stamps
    psArray *regions;                   ///< Regions for each stamp
    psArray *x, *y;                     ///< Coordinates for possible stamps (or NULL)
    psArray *flux;                      ///< Fluxes for possible stamps (or NULL)
    int footprint;                      ///< Half-size of stamps
} pmSubtractionStampList;

/// Allocate a list of stamps
pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, // Number of columns in image
                                                    int numRows, // Number of rows in image
                                                    const psRegion *region, // Region for stamps, or NULL
                                                    int footprint, // Half-size of stamps
                                                    float spacing // Rough average spacing between stamps
    );

/// Assertion for stamp list to be valid
#define PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(LIST, RETURNVALUE) { \
    PS_ASSERT_PTR_NON_NULL(LIST, RETURNVALUE); \
    PS_ASSERT_ARRAY_NON_NULL((LIST)->stamps, RETURNVALUE); \
    PS_ASSERT_ARRAY_NON_NULL((LIST)->regions, RETURNVALUE); \
    PS_ASSERT_INT_POSITIVE((LIST)->num, RETURNVALUE); \
    PS_ASSERT_ARRAY_SIZE((LIST)->stamps, (LIST)->num, RETURNVALUE); \
    PS_ASSERT_ARRAY_SIZE((LIST)->regions, (LIST)->num, RETURNVALUE); \
    PS_ASSERT_INT_NONNEGATIVE((LIST)->footprint, RETURNVALUE); \
    if ((LIST)->x || (LIST)->y || (LIST)->flux) { \
        PS_ASSERT_ARRAY_NON_NULL((LIST)->x, RETURNVALUE); \
        PS_ASSERT_ARRAY_NON_NULL((LIST)->y, RETURNVALUE); \
        PS_ASSERT_ARRAY_NON_NULL((LIST)->flux, RETURNVALUE); \
        PS_ASSERT_ARRAY_SIZE((LIST)->x, (LIST)->num, RETURNVALUE); \
        PS_ASSERT_ARRAY_SIZE((LIST)->y, (LIST)->num, RETURNVALUE); \
        PS_ASSERT_ARRAY_SIZE((LIST)->flux, (LIST)->num, RETURNVALUE); \
    } \
}

/// A stamp for image subtraction
typedef struct {
    float x, y;                         ///< Position
    float flux;                         ///< Flux
    float xNorm, yNorm;                 ///< Normalised position
    psKernel *image1;                   ///< Reference image postage stamp
    psKernel *image2;                   ///< Input image postage stamp
    psKernel *variance;                 ///< Variance image postage stamp, or NULL
    psArray *convolutions1;             ///< Convolutions of image 1 for each kernel component, or NULL
    psArray *convolutions2;             ///< Convolutions of image 2 for each kernel component, or NULL
    psImage *matrix1, *matrix2;         ///< Least-squares matrices for each image, or NULL
    psImage *matrixX;                   ///< Cross-matrix (for mode DUAL), or NULL
    psVector *vector1, *vector2;        ///< Least-squares vectors for each image, or NULL
    pmSubtractionStampStatus status;    ///< Status of stamp
} pmSubtractionStamp;

/// Allocate a stamp
pmSubtractionStamp *pmSubtractionStampAlloc(void);

/// Find stamps on an image
pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, ///< Output stamps, or NULL
                                                const psImage *image, ///< Image for which to find stamps
                                                const psImage *mask, ///< Mask, or NULL
                                                const psRegion *region, ///< Region to search, or NULL
                                                float threshold, ///< Threshold for stamps in the image
                                                int footprint, ///< Half-size for stamps
                                                float spacing, ///< Rough spacing for stamps
                                                pmSubtractionMode mode ///< Mode for subtraction
    );

/// Set stamps based on a list of x,y
pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, ///< x coordinates for each stamp
                                               const psVector *y, ///< y coordinates for each stamp
                                               const psImage *image, ///< Image for flux of stamp
                                               const psImage *mask, ///< Mask, or NULL
                                               const psRegion *region, ///< Region to search, or NULL
                                               int footprint, ///< Half-size for stamps
                                               float spacing, ///< Rough spacing for stamps
                                               pmSubtractionMode mode ///< Mode for subtraction
    );

/// Set stamps based on a list of sources
pmSubtractionStampList *pmSubtractionStampsSetFromSources(
    const psArray *sources,             ///< Sources for each stamp
    const psImage *image,               ///< Image for flux of stamp
    const psImage *subMask,             ///< Mask, or NULL
    const psRegion *region,             ///< Region to search, or NULL
    int footprint,                      ///< Half-size for stamps
    float spacing,                      ///< Rough spacing for stamps
    pmSubtractionMode mode              ///< Mode for subtraction
    );

/// Set stamps based on values in a file
pmSubtractionStampList *pmSubtractionStampsSetFromFile(
    const char *filename,               ///< Filename of file containing x,y (or x,y,flux) on each line
    const psImage *image,               ///< Image for flux of stamp
    const psImage *subMask,             ///< Mask, or NULL
    const psRegion *region,             ///< Region to search, or NULL
    int footprint,                      ///< Half-size for stamps
    float spacing,                      ///< Rough spacing for stamps
    pmSubtractionMode mode              ///< Mode for subtraction
    );

/// Extract stamps from the images
bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, ///< Stamps
                                psImage *image1, ///< Reference image
                                psImage *image2, ///< Input image (or NULL)
                                psImage *variance, ///< Variance map
                                int kernelSize ///< Kernel half-size
    );

#endif
