Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c	(revision 26572)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c	(revision 26573)
@@ -238,23 +238,23 @@
 
     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);
-    }
-
+	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);
@@ -276,20 +276,20 @@
     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;
-        }
+	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");
+	fprintf (stderr, "no valid stamps?\n");
     }
 
@@ -297,71 +297,71 @@
 
     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++) {
+	// 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);
-    }
-
+	    
+	    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++) {
+	// 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);
-    }
-
+	    
+	    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;
@@ -390,12 +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;
-        }
+	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;
@@ -425,5 +425,5 @@
 bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {
 
-    // if (!pmVisualIsVisual()) return true;
+    if (!pmVisualIsVisual()) return true;
 
     // generate 4 storage images large enough to hold the N stamps:
@@ -433,5 +433,5 @@
     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;
@@ -449,5 +449,5 @@
     differenceImage  = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-
+    
     psImageInit (sourceImage,      0.0);
     psImageInit (targetImage,      0.0);
@@ -462,5 +462,5 @@
 bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {
 
-    // if (!pmVisualIsVisual()) return true;
+    if (!pmVisualIsVisual()) return true;
 
     double sum;
@@ -475,8 +475,8 @@
     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];
-        }
+	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;
@@ -485,8 +485,8 @@
     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];
-        }
+	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;
@@ -495,18 +495,18 @@
     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];
-        }
+	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];
-        }
+	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;
@@ -515,8 +515,8 @@
     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];
-        }
+	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;
@@ -524,7 +524,7 @@
     // 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));
-        }
+	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;
@@ -533,16 +533,15 @@
 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;
-
-    // XXX a dumb test : dump the residual image and exit
-    {
-        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 (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false;
     if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;
@@ -558,5 +557,27 @@
     pmVisualScaleImage(kapa2, fresidualImage, "Frac Residual Stamps", 2, true);
     pmVisualScaleImage(kapa2, differenceImage, "Difference Stamps", 0, true);
-    pmVisualScaleImage(kapa2, residualImage, "Residual Stamps", 1, 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);
 
@@ -603,5 +624,5 @@
     for (int i = 0; i < kernels->num; i++) {
         x->data.F32[i] = i;
-        y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false);
+	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]);