/** Diagnostic plots for pmSubtraction
 * @author Chris Beaumont, IfA
 */

/* Include Files   */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <strings.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <pslib.h>

#include "pmKapaPlots.h"
#include "pmSubtraction.h"
#include "pmSubtractionStamps.h"
#include "pmSubtractionEquation.h"

#include "pmVisual.h"

#include "pmHDU.h"
#include "pmFPA.h"
#include "pmAstrometryObjects.h"

# if (HAVE_KAPA)
# include <kapa.h>

//variables to determine when things are plotted
static bool plotConvKernels      = true;
static bool plotStamps           = true;
static bool plotLeastSquares     = true;
static bool plotImage            = true;

// variables to store plotting window indices
static int kapa1 = -1;
static int kapa2 = -1;
static int kapa3 = -1;

/** function prototypes*/
static bool plotStampLocations(pmSubtractionStampList *stamps, pmReadout *ro);

// Initialization Routines



/** destroy windows at the end of a run*/
bool pmSubtractionVisualClose(void)
{
    if(kapa1 != -1) KiiClose(kapa1);
    if(kapa2 != -1) KiiClose(kapa2);
    return true;
}

// Plotting Routines

/** Display images of the convolution kernels
 *  @param convKernels the kernels to plot
 *    @return true for success */
bool pmSubtractionVisualPlotConvKernels(psImage *convKernels) {
    if (!pmVisualIsVisual() || !plotConvKernels) return true;
    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) {
        return false;
    }
    pmVisualScaleImage(kapa1, convKernels, "Convolution_Kernels", 0, true);
    pmVisualAskUser(&plotConvKernels);
    return true;
}


/** Display the postage stamps used to determine the convolution kernels
    @param stamps to display
    @return true for success */
bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro) {
    if (!pmVisualIsVisual() || !plotStamps) return true;
    if (!pmVisualInitWindow (&kapa1, "ppSub:Images")) {
        return false;
    }
    if (!pmVisualInitWindow (&kapa2, "ppSub:StampMasterImage")) {
        return false;
    }

    //Plot Location of stamps
    plotStampLocations(stamps, ro);

    //Find the stamp size
    int imageMax = -1;
    int numStamps = 0;
    psElemType type = PS_TYPE_F32; // will be overwritten
    for(int i = 0; i < stamps->num; i++) {
        pmSubtractionStamp *stamp = stamps->stamps->data[i];
        if (stamp == NULL) continue;

        psKernel *k1 = stamp->image1;
        if (k1 == NULL) continue;

        psImage *i1 = k1->image;
        if (i1 == NULL) continue;

        imageMax = PS_MAX(imageMax, i1->numCols);
        imageMax = PS_MAX(imageMax, i1->numRows);
        numStamps++;
        type = i1->type.type;
    }
    if (imageMax == -1) return false;

    int border = 30;
    int tileRowCount = (int) ceil(sqrt(numStamps));
    int canvasX = tileRowCount * (imageMax * 2 + border);
    int canvasY = tileRowCount * (imageMax + border);
    psImage *canvas = psImageAlloc (canvasX, canvasY, type);
    psImageInit (canvas, NAN);

    //overlay the images
    int stampNum = 0;
    int stampListNum = 0;
    while (stampNum < numStamps) {
        int x0 = (2 * imageMax + border) * (stampNum % tileRowCount);
        int y0 = (imageMax + border) * (stampNum / tileRowCount);

        pmSubtractionStamp *stamp = stamps->stamps->data[stampListNum++];
        if (stamp == NULL) continue;

        psKernel *k1 = stamp->image1;
        psKernel *k2 = stamp->image2;

        if (k1 == NULL) continue;
        psImage *i1 = k1->image;
        psImage *i2 = k2->image;

        if (i1 == NULL) continue;
        psImageOverlaySection(canvas, i1, x0 + 0 * imageMax, y0, "=");
        psImageOverlaySection(canvas, i2, x0 + 1 * imageMax + border / 4, y0, "=");
        stampNum++;
    }
    pmVisualScaleImage(kapa1, canvas, "Subtraction_Stamps", 0, true);

    pmVisualAskUser(&plotStamps);
    return true;
}

/** Plot the least-squares matrix of each stamp */
bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps, bool dual) {

    if (!pmVisualIsVisual() || !plotLeastSquares) return true;
    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
        return false;
    }

    //Find the stamp size
    int imageMax = -1;
    int numStamps = 0;
    psElemType type = PS_TYPE_F64;

    for(int i = 0; i < stamps->num; i++) {
        pmSubtractionStamp *stamp = stamps->stamps->data[i];
        if (stamp == NULL) continue;

        psImage *im = stamp->matrix;
        if (im == NULL) continue;

        imageMax = PS_MAX(imageMax, im->numCols);
        imageMax = PS_MAX(imageMax, im->numRows);
        numStamps++;
        type = im->type.type;
    }
    if (imageMax == -1) return false;

    int border = 15;
    imageMax += border;
    int tileRowCount = (int) ceil(sqrt(numStamps));
    int canvasX = tileRowCount * (imageMax);
    int canvasY = tileRowCount * (imageMax);
    psImage *canvas = psImageAlloc (canvasX, canvasY, type);
    psImageInit (canvas, NAN);

    //overlay the images
    int stampNum = 0;
    int stampListNum = 0;
    while (stampNum < numStamps) {
        int x0 = (imageMax) * (stampNum % tileRowCount);
        int y0 = (imageMax) * (stampNum / tileRowCount);

        pmSubtractionStamp *stamp = stamps->stamps->data[stampListNum++];
        if (stamp == NULL) continue;

        psImage *im = stamp->matrix;
        if (im == NULL) continue;

        psImageOverlaySection(canvas, im, x0, y0, "=");
        stampNum++;
    }

    psImage *canvas32 = pmVisualImageToFloat(canvas);
    pmVisualScaleImage(kapa1, canvas32, "Least_Squares", 0, true);

    pmVisualAskUser(&plotLeastSquares);
    psFree(canvas);
    psFree(canvas32);
    return true;
}

bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {
    if (!pmVisualIsVisual() || !plotImage) return true;
    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
        return false;
    }

    pmVisualScaleImage(kapa1, image, "Image", 0, true);
    pmVisualScaleImage(kapa1, ref, "Reference", 1, true);
    pmVisualScaleImage(kapa1, sub, "Subtraction", 2, true);
    pmVisualAskUser(&plotImage);
    return true;
}

bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels) {

    if (!pmVisualIsVisual()) return true;
    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
        return false;
    }

    // get the kernel sizes
    int footprint = kernels->size;

    // output image is a grid of NXsub by NYsub sub-images
    int NXsub = sqrt(kernels->num);
    int NYsub = kernels->num / NXsub;
    if (kernels->num % NXsub) NYsub++;

    int NXpix = NXsub * (2*footprint + 1 + 3);
    int NYpix = NYsub * (2*footprint + 1 + 3);

    psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
    psImageInit (output, 0.0);

    for (int i = 0; i < kernels->num; i++) {
	pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i];
	psKernel *kernel = preCalc->kernel;

	int xSub = i % NXsub;
	int ySub = i / NXsub;

	int xPix = xSub * (2*footprint + 1 + 3) + footprint;
	int yPix = ySub * (2*footprint + 1 + 3) + footprint;

	double sum = 0.0;
	for (int y = -footprint; y <= footprint; y++) {
	    for (int x = -footprint; x <= footprint; x++) {
		output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
		sum += kernel->kernel[y][x];
	    }
	}
	fprintf (stderr, "kernel %d, sum %f\n", i, sum);
    }							 
	
    pmVisualScaleImage(kapa1, output, "Image", 0, true);
    pmVisualAskUser(&plotImage);
    return true;
}

bool pmSubtractionVisualShowBasis(pmSubtractionStampList *stamps) {

    if (!pmVisualIsVisual()) return true;
    if (!pmVisualInitWindow (&kapa2, "ppSub:StampMasterImage")) {
        return false;
    }

    // get the kernel sizes
    int footprint = stamps->footprint;

    // choose the brightest stamp
    pmSubtractionStamp *maxStamp = NULL;
    float maxFlux = NAN;
    for (int i = 0; i < stamps->num; i++) {
	pmSubtractionStamp *stamp = stamps->stamps->data[i];
	if (!isfinite(stamp->flux)) continue;
	if (!stamp->convolutions1 && !stamp->convolutions2) continue;
	if (!maxStamp) {
	    maxFlux = stamp->flux;
	    maxStamp = stamp;
	    continue;
	}
	if (stamp->flux > maxFlux) {
	    maxFlux = stamp->flux;
	    maxStamp = stamp;
	}
    }

    if (!isfinite(maxStamp->flux)) {
	fprintf (stderr, "no valid stamps?\n");
    }

    int nKernels = 0;

    if (maxStamp->convolutions1) {
	// output image is a grid of NXsub by NYsub sub-images
	nKernels = maxStamp->convolutions1->n;
	int NXsub = sqrt(nKernels);
	int NYsub = nKernels / NXsub;
	if (nKernels % NXsub) NYsub++;

	int NXpix = NXsub * (2*footprint + 1 + 3);
	int NYpix = NYsub * (2*footprint + 1 + 3);

	psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
	psImageInit (output, 0.0);

	for (int i = 0; i < nKernels; i++) {
            psKernel *kernel = maxStamp->convolutions1->data[i];
	    
	    int xSub = i % NXsub;
	    int ySub = i / NXsub;
	    
	    int xPix = xSub * (2*footprint + 1 + 3) + footprint;
	    int yPix = ySub * (2*footprint + 1 + 3) + footprint;
	    
	    double sum = 0.0;
	    for (int y = -footprint; y <= footprint; y++) {
		for (int x = -footprint; x <= footprint; x++) {
		    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
		    sum += kernel->kernel[y][x];
		}
	    }
	    fprintf (stderr, "kernel %d, sum %f\n", i, sum);
	}		
	pmVisualScaleImage(kapa2, output, "Image", 0, true);
    }					 
	
    if (maxStamp->convolutions2) {
	// output image is a grid of NXsub by NYsub sub-images
	nKernels = maxStamp->convolutions2->n;
	int NXsub = sqrt(nKernels);
	int NYsub = nKernels / NXsub;
	if (nKernels % NXsub) NYsub++;

	int NXpix = NXsub * (2*footprint + 1 + 3);
	int NYpix = NYsub * (2*footprint + 1 + 3);

	psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
	psImageInit (output, 0.0);

	for (int i = 0; i < nKernels; i++) {
            psKernel *kernel = maxStamp->convolutions2->data[i];
	    
	    int xSub = i % NXsub;
	    int ySub = i / NXsub;
	    
	    int xPix = xSub * (2*footprint + 1 + 3) + footprint;
	    int yPix = ySub * (2*footprint + 1 + 3) + footprint;
	    
	    double sum = 0.0;
	    for (int y = -footprint; y <= footprint; y++) {
		for (int x = -footprint; x <= footprint; x++) {
		    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
		    sum += kernel->kernel[y][x];
		}
	    }
	    fprintf (stderr, "kernel %d, sum %f\n", i, sum);
	}		
	pmVisualScaleImage(kapa2, output, "Image", 1, true);
    }					 
	
    pmVisualAskUser(&plotImage);
    return true;
}

static bool plotStampLocations(pmSubtractionStampList *stamps, pmReadout *ro) {

    if (!pmVisualScaleImage(kapa2, ro->image, "Stamp_master_image", 0, true)) {
        fprintf(stderr, "Cannot display postage stamp master image. Skipping \n");
        return false;
    }

    int Noverlay;
    KiiOverlay *overlay;

    // note: this uses the Ohana allocation tools:
    // ALLOCATE (overlay, KiiOverlay, 3*peaks->n + 1);
    ALLOCATE (overlay, KiiOverlay, stamps->num);

    Noverlay = 0;
    char boxID[10];
    for (int i = 0; i < stamps->num; i++) {

        pmSubtractionStamp *stamp = stamps->stamps->data[i];
        if (stamp == NULL) continue;

        overlay[Noverlay].type = KII_OVERLAY_BOX;
	if ((stamp->x < 1.0) && (stamp->y < 1.0)) {
	    // fprintf (stderr, "stamp zero: %f %f\n", stamp->x, stamp->y);
	    continue;
	}
	if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
	    // fprintf (stderr, "stamp nan: %f %f\n", stamp->x, stamp->y);
	    continue;
	}
        overlay[Noverlay].x = stamp->x;
        overlay[Noverlay].y = stamp->y;
        overlay[Noverlay].dx = 40.0;
        overlay[Noverlay].dy = 40.0;
        overlay[Noverlay].angle = 0.0;
        sprintf(boxID, "%d", i);
        overlay[Noverlay].text = boxID;
        Noverlay ++;
    }

    KiiLoadOverlay (kapa2, overlay, Noverlay, "red");
    FREE (overlay);
    return true;
}

static int footprint = 0;
static int NX = 0;
static int NY = 0;
static psImage *sourceImage      = NULL;
static psImage *targetImage      = NULL;
static psImage *residualImage    = NULL;
static psImage *fresidualImage   = NULL;
static psImage *differenceImage  = NULL;
static psImage *convolutionImage = NULL;

bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {

    if (!pmVisualIsVisual()) return true;

    // generate 4 storage images large enough to hold the N stamps:

    footprint = stamps->footprint;

    float NXf = sqrt(stamps->num);
    NX = (int) NXf == NXf ? NXf : NXf + 1.0;
    
    float NYf = stamps->num / NX;
    NY = (int) NYf == NY ? NYf : NYf + 1.0;

    int NXpix = (2*footprint + 1) * NX;
    NXpix += (NX > 1) ? 3 * NX : 0;

    int NYpix = (2*footprint + 1) * NY;
    NYpix += (NY > 1) ? 3 * NY : 0;

    sourceImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    targetImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    residualImage    = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    fresidualImage   = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    differenceImage  = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    
    psImageInit (sourceImage,      0.0);
    psImageInit (targetImage,      0.0);
    psImageInit (residualImage,    0.0);
    psImageInit (fresidualImage,   0.0);
    psImageInit (differenceImage,  0.0);
    psImageInit (convolutionImage, 0.0);

    return true;
}

bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {

    if (!pmVisualIsVisual()) return true;

    double sum;

    int NXoff = index % NX;
    int NYoff = index / NX;

    int NXpix = NXoff * (2*footprint + 1 + 3) + footprint;
    int NYpix = NYoff * (2*footprint + 1 + 3) + footprint;

    // insert the (target) kernel into the (target) image:
    sum = 0.0;
    for (int y = -footprint; y <= footprint; y++) {
	for (int x = -footprint; x <= footprint; x++) {
	    targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x];
	    sum += targetImage->data.F32[y + NYpix][x + NXpix];
	}
    }
    targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;

    // insert the (source) kernel into the (source) image:
    sum = 0.0;
    for (int y = -footprint; y <= footprint; y++) {
	for (int x = -footprint; x <= footprint; x++) {
	    sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x];
	    sum += sourceImage->data.F32[y + NYpix][x + NXpix];
	}
    }
    sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;

    // insert the (convolution) kernel into the (convolution) image:
    sum = 0.0;
    for (int y = -footprint; y <= footprint; y++) {
	for (int x = -footprint; x <= footprint; x++) {
	    convolutionImage->data.F32[y + NYpix][x + NXpix] = -convolution->kernel[y][x];
	    sum += convolutionImage->data.F32[y + NYpix][x + NXpix];
	}
    }
    convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    
    // insert the (difference) kernel into the (difference) image:
    sum = 0.0;
    for (int y = -footprint; y <= footprint; y++) {
	for (int x = -footprint; x <= footprint; x++) {
	    differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm;
	    sum += differenceImage->data.F32[y + NYpix][x + NXpix];
	}
    }
    differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;

    // insert the (residual) kernel into the (residual) image:
    sum = 0.0;
    for (int y = -footprint; y <= footprint; y++) {
	for (int x = -footprint; x <= footprint; x++) {
	    residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm - convolution->kernel[y][x];
	    sum += residualImage->data.F32[y + NYpix][x + NXpix];
	}
    }
    residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;

    // insert the (fresidual) kernel into the (fresidual) image:
    for (int y = -footprint; y <= footprint; y++) {
	for (int x = -footprint; x <= footprint; x++) {
	    fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0));
	}
    }
    return true;
}

bool pmSubtractionVisualShowFit(double norm) {

    // for testing, dump the residual image and exit
    if (0) {
	psMetadata *header = psMetadataAlloc();
        psMetadataAddF32 (header, PS_LIST_TAIL, "NORM", 0, "Normalization", norm);
	psFits *fits = psFitsOpen("resid.fits", "w");
	psFitsWriteImage(fits, header, residualImage, 0, NULL);
	psFitsClose(fits);
	// exit (0);
    }

    if (!pmVisualIsVisual()) return true;
    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false;
    if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;

    KiiEraseOverlay (kapa1, "red");
    KiiEraseOverlay (kapa2, "red");

    pmVisualScaleImage(kapa1, targetImage, "Target Stamps", 0, true);
    pmVisualScaleImage(kapa1, sourceImage, "Source Stamps", 1, true);
    pmVisualScaleImage(kapa1, convolutionImage, "Convolution Stamps", 2, true);
    KiiCenter (kapa1, 0.5*targetImage->numCols, 0.5*targetImage->numRows, 1);

    pmVisualScaleImage(kapa2, fresidualImage, "Frac Residual Stamps", 2, true);
    pmVisualScaleImage(kapa2, differenceImage, "Difference Stamps", 0, true);

    if (1) {
	KiiImage image;
	KapaImageData data;
	Coords coords;
	strcpy (coords.ctype, "RA---TAN");

	image.data2d = residualImage->data.F32;
	image.Nx = residualImage->numCols;
	image.Ny = residualImage->numRows;
	strcpy (data.name, "Residual Stamps");
	strcpy (data.file, "Residual Stamps");

	data.zero  = -300.0;
	data.range = +600.0;
	data.logflux = 0;

	KiiSetChannel (kapa2, 1);
	KiiNewPicture2D (kapa2, &image, &data, &coords);
    } else {
	pmVisualScaleImage(kapa2, residualImage, "Residual Stamps", 1, true);
    }

    KiiCenter (kapa2, 0.5*residualImage->numCols, 0.5*residualImage->numRows, 1);

    pmVisualAskUser(NULL);

    psFree(targetImage);
    psFree(sourceImage);
    psFree(convolutionImage);
    psFree(differenceImage);
    psFree(residualImage);
    psFree(fresidualImage);

    targetImage = NULL;
    sourceImage = NULL;
    convolutionImage = NULL;
    differenceImage = NULL;
    residualImage = NULL;
    fresidualImage = NULL;

    return true;
}

bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels) {

    Graphdata graphdata;

    if (!pmVisualIsVisual()) return true;
    if (!pmVisualInitWindow(&kapa3, "ppSub:plots")) return false;

    KapaClearSections (kapa3);
    KapaInitGraph (&graphdata);

    psVector *x = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
    psVector *y = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);

    graphdata.xmin = -1.0;
    graphdata.xmax = kernels->num + 1.0;
    graphdata.ymin = +32.0;
    graphdata.ymax = -32.0;

    psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, 0.0, 0.0);

    // construct the plot vectors
    for (int i = 0; i < kernels->num; i++) {
        x->data.F32[i] = i;
	y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false);
        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[i]);
        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[i]);
    }
    x->n = y->n = kernels->num;

    float range;
    range = graphdata.xmax - graphdata.xmin;
    graphdata.xmax += 0.05*range;
    graphdata.xmin -= 0.05*range;
    range = graphdata.ymax - graphdata.ymin;
    graphdata.ymax += 0.05*range;
    graphdata.ymin -= 0.05*range;

    KapaSetLimits (kapa3, &graphdata);

    KapaSetFont (kapa3, "helvetica", 14);
    KapaBox (kapa3, &graphdata);
    KapaSendLabel (kapa3, "kernel number", KAPA_LABEL_XM);
    KapaSendLabel (kapa3, "coeff", KAPA_LABEL_YM);

    graphdata.color = KapaColorByName ("black");
    graphdata.ptype = 2;
    graphdata.size = 0.5;
    graphdata.style = 2;

    KapaPrepPlot   (kapa3, x->n, &graphdata);
    KapaPlotVector (kapa3, x->n, x->data.F32, "x");
    KapaPlotVector (kapa3, x->n, y->data.F32, "y");

    psFree (x);
    psFree (y);

    pmVisualAskUser(NULL);
    return true;
}

#else
bool pmSubtractionVisualClose(void) {return true;}
bool pmSubtractionVisualPlotConvKernels(psImage *convKernels) {return true;}
bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro) {return true;}
bool pmSubtractionVisualPlotLeastSquares(pmSubtractionStampList *stamps) {return true;}
bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {return true;}
bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {return true;}
bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {return true;}
bool pmSubtractionVisualShowFit() {return true;}
bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels);
#endif
