Index: trunk/psModules/src/imcombine/pmSubtractionSimple.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionSimple.c	(revision 35717)
+++ trunk/psModules/src/imcombine/pmSubtractionSimple.c	(revision 35726)
@@ -30,5 +30,8 @@
 			      const pmReadout *ro2,
 			      const psArray *sources,
-			      int size
+			      int size,
+			      psImageMaskType maskVal,
+			      psImageMaskType maskBad,
+			      psImageMaskType maskPoor
 			      ) {
   //
@@ -38,58 +41,78 @@
   float sigma1,sigma2,sigmaKern;
   int convolution_direction = 0;
+  psImage *image1;
+  psImage *mask1;
+  psImage *var1;
+
+  psImage *image2;
+  psImage *mask2;
+  psImage *var2;
+
+  psImage *imageC1;
+  psImage *maskC1;
+  psImage *varC1;
+
+  psImage *imageC2;
+  psImage *maskC2;
+  psImage *varC2;
 
   // Allocate images, as this is usually done by subtractionMatchAlloc after this function is called.  
-  psImageMaskType maskBad = 1;
   int numCols = ro1->image->numCols;
   int numRows = ro1->image->numRows;
-  
-  if (!conv1->image) {
-    conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-  }
-  psImageInit(conv1->image, NAN);
-  if (ro1->variance) {
-    if (!conv1->variance) {
-      conv1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-    }
-    psImageInit(conv1->variance, NAN);
-  }
-  if (!conv1->mask) {
-    conv1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
-  }
-  psImageInit(conv1->mask, maskBad);
-  if (!conv2->image) {
-    conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-  }
-  psImageInit(conv2->image, NAN);
-  if (ro2->variance) {
-    if (!conv2->variance) {
-      conv2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-    }
-    psImageInit(conv2->variance, NAN);
-  }
-  if (!conv2->mask) {
-    conv2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
-  }
-  psImageInit(conv2->mask, maskBad);
-
-  
-  
-  
-  psImage *image1 = ro1->image;
-  psImage *mask1  = ro1->mask;
-  psImage *var1   = ro1->variance;
-
-  psImage *image2 = ro2->image;
-  psImage *mask2  = ro2->mask;
-  psImage *var2   = ro2->variance;
-
-  psImage *imageC1 = conv1->image;
-  psImage *maskC1  = conv1->mask;
-  psImage *varC1   = conv1->variance;
-
-  psImage *imageC2 = conv2->image;
-  psImage *maskC2  = conv2->mask;
-  psImage *varC2   = conv2->variance;
-  
+
+  if (conv1) {
+    if (!conv1->image) {
+      conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+    }
+    psImageInit(conv1->image, NAN);
+    if (ro1->variance) {
+      if (!conv1->variance) {
+	conv1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+      }
+      psImageInit(conv1->variance, NAN);
+    }
+    if (!conv1->mask) {
+      conv1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+    }
+    psImageInit(conv1->mask, maskBad);
+  }
+  if (conv2) {
+    if (!conv2->image) {
+      conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+    }
+    psImageInit(conv2->image, NAN);
+    if (ro2->variance) {
+      if (!conv2->variance) {
+	conv2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+      }
+      psImageInit(conv2->variance, NAN);
+    }
+    if (!conv2->mask) {
+      conv2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+    }
+    psImageInit(conv2->mask, maskBad);
+  }
+  
+  // Assign local aliases to images
+  image1 = ro1->image;
+  mask1  = ro1->mask;
+  var1   = ro1->variance;
+
+  image2 = ro2->image;
+  mask2  = ro2->mask;
+  var2   = ro2->variance;
+
+  if (conv1) {
+    imageC1 = conv1->image;
+    maskC1  = conv1->mask;
+    varC1   = conv1->variance;
+  }
+  
+  if (conv2) {
+    imageC2 = conv2->image;
+    maskC2  = conv2->mask;
+    varC2   = conv2->variance;
+  }
+ 
   //
   // Determine Gaussian widths
@@ -105,4 +128,10 @@
     sigmaKern = sqrt(PS_SQR(sigma2) - PS_SQR(sigma1));
   }
+  if (!conv1) {
+    convolution_direction = 2;
+  }
+  if (!conv2) {
+    convolution_direction = 1;
+  }
 
   //
@@ -110,43 +139,55 @@
   psVector *kernelVec = pmSubtractionKernelSIMPLE(sigmaKern,0,size); // This is normalized to unity.
 
-  float normalization = 1.0;
-
-  for (int i = 0; i < 2 * size + 1; i++) {
-    kernelVec->data.F32[i] *= sqrt(normalization); // Needs to be this because we apply this vector twice.
-  }
   
   //
   // Do convolutions
   psImage *temp = psImageAlloc(numCols,numRows,PS_TYPE_F32);
-  for (int y = 0; y < image1->numRows; y++) {
-    for (int x = 0; x < image1->numCols; x++) {
-
+  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) {
-	imageC1->data.F32[y][x] = 0.0;
-	maskC1->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
-	varC1->data.F32[y][x] = 0.0;
-
-	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];
+	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) {
-	imageC2->data.F32[y][x] = 0.0;
-	maskC2->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
-	varC2->data.F32[y][x] = 0.0;
-
-	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];
-      }
+	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];
 	  }
 	}
@@ -154,13 +195,19 @@
     }
   }
-  for (int y = 0; y < image1->numRows; y++) {
-    for (int x = 0; x < image1->numCols; 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];
 	    }
 	  }
@@ -171,4 +218,26 @@
   psFree(kernelVec);
 
+  //
+  // Do normalization
+  float normalization = 1.0;
+
+  // Something with the source list here
+
+  // 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;
+	}
+      }
+    }
+  }
+  
+
+  //
   // Make a fake pmSubtractionKernels element so we can add it appropriately.
   // I call it fake because we've successfully done everything at this point
@@ -199,18 +268,11 @@
   kernels->solution1->data.F32[1] = 0.0;
   kernels->solution1->data.F32[2] = 0.0;
-/*   kernels->solution2 = psVectorAlloc(3,PS_TYPE_F64); */
-/*   kernels->solution2->data.F32[0] = 1.0; */
   kernels->solution1err = psVectorAlloc(3,PS_TYPE_F32);
   kernels->solution1err->data.F32[0] = 0.0;
   kernels->solution1err->data.F32[1] = 0.0;
   kernels->solution1err->data.F32[2] = 0.0;
-/*   kernels->solution2err = psVectorAlloc(3,PS_TYPE_F32); */
-/*   kernels->solution2err->data.F32[0] = 0.0; */
-/*   kernels->solution2err->data.F32[1] = 0.0; */
-/*   kernels->solution2err->data.F32[2] = 0.0; */
-  
-
+  
+  //
   // Actually add it to the headers
-  
   psMetadata *outAnalysis = psMetadataAlloc();
   psMetadata *outHeader   = psMetadataAlloc();
@@ -248,6 +310,7 @@
   psFree(outAnalysis);
   psFree(outHeader);
+
+  //
   // Return
-
   return(true);
 }