#ifndef PM_SUBTRACTION_KERNELS_H
#define PM_SUBTRACTION_KERNELS_H

#include <pslib.h>

/// Type of subtraction kernel
typedef enum {
    PM_SUBTRACTION_KERNEL_POIS,         ///< Pan-STARRS Optimal Image Subtraction --- delta functions
    PM_SUBTRACTION_KERNEL_ISIS,         ///< Traditional kernel --- gaussians modified by polynomials
    PM_SUBTRACTION_KERNEL_SPAM,         ///< Summed Pixels for Advanced Matching --- summed delta functions
    PM_SUBTRACTION_KERNEL_FRIES,        ///< Fibonacci Radius Increases Excellence of Subtraction
    PM_SUBTRACTION_KERNEL_GUNK,         ///< Grid United with Normal Kernel --- POIS and ISIS hybrid
    PM_SUBTRACTION_KERNEL_RINGS,        ///< Rings Instead of the Normal Gaussian Subtraction
} pmSubtractionKernelsType;

/// Kernels specification
typedef struct {
    pmSubtractionKernelsType type;      ///< Type of kernels --- allowing the use of multiple kernels
    psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS)
    psVector *widths;                   ///< Gaussian widths (ISIS) or ring radius (RINGS)
    psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
    psVector *xOrder, *yOrder;          ///< Spatial Polynomial order (for all)
    int subIndex;                       ///< Index of kernel to be subtracted (to maintain flux conservation)
    psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS)
    int size;                           ///< The half-size of the kernel
    int inner;                          ///< The size of an inner region
    int spatialOrder;                   ///< The spatial order of the kernels
} pmSubtractionKernels;

/// General allocator for pmSubtractionKernels
///
/// Unlike the functions for the specific kernel type, this function does not set up the basis functions, but
/// merely allocates space for their storage.
pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, ///< Number of basis functions
                                                pmSubtractionKernelsType type, ///< Kernel type
                                                int size, ///< Half-size of kernel
                                                int spatialOrder ///< Order of spatial variations
    );

/// Generate POIS kernels
pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims)
                                               int spatialOrder ///< Order of spatial variations
    );

/// Generate ISIS kernels
pmSubtractionKernels *pmSubtractionKernelsISIS(int size, ///< Half-size of the kernel
                                               int spatialOrder, ///< Order of spatial variations
                                               const psVector *sigmas, ///< Gaussian widths
                                               const psVector *orders ///< Polynomial order of gaussians
                                               );

/// Generate SPAM kernels
pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel
                                               int spatialOrder, ///< Order of spatial variations
                                               int inner, ///< Inner radius to preserve unbinned
                                               int binning ///< Kernel binning factor
    );

/// Generate FRIES kernels
pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel
                                                int spatialOrder, ///< Order of spatial variations
                                                int inner ///< Inner radius to preserve unbinned
    );

/// Generate GUNK kernels
pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, ///< Half-size of the kernel
                                               int spatialOrder, ///< Order of spatial variations
                                               const psVector *sigmas, ///< Gaussian widths
                                               const psVector *orders, ///< Polynomial order of gaussians
                                               int inner ///< Inner radius containing grid of delta functions
    );

/// Generate RINGS kernels
pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, ///< Half-size of the kernel
                                                int spatialOrder, ///< Order of spatial variations
                                                int inner, ///< Inner radius to preserve unbinned
                                                int ringsOrder ///< Polynomial order
    );


/// Generate a kernel of a specified type
pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, ///< Kernel type
                                                   int size, ///< Half-size of the kernel
                                                   int spatialOrder, ///< Order of spatial variations
                                                   const psVector *sigmas, ///< Gaussian widths
                                                   const psVector *orders, ///< Polynomial order of gaussians
                                                   int inner, ///< Inner radius to preserve unbinned
                                                   int binning, ///< Kernel binning factor
                                                   int ringsOrder ///< Polynomial order for RINGS
    );

#endif
