Index: /trunk/pswarp/src/pswarp.h
===================================================================
--- /trunk/pswarp/src/pswarp.h	(revision 19394)
+++ /trunk/pswarp/src/pswarp.h	(revision 19395)
@@ -35,5 +35,5 @@
     int nXpts, nYpts;                   // number of x,y samples in the grid
     int nXpix, nYpix;                   // x,y spacing in src image pixels of grid samples
-    double xMin,  yMin;			// coordinate of first grid sample
+    double xMin,  yMin;                 // coordinate of first grid sample
 } pswarpMapGrid;
 
@@ -51,5 +51,6 @@
 
     // output values for this tile
-    long goodPixels;
+    long goodPixels;                    // Number of good pixels
+    int xMin, xMax, yMin, yMax;         // Bounds of tile
 } pswarpTransformTileArgs;
 
Index: /trunk/pswarp/src/pswarpTransformReadout.c
===================================================================
--- /trunk/pswarp/src/pswarpTransformReadout.c	(revision 19394)
+++ /trunk/pswarp/src/pswarpTransformReadout.c	(revision 19395)
@@ -1,3 +1,3 @@
-# include "pswarp.h"
+#include "pswarp.h"
 
 // NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT
@@ -5,16 +5,17 @@
 {
     // XXX this implementation currently ignores the use of the region
-    psImage *region = NULL;
+    psImage *region = NULL;             // Region to transform
 
-    psTimerStart ("warp");
+    psTimerStart("warp");
 
     // Get warp parameters
-    bool mdok;
-    int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX");
-    int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY");
-    psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.MODE");
+    bool mdok;                          // Status of MD lookup
+    int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX"); // Number of grid points in x
+    int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY"); // Number of grid points in y
+    psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments,
+                                                                   "INTERPOLATION.MODE"); // Mode for interp
 
     // load the recipe
-    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
+    psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PSWARP_RECIPE);
     psAssert (recipe, "missing recipe %s", PSWARP_RECIPE);
 
@@ -23,22 +24,47 @@
     psMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config);
     psMaskType maskBad  = pmConfigMaskGet("BAD.WARP", config);
-    psAssert (mdok, "MASK.INPUT was not defined");
+    psAssert(mdok, "MASK.INPUT was not defined");
 
-    int nThreads = psMetadataLookupS32 (&mdok, config->arguments, "NTHREADS");
-    if (!mdok) nThreads = 0;
-
-    // Flux fraction for "poor"
-    float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC");
+    int nThreads = psMetadataLookupS32(&mdok, config->arguments, "NTHREADS"); // Number of threads
+    if (!mdok) {
+        nThreads = 0;
+    }
+    float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); // Flux fraction for "poor"
 
     // pswarpMapGridFromImage builds a set of locally-linear maps which convert the
     // output coordinates to input coordinates
-    pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, nGridX, nGridY);
+    pswarpMapGrid *grid = pswarpMapGridFromImage(input, output, nGridX, nGridY);
 
     // XXX optionally modify the grid based on this result and force the maxError < XXX
-    double maxError = pswarpMapGridMaxError (grid);
-    psLogMsg ("pswarp", 3, "maximum error using this grid sampling: %f\n", maxError);
+    double maxError = pswarpMapGridMaxError(grid); // Maximum (positional) error from using grid
+    psLogMsg("pswarp", 3, "maximum error using this grid sampling: %lf\n", maxError);
+
+    // Get range of interest
+    int xOutMin, xOutMax, yOutMin, yOutMax; // Output pixel range
+    pswarpMatchRange(&xOutMin, &yOutMin, &xOutMax, &yOutMax, input, output);
+
+    // Check the range of the grid coordinates
+#define CHECK_GRID_RANGE() { \
+        xGridMin = PS_MIN(xGridMin, xGrid); \
+        xGridMax = PS_MAX(xGridMax, xGrid); \
+        yGridMin = PS_MIN(yGridMin, yGrid); \
+        yGridMax = PS_MAX(yGridMax, yGrid); \
+    }
+
+    int xGridMin = nGridX, xGridMax = 0, yGridMin = nGridY, yGridMax = 0; // Grid range
+    int xGrid, yGrid;                   // Grid coordinates
+    pswarpMapGridSetGrid(grid, xOutMin, yOutMin, &xGrid, &yGrid);
+    CHECK_GRID_RANGE();
+    pswarpMapGridSetGrid(grid, xOutMin, yOutMax, &xGrid, &yGrid);
+    CHECK_GRID_RANGE();
+    pswarpMapGridSetGrid(grid, xOutMax, yOutMin, &xGrid, &yGrid);
+    CHECK_GRID_RANGE();
+    pswarpMapGridSetGrid(grid, xOutMax, yOutMax, &xGrid, &yGrid);
+    CHECK_GRID_RANGE();
 
     // Interpolation options : move these from the arguments to explicit assignments
-    psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(interpolationMode, input->image, input->weight, input->mask, maskIn, NAN, NAN, maskBad, maskPoor, poorFrac);
+    psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(interpolationMode, input->image,
+                                                                       input->weight, input->mask, maskIn,
+                                                                       NAN, NAN, maskBad, maskPoor, poorFrac);
 
     if (input->weight && !output->weight) {
@@ -51,29 +77,13 @@
     }
 
-    static int maskNum = 0;
-    {
-        psString name = NULL;
-        psStringAppend(&name, "mask_%d.fits", maskNum++);
-        psFits *fits = psFitsOpen(name, "w");
-        psFree(name);
-        psFitsWriteImage(fits, NULL, output->mask, 0, NULL);
-        psFitsClose(fits);
-    }
-
-    // total number of good pixels across all tiles (summed below)
-    int goodPixels = 0;
-
     // create jobs and supply them to the threads
-    for (int gridY = 0; gridY < grid->nYpts; gridY++) {
-        for (int gridX = 0; gridX < grid->nXpts; gridX++) {
-
+    for (int gridY = yGridMin; gridY < yGridMax; gridY++) {
+        for (int gridX = xGridMin; gridX < xGridMax; gridX++) {
             pswarpTransformTileArgs *args = pswarpTransformTileArgsAlloc();
-
-            // these items are just views to the data; they are not freed with args
-            args->input = input;
-            args->output = output;
-            args->grid = grid;
-            args->interp = interp;
-            args->region = region;
+            args->input = psMemIncrRefCounter(input);
+            args->output = psMemIncrRefCounter(output);
+            args->grid = psMemIncrRefCounter(grid);
+            args->interp = psMemIncrRefCounter(interp);
+            args->region = psMemIncrRefCounter(region);
 
             args->gridX = gridX;
@@ -83,12 +93,6 @@
             // allocate a job
             psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE");
-
-            // construct the arguments for this job
-            // job is pswarpTransformTile (gridX, gridY);
-            psArrayAdd (job->args, 1, args);
-            // fprintf (stderr, "adding job %d,%d, Nargs: %ld\n", gridX, gridY, job->args->n);
-
-            // call: pswarpTransformTile (args);
-            if (!psThreadJobAddPending (job)) {
+            psArrayAdd(job->args, 1, args);
+            if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to warp image.");
                 return false;
@@ -109,4 +113,6 @@
     // we have only supplied one type of job, so we can assume the types here
     psThreadJob *job = NULL;
+    int xMin = output->image->numCols, xMax = 0, yMin = output->image->numRows, yMax = 0; // Bounds
+    int goodPixels = 0;                 // total number of good pixels across all tiles
     while ((job = psThreadJobGetDone()) != NULL) {
         if (job->args->n < 1) {
@@ -116,15 +122,26 @@
             // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n);
             goodPixels += args->goodPixels;
+            xMin = PS_MIN(args->xMin, xMin);
+            xMax = PS_MAX(args->xMax, xMax);
+            yMin = PS_MIN(args->yMin, yMin);
+            yMax = PS_MAX(args->yMax, yMax);
         }
-        psFree (job);
+        psFree(job);
     }
     psFree(interp);
     psFree(grid);
 
+    if (xMin < xMax && yMin < yMax) {
+        psTrace("pswarp.transform", 1, "Bounds [%d:%d,%d:%d]\n", xMin, xMax, yMin, yMax);
+    } else {
+        psTrace("pswarp.transform", 1, "No overlap\n");
+    }
 
     // Store the variance factor and number of good pixels
     if (goodPixels > 0) {
         // Variance factor: large factor --> small scale
-        float varFactor = psImageInterpolateVarianceFactor(input->image->numCols / 2.0 + input->image->col0, input->image->numRows / 2.0 + input->image->row0, interp);
+        float varFactor = psImageInterpolateVarianceFactor(input->image->numCols / 2.0 + input->image->col0,
+                                                           input->image->numRows / 2.0 + input->image->row0,
+                                                           interp);
         psMetadataItem *vfItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_VARFACTOR);
         if (vfItem) {
@@ -145,25 +162,16 @@
 
     if (goodPixels > 0) {
-        if (!pswarpTransformSources (output, input, config)) {
+        if (!pswarpTransformSources(output, input, config)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
             return false;
         }
+
+        // Data is only written out if there are good pixels
+        output->data_exists = true;
+        output->parent->data_exists = true;
+        output->parent->parent->data_exists = true;
     }
 
-    // XXX should we not write anything out if there are no good pixels?
-    output->data_exists = true;
-    output->parent->data_exists = true;
-    output->parent->parent->data_exists = true;
-
-    {
-        psString name = NULL;
-        psStringAppend(&name, "mask_%d.fits", maskNum++);
-        psFits *fits = psFitsOpen(name, "w");
-        psFree(name);
-        psFitsWriteImage(fits, NULL, output->mask, 0, NULL);
-        psFitsClose(fits);
-    }
-
-    psLogMsg ("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp"));
+    psLogMsg("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp"));
 
     return true;
Index: /trunk/pswarp/src/pswarpTransformTile.c
===================================================================
--- /trunk/pswarp/src/pswarpTransformTile.c	(revision 19394)
+++ /trunk/pswarp/src/pswarpTransformTile.c	(revision 19395)
@@ -1,13 +1,18 @@
-# include "pswarp.h"
+#include "pswarp.h"
 
-void pswarpTransformTileArgsFree (pswarpTransformTileArgs *args) {
-    if (!args) return;
+static void transformTileArgsFree(pswarpTransformTileArgs *args)
+{
+    psFree(args->input);
+    psFree(args->output);
+    psFree(args->grid);
+    psFree(args->interp);
+    psFree(args->region);
     return;
 }
 
-pswarpTransformTileArgs *pswarpTransformTileArgsAlloc() {
-
-    pswarpTransformTileArgs *args = (pswarpTransformTileArgs *)psAlloc(sizeof(pswarpTransformTileArgs));
-    psMemSetDeallocator(args, (psFreeFunc)pswarpTransformTileArgsFree);
+pswarpTransformTileArgs *pswarpTransformTileArgsAlloc()
+{
+    pswarpTransformTileArgs *args = psAlloc(sizeof(pswarpTransformTileArgs));
+    psMemSetDeallocator(args, (psFreeFunc)transformTileArgsFree);
 
     args->input = NULL;
@@ -21,46 +26,38 @@
 
     args->goodPixels = 0;
+    args->xMin = PS_MAX_S32;
+    args->xMax = PS_MIN_S32;
+    args->yMin = PS_MAX_S32;
+    args->yMax = PS_MIN_S32;
 
     return args;
 }
 
-bool pswarpTransformTile (pswarpTransformTileArgs *args) {
+bool pswarpTransformTile(pswarpTransformTileArgs *args)
+{
+    psImage *inImage = args->input->image; // Input image
+    psImage *outImage = args->output->image; // Output image
 
-    // int inCol0 = args->input->image->col0;
-    // int inRow0 = args->input->image->row0;
-    int inNcol = args->input->image->numCols;
-    int inNrow = args->input->image->numRows;
+    int inNumCols = inImage->numCols, inNumRows = inImage->numRows; // Size of input image
+    int outNumCols = outImage->numCols, outNumRows = outImage->numRows; // Size of output image
+    int outCol0 = outImage->col0, outRow0 = outImage->row0; // Offset of output image
 
-    int outCol0 = args->output->image->col0;
-    int outRow0 = args->output->image->row0;
-    int outNcol = args->output->image->numCols;
-    int outNrow = args->output->image->numRows;
+    psPlane minPt, maxPt;               // Minimum and maximum points for this tile
+    pswarpMapGridCoordRange(args->grid, args->gridX, args->gridY, &minPt, &maxPt);
 
-    // get the coordinate range for this grid tile
-    psPlane minPt, maxPt;
-    pswarpMapGridCoordRange (args->grid, args->gridX, args->gridY, &minPt, &maxPt);
-
-    psF32 **outImageData     = (args->output->image)  ? args->output->image->data.F32 : NULL;
+    // Dereference images for convenience
+    psF32 **outImageData     = args->output->image->data.F32;
     psF32 **outVarData       = (args->output->weight) ? args->output->weight->data.F32 : NULL;
     psMaskType **outMaskData = (args->output->mask)   ? args->output->mask->data.PS_TYPE_MASK_DATA : NULL;
     psMaskType **inMaskData  = (args->input->mask)    ? args->input->mask->data.PS_TYPE_MASK_DATA : NULL;
 
-    pswarpMap *map = args->grid->maps[args->gridX][args->gridY];
-
-    psImage *region = args->region;
-
-    double xInRaw, yInRaw;
-
-    // output values for this pixel
-    double imageValue;
-    double varValue;
-    psMaskType maskValue;
+    pswarpMap *map = args->grid->maps[args->gridX][args->gridY]; // Map for this tile
+    psImage *region = args->region;     // Region to transform
 
     // Bounds for iteration
     int xMin = PS_MAX(minPt.x, 0);
-    int xMax = PS_MIN(maxPt.x, outNcol);
+    int xMax = PS_MIN(maxPt.x, outNumCols);
     int yMin = PS_MAX(minPt.y, 0);
-    int yMax = PS_MIN(maxPt.y, outNrow);
-
+    int yMax = PS_MIN(maxPt.y, outNumRows);
 
     // Iterate over the output image pixels (parent frame)
@@ -70,22 +67,20 @@
 
             // Only transform those pixels requested
-            if (region && region->data.U8[y][x]) continue;
+            if (region && region->data.U8[y][x]) {
+                continue;
+            }
 
             // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
             // both are in the parent frames of the input and output images.
-            pswarpMapApply (&xInRaw, &yInRaw, map, x + 0.5, y + 0.5);
-
-            double xIn = xInRaw - outCol0;      // Position on input image
-            double yIn = yInRaw - outRow0;      // Position on input image
-
-            if (xIn < 0) continue;
-            if (yIn >= inNcol) continue;
-            if (yIn < 0) continue;
-            if (yIn >= inNrow) continue;
-
-            goodPixels++;
+            double xInRaw, yInRaw;      // Input raw pixel coordinates
+            pswarpMapApply(&xInRaw, &yInRaw, map, x + 0.5, y + 0.5);
+            double xIn = xInRaw - outCol0, yIn = yInRaw - outRow0; // Position on input image
+            if (xIn < 0 || xIn >= inNumCols || yIn < 0 || yIn >= inNumRows) {
+                continue;
+            }
 
             // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates
-            maskValue = inMaskData ? inMaskData[(int)yIn][(int)xIn] : 0;
+            double imageValue, varValue; // Value of image and variance map
+            psMaskType maskValue = inMaskData ? inMaskData[(int)yIn][(int)xIn] : 0; // Value of mask
             if (!psImageInterpolate(&imageValue, &varValue, &maskValue, xIn, yIn, args->interp)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
@@ -93,8 +88,6 @@
             }
 
-            int xOut = x - outCol0;     // Position on output image
-            int yOut = y - outRow0;     // Position on output image
+            int xOut = x - outCol0, yOut = y - outRow0; // Position on output image
 
-            // not all images need be transformed
             if (outImageData) {
                 outImageData[yOut][xOut] = imageValue;
@@ -106,7 +99,15 @@
                 outMaskData[yOut][xOut] = maskValue;
             }
+
+            goodPixels++;
         }
     }
+
     args->goodPixels = goodPixels;
+    args->xMin = xMin;
+    args->xMax = xMax;
+    args->yMin = yMin;
+    args->yMax = yMax;
+
     return true;
 }