Index: trunk/psModules/src/imcombine/pmSubtractionVisual.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionVisual.c	(revision 23487)
+++ trunk/psModules/src/imcombine/pmSubtractionVisual.c	(revision 26893)
@@ -16,5 +16,7 @@
 
 #include "pmKapaPlots.h"
+#include "pmSubtraction.h"
 #include "pmSubtractionStamps.h"
+#include "pmSubtractionEquation.h"
 
 #include "pmVisual.h"
@@ -32,7 +34,9 @@
 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,8 +50,6 @@
 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,8 +62,8 @@
 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,5 +76,5 @@
 bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro) {
     if (!pmVisualIsVisual() || !plotStamps) return true;
-    if (!pmVisualInitWindow (&kapa, "ppSub:Images")) {
+    if (!pmVisualInitWindow (&kapa1, "ppSub:Images")) {
         return false;
     }
@@ -134,5 +136,5 @@
         stampNum++;
     }
-    pmVisualScaleImage(kapa, canvas, "Subtraction_Stamps", 0, true);
+    pmVisualScaleImage(kapa1, canvas, "Subtraction_Stamps", 0, true);
 
     pmVisualAskUser(&plotStamps);
@@ -144,5 +146,5 @@
 
     if (!pmVisualIsVisual() || !plotLeastSquares) return true;
-    if (!pmVisualInitWindow (&kapa, "PPSub:Images")) {
+    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
         return false;
     }
@@ -157,7 +159,5 @@
         if (stamp == NULL) continue;
 
-        psImage *im = stamp->matrix1;
-        if (im == NULL) im = stamp->matrix2;
-        if (im == NULL) im = stamp->matrixX;
+        psImage *im = stamp->matrix;
         if (im == NULL) continue;
 
@@ -187,7 +187,5 @@
         if (stamp == NULL) continue;
 
-        psImage *im = stamp->matrix1;
-        if (im == NULL) im = stamp->matrix2;
-        if (im == NULL) im = stamp->matrixX;
+        psImage *im = stamp->matrix;
         if (im == NULL) continue;
 
@@ -197,5 +195,5 @@
 
     psImage *canvas32 = pmVisualImageToFloat(canvas);
-    pmVisualScaleImage(kapa, canvas32, "Least_Squares", 0, true);
+    pmVisualScaleImage(kapa1, canvas32, "Least_Squares", 0, true);
 
     pmVisualAskUser(&plotLeastSquares);
@@ -207,11 +205,163 @@
 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;
+}
+
+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;
@@ -240,4 +390,12 @@
 
         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;
@@ -255,4 +413,252 @@
 }
 
+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;}
@@ -261,3 +667,7 @@
 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