Index: trunk/pois/src/poisConvolveImage.c
===================================================================
--- trunk/pois/src/poisConvolveImage.c	(revision 3801)
+++ trunk/pois/src/poisConvolveImage.c	(revision 5717)
@@ -7,91 +7,90 @@
 /* Convolve an image by the kernel */
 psImage *poisConvolveImage(const psImage *input, // Input image, to be convolved
-			   const psImage *mask,	// Mask image
-			   const psVector *solution, // The solution vector, containing the coefficients
-			   const psArray *kernelParams, // The kernel parameters
-			   const poisConfig *config // The configuration
+                           const psImage *mask, // Mask image
+                           const psVector *solution, // The solution vector, containing the coefficients
+                           const psArray *kernelParams, // The kernel parameters
+                           const poisConfig *config // The configuration
     )
 {
-    assert(solution->n == kernelParams->n + 1);	// Extra parameter for background
+    assert(solution->n == kernelParams->n + 1); // Extra parameter for background
     assert(solution->type.type == PS_TYPE_F64);
     assert(input->numCols == mask->numCols && input->numRows == mask->numRows);
     assert(mask->type.type == PS_TYPE_U8);
     assert(input->type.type == PS_TYPE_F32);
-    
-    int nx = input->numCols;		// The size of the image in x
-    int ny = input->numRows;		// The size of the image in y
-    int xKernel = config->xKernel;	// The half-size of the kernel in x
-    int yKernel = config->yKernel;	// The half-size of the kernel in y
-    float background = solution->data.F64[solution->n - 1]; // The difference in background
-    int numParams = kernelParams->n;	// Number of kernel parameters
+
+    int nx = input->numCols;            // The size of the image in x
+    int ny = input->numRows;            // The size of the image in y
+    float nxHalf = 0.5 * (float)nx;     // Half the size of the image
+    float nyHalf = 0.5 * (float)ny;     // Half the size of the image
+    int xKernel = config->xKernel;      // The half-size of the kernel in x
+    int yKernel = config->yKernel;      // The half-size of the kernel in y
+    float background = solution->data.F64[solution->n-1]; // The difference in background
+    int numParams = kernelParams->n;    // Number of kernel parameters
 
     psImage *convolved = psImageAlloc(nx, ny, PS_TYPE_F32); // The convolved image, to be returned
 
-    psDPolynomial2D *poly = psDPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
-						 PS_POLYNOMIAL_ORD); // Polynomial
+    psPolynomial2D *poly = psPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
+                                               PS_POLYNOMIAL_ORD ); // Polynomial
     for (int i = 0; i < config->spatialOrder + 1; i++) {
-	for (int j = 0; j < config->spatialOrder + 1; j++) {
-	    poly->coeff[i][j] = 1.0;
-	    poly->mask[i][j] = 1;	// Mask all coefficients; unmask to evaluate
-	}
+        for (int j = 0; j < config->spatialOrder + 1; j++) {
+            poly->coeff[i][j] = 1.0;
+            poly->mask[i][j] = 1;       // Mask all coefficients; unmask to evaluate
+        }
     }
 
     // Convolve only the middle part
     for (int y = yKernel; y < ny - yKernel; y++) {
-	for (int x = xKernel; x < nx - xKernel; x++) {
-	    if (! mask->data.U8[y][x] & (POIS_MASK_BAD | POIS_MASK_NEAR_BAD)) {
-		double conv = background; // Convolved value
+        for (int x = xKernel; x < nx - xKernel; x++) {
+            if (! mask->data.U8[y][x] & (POIS_MASK_BAD | POIS_MASK_NEAR_BAD)) {
+                double conv = background; // Convolved value
 
-#if 0
-		float imageX = (x - nxHalf) / nxHalf; // Normalised position in x
-		float imageY = (y - nyHalf) / nyHalf; // Normalised position in y
-#endif
-		
-		// Iterate over the kernel basis functions
-		for (int k = 0; k < numParams; k++) {
-		    poisKernelBasis *kernel = kernelParams->data[k]; // The kernel basis function
-		    int u = kernel->u;
-		    int v = kernel->v;
-		    
-#if 0
-		    int xOrder = kernel->xOrder;
-		    int yOrder = kernel->yOrder;
-		    // Evaluate the polynomial
-		    poly->mask[xOrder][yOrder] = 0;
-		    double polyVal = psDPolynomial2DEval(poly, imageX, imageY);
-		    poly->mask[xOrder][yOrder] = 1;
-#else
-		    double polyVal = 1.0;
-#endif
+                float imageX = (x - nxHalf) / nxHalf; // Normalised position in x
+                float imageY = (y - nyHalf) / nyHalf; // Normalised position in y
 
-		    if (k == 0) {
-			conv += solution->data.F64[k] * input->data.F32[y - v][x - u] * polyVal;
-		    } else {
-			conv += solution->data.F64[k] *
-			    (input->data.F32[y - v][x - u] * polyVal - input->data.F32[y][x]);
-		    }
-		}
-		convolved->data.F32[y][x] = (float)conv;
-	    }
-	}
-	    
-	/* Pad the rest with zeros */
-	for (int x = 0; x < xKernel; x++) {
-	    convolved->data.F32[y][x] = 0.0;
-	}
-	for (int x = nx - xKernel; x < nx; x++) {
-	    convolved->data.F32[y][x] = 0.0;
-	}
+                // Iterate over the kernel basis functions
+                for (int k = 0; k < numParams; k++) {
+                    poisKernelBasis *kernel = kernelParams->data[k]; // The kernel basis function
+                    int u = kernel->u;
+                    int v = kernel->v;
+                    int xOrder = kernel->xOrder;
+                    int yOrder = kernel->yOrder;
+
+                    double polyVal = 1.0;
+                    if (xOrder != 0 || yOrder != 0) {
+                        // Evaluate the polynomial
+                        poly->mask[xOrder][yOrder] = 0;
+                        polyVal = psPolynomial2DEval(poly, imageX, imageY);
+                        poly->mask[xOrder][yOrder] = 1;
+                    }
+
+                    if (k == 0) {
+                        conv += solution->data.F64[k] * input->data.F32[y - v][x - u] * polyVal;
+                    } else {
+                        conv += solution->data.F64[k] *
+                            (input->data.F32[y - v][x - u] * polyVal - input->data.F32[y][x]);
+                    }
+                }
+                convolved->data.F32[y][x] = (float)conv;
+            }
+        }
+
+        /* Pad the rest with zeros */
+        for (int x = 0; x < xKernel; x++) {
+            convolved->data.F32[y][x] = 0.0;
+        }
+        for (int x = nx - xKernel; x < nx; x++) {
+            convolved->data.F32[y][x] = 0.0;
+        }
     }
 
     for (int y = 0; y < yKernel; y++) {
-	for (int x = 0; x < nx; x++) {
-	    convolved->data.F32[y][x] = 0.0;
-	}
+        for (int x = 0; x < nx; x++) {
+            convolved->data.F32[y][x] = 0.0;
+        }
     }
     for (int y = ny - yKernel; y < ny; y++) {
-	for (int x = 0; x < nx; x++) {
-	    convolved->data.F32[y][x] = 0.0;
-	}
+        for (int x = 0; x < nx; x++) {
+            convolved->data.F32[y][x] = 0.0;
+        }
     }