#ifndef PM_SUBTRACTION_EQUATION_H
#define PM_SUBTRACTION_EQUATION_H

#include "pmSubtractionStamps.h"
#include "pmSubtractionKernels.h"

typedef enum {
    PM_SUBTRACTION_EQUATION_NONE    = 0x00,
    PM_SUBTRACTION_EQUATION_NORM    = 0x01,
    PM_SUBTRACTION_EQUATION_BG      = 0x02,
    PM_SUBTRACTION_EQUATION_KERNELS = 0x04,
    PM_SUBTRACTION_EQUATION_ALL     = 0x07, // value should be NORM | BG | KERNELS
} pmSubtractionEquationCalculationMode;

/// Execute a thread job to calculate the least-squares equation for a stamp
bool pmSubtractionCalculateEquationThread(psThreadJob *job ///< Job to execute
    );

/// Calculate the least-squares equation to match the image quality for a single stamp
bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, ///< Stamps
                                         const pmSubtractionKernels *kernels, ///< Kernel parameters
                                         int index, ///< Index of stamp
                                         const pmSubtractionEquationCalculationMode mode
    );

/// Calculate the least-squares equation to match the image quality
bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps
                                    const pmSubtractionKernels *kernels, ///< Kernel parameters
                                    const pmSubtractionEquationCalculationMode mode
    );

/// Solve the least-squares equation to match the image quality
bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters
                                const pmSubtractionStampList *stamps, ///< Stamps
                                const pmSubtractionEquationCalculationMode mode
    );

/// Calculate deviations
psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, ///< Stamps
                                           pmSubtractionKernels *kernels ///< Kernel parameters
    );

/// Calculate the value of a polynomial, specified by coefficients and polynomial values
double p_pmSubtractionCalculatePolynomial(const psVector *coeff, ///< Coefficients
                                          const psImage *polyValues, ///< Polynomial values
                                          int order, ///< Order of polynomials
                                          int index, ///< Index at which to begin
                                          int step ///< Step between subsequent indices
    );

/// Return the specified coefficient in the solution
double p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, ///< Kernel parameters
                                    const psImage *polyValues, ///< Polynomial values
                                    int index, ///< Coefficient index to calculate
                                    bool wantDual ///< Calculate the coefficient for the dual solution?
    );

/// Return the normalisation in the solution
double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels ///< Kernel parameters
    );

/// Return the background (difference) in the solution
double p_pmSubtractionSolutionBackground(const pmSubtractionKernels *kernels, ///< Kernel parameters
                                         const psImage *polyValues ///< Polynomial values
    );

#endif
