Index: /trunk/pois/src/poisCalculateDeviations.c
===================================================================
--- /trunk/pois/src/poisCalculateDeviations.c	(revision 6890)
+++ /trunk/pois/src/poisCalculateDeviations.c	(revision 6891)
@@ -33,4 +33,5 @@
     if (!deviations) {
         deviations = psVectorAlloc(stamps->n, PS_TYPE_F32);
+        deviations->n = stamps->n;
     }
 
@@ -64,9 +65,12 @@
             psRegion stampTrim = { config->xKernel, config->xKernel + 2 * footprint, config->yKernel,
                                    config->yKernel + 2 * footprint };
+#if 0
             psRegion maskTrim = { x - xSize + config->xKernel, x + xSize - config->xKernel,
                                   y - ySize + config->yKernel, y + ySize - config->yKernel };
-
+            psImage *maskStampTrim = psImageSubset(maskStamp, maskTrim);
+#else
+            psImage *maskStampTrim = psImageSubset(maskStamp, stampTrim);
+#endif
             psImage *subStampTrim = psImageSubset(subStamp, stampTrim);
-            psImage *maskStampTrim = psImageSubset(maskStamp, maskTrim);
             (void)psImageStats(stats, subStampTrim, maskStampTrim, POIS_MASK_BAD | POIS_MASK_NEAR_BAD);
 
Index: /trunk/pois/src/poisCalculateEquation.c
===================================================================
--- /trunk/pois/src/poisCalculateEquation.c	(revision 6890)
+++ /trunk/pois/src/poisCalculateEquation.c	(revision 6891)
@@ -60,4 +60,5 @@
                 psTrace("pois.calculateEquation", 3, "Allocating vector: %d\n", numSolveParams);
                 stamp->vector = psVectorAlloc(numSolveParams, PS_TYPE_F64);
+                stamp->vector->n = numSolveParams;
             }
 
Index: /trunk/pois/src/poisCheckKernel.c
===================================================================
--- /trunk/pois/src/poisCheckKernel.c	(revision 6890)
+++ /trunk/pois/src/poisCheckKernel.c	(revision 6891)
@@ -4,6 +4,6 @@
 
 bool poisCheckKernel(const psVector *solution,// Kernel solution
-		     const psArray *kernels, // Kernel basis functions
-		     const poisConfig *config // Configuration
+                     const psArray *kernels, // Kernel basis functions
+                     const poisConfig *config // Configuration
     )
 {
@@ -13,38 +13,39 @@
     psVector *radKernel = psVectorAlloc((2*config->xKernel + 1) * (2*config->yKernel + 1), PS_TYPE_F32);
     psVector *radii = psVectorAlloc((2*config->xKernel + 1) * (2*config->yKernel + 1), PS_TYPE_F32); // Radius
+    radKernel->n = radii->n = (2*config->xKernel + 1) * (2*config->yKernel + 1);
     for (int i = 0; i < solution->n - 1; i++) {
-	poisKernelBasis *kernel = kernels->data[i]; // The kernel basis function
-	if (kernel->xOrder == 0 && kernel->yOrder == 0) {
-	    int u = kernel->u;		// x offset
-	    int v = kernel->v;		// y offset
-	    float distance = sqrtf((float)(u*u) + (float)(v*v)); // Distance from the centre
-	    radKernel->data.F32[i] = solution->data.F64[i];
-	    radii->data.F32[i] = distance;
-	}
+        poisKernelBasis *kernel = kernels->data[i]; // The kernel basis function
+        if (kernel->xOrder == 0 && kernel->yOrder == 0) {
+            int u = kernel->u;          // x offset
+            int v = kernel->v;          // y offset
+            float distance = sqrtf((float)(u*u) + (float)(v*v)); // Distance from the centre
+            radKernel->data.F32[i] = solution->data.F64[i];
+            radii->data.F32[i] = distance;
+        }
     }
 
     psVector *sortIndex = psVectorSortIndex(NULL, radii); // Indices from sort
-    float distance = 0.0;		// Distance from the centre that's being examined
-    float sum = 0.0;			// Sum of kernel for that distance
-    int num = 0;			// Number of pixels at that distance
+    float distance = 0.0;               // Distance from the centre that's being examined
+    float sum = 0.0;                    // Sum of kernel for that distance
+    int num = 0;                        // Number of pixels at that distance
     int centreIndex = sortIndex->data.U32[0]; // Index of the centre pixel
-    int numSuspect = 0;			// Number of radii considered suspect
+    int numSuspect = 0;                 // Number of radii considered suspect
     for (int i = 0; i < (2*config->xKernel + 1) * (2*config->yKernel + 1); i++) {
-	unsigned int index = sortIndex->data.U32[i];
-	if (radii->data.F32[index] != distance) {
-	    psTrace("pois.checkKernel", 7, "Radius: %f\tMean: %f\n", distance, sum/(float)num);
-	    if (fabs(distance*sum/(float)num) > radKernel->data.F32[centreIndex]) {
-		numSuspect++;
-	    }
-	    sum = 0.0;
-	    num = 0;
-	    distance = radii->data.F32[index];
-	}
-	sum += radKernel->data.F32[index];
-	num++;
+        unsigned int index = sortIndex->data.U32[i];
+        if (radii->data.F32[index] != distance) {
+            psTrace("pois.checkKernel", 7, "Radius: %f\tMean: %f\n", distance, sum/(float)num);
+            if (fabs(distance*sum/(float)num) > radKernel->data.F32[centreIndex]) {
+                numSuspect++;
+            }
+            sum = 0.0;
+            num = 0;
+            distance = radii->data.F32[index];
+        }
+        sum += radKernel->data.F32[index];
+        num++;
     }
     psTrace("pois.checkKernel", 7, "Radius: %f\tMean: %f\n", distance, sum/(float)num);
     if (fabs(distance*sum/(float)num) > radKernel->data.F32[centreIndex]) {
-	numSuspect++;
+        numSuspect++;
     }
 
@@ -56,6 +57,6 @@
 
     if (numSuspect > 0) {
-	psWarning("There is significant power at large radii.  The kernel may be bad.\n");
-	return false;
+        psWarning("There is significant power at large radii.  The kernel may be bad.\n");
+        return false;
     }
 
Index: /trunk/pois/src/poisFindStamps.c
===================================================================
--- /trunk/pois/src/poisFindStamps.c	(revision 6890)
+++ /trunk/pois/src/poisFindStamps.c	(revision 6891)
@@ -31,4 +31,5 @@
     if (stamps == NULL) {
         stamps = psArrayAlloc(nsx * nsy);
+        stamps->n = nsx * nsy;
         // Initialise
         for (int i = 0; i < nsx * nsy; i++) {
Index: /trunk/pois/src/poisKernelBasisFunctions.c
===================================================================
--- /trunk/pois/src/poisKernelBasisFunctions.c	(revision 6890)
+++ /trunk/pois/src/poisKernelBasisFunctions.c	(revision 6891)
@@ -14,61 +14,62 @@
     int yKernelHalfSize = config->yKernel; // Half size of kernel in y
     int spatialOrder = config->spatialOrder; // Order of spatial variations
-	
+
     // Number of basis functions
     int nBF = (2 * xKernelHalfSize + 1) * (2 * yKernelHalfSize + 1) * (spatialOrder + 1) *
-	(spatialOrder + 2) / 2;
+        (spatialOrder + 2) / 2;
 
-    psArray *array = psArrayAlloc(nBF);	// Array to hold the basis functions
+    psArray *array = psArrayAlloc(nBF); // Array to hold the basis functions
+    array->n = nBF;
 
-    int num = 0;			// Kernel parameter number
+    int num = 0;                        // Kernel parameter number
 
     // Put the (u,v) = (0,0) component right at the start of the list for convenience
     for (int order = 0; order <= spatialOrder; order++) {
-	for (int i = 0; i <= order; i++) {
-	    int j = order - i;
-	    poisKernelBasis *kernel = psAlloc(sizeof(poisKernelBasis));
-	    kernel->u = 0;
-	    kernel->v = 0;
-	    kernel->xOrder = i;
-	    kernel->yOrder = j;
+        for (int i = 0; i <= order; i++) {
+            int j = order - i;
+            poisKernelBasis *kernel = psAlloc(sizeof(poisKernelBasis));
+            kernel->u = 0;
+            kernel->v = 0;
+            kernel->xOrder = i;
+            kernel->yOrder = j;
 
-	    // Stuff into the array
-	    array->data[num] = kernel;
+            // Stuff into the array
+            array->data[num] = kernel;
 
-	    psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", i, kernel->u,
-		    kernel->v, kernel->xOrder, kernel->yOrder);
-	    num++;
-	}
+            psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", i, kernel->u,
+                    kernel->v, kernel->xOrder, kernel->yOrder);
+            num++;
+        }
     }
 
     // Iterate over (u,v)
     for (int u = -xKernelHalfSize; u <= xKernelHalfSize; u++) {
-	for (int v = -yKernelHalfSize; v <= yKernelHalfSize; v++) {
+        for (int v = -yKernelHalfSize; v <= yKernelHalfSize; v++) {
 
-	    // Already did (u,v) = (0,0): it's at the start, so skip it now.
-	    if ((u != 0) || (v != 0)) {
+            // Already did (u,v) = (0,0): it's at the start, so skip it now.
+            if ((u != 0) || (v != 0)) {
 
-		// Iterate over spatial order
-		for (int order = 0; order <= spatialOrder; order++) {
+                // Iterate over spatial order
+                for (int order = 0; order <= spatialOrder; order++) {
 
-		    for (int i = 0; i <= order; i++) {
-			int j = order - i;
-			poisKernelBasis *kernel = psAlloc(sizeof(poisKernelBasis));
+                    for (int i = 0; i <= order; i++) {
+                        int j = order - i;
+                        poisKernelBasis *kernel = psAlloc(sizeof(poisKernelBasis));
 
-			kernel->u = u;
-			kernel->v = v;
-			kernel->xOrder = i;
-			kernel->yOrder = j;
+                        kernel->u = u;
+                        kernel->v = v;
+                        kernel->xOrder = i;
+                        kernel->yOrder = j;
 
-			// Stuff into the array
-			array->data[num] = kernel;
+                        // Stuff into the array
+                        array->data[num] = kernel;
 
-			psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", num,
-				kernel->u, kernel->v, kernel->xOrder, kernel->yOrder);
-			num++;
-		    }
-		} // Iterating over spatial order
-	    } // Except (u,v) = (0,0)
-	}
+                        psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", num,
+                                kernel->u, kernel->v, kernel->xOrder, kernel->yOrder);
+                        num++;
+                    }
+                } // Iterating over spatial order
+            } // Except (u,v) = (0,0)
+        }
     } // Iterating over (u,v)
 
Index: /trunk/pois/src/poisReadStamps.c
===================================================================
--- /trunk/pois/src/poisReadStamps.c	(revision 6890)
+++ /trunk/pois/src/poisReadStamps.c	(revision 6891)
@@ -17,4 +17,5 @@
     if (stamps == NULL) {
         stamps = psArrayAlloc(config->nsx);
+        stamps->n = config->nsx;
         // Initialise
         for (int i = 0; i < config->nsx; i++) {
Index: /trunk/pois/src/poisSolveEquation.c
===================================================================
--- /trunk/pois/src/poisSolveEquation.c	(revision 6890)
+++ /trunk/pois/src/poisSolveEquation.c	(revision 6891)
@@ -4,7 +4,7 @@
 #include "pois.h"
 
-psVector *poisSolveEquation(psVector *solution,	// Solution vector, or NULL
-			    const psArray *stamps, // Array of stamps
-			    const poisConfig *config
+psVector *poisSolveEquation(psVector *solution, // Solution vector, or NULL
+                            const psArray *stamps, // Array of stamps
+                            const poisConfig *config
     )
 {
@@ -14,5 +14,5 @@
 
     int size = (2 * config->xKernel + 1) * (2 * config->yKernel + 1) * (config->spatialOrder + 1) *
-	(config->spatialOrder + 2) / 2 + 1; // Size of matrix and vector
+        (config->spatialOrder + 2) / 2 + 1; // Size of matrix and vector
 
     assert(!solution || solution->n == size);
@@ -20,26 +20,27 @@
     psImage *matrix = psImageAlloc(size, size, PS_TYPE_F64);
     psVector *vector = psVectorAlloc(size, PS_TYPE_F64);
+    vector->n = size;
     for (int j = 0; j < size; j++) {
-	for (int i = 0; i < size; i++) {
-	    matrix->data.F64[j][i] = 0.0;
-	}
-	vector->data.F64[j] = 0.0;
+        for (int i = 0; i < size; i++) {
+            matrix->data.F64[j][i] = 0.0;
+        }
+        vector->data.F64[j] = 0.0;
     }
 
     for (int s = 0; s < stamps->n; s++) {
-	poisStamp *stamp = stamps->data[s]; // Stamp of interest
-	psImage *stampMatrix = stamp->matrix; // Corresponding matrix
-	psVector *stampVector = stamp->vector; // Corresponding vector
+        poisStamp *stamp = stamps->data[s]; // Stamp of interest
+        psImage *stampMatrix = stamp->matrix; // Corresponding matrix
+        psVector *stampVector = stamp->vector; // Corresponding vector
 
-	// Check inputs
-	assert(matrix->type.type == PS_TYPE_F64);
-	assert(vector->type.type == PS_TYPE_F64);
-	assert(matrix->numCols == size && matrix->numRows == size);
-	assert(vector->n == size);
+        // Check inputs
+        assert(matrix->type.type == PS_TYPE_F64);
+        assert(vector->type.type == PS_TYPE_F64);
+        assert(matrix->numCols == size && matrix->numRows == size);
+        assert(vector->n == size);
 
-	if (stamp->status == POIS_STAMP_USED) {
-	    (void)psBinaryOp(matrix, matrix, "+", stampMatrix);
-	    (void)psBinaryOp(vector, vector, "+", stampVector);
-	}
+        if (stamp->status == POIS_STAMP_USED) {
+            (void)psBinaryOp(matrix, matrix, "+", stampMatrix);
+            (void)psBinaryOp(vector, vector, "+", stampVector);
+        }
     }
 
@@ -53,8 +54,8 @@
     printf("Matrix:\n");
     for (int i = 0; i < matrix->numCols; i++) {
-	for (int j = 0; j < matrix->numRows; j++) {
-	    printf("%f ", matrix->data.F64[i][j]);
-	}
-	printf("\n");
+        for (int j = 0; j < matrix->numRows; j++) {
+            printf("%f ", matrix->data.F64[i][j]);
+        }
+        printf("\n");
     }
 #endif