Index: trunk/psModules/src/imcombine/pmSubtractionSimple.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionSimple.c	(revision 35817)
+++ trunk/psModules/src/imcombine/pmSubtractionSimple.c	(revision 35821)
@@ -63,4 +63,23 @@
 }  
 
+// Copied from pmSubtraction
+static void solvedKernelPreCalc(psKernel *kernel, // Kernel, updated
+				const pmSubtractionKernels *kernels, // Kernel basis functions
+				float value,                         // Normalisation value for basis function
+				int index                  // Index of basis function of interest
+				)
+{
+  int size = kernels->size;           // Kernel half-size
+  pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated values
+  for (int v = -size; v <= size; v++) {
+    for (int u = -size; u <= size; u++) {
+      kernel->kernel[v][u] +=  value * preCalc->kernel->kernel[v][u];
+    }
+  }
+
+  return;
+}
+//End copy
+
 bool pmSubtractionSimpleMatch(pmReadout *conv1,
 			      pmReadout *conv2,
@@ -101,5 +120,5 @@
 
   if (conv1) {
-    conv1->covariance = psMemIncrRefCounter(ro1->covariance);
+    //    conv1->covariance = psMemIncrRefCounter(ro1->covariance);
     if (!conv1->image) {
       conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
@@ -118,5 +137,5 @@
   }
   if (conv2) {
-    conv2->covariance = psMemIncrRefCounter(ro2->covariance);
+    //    
     if (!conv2->image) {
       conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
@@ -188,112 +207,7 @@
 
   //
-  // 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
-  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);
-  }    
-
-  //
-  // Do normalization
-  float normalization = 1.0;
-
-  // 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) {
-    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));
-    }
-  }
-  
-
-  //
   // Make a fake pmSubtractionKernels element so we can add it appropriately.
-  // I call it fake because we've successfully done everything at this point
-  // without having to define these things.
+  // Defining everything here is a bit clunky, but it's necessary to get the covariance
+  // correct.
   psVector *fwhms = psVectorAlloc(1,PS_TYPE_F32);
   fwhms->data.F32[0] = sigmaKern * sig2fwhm;
@@ -329,4 +243,115 @@
   kernels->numStamps = sources->n;
   
+  psKernel *kernel = psKernelAlloc(-size,size,-size,size);
+  solvedKernelPreCalc(kernel,kernels,1.0,0);
+  
+  //
+  // Do convolutions
+  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);
+    conv1->covariance = psImageCovarianceCalculate(kernel,ro1->covariance);
+    if (conv2) {
+      imageC2 = psImageCopy(imageC2,image2,PS_TYPE_F32);
+      varC2   = psImageCopy(varC2,var2,PS_TYPE_F32);
+      maskC2  = psImageCopy(maskC2,mask2,PS_TYPE_IMAGE_MASK);
+      conv2->covariance = psMemIncrRefCounter(ro2->covariance);
+    }
+    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);
+    conv2->covariance = psImageCovarianceCalculate(kernel,ro2->covariance);
+    if (conv1) {
+      imageC1 = psImageCopy(imageC1,image1,PS_TYPE_F32);
+      varC1   = psImageCopy(varC1,var1,PS_TYPE_F32);
+      maskC1  = psImageCopy(maskC1,mask1,PS_TYPE_IMAGE_MASK);
+      conv1->covariance = psMemIncrRefCounter(ro1->covariance);
+    }
+    pmSubtractionBorder(imageC2,varC2,maskC2,maskBox,maskBlank);
+    pmSubtractionMaskApply(imageC2,varC2,maskC2,PM_SUBTRACTION_MODE_2);
+  }    
+
+  psFree(kernel); // No longer needed after doing covariance calculation
+
+  //
+  // Do normalization
+  float normalization = 1.0;
+
+  // 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) {
+    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));
+    }
+  }
+  
+
+  //
+  
   //
   // Actually add it to the headers