Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14739)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14752)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.56 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-05 00:15:28 $
+ *  @version $Revision: 1.57 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-05 21:17:59 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -75,12 +75,20 @@
 }
 
-// Apply the appropriate weighting to the kernel value
-static inline float kernelWeighting(float value, // Value to apply weighting
-                                    bool variance // Do variance weighting?
-                                    )
-{
-    return variance ? PS_SQR(value) : value;
-}
-
+// Generate the kernel to apply to the variance from the normal kernel
+static psKernel *varianceKernel(psKernel *normalKernel // Normal kernel
+                                )
+{
+    psKernel *kernel = psKernelAlloc(normalKernel->xMin, normalKernel->xMax,
+                                     normalKernel->yMin, normalKernel->yMax); // Kernel to return
+
+    // Take the square of the normal kernel
+    for (int v = -size; v <= size; v++) {
+        for (int u = -size; u <= size; u++) {
+            kernel->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]);
+        }
+    }
+
+    return kernel;
+}
 
 // Generate an image of the solved kernel
@@ -88,6 +96,5 @@
                                      const psVector *solution, // Solution to the least-squares problem
                                      const pmSubtractionKernels *kernels, // Kernel basis functions
-                                     const psImage *polyValues, // Spatial polynomial values
-                                     bool varianceWeighting // Do variance weighting?
+                                     const psImage *polyValues // Spatial polynomial values
                                      )
 {
@@ -113,12 +120,10 @@
 
     for (int i = 0; i < numKernels; i++) {
-        double sum = 0.0;               // Sum of multiple polynomial terms
+        double value = 0.0;             // Value of this component, from the sum of multiple polynomial terms
         for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
             for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
-                sum += polyValues->data.F64[yOrder][xOrder] * solution->data.F64[index];
-            }
-        }
-        float value = kernelWeighting(sum, varianceWeighting); // Value to be added to kernel
-        float subValue = kernelWeighting(- sum, varianceWeighting); // Value to be added for subtraction
+                value += polyValues->data.F64[yOrder][xOrder] * solution->data.F64[index];
+            }
+        }
 
         switch (kernels->type) {
@@ -129,5 +134,5 @@
               if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
                   // The (0,0) element is subtracted from most kernels to preserve photometric scaling
-                  kernel->kernel[0][0] += subValue;
+                  kernel->kernel[0][0] -= value;
               }
               break;
@@ -142,6 +147,5 @@
 
               // Normalising sum of kernel component to unity
-              float norm = kernelWeighting(1.0 / (float)((uStop - uStart + 1) * (vStop - vStart + 1)),
-                                           varianceWeighting);
+              float norm = 1.0 / (float)((uStop - uStart + 1) * (vStop - vStart + 1));
 
               for (int v = vStart; v <= vStop; v++) {
@@ -150,5 +154,5 @@
                       if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
                           // The (0,0) element is subtracted from most kernels to preserve photometric scaling
-                          kernel->kernel[0][0] += subValue;
+                          kernel->kernel[0][0] -= value;
                       }
                   }
@@ -163,6 +167,5 @@
                   for (int v = -size; v <= size; v++) {
                       for (int u = -size; u <= size; u++) {
-                          kernel->kernel[v][u] += kernelWeighting(preCalc->kernel[v][u],
-                                                                  varianceWeighting) * value;
+                          kernel->kernel[v][u] += preCalc->kernel[v][u] * value;
                           // Photometric scaling is built into the preCalc kernel --- no subtraction!
                       }
@@ -176,5 +179,5 @@
                   if (subtract) {
                       // The (0,0) element is subtracted from most kernels to preserve photometric scaling
-                      kernel->kernel[0][0] += subValue;
+                      kernel->kernel[0][0] -= value;
                   }
               }
@@ -186,6 +189,5 @@
               for (int v = -size; v <= size; v++) {
                   for (int u = -size; u <= size; u++) {
-                      kernel->kernel[v][u] += kernelWeighting(preCalc->kernel[v][u],
-                                                              varianceWeighting) * value;
+                      kernel->kernel[v][u] += preCalc->kernel[v][u] * value;
                       // Photometric scaling is built into the preCalc kernel --- no subtraction!
                   }
@@ -206,5 +208,5 @@
               for (int j = 0; j < num; j++) {
                   int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
-                  kernel->kernel[v][u] += kernelWeighting(poly->data.F32[j], varianceWeighting) * value;
+                  kernel->kernel[v][u] += poly->data.F32[j] * value;
               }
               // Photometric scaling is built into the kernel --- no subtraction!
@@ -939,5 +941,5 @@
 
     // The appropriate kernel
-    psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues, false);
+    psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues);
 
     psFree(polyValues);
@@ -1091,7 +1093,7 @@
                                            2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows);
 
-            kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, false);
+            kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues);
             if (inWeight) {
-                kernelWeight = solvedKernel(kernelWeight, solution, kernels, polyValues, true);
+                kernelWeight = varianceKernel(kernelImage);
             }