Changeset 26321 for branches/pap/psModules/src/imcombine

Timestamp:

Dec 2, 2009, 12:04:51 PM (16 years ago)

Author:

Paul Price

Message:

Merging Gene's work on adding a window function to the stamps. I'd been thinking about something like this for a while, so glad someone's done it. It's got visualisation turned on though.

Location:

branches/pap/psModules/src/imcombine

Files:

• . (modified) (1 prop)
• pmSubtractionEquation.c (modified) (11 diffs)
• pmSubtractionKernels.c (modified) (1 diff)
• pmSubtractionMatch.c (modified) (3 diffs)
• pmSubtractionStamps.c (modified) (4 diffs)
• pmSubtractionStamps.h (modified) (2 diffs)
• pmSubtractionVisual.c (modified) (11 diffs)
• pmSubtractionVisual.h (modified) (1 diff)

branches/pap/psModules/src/imcombine/pmSubtractionEquation.c

@@ -18 +18 @@
 
 #define USE_WEIGHT                      // Include weight (1/variance) in equation?
+#define USE_WINDOW                      // Include weight (1/variance) in equation?
 
 
@@ -30 +31 @@
                                   const psKernel *reference, // Reference image (convolution source)
                                   const psKernel *weight,  // Weight image
+                                  const psKernel *window,  // Window image
                                   const psArray *convolutions,         // Convolutions for each kernel
                                   const pmSubtractionKernels *kernels, // Kernels
@@ -76 +78 @@
 #ifdef USE_WEIGHT
                     cc *= weight->kernel[y][x];
+#endif
+#ifdef USE_WINDOW
+                    if (window) {
+                        cc *= window->kernel[y][x];
+                    }
 #endif
                     sumCC += cc;
@@ -107 +114 @@
                 rc *= wtVal;
                 c *= wtVal;
+#endif
+#ifdef USE_WINDOW
+                if (window) {
+                    float winVal = window->kernel[y][x];
+                    ic *= winVal;
+                    rc *= winVal;
+                    c  *= winVal;
+                }
 #endif
                 sumIC += ic;
@@ -144 +159 @@
             ref *= wtVal;
             one *= wtVal;
+#endif
+#ifdef USE_WINDOW
+            if (window) {
+                float  winVal = window->kernel[y][x];
+                rr      *= winVal;
+                ir      *= winVal;
+                in      *= winVal;
+                ref *= winVal;
+                one *= winVal;
+            }
 #endif
             sumRR += rr;
@@ -603 +628 @@
 #endif
 
+    // XXX visualize the set of convolved stamps
+
     psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder,
                                                     stamp->xNorm, stamp->yNorm); // Polynomial terms
@@ -620 +647 @@
       case PM_SUBTRACTION_MODE_1:
         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1,
-                                       stamp->weight, stamp->convolutions1, kernels, polyValues,
-                                       footprint);
+                                       stamp->weight, stamps->window, stamp->convolutions1, kernels,
+                                       polyValues, footprint);
         break;
       case PM_SUBTRACTION_MODE_2:
         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2,
-                                       stamp->weight, stamp->convolutions2, kernels, polyValues,
-                                       footprint);
+                                       stamp->weight, stamps->window, stamp->convolutions2, kernels,
+                                       polyValues, footprint);
         break;
       case PM_SUBTRACTION_MODE_DUAL:
@@ -1147 +1174 @@
     double devNorm = 1.0 / (double)numPixels; // Normalisation for deviations
     int numKernels = kernels->num;      // Number of kernels
+    double norm = NAN;
 
     psImage *polyValues = NULL;         // Polynomial values
     psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+
+    // set up holding images for the visualization
+    pmSubtractionVisualShowFitInit (stamps);
 
     for (int i = 0; i < stamps->num; i++) {
@@ -1160 +1191 @@
         // Calculate coefficients of the kernel basis functions
         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
-        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+        norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
 
@@ -1214 +1245 @@
                 }
             }
+
+            // XXX visualize the target, source, convolution and residual
+            pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
+
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
@@ -1300 +1335 @@
         }
 #endif
-
+    
     }
     psFree(residual);
     psFree(polyValues);
 
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f", norm);
+    pmSubtractionVisualShowFit();
+    pmSubtractionVisualPlotFit(kernels);
+
     return deviations;
 }

branches/pap/psModules/src/imcombine/pmSubtractionKernels.c

@@ -111 +111 @@
 
 pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder,
-                                                    const psVector *fwhms, const psVector *orders,
+                                                    const psVector *fwhmsIN, const psVector *ordersIN,
                                                     float penalty, pmSubtractionMode mode)
 {
-    PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL);
-    PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(orders, NULL);
-    PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(fwhms, orders, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(fwhmsIN, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(ordersIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(ordersIN, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(fwhmsIN, ordersIN, NULL);
     PS_ASSERT_INT_POSITIVE(size, NULL);
     PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
+
+    // check the requested fwhm values: any values <= 0.0 should be dropped
+    psVector *fwhms  = psVectorAllocEmpty (fwhmsIN->n, PS_TYPE_F32);
+    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
+    for (int i = 0; i < fwhmsIN->n; i++) {
+        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
+        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
+        psVectorAppend(orders, ordersIN->data.S32[i]);
+    }
 
     int numGaussians = fwhms->n;       // Number of Gaussians

branches/pap/psModules/src/imcombine/pmSubtractionMatch.c

@@ -483 +483 @@
             }
 
+
+            // generate the window function from the set of stamps
+            if (!pmSubtractionStampsGetWindow(stamps, size)) {
+                goto MATCH_ERROR;
+            }
 
             // Define kernel basis functions
@@ -550 +555 @@
                 }
 
+                // generate the window function from the set of stamps
+                if (!pmSubtractionStampsGetWindow(stamps, size)) {
+                    goto MATCH_ERROR;
+                }
+
                 psTrace("psModules.imcombine", 3, "Calculating equation...\n");
                 if (!pmSubtractionCalculateEquation(stamps, kernels)) {
@@ -587 +597 @@
             }
 
+            // if we hit the max number of iterations and we have rejected stamps, re-solve
             if (numRejected > 0) {
                 psTrace("psModules.imcombine", 3, "Solving equation...\n");

branches/pap/psModules/src/imcombine/pmSubtractionStamps.c

@@ -46 +46 @@
     psFree(list->y);
     psFree(list->flux);
+    psFree(list->window);
 }
 
@@ -220 +221 @@
     list->y = NULL;
     list->flux = NULL;
+    list->window = NULL;
     list->footprint = footprint;
     list->sysErr = sysErr;
@@ -239 +241 @@
     out->y = NULL;
     out->flux = NULL;
+    out->window = psMemIncrRefCounter(in->window);
     out->footprint = in->footprint;
 
@@ -606 +609 @@
 }
 
+
+bool pmSubtractionStampsGetWindow(pmSubtractionStampList *stamps, int kernelSize)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
+    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
+
+    int size = kernelSize + stamps->footprint; // Size of postage stamps
+
+    psFree (stamps->window);
+    stamps->window = psKernelAlloc(-size, size, -size, size);
+    psImageInit(stamps->window->image, 0.0);
+
+    // generate normalizations for each stamp
+    psVector *norm1 = psVectorAlloc(stamps->num, PS_TYPE_F32);
+    psVector *norm2 = psVectorAlloc(stamps->num, PS_TYPE_F32);
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (!stamp) continue;
+        if (!stamp->image1) continue;
+        if (!stamp->image2) continue;
+
+        float sum1 = 0.0;
+        float sum2 = 0.0;
+        for (int y = -size; y <= size; y++) {
+            for (int x = -size; x <= size; x++) {
+                sum1 += stamp->image1->kernel[y][x];
+                sum2 += stamp->image2->kernel[y][x];
+            }
+        }
+        norm1->data.F32[i] = sum1;
+        norm2->data.F32[i] = sum2;
+    }
+
+    // storage vector for flux data
+    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psVector *flux = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
+
+    // generate the window pixels
+    for (int y = -size; y <= size; y++) {
+        for (int x = -size; x <= size; x++) {
+
+            flux->n = 0;
+            for (int i = 0; i < stamps->num; i++) {
+                pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+                if (!stamp) continue;
+                if (!stamp->image1) continue;
+                if (!stamp->image2) continue;
+                
+                psVectorAppend (flux, stamp->image1->kernel[y][x] / norm1->data.F32[i]);
+                psVectorAppend (flux, stamp->image2->kernel[y][x] / norm2->data.F32[i]);
+            }
+
+            psStatsInit (stats);
+            if (!psVectorStats (stats, flux, NULL, NULL, 0)) {
+                psAbort ("failed to generate stats");
+            }
+            stamps->window->kernel[y][x] = stats->robustMedian;
+        }
+    }
+
+    // XXX TESTING: 
+    {
+        psFits *fits = psFitsOpen ("window.fits", "w");
+        psFitsWriteImage (fits, NULL, stamps->window->image, 0, NULL);
+        psFitsClose (fits);
+    }
+
+    psFree (stats);
+    psFree (flux);
+    psFree (norm1);
+    psFree (norm2);
+    return true;
+}
 
 bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,

branches/pap/psModules/src/imcombine/pmSubtractionStamps.h

@@ -24 +24 @@
     psArray *flux;                      ///< Fluxes for possible stamps (or NULL)
     int footprint;                      ///< Half-size of stamps
+    psKernel *window;                   ///< window function generated from ensemble of stamps
     float sysErr;                       ///< Systematic error
 } pmSubtractionStampList;
@@ -136 +137 @@
     );
 
+bool pmSubtractionStampsGetWindow(pmSubtractionStampList *stamps, int kernelSize);
+
 /// Extract stamps from the images
 bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, ///< Stamps

branches/pap/psModules/src/imcombine/pmSubtractionVisual.c

@@ -16 +16 @@
 
 #include "pmKapaPlots.h"
+#include "pmSubtraction.h"
 #include "pmSubtractionStamps.h"
+#include "pmSubtractionEquation.h"
 
 #include "pmVisual.h"
@@ -32 +34 @@
 static bool plotLeastSquares     = true;
 static bool plotImage            = true;
+
 // variables to store plotting window indices
-static int kapa  = -1;
+static int kapa1 = -1;
 static int kapa2 = -1;
+static int kapa3 = -1;
 
 /** function prototypes*/
@@ -46 +50 @@
 bool pmSubtractionVisualClose(void)
 {
-    if(kapa != -1)
-        KiiClose(kapa);
-    if(kapa2 != -1)
-        KiiClose(kapa2);
+    if(kapa1 != -1) KiiClose(kapa1);
+    if(kapa2 != -1) KiiClose(kapa2);
     return true;
 }
@@ -60 +62 @@
 bool pmSubtractionVisualPlotConvKernels(psImage *convKernels) {
     if (!pmVisualIsVisual() || !plotConvKernels) return true;
-    if (!pmVisualInitWindow(&kapa, "ppSub:Images")) {
-        return false;
-    }
-    pmVisualScaleImage(kapa, convKernels, "Convolution_Kernels", 0, true);
+    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) {
+        return false;
+    }
+    pmVisualScaleImage(kapa1, convKernels, "Convolution_Kernels", 0, true);
     pmVisualAskUser(&plotConvKernels);
     return true;
@@ -74 +76 @@
 bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro) {
     if (!pmVisualIsVisual() || !plotStamps) return true;
-    if (!pmVisualInitWindow (&kapa, "ppSub:Images")) {
+    if (!pmVisualInitWindow (&kapa1, "ppSub:Images")) {
         return false;
     }
@@ -134 +136 @@
         stampNum++;
     }
-    pmVisualScaleImage(kapa, canvas, "Subtraction_Stamps", 0, true);
+    pmVisualScaleImage(kapa1, canvas, "Subtraction_Stamps", 0, true);
 
     pmVisualAskUser(&plotStamps);
@@ -144 +146 @@
 
     if (!pmVisualIsVisual() || !plotLeastSquares) return true;
-    if (!pmVisualInitWindow (&kapa, "PPSub:Images")) {
+    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
         return false;
     }
@@ -197 +199 @@
 
     psImage *canvas32 = pmVisualImageToFloat(canvas);
-    pmVisualScaleImage(kapa, canvas32, "Least_Squares", 0, true);
+    pmVisualScaleImage(kapa1, canvas32, "Least_Squares", 0, true);
 
     pmVisualAskUser(&plotLeastSquares);
@@ -207 +209 @@
 bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {
     if (!pmVisualIsVisual() || !plotImage) return true;
-    if (!pmVisualInitWindow (&kapa, "PPSub:Images")) {
-        return false;
-    }
-
-    pmVisualScaleImage(kapa, image, "Image", 0, true);
-    pmVisualScaleImage(kapa, ref, "Reference", 1, true);
-    pmVisualScaleImage(kapa, sub, "Subtraction", 2, 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;
@@ -255 +257 @@
 }
 
+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 *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);
+    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 (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;
+
+    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:
+    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];
+        }
+    }
+
+    // insert the (source) kernel into the (source) image:
+    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];
+        }
+    }
+    
+    // insert the (convolution) kernel into the (convolution) image:
+    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];
+        }
+    }
+    
+    // insert the (difference) kernel into the (difference) image:
+    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;
+        }
+    }
+
+    // insert the (residual) kernel into the (residual) image:
+    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];
+        }
+    }
+    return true;
+}
+
+bool pmSubtractionVisualShowFit() {
+
+    // if (!pmVisualIsVisual()) return true;
+    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false;
+    if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;
+
+    pmVisualScaleImage(kapa1, targetImage, "Target Stamps", 0, true);
+    pmVisualScaleImage(kapa1, sourceImage, "Source Stamps", 1, true);
+    pmVisualScaleImage(kapa1, convolutionImage, "Convolution Stamps", 2, true);
+
+    pmVisualScaleImage(kapa2, differenceImage, "Difference Stamps", 0, true);
+    pmVisualScaleImage(kapa2, residualImage, "Residual Stamps", 1, true);
+    pmVisualAskUser(NULL);
+
+    psFree(targetImage);
+    psFree(sourceImage);
+    psFree(convolutionImage);
+    psFree(differenceImage);
+    psFree(residualImage);
+
+    targetImage = NULL;
+    sourceImage = NULL;
+    convolutionImage = NULL;
+    differenceImage = NULL;
+    residualImage = 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;}
@@ -261 +443 @@
 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

branches/pap/psModules/src/imcombine/pmSubtractionVisual.h

@@ -7 +7 @@
 bool pmSubtractionVisualPlotLeastSquares(pmSubtractionStampList *stamps);
 bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub);
+bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps);
+bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index);
+bool pmSubtractionVisualShowFit();
+bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels);
 
 #endif