/** @file  ImageSubtract.h
 *
 *  This file will contain code which creates a set of kernel basis
 *  functions, solves for their solution, and applies them to an image.
 *
 *  @author Paul Price, IfA (original prototype)
 *  @author GLG, MHPCC
 *
 *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-09-28 20:43:52 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

#if !defined(PM_IMAGE_COMBINE_H)
#define PM_IMAGE_COMBINE_H

#if HAVE_CONFIG_H
#include <config.h>
#endif

#include<stdio.h>
#include<math.h>
#include "pslib.h"
#include "psConstants.h"

typedef enum {
    PM_SUBTRACTION_KERNEL_POIS,         ///< POIS kernel --- delta functions
    PM_SUBTRACTION_KERNEL_ISIS          ///< ISIS kernel --- gaussians modified by polynomials
} pmSubtractionKernelsType;

typedef struct
{
    pmSubtractionKernelsType type;      ///< Type ofKernels --- 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 (to
    ///< accelerate ISIS; don't use for POIS)
    psS32 size;                         ///< The halfsize of the kernel
    psS32 spatialOrder;                 ///< The spatial order of the kernels
}
psSubtractionKernels;

psSubtractionKernels *pmSubtractionKernelsAllocPOIS(
    int size,
    int SpatialOrder
);

psSubtractionKernels *pmSubtractionKernelsAllocISIS(
    const psVector *sigmas,
    const psVector *orders,
    int size,
    int SpatialOrder
);

typedef enum {
    PM_STAMP_INIT,                      ///< Initial state
    PM_STAMP_USED,                      ///< Use this stamp
    PM_STAMP_REJECTED,                  ///< This stamp has been rejected
    PM_STAMP_RECALC,                    ///< Having been reset, this stamp is to be recalculated
    PM_STAMP_NONE                       ///< No stamp in this region
} pmStampStatus;

typedef struct
{
    int x, y;                           ///< Position
    int p_xSize;
    int p_ySize;
    int p_xMin;
    int p_xMax;
    int p_yMin;
    int p_yMax;
    psImage *matrix;                    ///< Associated matrix
    psVector *vector;                   ///< Assoicated vector
    pmStampStatus status;               ///< Status ofstamp
}
pmStamp;

psArray *pmSubtractionFindStamps(
    psArray *stamps,                    ///< Output stamps, or NULL
    const psImage *image,               ///< Image for which to find stamps
    const psImage *mask,                ///< Mask
    psU32 maskVal,                      ///< Value for mask
    psF32 threshold,                    ///< Threshold for stamps in the image
    psS32 xNum,                         ///< Number of stamps in x
    psS32 yNum,                         ///< Number of stamps in y
    psS32 border                        ///< Border around image to ignore (should be size of kernel)
);

bool pmSubtractionCalculateEquation(
    psArray *stamps,                    ///< The stamps for which to calculate the equation,
    const psImage *reference,           ///< Reference image
    const psImage *input,               ///< Input image
    const psSubtractionKernels *kernels,///< The kernel basis functions
    psS32 footprint                     ///< Half-size of region over which to calculate equation
);


psVector *pmSubtractionSolveEquation(
    psVector *solution,                 ///< Solution vector, or NULL
    const psArray *stamps               ///< Array of stamps
);

bool pmSubtractionRejectStamps(
    psArray *stamps,                    ///< Array of stamps to check for rejection
    psImage *mask,                      ///< Mask image
    psU32 badStampMaskVal,              ///< Value to use in mask for bad stamp
    psS32 footprint,                    ///< Region to mask if stamp is bad
    psF32 sigmaRej,                     ///< Number of RMS deviations above zero at which to reject
    const psImage *refImage,            ///< Reference image
    const psImage *inImage,             ///< Input image
    const psVector *solution,           ///< Solution vector
    const psSubtractionKernels *kernels ///< Array of kernel parameters
);

psImage *pmSubtractionKernelImage(
    psImage *out,
    const psVector *solution,
    const psSubtractionKernels *kernels,
    psF32 x,
    psF32 y
);

#endif
