#ifndef PM_SUBTRACTION_KERNELS_H
#define PM_SUBTRACTION_KERNELS_H

#include <pslib.h>

/// Type of subtraction kernel
typedef enum {
    PM_SUBTRACTION_KERNEL_POIS,         ///< POIS kernel --- delta functions
    PM_SUBTRACTION_KERNEL_ISIS          ///< ISIS kernel --- gaussians modified by polynomials
} 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 *sigma;                    ///< Width of Gaussian (for ISIS)
    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 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 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
                                                   );

#endif
