Index: /trunk/stac/src/Makefile
===================================================================
--- /trunk/stac/src/Makefile	(revision 3609)
+++ /trunk/stac/src/Makefile	(revision 3610)
@@ -1,10 +1,10 @@
 SHELL = /bin/sh
 CC = gcc
-CFLAGS += -g -std=c99 -I/home/mithrandir/price/psLib/include -DCRFLUX -DTESTING
-PSLIB += -L/home/mithrandir/price/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm -lxml2
+CFLAGS += -g -std=c99 -I/home/mithrandir/price/psLib/include -D_GNU_SOURCE # -DCRFLUX -DTESTING
+PSLIB += -L/home/mithrandir/price/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm -lxml2 -lmysqlclient
 LDFLAGS += $(PSLIB)
 
 STAC = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
-	stacCombine.o stacInvertMaps.o stacRejection.o stacAreaOfInterest.o stacSize.o
+	stacCombine.o stacInvertMaps.o stacRejection.o stacAreaOfInterest.o stacSize.o stacScales.o
 
 SHIFT = shift.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
@@ -15,4 +15,6 @@
 TEST_SHIFT = testShift.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
 	stacInvertMaps.o
+
+PHOTSCALE = photScale.o
 
 .PHONY: tags clean empty test profile optimise
@@ -26,4 +28,7 @@
 shift:		$(SHIFT)
 		$(CC) $(CFLAGS) -o $@ $(SHIFT) $(LDFLAGS) $(OPTFLAGS)
+
+photScale:	$(PHOTSCALE)
+		$(CC) $(CFLAGS) -o $@ $(PHOTSCALE) $(LDFLAGS) $(OPTFLAGS)
 
 calcGrad:	calcGrad.o
Index: /trunk/stac/src/stac.c
===================================================================
--- /trunk/stac/src/stac.c	(revision 3609)
+++ /trunk/stac/src/stac.c	(revision 3610)
@@ -27,4 +27,5 @@
     // Set trace levels
     psTraceSetLevel(".",0);
+    psTraceSetLevel(".psLib", 0);
     psTraceSetLevel("stac.checkMemory",10);
     psTraceSetLevel("stac.config",10);
@@ -38,4 +39,5 @@
     psTraceSetLevel("stac.area",10);
     psTraceSetLevel("stac.size",10);
+    psTraceSetLevel("stac.scales",7);
 
     // Set logging level
@@ -65,7 +67,14 @@
     psArray *transformed = NULL;
     psArray *transformedErrors = NULL;
-    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, config);
+    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, NULL, NULL,
+			config);
 
     psTrace("stac.time",1,"Transformation completed at %f seconds\n", getTime() - startTime);
+
+    psVector *scales = NULL;		// Relative scales between images
+    psVector *offsets = NULL;		// Offsets between images
+    (void)stacScales(&scales, &offsets, transformed, config);
+    // Rescale the images
+    (void)stacRescale(transformed, transformedErrors, NULL, scales, offsets);
 
     // Combine with rejection
@@ -122,7 +131,7 @@
     }
 
-    // Redo transformation with the masks
+    // Redo transformation with the masks and scales/offsets
     (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, rejectedSource,
-			combineRegion, config);
+			combineRegion, scales, offsets, config);
 
     // Combine the newly-transformed CR-free images, no rejection
@@ -139,6 +148,22 @@
     psFree(outFile);
 
+    // Save shifted images
+    if (config->saveShifts) {
+	for (int i = 0; i < transformed->n; i++) {
+	    char shiftName[MAXCHAR];	// Filename of shifted image
+	    sprintf(shiftName, "%s.shift", config->inputs->data[i]);
+	    psFits *shiftFile = psFitsAlloc(shiftName);
+	    if (!psFitsWriteImage(shiftFile, NULL, transformed->data[i], 0, NULL)) {
+		psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
+	    }
+	    psTrace("stac", 1, "Shifted image %d written to %s\n", i, shiftName);
+	    psFree(shiftFile);
+	}
+    }
+
     // Free everything I've used
     stacConfigFree(config);
+    psFree(scales);
+    psFree(offsets);
     psFree(combineRegion);
     psFree(regions);
Index: /trunk/stac/src/stac.h
===================================================================
--- /trunk/stac/src/stac.h	(revision 3609)
+++ /trunk/stac/src/stac.h	(revision 3610)
@@ -8,6 +8,6 @@
 #define MAXCHAR 80
 
-// Temporary implementation
-psPlaneTransform *psPlaneTransformAlloc(psS32 n1, psS32 n2);
+// Get the current time
+double getTime(void);
 
 /************************************************************************************************************/
@@ -21,8 +21,11 @@
     psArray *inputs;			// Input file names
     const char *output;			// Output file name
+    bool saveShifts;			// Save shifted images?
+    const char *starFile;		// File with star coordinates
+    const char *starMapFile;		// File with map for stars
     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
+    psVector *saturated;		// Saturation level for each image
+    psVector *bad;			// Bad level for each image
     float reject;			// Rejection level
     float frac;				// Fraction of input pixel that must be masked before the pixel is
@@ -30,4 +33,6 @@
     float grad;				// Multiplier of the gradient
     int nReject;			// Number of rejection iterations
+    float xMapDiff, yMapDiff;		// Difference between pure map and output image coordinates
+    float aper;				// Aperture size (pixels)
 } stacConfig;
 
@@ -52,4 +57,10 @@
 // Read the input files and return an array of images
 psArray *stacReadImages(stacConfig *config);
+
+// Read a file of coordinates (x,y)
+psArray *stacReadCoords(const char *filename);
+
+// Read a single map file and return the transformation
+psPlaneTransform *stacReadMap(const char *filename);
 
 // Read the map files and return an array of transformations
@@ -77,4 +88,6 @@
 		   const psArray *masks, // Masks of input images
 		   const psImage *region, // Region of interest for transformation
+		   const psVector *scales, // Relative scales
+		   const psVector *offsets, // Relative offsets
 		   const stacConfig *config // Configuration
     );
@@ -181,3 +194,26 @@
 /************************************************************************************************************/
 
+// stacScale.c
+
+// Calculate the background in an image
+float stacBackground(const psImage *image, // Image for which to get the background
+		     int sample		// Sample in increments of this value
+    );
+
+
+// Calculate the relative scales and offsets between the images
+bool stacScales(psVector **scalesPtr,	// Scales to return
+		psVector **offsetsPtr,	// Offsets to return
+		const psArray *images,	// Images on which to measure the scales and offsets
+		const stacConfig *config // Configuration
+    );
+
+// Rescale images
+bool stacRescale(psArray *images,	// Images to rescale
+		 psArray *errImages,	// Variance images to rescale
+		 const psImage *mask,	// Mask indicating which pixels to scale
+		 const psVector *scales,// Scales for images
+		 const psVector *offsets // Offsets for images
+    );
+
 #endif
Index: /trunk/stac/src/stacCombine.c
===================================================================
--- /trunk/stac/src/stacCombine.c	(revision 3609)
+++ /trunk/stac/src/stacCombine.c	(revision 3610)
@@ -68,6 +68,6 @@
     int numRows = ((psImage*)images->data[0])->numRows;	// Image size
     int numCols = ((psImage*)images->data[0])->numCols; // Image size
-    float saturated = config->saturated;// Saturation limit
-    float bad = config->bad;		// Bad pixel limit
+    psVector *saturated = config->saturated;// Saturation limits
+    psVector *bad = config->bad;	// Bad pixel limits
     float reject = config->reject;	// Rejection (k-sigma)
     
@@ -134,5 +134,6 @@
 		    pixels->data.F32[i] = pixel;
 		    deltas->data.F32[i] = delta;
-		    if ((pixel >= saturated) || (pixel <= bad) || (! isfinite(pixel)) || (! isfinite(delta))) {
+		    if ((pixel >= saturated->data.F32[i]) || (pixel <= bad->data.F32[i]) ||
+			(! isfinite(pixel)) || (! isfinite(delta))) {
 			mask->data.U8[i] = (psU8)1; // Don't use!
 		    } else {
Index: /trunk/stac/src/stacConfig.c
===================================================================
--- /trunk/stac/src/stacConfig.c	(revision 3609)
+++ /trunk/stac/src/stacConfig.c	(revision 3610)
@@ -9,5 +9,5 @@
 {
     fprintf (stderr, "STAC: Simultaneous Telescope Array Combination\n"
-	     "Usage: %s [-h] [-v] [-g GAIN] [-r RN] [-o NX NY] [-s SAT] [-b BAD] [-k REJ] [-n NREJECT] [-k FRAC] [-G GRAD] OUT IN1 IN2...\n"
+	     "Usage: %s [-h] [-v] [-g GAIN] [-r RN] [-o NX NY] [-S] [-s SAT] [-b BAD] [-p FILE MAP] [-a APER] [-k REJ] [-n NREJECT] [-k FRAC] [-G GRAD] OUT IN1 IN2...\n"
 	     "where\n"
 	     "\t-h           Help (this info)\n"
@@ -16,6 +16,9 @@
 	     "\t-r           Read noise (e) for detectors\n"
 	     "\t-o NX NY     Size of output image (512, 512)\n"
-	     "\t-s SAT       Saturation level (65536)\n"
+	     "\t-S           Save shifted images (false)\n"
+	     "\t-s SAT       Saturation level (65535)\n"
 	     "\t-b BAD       Bad level (0)\n"
+	     "\t-p FILE MAP  Specify file containing star coordinates, with map\n"
+	     "\t-a APER      Aperture size to use for photometry (3 pix)\n"
 	     "\t-k REJ       Rejection level (k-sigma; 3.5)\n"
 	     "\t-n NREJECT   Number of rejection iterations (2)\n"
@@ -39,8 +42,12 @@
     config->inputs = NULL;
     config->output = NULL;
+    config->saveShifts = false;
+    config->starFile = NULL;
+    config->starMapFile = NULL;
+    config->aper = 3.0;
     config->outnx = 512;
     config->outny = 512;
-    config->saturated = 65535.0;
-    config->bad = 0.0;
+    config->saturated = NULL;		// Saturations; will fill this in later, when we know how many images
+    config->bad = NULL;			// Bad levels; will fill this in later, when we know how many images
     config->reject = 3.0;
     config->frac = 0.5;
@@ -54,7 +61,13 @@
 void stacConfigFree(stacConfig *config)
 {
-    // Free the vector, if necessary
+    // Free the vectors, if necessary
     if (config->inputs) {
 	psFree(config->inputs);
+    }
+    if (config->saturated) {
+	psFree(config->saturated);
+    }
+    if (config->bad) {
+	psFree(config->bad);
     }
     // Free everything
@@ -69,4 +82,7 @@
     stacConfig *config = stacConfigAlloc(); // Configuration values
     const char *programName = argv[0];	// Program name
+
+    float saturated = 65535.0;		// Saturation level
+    float bad = 0.0;			// Bad level
 
     /* Variables for getopt */
@@ -76,5 +92,5 @@
 
     /* Parse command-line arguments using getopt */
-    while ((opt = getopt(argc, argv, "hvg:r:o:s:b:k:n:f:G:")) != -1) {
+    while ((opt = getopt(argc, argv, "hvSg:r:o:s:b:k:n:f:G:p:a:")) != -1) {
         switch (opt) {
 	  case 'h':
@@ -87,4 +103,5 @@
 	    if (sscanf(optarg, "%f", &config->gain) != 1) {
 		help(programName);
+		exit(EXIT_FAILURE);
 	    }
 	    break;
@@ -92,4 +109,5 @@
 	    if (sscanf(optarg, "%f", &config->readnoise) != 1) {
 		help(programName);
+		exit(EXIT_FAILURE);
 	    }
 	    break;
@@ -97,8 +115,5 @@
             if ((sscanf(argv[optind-1], "%d", &config->outnx) != 1) ||
                 (sscanf(argv[optind++], "%d", &config->outny) != 1)) {
-                /*
-                  Note: incrementing optind, so I can read more than
-                  one parameter.
-                */
+                // Note: incrementing optind, so I can read more than one parameter.
                 help(programName);
                 exit(EXIT_FAILURE);
@@ -106,11 +121,28 @@
 	    break;
 	  case 's':
-	    if (sscanf(optarg, "%f", &config->saturated) != 1) {
-		help(programName);
+	    if (sscanf(optarg, "%f", &saturated) != 1) {
+		help(programName);
+		exit(EXIT_FAILURE);
 	    }
 	    break;
 	  case 'b':
-	    if (sscanf(optarg, "%f", &config->bad) != 1) {
-		help(programName);
+	    if (sscanf(optarg, "%f", &bad) != 1) {
+		help(programName);
+		exit(EXIT_FAILURE);
+	    }
+	    break;
+	  case 'p':
+	    if (argc < optind+1) {
+		help(programName);
+		exit(EXIT_FAILURE);
+	    }
+	    config->starFile = argv[optind-1];
+	    config->starMapFile = argv[optind++];
+	    // Note: incrementing optind, so I can read more than one parameter.
+	    break;
+	  case 'a':
+	    if (sscanf(optarg, "%f", &config->aper) != 1) {
+		help(programName);
+		exit(EXIT_FAILURE);
 	    }
 	    break;
@@ -118,4 +150,5 @@
 	    if (sscanf(optarg, "%f", &config->reject) != 1) {
 		help(programName);
+		exit(EXIT_FAILURE);
 	    }
 	    break;
@@ -123,4 +156,5 @@
 	    if (sscanf(optarg, "%d", &config->nReject) != 1) {
 		help(programName);
+		exit(EXIT_FAILURE);
 	    }
 	    break;
@@ -128,4 +162,5 @@
 	    if (sscanf(optarg, "%f", &config->frac) != 1) {
 		help(programName);
+		exit(EXIT_FAILURE);
 	    }
 	    break;
@@ -133,5 +168,9 @@
 	    if (sscanf(optarg, "%f", &config->grad) != 1) {
 		help(programName);
-	    }
+		exit(EXIT_FAILURE);
+	    }
+	    break;
+	  case 'S':
+	    config->saveShifts = true;
 	    break;
 	  default:
@@ -156,4 +195,12 @@
     }
 
+    // Create the saturation and bad vectors
+    config->saturated = psVectorAlloc(argc-1, PS_TYPE_F32);
+    config->bad = psVectorAlloc(argc-1, PS_TYPE_F32);
+    for (int i = 0; i < argc - 1; i++) {
+	((psVector*)config->saturated)->data.F32[i] = saturated;
+	((psVector*)config->bad)->data.F32[i] = bad;
+    }
+
     // Debugging output
     psTrace("stac.config", 8, "Parsed command line for configuration\n");
Index: /trunk/stac/src/stacRead.c
===================================================================
--- /trunk/stac/src/stacRead.c	(revision 3609)
+++ /trunk/stac/src/stacRead.c	(revision 3610)
@@ -3,4 +3,6 @@
 #include "pslib.h"
 #include "stac.h"
+
+#define BUFFER 100			// Size of buffer for incrementally reading coordinates
 
 psArray *stacReadImages(stacConfig *config)
@@ -40,4 +42,123 @@
 }
 
+psArray *stacReadCoords(const char *filename)
+{
+    FILE *file = fopen(filename, "r");
+    if (file == NULL) {
+	psLogMsg("stac.read.coords", PS_LOG_ERROR, "Cannot open coordinate file, %s\n", filename);
+	return NULL;
+    }
+
+    psTrace("stac.read.coords", 5, "Reading coordinate file, %s\n", filename);
+
+    psArray *coords = psArrayAlloc(BUFFER); // The array of coordinates to be returned
+    float x, y;				// Coordinates to read
+    int num = 0;			// Number of coordinates read
+    while (fscanf(file, "%f %f\n", &x, &y) == 2) {
+	psPlane *coord = psPlaneAlloc();// A coordinate
+	coord->x = x;
+	coord->y = y;
+	coords->data[num] = coord;
+	num++;
+	if (num % BUFFER) {
+	    coords = psArrayRealloc(coords, num + BUFFER);
+	}
+    }
+    coords->n = num;
+
+    psTrace("stac.read.coords", 5, "%d coordinates read.\n", num);
+
+    fclose(file);
+    return coords;
+}
+
+
+psPlaneTransform *stacReadMap(const char *filename)
+{
+    FILE *mapfp = fopen(filename, "r");
+    if (mapfp == NULL) {
+	psLogMsg("stac.read.map", PS_LOG_ERROR, "Cannot open map file, %s\n", filename);
+	return NULL;
+    }
+    // Read the file
+    psTrace("stac.read.map", 5, "Reading map file %s....\n", filename);
+    
+    // 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.
+
+    int order = 0;			// Polynomial order
+    if (fscanf(mapfp, "%d", &order) != 1) {
+	psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
+	fclose(mapfp);
+	return NULL;
+    }
+    
+    psTrace("stac.read.map", 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.map", 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.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
+		fclose(mapfp);
+		psFree(map);
+		return NULL;
+	    }
+	    psTrace("stac.read.map", 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.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
+		fclose(mapfp);
+		psFree(map);
+		return NULL;
+	    }
+	    psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
+	}
+    }
+    
+    fclose(mapfp);
+
+    return map;
+}
+    
+
 
 psArray *stacReadMaps(stacConfig *config)
@@ -54,84 +175,10 @@
 	    exit(EXIT_FAILURE);
 	}
-	// Open the file
+	// Read the file
 	sprintf(mapfile,"%s.map",filenames->data[i]);
-	FILE *mapfp = fopen(mapfile,"r");
-	if (mapfp == NULL) {
-	    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Cannot open map file, %s\n",mapfile);
-	    exit(EXIT_FAILURE);
+	maps->data[i] = stacReadMap(mapfile);
+	if (maps->data[i] == NULL) {
+	    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map: %s\n", mapfile);
 	}
-	// Read the file
-	psTrace("stac.read.maps",5,"Reading map file %s....\n",mapfile);
-
-	// 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.
-
-	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);
-	}
-
-	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;
 
     }
Index: /trunk/stac/src/stacSize.c
===================================================================
--- /trunk/stac/src/stacSize.c	(revision 3609)
+++ /trunk/stac/src/stacSize.c	(revision 3610)
@@ -91,4 +91,6 @@
 
     // Tweak the maps to account for the offset
+    config->xMapDiff = floor(xMin);
+    config->yMapDiff = floor(yMin);
     for (int i = 0; i < nImages; i++) {
 	psPlaneTransform *map = maps->data[i]; // The map of interest
Index: /trunk/stac/src/stacTransform.c
===================================================================
--- /trunk/stac/src/stacTransform.c	(revision 3609)
+++ /trunk/stac/src/stacTransform.c	(revision 3610)
@@ -113,4 +113,6 @@
 		   const psArray *masks, // Masks of input images
 		   const psImage *region, // Region of interest for transformation
+		   const psVector *scales, // Relative scales
+		   const psVector *offsets, // Relative offsets
 		   const stacConfig *config // Configuration
     )
@@ -123,4 +125,8 @@
     assert(images->n == maps->n);
     assert(!errors || (images->n == errors->n));
+    assert(!scales || scales->n == images->n);
+    assert(!offsets || offsets->n == images->n);
+    assert(!scales || scales->type.type == PS_TYPE_F32);
+    assert(!offsets || offsets->type.type == PS_TYPE_F32);
 
     // Allocate the output images if required, otherwise check the number
@@ -170,4 +176,12 @@
 	psImage *outImage = (*outputs)->data[n]; // The output image
 	psImage *outError = (*outErrors)->data[n]; // The output error image
+	float offset = 0.0;		// Relative offset
+	float scale = 1.0;		// Relative scale
+	if (offsets) {
+	    offset = offsets->data.F32[n];
+	}
+	if (scales) {
+	    scale = scales->data.F32[n];
+	}
 
 #if 0
@@ -206,4 +220,8 @@
 											    detector->y,
 											    mask, 1, 0.0);
+
+		    outImage->data.F32[y][x] = (outImage->data.F32[y][x] - offset) / scale;
+		    outImage->data.F32[y][x] = outImage->data.F32[y][x] / SQUARE(scale);
+
 		} // Pixels of interest