Index: trunk/psModules/src/imcombine/Makefile.am
===================================================================
--- trunk/psModules/src/imcombine/Makefile.am	(revision 35768)
+++ trunk/psModules/src/imcombine/Makefile.am	(revision 35771)
@@ -20,4 +20,5 @@
 	pmSubtractionStamps.c	\
 	pmSubtractionThreads.c	\
+	pmSubtractionSimple.c   \
 	pmPSFEnvelope.c         \
 	pmSubtractionVisual.c   \
@@ -42,4 +43,5 @@
 	pmSubtractionStamps.h	\
 	pmSubtractionThreads.h	\
+	pmSubtractionSimple.h   \
 	pmPSFEnvelope.h         \
 	pmSubtractionVisual.h   \
Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 35768)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 35771)
@@ -167,4 +167,9 @@
               break;
           }
+	case PM_SUBTRACTION_KERNEL_SIMPLE: {
+              solvedKernelPreCalc(kernel, kernels, 1.0, i);
+              break;
+	}	  
+	  
           case PM_SUBTRACTION_KERNEL_ISIS:
           case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 35768)
+++ trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 35771)
@@ -65,4 +65,25 @@
     return result;
 }
+
+// Generate 1D convolution kernel for SIMPLE
+psVector *pmSubtractionKernelSIMPLE(float sigma, // Gaussian width
+				    int order,   // Unused polynomial order
+				    int size     // Kernel half-size
+				    )
+{
+  int fullSize = 2 * size + 1;    // Full size of the kernel
+  psVector *kernel = psVectorAlloc(fullSize, PS_TYPE_F32); // Kernel to return
+  float expNorm = -0.5 / PS_SQR(sigma); // Normalization for exponential
+  float norm    = 1.0 / sqrtf(2.0 * M_PI * sigma * sigma); // Correct Normalization for Gaussian
+  if (sigma < 0.1) {
+    kernel->data.F32[size] = 1.0;
+    return(kernel);
+  }
+  for (int i = 0, x = -size; x <= size; i++,x++) {
+    kernel->data.F32[i] = norm * expf(expNorm * PS_SQR(x));
+  }
+  return(kernel);
+}
+
 
 // Generate 1D convolution kernel for ISIS
@@ -771,4 +792,11 @@
         preCalc->poly    = NULL;
         break;
+    case PM_SUBTRACTION_KERNEL_SIMPLE:
+      preCalc->xKernel = pmSubtractionKernelSIMPLE(sigma,uOrder,size);
+      preCalc->yKernel = pmSubtractionKernelSIMPLE(sigma,vOrder,size);
+      preCalc->uCoords = NULL;
+      preCalc->vCoords = NULL;
+      preCalc->poly    = NULL;
+      break;
       case PM_SUBTRACTION_KERNEL_RINGS:
         // the RINGS kernel uses the uCoords, vCoords, and poly elements of the structure
@@ -1230,4 +1258,6 @@
       case PM_SUBTRACTION_KERNEL_RINGS:
         return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, bounds, mode);
+    case PM_SUBTRACTION_KERNEL_SIMPLE:
+      return pmSubtractionKernelsISIS(size,spatialOrder,fwhms,orders,penalty,bounds,mode);
       default:
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);
@@ -1298,4 +1328,5 @@
       case PM_SUBTRACTION_KERNEL_HERM:
       case PM_SUBTRACTION_KERNEL_DECONV_HERM:
+    case PM_SUBTRACTION_KERNEL_SIMPLE:
         PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
 
@@ -1379,5 +1410,7 @@
         return PM_SUBTRACTION_KERNEL_RINGS;
     }
-
+    if (strncasecmp(type, "SIMPLE", nameLength) == 0) {
+      return PM_SUBTRACTION_KERNEL_SIMPLE;
+    }
     psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", type);
     return PM_SUBTRACTION_KERNEL_NONE;
@@ -1436,4 +1469,7 @@
 	break;
 
+    case PM_SUBTRACTION_KERNEL_SIMPLE:
+      psStringAppend(&kernels->description, "SIMPLE(%d,%s)",kernels->size,params);
+      break;
       default:
         psAbort("unknown kernel");
Index: trunk/psModules/src/imcombine/pmSubtractionKernels.h
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionKernels.h	(revision 35768)
+++ trunk/psModules/src/imcombine/pmSubtractionKernels.h	(revision 35771)
@@ -22,4 +22,9 @@
     } \
     if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_ISIS_RADIAL) { \
+        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \
+        PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \
+        PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \
+    } \
+    if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_SIMPLE) { \
         PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \
         PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \
@@ -71,4 +76,10 @@
 }
 
+// Generate 1D convolution kernel for SIMPLE
+psVector *pmSubtractionKernelSIMPLE(float sigma, // Gaussian width
+				    int order,   // Unused polynomial order
+				    int size     // Kernel half-size
+				    );
+
 // Generate 1D convolution kernel for ISIS
 psVector *pmSubtractionKernelISIS(float sigma, // Gaussian width
Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionMatch.c	(revision 35768)
+++ trunk/psModules/src/imcombine/pmSubtractionMatch.c	(revision 35771)
@@ -21,4 +21,5 @@
 #include "pmSubtractionThreads.h"
 #include "pmSubtractionVisual.h"
+#include "pmSubtractionSimple.h"
 #include "pmErrorCodes.h"
 
@@ -537,4 +538,12 @@
     }
 
+    // Bail here if we're doing the simple matching
+    if (type == PM_SUBTRACTION_KERNEL_SIMPLE) {
+      if (!pmSubtractionSimpleMatch(conv1,conv2,ro1,ro2,sources,size,maskVal,maskBad,maskPoor)) {
+	return false;
+      }
+      return(true);
+    }
+
     // Where does our variance map come from?
     // Getting the variance exactly right is not necessary --- it's just used for weighting.
@@ -593,4 +602,5 @@
 	goto MATCH_ERROR;
     }
+
     
     memCheck("start");
@@ -648,7 +658,7 @@
         psFree(bg);
     }
-
+    
     subtractionMatchAlloc(conv1, conv2, ro1, ro2, subMask, maskBad, subMode, numCols, numRows);
-
+    
     // Iterate over iso-kernel regions
     for (int j = 0; j < yRegions; j++) {
Index: trunk/psModules/src/imcombine/pmSubtractionSimple.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionSimple.c	(revision 35768)
+++ trunk/psModules/src/imcombine/pmSubtractionSimple.c	(revision 35771)
@@ -22,6 +22,44 @@
 #include "pmSubtractionVisual.h"
 #include "pmErrorCodes.h"
+#include "pmSpan.h"
+#include "pmFootprintSpans.h"
+#include "pmFootprint.h"
+#include "pmPeaks.h"
+#include "pmMoments.h"
+#include "pmModelFuncs.h"
+
+#include "pmSourceMasks.h"
+#include "pmSourceExtendedPars.h"
+#include "pmSourceDiffStats.h"
+#include "pmSourceSatstar.h"
+
+#include "pmSource.h"
 
 #include "pmSubtractionSimple.h"
+
+
+bool simple_do_boxphot(int *nPix,
+		       float *flux,
+		       pmSource *source,
+		       psImage *image,
+		       psImage *mask,
+		       psImageMaskType maskVal,
+		       int size) {
+  *nPix = 0;
+  *flux = 0.0;
+  for (int y = source->peak->yf - size;y <= source->peak->yf + size;y++) {
+    for (int x = source->peak->xf - size; x <= source->peak->xf + size; x++) {
+      if ((y > 0)&&(y < image->numRows)&&
+	  (x > 0)&&(x < image->numCols)) {
+	if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) {
+	  *nPix += 1;
+	  *flux += image->data.F32[y][x];
+	}
+      }
+    }
+  }
+  *flux = log10(*flux);
+  return(true);
+}  
 
 bool pmSubtractionSimpleMatch(pmReadout *conv1,
@@ -40,4 +78,5 @@
   float fwhm1,fwhm2;
   float sigma1,sigma2,sigmaKern;
+  float chisq = 1.0;
   int convolution_direction = 0;
   psImage *image1;
@@ -62,4 +101,5 @@
 
   if (conv1) {
+    conv1->covariance = psMemIncrRefCounter(ro1->covariance);
     if (!conv1->image) {
       conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
@@ -78,4 +118,5 @@
   }
   if (conv2) {
+    conv2->covariance = psMemIncrRefCounter(ro2->covariance);
     if (!conv2->image) {
       conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
@@ -129,92 +170,53 @@
   }
   if (!conv1) {
+    if (convolution_direction == 1) {
+      chisq = 100;
+    }
     convolution_direction = 2;
   }
   if (!conv2) {
+    if (convolution_direction == 2) {
+      chisq = 100;
+    }
     convolution_direction = 1;
   }
-
-  //
-  // Determine Normalization scaling
-  psVector *kernelVec = pmSubtractionKernelSIMPLE(sigmaKern,0,size); // This is normalized to unity.
-
-  
+  
+
+  int maskBox = (int) ceil(sigmaKern * 1.1774); // diameter is 1/2 FWHM
+  int maskBlank = 8;  // I should be able to get this from a reference, right?
+
+  //
+  // Construct required kernel.  No longer needed as we can direct convolve
+  //  psVector *kernelVec = pmSubtractionKernelSIMPLE(sigmaKern,0,size); // This is normalized to unity.
+  //  psFree(kernelVec);
+
   //
   // Do convolutions
-  psImage *temp = psImageAlloc(numCols,numRows,PS_TYPE_F32);
-  psImage *tempMask = psImageAlloc(numCols,numRows,PS_TYPE_IMAGE_MASK);
-  psImage *tempVar  = psImageAlloc(numCols,numRows,PS_TYPE_F32);
-
-  for (int y = 0; y < numRows; y++) {
-    for (int x = 0; x < numCols; x++) {
-      temp->data.F32[y][x] = 0.0;
-      tempVar->data.F32[y][x] = 0.0;
-      tempMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
-
-      // Copy the image we're not convolving into the convolved output
-      if (convolution_direction == 1) {
-	if (conv1) {
-	  imageC1->data.F32[y][x] = 0.0;
-	  maskC1->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
-	  varC1->data.F32[y][x] = 0.0;
-	}
-	if (conv2) {
-	  imageC2->data.F32[y][x] = image2->data.F32[y][x];
-	  maskC2->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = mask2->data.PS_TYPE_IMAGE_MASK_DATA[y][x];
-	  varC2->data.F32[y][x] = var2->data.F32[y][x];
-	}
-      }
-      else if (convolution_direction == 2) {
-	if (conv2) {
-	  imageC2->data.F32[y][x] = 0.0;
-	  maskC2->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
-	  varC2->data.F32[y][x] = 0.0;
-	}
-	if (conv1) {
-	  imageC1->data.F32[y][x] = image1->data.F32[y][x];
-	  maskC1->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = mask1->data.PS_TYPE_IMAGE_MASK_DATA[y][x];
-	  varC1->data.F32[y][x] = var1->data.F32[y][x];
-	}
-      }
-
-      // Do y-direction convolution
-      for (int v = -size; v <= size; v++) {
-	if ((y-v >= 0)&&(y-v < numRows)) {
-	  if (convolution_direction == 1) {
-	    // Insert mask/finite check here.
-	    temp->data.F32[y][x] += kernelVec->data.F32[v + size] * image1->data.F32[y - v][x];
-	    tempVar->data.F32[y][x] += kernelVec->data.F32[v + size] * var1->data.F32[y - v][x];
-	  }
-	  else if (convolution_direction == 2) {
-	    // And here.
-	    temp->data.F32[y][x] += kernelVec->data.F32[v + size] * image2->data.F32[y - v][x];
-	    tempVar->data.F32[y][x] += kernelVec->data.F32[v + size] * var2->data.F32[y - v][x];
-	  }
-	}
-      }
-    }
-  }
-
-  // Do x-direction convolution
-  for (int y = 0; y < numRows; y++) {
-    for (int x = 0; x < numCols; x++) {
-	for (int u = -size; u <= size; u++) {
-	  if ((x-u >= 0)&&(x-u < numCols)) {
-	    if (convolution_direction == 1) {
-	      // And here
-	      imageC1->data.F32[y][x] += kernelVec->data.F32[u + size] * temp->data.F32[y][x - u];
-	      varC1->data.F32[y][x] += kernelVec->data.F32[u + size] * temp->data.F32[y][x - u];
-	    }
-	    else if (convolution_direction == 2) {
-	      // And here
-	      imageC2->data.F32[y][x] += kernelVec->data.F32[u + size] * temp->data.F32[y][x - u];
-	      varC2->data.F32[y][x] += kernelVec->data.F32[u + size] * temp->data.F32[y][x - u];
-	    }
-	  }
-	}
-    }
-  }
-  psFree(temp);
-  psFree(kernelVec);
+  if (convolution_direction == 1) {
+    psImageSmoothMask_Threaded(imageC1,image1,mask1,maskVal,sigmaKern,6,1e-6);
+    psImageSmoothMask_Threaded(varC1,var1,mask1,maskVal,sigmaKern * M_SQRT1_2,6,1e-6);
+    maskC1 = psImageConvolveMask(maskC1,mask1,maskVal,maskBad,
+				 -maskBox,maskBox,-maskBox,maskBox);
+    if (conv2) {
+      imageC2 = psImageCopy(imageC2,image2,PS_TYPE_F32);
+      varC2   = psImageCopy(varC2,var2,PS_TYPE_F32);
+      maskC2  = psImageCopy(maskC2,mask2,PS_TYPE_IMAGE_MASK);
+    }
+    pmSubtractionBorder(imageC1,varC1,maskC1,maskBox,maskBlank);
+    pmSubtractionMaskApply(imageC1,varC1,maskC1,PM_SUBTRACTION_MODE_1);
+  }
+  else if (convolution_direction == 2) {
+    psImageSmoothMask_Threaded(imageC2,image2,mask2,maskVal,sigmaKern,6,1e-6);
+    psImageSmoothMask_Threaded(varC2,var2,mask2,maskVal,sigmaKern * M_SQRT1_2,6,1e-6);
+    maskC2 = psImageConvolveMask(maskC2,mask2,maskVal,maskBad,
+				 -maskBox,maskBox,-maskBox,maskBox);
+    if (conv1) {
+      imageC1 = psImageCopy(imageC1,image1,PS_TYPE_F32);
+      varC1   = psImageCopy(varC1,var1,PS_TYPE_F32);
+      maskC1  = psImageCopy(maskC1,mask1,PS_TYPE_IMAGE_MASK);
+    }
+    pmSubtractionBorder(imageC2,varC2,maskC2,maskBox,maskBlank);
+    pmSubtractionMaskApply(imageC2,varC2,maskC2,PM_SUBTRACTION_MODE_2);
+  }    
 
   //
@@ -222,17 +224,67 @@
   float normalization = 1.0;
 
-  // Something with the source list here
+  // Scan source list, do box photometry on peaks, and then solve the linear relation.
+  int photRadius = (int) floor(PS_MAX(sigma1,sigma2) * 2.0 * sqrt(2.0 * log(2.0))); // Go out a FWHM diameter from the center.
+  psVector *logFluxDifferences = psVectorAlloc(sources->n,PS_TYPE_F32);
+  psVector *fitMask = psVectorAlloc(sources->n,PS_TYPE_VECTOR_MASK);
+  for (int i = 0; i < sources->n; i++) {
+    pmSource *source = sources->data[i];
+    int nPix1,nPix2;
+    float flux1,flux2;
+
+    if (convolution_direction == 1) {
+      simple_do_boxphot(&nPix1,&flux1,source,imageC1,maskC1,maskBad,photRadius);
+      if (conv2) {
+	simple_do_boxphot(&nPix2,&flux2,source,imageC2,maskC2,maskBad,photRadius);
+      }
+      else {
+	simple_do_boxphot(&nPix2,&flux2,source,image2,mask2,maskBad,photRadius);
+      }
+    }
+    else if (convolution_direction == 2) {
+      simple_do_boxphot(&nPix2,&flux2,source,imageC2,maskC2,maskBad,photRadius);
+      if (conv1) {
+	simple_do_boxphot(&nPix1,&flux1,source,imageC1,maskC1,maskBad,photRadius);
+      }
+      else {
+	simple_do_boxphot(&nPix1,&flux1,source,image1,mask1,maskBad,photRadius);
+      }
+    }
+    logFluxDifferences->data.F32[i] = flux2 - flux1;
+    fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
+    if ((PS_MIN(nPix1,nPix2) <= 0.75 * PS_MAX(nPix1,nPix2))||
+	(!isfinite(flux1))||(!isfinite(flux2))) {
+      fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
+    }
+
+    //    fprintf(stderr,"SOURCES: %d %g %g %g -> %d %d %g %g %d %g\n",i,source->peak->xf,source->peak->yf,source->psfMag,
+    //	    nPix1,nPix2,flux1,flux2,fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i],logFluxDifferences->data.F32[i]);
+    
+  }
+
+  // Given the differences in log-flux space, the normalization factor is just the exponential of the median difference
+  psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+  if (!psVectorStats(stats,logFluxDifferences,NULL,fitMask,0xff)) {
+    // This should complain.
+    normalization = 1.0;
+  }
+
+  normalization = pow(10,stats->robustMedian);
+  // fprintf(stderr,"NORM: %g+/-%g\n",stats->robustMedian,stats->robustStdev);
+  
+  psFree(stats);
+  psFree(logFluxDifferences);
+  psFree(fitMask);
 
   // Apply normalization
   if (normalization != 1.0) {
-    for (int y = 0; y < numRows; y++) {
-      for (int x = 0; x < numCols; x++) {
-	if ((conv1)&&((convolution_direction == 1))) {
-	  imageC1->data.F32[y][x] *= normalization;
-	}
-	else if ((conv2)&&(convolution_direction == 2)) {
-	  imageC2->data.F32[y][x] *= normalization;
-	}
-      }
+    if ((conv1)&&((convolution_direction == 1))) {
+      psBinaryOp(imageC1,imageC1,"*",psScalarAlloc((float) normalization, PS_TYPE_F32));
+      psBinaryOp(varC1,varC1,"*",psScalarAlloc((float) PS_SQR(normalization), PS_TYPE_F32));
+    }
+    else if ((conv2)&&(convolution_direction == 2)) {
+      normalization = 1.0 / normalization; // Because we fit one way, but are using it in the other.
+      psBinaryOp(imageC2,imageC2,"*",psScalarAlloc((float) normalization, PS_TYPE_F32));
+      psBinaryOp(varC2,varC2,"*",psScalarAlloc((float) PS_SQR(normalization), PS_TYPE_F32));
     }
   }
@@ -263,5 +315,5 @@
     psFree(kernels->preCalc->data[0]);
   }
-  kernels->preCalc->data[0] = preCalc;
+  kernels->preCalc->data[0] = psMemIncrRefCounter(preCalc);
   kernels->solution1 = psVectorAlloc(3,PS_TYPE_F64);
   kernels->solution1->data.F32[0] = 1.0;
@@ -272,4 +324,7 @@
   kernels->solution1err->data.F32[1] = 0.0;
   kernels->solution1err->data.F32[2] = 0.0;
+  kernels->mean = 0.0;
+  kernels->rms = chisq; // This is the chi^2 value that's passed to ppStack
+  kernels->numStamps = sources->n;
   
   //
@@ -284,4 +339,10 @@
     psFree(kernels);
   }
+  // This is a hack.  Yes, I know.  pmSAnalysis doesn't get the normalization correct, so I just do it here.
+  psMetadataRemoveKey(outAnalysis,PM_SUBTRACTION_ANALYSIS_NORM);
+  psMetadataRemoveKey(outHeader,PM_SUBTRACTION_ANALYSIS_NORM);
+  psMetadataAddF32(outAnalysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM, PS_META_REPLACE, "Normalisation", normalization);
+  psMetadataAddF32(outHeader,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM, PS_META_REPLACE, "Normalisation", normalization);
+
   psFree(fwhms);
   psFree(orders);
@@ -291,7 +352,13 @@
   if (conv1) {
     conv1->analysis = psMetadataCopy(conv1->analysis, outAnalysis);
+    conv1->data_exists = true;
+    conv1->parent->data_exists = true;
+    conv1->parent->parent->data_exists = true;
   }
   if (conv2) {
     conv2->analysis = psMetadataCopy(conv2->analysis, outAnalysis);
+    conv2->data_exists = true;
+    conv2->parent->data_exists = true;
+    conv2->parent->parent->data_exists = true;
   }
 
Index: trunk/psModules/src/imcombine/pmSubtractionSimple.h
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionSimple.h	(revision 35768)
+++ trunk/psModules/src/imcombine/pmSubtractionSimple.h	(revision 35771)
@@ -13,5 +13,8 @@
 			      const pmReadout *ro2,
 			      const psArray *sources,
-			      int size
+			      int size,
+			      psImageMaskType maskVal,
+			      psImageMaskType maskBad,
+			      psImageMaskType maskPoor
 			      );
 
Index: trunk/psModules/src/imcombine/pmSubtractionTypes.h
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionTypes.h	(revision 35768)
+++ trunk/psModules/src/imcombine/pmSubtractionTypes.h	(revision 35771)
@@ -42,4 +42,5 @@
     PM_SUBTRACTION_KERNEL_GUNK,         ///< Grid United with Normal Kernel --- POIS and ISIS hybrid
     PM_SUBTRACTION_KERNEL_RINGS,        ///< Rings Instead of the Normal Gaussian Subtraction
+    PM_SUBTRACTION_KERNEL_SIMPLE,       ///< Simple Gaussian kernel to avoid complications
 } pmSubtractionKernelsType;