Index: /trunk/stac/src/Makefile
===================================================================
--- /trunk/stac/src/Makefile	(revision 2896)
+++ /trunk/stac/src/Makefile	(revision 2897)
@@ -5,8 +5,8 @@
 LDFLAGS += $(PSLIB)
 
-OBJECTS = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
-	stacCombine.o stacInvertMaps.o stacRejection.o psPlaneTransform.o stacAreaOfInterest.o
+STAC = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
+	stacCombine.o stacInvertMaps.o stacRejection.o psPlaneTransform.o stacAreaOfInterest.o stacSize.o
 
-TARGET = stac
+FITTRANS = fitTrans.o fitTransConfig.o fitTransRead.o stacCheckMemory.o fitTransSolve.o psPlaneTransform.o
 
 .PHONY: tags clean empty test profile optimise
@@ -15,6 +15,6 @@
 		$(CC) -c $(CFLAGS) $(OPTFLAGS) $<
 
-$(TARGET):	$(OBJECTS)
-		$(CC) $(CFLAGS) -o $@ $(OBJECTS) $(LDFLAGS) $(OPTFLAGS)
+stac:		$(STAC)
+		$(CC) $(CFLAGS) -o $@ $(STAC) $(LDFLAGS) $(OPTFLAGS)
 
 calcGrad:	calcGrad.o
@@ -27,4 +27,6 @@
 		$(CC) $(CFLAGS) -o $@ getPixels.o $(LDFLAGS) $(OPTFLAGS)
 
+fitTrans:	$(FITTRANS)
+		$(CC) $(CFLAGS) -o $@ $(FITTRANS) $(LDFLAGS) $(OPTFLAGS)
 
 clean:
Index: /trunk/stac/src/stac.c
===================================================================
--- /trunk/stac/src/stac.c	(revision 2896)
+++ /trunk/stac/src/stac.c	(revision 2897)
@@ -37,4 +37,5 @@
     psTraceSetLevel("stac.time",10);
     psTraceSetLevel("stac.area",10);
+    psTraceSetLevel("stac.size",10);
 
     // Set logging level
@@ -52,6 +53,9 @@
     psTrace("stac.time", 1, "I/O completed at %f seconds\n", getTime() - startTime);
 
+    // Get size
+    stacSize(config, inputs, maps);
+
     // Invert maps
-    psArray *inverseMaps = stacInvertMaps(maps);
+    psArray *inverseMaps = stacInvertMaps(maps, config);
 
     // Generate errors
Index: /trunk/stac/src/stac.h
===================================================================
--- /trunk/stac/src/stac.h	(revision 2896)
+++ /trunk/stac/src/stac.h	(revision 2897)
@@ -22,4 +22,5 @@
     const char *output;			// Output file name
     int outnx, outny;			// Size of output image
+    float xOffset, yOffset;		// Offset to get lower-left corner at 0,0
     float saturated;			// Saturation level
     float bad;				// Bad level
@@ -127,5 +128,6 @@
 
 // Invert an array of maps
-psArray *stacInvertMaps(const psArray *maps // Array of maps to invert
+psArray *stacInvertMaps(const psArray *maps, // Array of maps to invert
+			const stacConfig *config // Configuration
     );
 
@@ -168,3 +170,14 @@
 /************************************************************************************************************/
 
+// stacSize.c
+
+// Calculate the size of the output image
+bool stacSize(stacConfig *config,	// Configuration, containing the output size
+	      const psArray *images,	// Input images
+	      psArray *maps		// Transformation maps
+    );
+
+
+/************************************************************************************************************/
+
 #endif
Index: /trunk/stac/src/stacInvertMaps.c
===================================================================
--- /trunk/stac/src/stacInvertMaps.c	(revision 2896)
+++ /trunk/stac/src/stacInvertMaps.c	(revision 2897)
@@ -5,6 +5,8 @@
 
 #define MAX(x,y) ((x) > (y) ? (x) : (y))
-
-psArray *stacInvertMaps(const psArray *maps // Array of maps to invert
+#define NUM_GRID 20
+
+psArray *stacInvertMaps(const psArray *maps, // Array of maps to invert
+			const stacConfig *config // Configuration
     )
 {
@@ -14,38 +16,130 @@
     psTrace("stac.invertMaps", 1, "Inverting maps....\n");
 
-    for (int i = 0; i < nMaps; i++) {
-	// Can't handle higher order than linear yet
-	if (((psPlaneTransform*)maps->data[i])->x->nX != 2 ||
-	    ((psPlaneTransform*)maps->data[i])->x->nY != 2 ||
-	    ((psPlaneTransform*)maps->data[i])->y->nX != 2 ||
-	    ((psPlaneTransform*)maps->data[i])->y->nY != 2) {
-	    psError("stac.invertMaps",
-		    "STAC cannot currently support orders other than linear.\n");
-	    psFree(inverted);
-	    return NULL;
-	}
-
-	psPlaneTransform *newMap = psPlaneTransformAlloc(2, 2); // Inverted map
-	psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[i]; // Uninverted map
-
-	// Now, simply do a 2x2 matrix inversion
-
-	double a = oldMap->x->coeff[1][0];
-	double b = oldMap->x->coeff[0][1];
-	double c = oldMap->y->coeff[1][0];
-	double d = oldMap->y->coeff[0][1];
-	double e = oldMap->x->coeff[0][0];
-	double f = oldMap->y->coeff[0][0];
-
-	double invDet = 1.0 / (a * d - b * c); // Inverse of the determinant
-
-	// Not entirely sure why this works, but it appears to do so.......................................
-	newMap->x->coeff[1][0] = invDet * a;
-	newMap->x->coeff[0][1] = - invDet * b;
-	newMap->y->coeff[1][0] = - invDet * c;
-	newMap->y->coeff[0][1] = invDet * d;
-
-	newMap->x->coeff[0][0] = - invDet * (d * e + c * f);
-	newMap->y->coeff[0][0] = - invDet * (b * e + a * f);
+    // Coordinates for the transformations
+    psPlane *inCoord = psAlloc(sizeof(psPlane));
+    psPlane *outCoord = psAlloc(sizeof(psPlane));
+
+    for (int mapNum = 0; mapNum < nMaps; mapNum++) {
+
+	psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[mapNum]; // Uninverted map
+	// Check input
+	if (oldMap->x->nX != oldMap->x->nY) {
+	    psError("stac.invertMaps", "Polynomial order in x and y don't match for map %d\n", mapNum);
+	    return NULL;
+	}
+	if (oldMap->y->nX != oldMap->y->nY) {
+	    psError("stac.invertMaps", "Polynomial order in x and y don't match for map %d\n", mapNum);
+	    return NULL;
+	}
+	if (oldMap->x->nX != oldMap->y->nX) {
+	    psError("stac.invertMaps", "Polynomial order in x and y don't match for map %d\n", mapNum);
+	    return NULL;
+	}
+	int order = oldMap->x->nX;	// Polynomial order
+	psTrace("stac.invertMaps", 4, "Generating order %d polynomial inverse transformation.\n", order);
+	psPlaneTransform *newMap = psPlaneTransformAlloc(order, order);	// Inverted map
+
+	// Create fake polynomial to use in evaluation
+	psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
+	for (int i = 0; i < order; i++) {
+	    for (int j = 0; j < order; j++) {
+		fakePoly->coeff[i][j] = 1.0; // Set all coeffecients to 1
+		fakePoly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate
+	    }
+	}
+
+	// A grid of xin,yin --> xout,yout
+	psVector *xIn = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
+	psVector *yIn = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
+	psVector *xOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
+	psVector *yOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
+
+	// Create grid of points
+	for (int yint = 0; yint < NUM_GRID; yint++) {
+	    inCoord->y = (float)(yint * config->outny) / (float)(NUM_GRID - 1);
+	    for (int xint = 0; xint < NUM_GRID; xint++) {
+		inCoord->x = (float)(xint * config->outnx) / (float)(NUM_GRID - 1);
+
+		(void)psPlaneTransformApply(outCoord, oldMap, inCoord);
+
+		xOut->data.F32[yint*NUM_GRID + xint] = inCoord->x;
+		yOut->data.F32[yint*NUM_GRID + xint] = inCoord->y;
+		xIn->data.F32[yint*NUM_GRID + xint] = outCoord->x;
+		yIn->data.F32[yint*NUM_GRID + xint] = outCoord->y;
+	    }
+	}
+
+	// Initialise the matrix and vectors
+	int nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
+	psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
+	psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
+	psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
+	for (int i = 0; i < nCoeff; i++) {
+	    for (int j = 0; j < nCoeff; j++) {
+		matrix->data.F64[i][j] = 0.0;
+	    }
+	    xVector->data.F64[i] = 0.0;
+	    yVector->data.F64[i] = 0.0;
+	}
+
+	// Iterate over the grid points
+	for (int g = 0; g < NUM_GRID*NUM_GRID; g++) {
+
+	    // Iterate over the polynomial coefficients, accumulating the matrix and vectors
+	    for (int i = 0, ijIndex = 0; i < order; i++) {
+		for (int j = 0; j < order - i; j++, ijIndex++) {
+
+		    fakePoly->mask[i][j] = 0;
+		    double ijPoly = psDPolynomial2DEval(xIn->data.F32[g], yIn->data.F32[g], fakePoly);
+		    fakePoly->mask[i][j] = 1;
+
+		    for (int m = 0, mnIndex = 0; m < order; m++) {
+			for (int n = 0; n < order - m; n++, mnIndex++) {
+			    
+			    fakePoly->mask[m][n] = 0;
+			    double mnPoly = psDPolynomial2DEval(xIn->data.F32[g], yIn->data.F32[g], fakePoly);
+			    fakePoly->mask[m][n] = 1;
+
+			    matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
+			}
+		    }
+		    
+		    xVector->data.F64[ijIndex] += ijPoly * (double)xOut->data.F32[g];
+		    yVector->data.F64[ijIndex] += ijPoly * (double)yOut->data.F32[g];
+		}
+	    } // Iterating over coefficients
+	} // Iterating over grid points
+
+	// Solution via LU Decomposition
+	psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
+	psImage *luMatrix = psMatrixLUD(NULL, permutation, matrix); // LU decomposed matrix
+	psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
+	psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
+
+	// Stuff coefficients into transformation
+	for (int i = 0, ijIndex = 0; i < order; i++) {
+	    for (int j = 0; j < order - i; j++, ijIndex++) {
+		newMap->x->coeff[i][j] = xSolution->data.F64[ijIndex];
+		newMap->y->coeff[i][j] = ySolution->data.F64[ijIndex];
+	    }
+	}
+	inverted->data[mapNum] = newMap;
+
+#ifdef TESTING
+	// Print x coefficients
+	psTrace("stac.invertMaps", 7, "x' = \n");
+	for (int i = 0; i < order; i++) {
+	    for (int j = 0; j < order - i; j++) {
+		psTrace("stac.invertMaps", 7, "      %f x^%d y^%d\n", newMap->x->coeff[i][j], i, j);
+	    }
+	}
+	// Print y coefficients
+	psTrace("stac.invertMaps", 7, "y' = \n");
+	for (int i = 0; i < order; i++) {
+	    for (int j = 0; j < order - i; j++) {
+		psTrace("stac.invertMaps", 7, "      %f x^%d y^%d\n", newMap->y->coeff[i][j], i, j);
+	    }
+	}
+#endif
 
 #ifdef TESTING
@@ -64,6 +158,58 @@
 #endif
 
-	inverted->data[i] = newMap;	// Stuff into the array
+	psFree(permutation);
+	psFree(luMatrix);
+	psFree(matrix);
+	psFree(xVector);
+	psFree(yVector);
+	psFree(xSolution);
+	psFree(ySolution);
+	psFree(fakePoly);
+	psFree(xIn);
+	psFree(yIn);
+	psFree(xOut);
+	psFree(yOut);
     }
+
+    psFree(inCoord);
+    psFree(outCoord);
+		
+
+#if 0
+	// Can't handle higher order than linear yet
+	if (((psPlaneTransform*)maps->data[i])->x->nX != 2 ||
+	    ((psPlaneTransform*)maps->data[i])->x->nY != 2 ||
+	    ((psPlaneTransform*)maps->data[i])->y->nX != 2 ||
+	    ((psPlaneTransform*)maps->data[i])->y->nY != 2) {
+	    psError("stac.invertMaps",
+		    "STAC cannot currently support orders other than linear.\n");
+	    psFree(inverted);
+	    return NULL;
+	}
+
+	psPlaneTransform *newMap = psPlaneTransformAlloc(2, 2); // Inverted map
+	psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[i]; // Uninverted map
+
+	// Now, simply do a 2x2 matrix inversion
+
+	double a = oldMap->x->coeff[1][0];
+	double b = oldMap->x->coeff[0][1];
+	double c = oldMap->y->coeff[1][0];
+	double d = oldMap->y->coeff[0][1];
+	double e = oldMap->x->coeff[0][0];
+	double f = oldMap->y->coeff[0][0];
+
+	double invDet = 1.0 / (a * d - b * c); // Inverse of the determinant
+
+	// Not entirely sure why this works, but it appears to do so.......................................
+	newMap->x->coeff[1][0] = invDet * a;
+	newMap->x->coeff[0][1] = - invDet * b;
+	newMap->y->coeff[1][0] = - invDet * c;
+	newMap->y->coeff[0][1] = invDet * d;
+
+	newMap->x->coeff[0][0] = - invDet * (d * e + c * f);
+	newMap->y->coeff[0][0] = - invDet * (b * e + a * f);
+#endif
+
 
     return inverted;
Index: /trunk/stac/src/stacRead.c
===================================================================
--- /trunk/stac/src/stacRead.c	(revision 2896)
+++ /trunk/stac/src/stacRead.c	(revision 2897)
@@ -3,6 +3,4 @@
 #include "pslib.h"
 #include "stac.h"
-
-#define ORDER 2
 
 psArray *stacReadImages(stacConfig *config)
@@ -54,57 +52,75 @@
 	// Read the file
 	psTrace("stac.read.maps",5,"Reading map file %s....\n",mapfile);
-	psPlaneTransform *map = psPlaneTransformAlloc(ORDER,ORDER);
 
-	// Format is as following:
-	// A0  A1
-	// B00 B01
-	// B10 B11
-	//
-	// where:
-	// x_sky = A0 + B00*x + B10*y
-	// y_sky = A1 + B01*x + B11*y
+	// Format is now:
+	// order
+	// x coefficients
+	// y coefficients
+	// 
+	// where the order is 1 for linear, 2 for quadratic, 3 for cubic.
+	// and the coefficients are read by the following pseudo-code:
+	// 
+	// 	for (int k = 0; k < order + 1; k++)
+	// 	    for (int j = 0; j < k; j++)
+	// 		int i = k - j;
+	//              read coefficient of x^i y^j
+	// 
+	// This produces the ordering:
+	// x^0 y^0, x^1 y^0, x^0 y^1, x^2 y^0, x^1 y^1, x^0 y^2, x^3 y^0, x^2 y^1, x^1 y^2, x^0 y^3
+	// 
+	// This is, of course, for third order polynomials.
+	// For lower orders, the list is truncated at the appropriate level.
 
-	double a0, a1, b00, b10, b01, b11;
-
-	if (fscanf(mapfp,"%lf %lf %lf %lf %lf %lf", &a0, &a1, &b00, &b01, &b10, &b11) != 6) {
-	    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Couldn't read all 6 transformation values from %s\n",
-		     mapfile);
+	int order = 0;			// Polynomial order
+	if (fscanf(mapfp, "%d", &order) != 1) {
+	    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
 	    exit(EXIT_FAILURE);
 	}
-	fclose(mapfp);
-	map->x->coeff[0][0] = a0;
-	map->x->coeff[1][0] = b00;
-	map->x->coeff[0][1] = b10;
-	map->y->coeff[0][0] = a1;
-	map->y->coeff[1][0] = b01;
-	map->y->coeff[0][1] = b11;
 
-	// Mask out the cross-terms
-	map->x->mask[1][1] = 1;
-	map->y->mask[1][1] = 1;
-	// Set the cross-terms to zero, just in case...
-	map->x->coeff[1][1] = 0.0;
-	map->y->coeff[1][1] = 0.0;
+	psTrace("stac.read.maps", 5, "Polynomial order: %d\n", order);
+
+	psPlaneTransform *map = psPlaneTransformAlloc(order + 1, order + 1); // The transformation
+	// Set coefficient masks
+	for (int i = 0; i < order + 1; i++) {
+	    for (int j = 0; j < order + 1; j++) {
+		if (i + j > order + 1) {
+		    map->x->mask[i][j] = 1;
+		    map->y->mask[i][j] = 1;
+		} else {
+		    map->x->mask[i][j] = 0;
+		    map->y->mask[i][j] = 0;
+		}
+	    }
+	}
+
+	// Read x coefficients
+	psTrace("stac.read.maps", 7, "x' = \n");
+	for (int k = 0; k < order + 1; k++) {
+	    for (int j = 0; j < k + 1; j++) {
+		int i = k - j;
+		if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
+		    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
+		    exit(EXIT_FAILURE);
+		}
+		psTrace("stac.read.maps", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
+	    }
+	}
+	// Read y coefficients
+	psTrace("stac.read.maps", 7, "y' = \n");
+	for (int k = 0; k < order + 1; k++) {
+	    for (int j = 0; j < k + 1; j++) {
+		int i = k - j;
+		if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
+		    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
+		    exit(EXIT_FAILURE);
+		}
+		psTrace("stac.read.maps", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
+	    }
+	}
+
+   	fclose(mapfp);
 
 	// Plug it into the array
 	maps->data[i] = map;
-
-#ifdef TESTING
-	// Check psLib's answer is the same as I would get with old version
-	double detx = 256.0, dety = 256.0;
-	// My code
-	double myx = a0 + b00*detx + b10*dety;
-	double myy = a1 + b01*detx + b11*dety;
-	// psLib
-	psPlane *det = psAlloc(sizeof(psPlane));
-	det->x = 256.0;
-	det->y = 256.0;
-	psPlane *pslib = psPlaneTransformApply(NULL, map, det);
-
-	psTrace("stac.read.maps.test",0,"me: %lf,%lf\n",myx,myy);
-	psTrace("stac.read.maps.test",0,"ps: %lf,%lf\n",pslib->x,pslib->y);
-	psFree(det);
-	psFree(pslib);
-#endif
 
     }
Index: /trunk/stac/src/stacSize.c
===================================================================
--- /trunk/stac/src/stacSize.c	(revision 2897)
+++ /trunk/stac/src/stacSize.c	(revision 2897)
@@ -0,0 +1,103 @@
+#include <stdio.h>
+#include "pslib.h"
+#include "stac.h"
+
+bool stacSize(stacConfig *config,	// Configuration, containing the output size
+	      const psArray *images,	// Input images
+	      psArray *maps		// Transformation maps
+    )
+{
+    int nImages = images->n;		// Number of input images
+    if (nImages != maps->n) {
+	psError("stac.size", "Number of images (%d) and maps (%d) do not match.\n", nImages, maps->n);
+	return false;
+    }
+
+    psPlane *inCoord = psAlloc(sizeof(psPlane)); // Input coordinates
+    psPlane *outCoord = psAlloc(sizeof(psPlane)); // Output coordinates
+
+    // Initial "guess" at limits
+    float xMin = 0.0;
+    float xMax = 0.0;
+    float yMin = 0.0;
+    float yMax = 0.0;
+
+    for (int i = 0; i < nImages; i++) {
+	psImage *image = images->data[i]; // The image
+	psPlaneTransform *map = maps->data[i]; // The map
+
+	// Lower left corner
+	inCoord->x = 0;
+	inCoord->y = 0;
+	(void)psPlaneTransformApply(outCoord, map, inCoord);
+	if (outCoord->x < xMin) {
+	    xMin = outCoord->x;
+	} else if (outCoord->x > xMax) {
+	    xMax = outCoord->x;
+	}
+	if (outCoord->y < yMin) {
+	    yMin = outCoord->y;
+	} else if (outCoord->y > yMax) {
+	    yMax = outCoord->y;
+	}
+
+	// Lower right corner
+	inCoord->x = image->numCols;
+	(void)psPlaneTransformApply(outCoord, map, inCoord);
+	if (outCoord->x < xMin) {
+	    xMin = outCoord->x;
+	} else if (outCoord->x > xMax) {
+	    xMax = outCoord->x;
+	}
+	if (outCoord->y < yMin) {
+	    yMin = outCoord->y;
+	} else if (outCoord->y > yMax) {
+	    yMax = outCoord->y;
+	}
+
+	// Upper right corner
+	inCoord->y = image->numRows;
+	(void)psPlaneTransformApply(outCoord, map, inCoord);
+	if (outCoord->x < xMin) {
+	    xMin = outCoord->x;
+	} else if (outCoord->x > xMax) {
+	    xMax = outCoord->x;
+	}
+	if (outCoord->y < yMin) {
+	    yMin = outCoord->y;
+	} else if (outCoord->y > yMax) {
+	    yMax = outCoord->y;
+	}
+
+	// Upper left corner
+	inCoord->x = image->numCols;
+	(void)psPlaneTransformApply(outCoord, map, inCoord);
+	if (outCoord->x < xMin) {
+	    xMin = outCoord->x;
+	} else if (outCoord->x > xMax) {
+	    xMax = outCoord->x;
+	}
+	if (outCoord->y < yMin) {
+	    yMin = outCoord->y;
+	} else if (outCoord->y > yMax) {
+	    yMax = outCoord->y;
+	}
+    }
+
+    // Tweak the maps to account for the offset
+    for (int i = 0; i < nImages; i++) {
+	psPlaneTransform *map = maps->data[i]; // The map of interest
+	map->x->coeff[0][0] -= floor(xMin);
+	map->y->coeff[0][0] -= floor(yMin);
+    }
+
+    config->outnx = (int)(xMax + 0.5) - (int)xMin;
+    config->outny = (int)(yMax + 0.5) - (int)yMin;
+
+    psTrace("stac.size", 1, "Output size is to be %dx%d\n", config->outnx, config->outny);
+
+    psFree(inCoord);
+    psFree(outCoord);
+
+    return true;
+}