Index: branches/pap/pswarp/src/pswarpTransformTile.c
===================================================================
--- branches/pap/pswarp/src/pswarpTransformTile.c	(revision 27708)
+++ branches/pap/pswarp/src/pswarpTransformTile.c	(revision 27933)
@@ -44,4 +44,5 @@
     args->yMax = PS_MIN_S32;
     args->covariance = NULL;
+    args->jacobian = NAN;
 
     return args;
@@ -50,6 +51,8 @@
 bool pswarpTransformTile(pswarpTransformTileArgs *args)
 {
-    psImage *inImage = args->input->image; ///< Input image
-    psImage *outImage = args->output->image; ///< Output image
+    pmReadout *input = args->input;        // Input readout
+    psImage *inImage = input->image, *inMask = input->mask; // Input images
+    pmReadout *output = args->output;      // Output readout
+    psImage *outImage = output->image, *outVariance = output->variance, *outMask = output->mask; // Outputs
 
     int inNumCols = inImage->numCols, inNumRows = inImage->numRows; ///< Size of input image
@@ -61,8 +64,8 @@
 
     // Dereference images for convenience
-    psF32 **outImageData     = args->output->image->data.F32;
-    psF32 **outVarData       = (args->output->variance) ? args->output->variance->data.F32 : NULL;
-    psImageMaskType **outMaskData = (args->output->mask)   ? args->output->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
-    psImageMaskType **inMaskData  = (args->input->mask)    ? args->input->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
+    psF32 **outImageData     = outImage->data.F32;
+    psF32 **outVarData       = outVariance ? outVariance->data.F32 : NULL;
+    psImageMaskType **outMaskData = outMask ? outMask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
+    psImageMaskType **inMaskData  = inMask ? inMask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
 
     pswarpMap *map = args->grid->maps[args->gridX][args->gridY]; ///< Map for this tile
@@ -74,4 +77,7 @@
     int yMin = PS_MAX(minPt.y, 0);
     int yMax = PS_MIN(maxPt.y, outNumRows);
+
+    double jacobian = map->Xx * map->Yy - map->Yx * map->Xy; // Jacobian of transformation
+    double jacobian2 = PS_SQR(jacobian);                     // Square Jacobian
 
     // Iterate over the output image pixels (parent frame)
@@ -87,7 +93,6 @@
             // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
             // both are in the parent frames of the input and output images.
-            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
+            double xIn, yIn;            // Input pixel coordinates
+            pswarpMapApply(&xIn, &yIn, map, x + 0.5, y + 0.5);
             if (xIn < 0 || xIn >= inNumCols || yIn < 0 || yIn >= inNumRows) {
                 continue;
@@ -105,8 +110,8 @@
 
             if (outImageData) {
-                outImageData[yOut][xOut] = imageValue;
+                outImageData[yOut][xOut] = imageValue * jacobian;
             }
             if (outVarData) {
-                outVarData[yOut][xOut] = varValue;
+                outVarData[yOut][xOut] = varValue * jacobian2;
             }
             if (outMaskData) {
@@ -122,10 +127,7 @@
         double xIn, yIn;                // Position of interest on input
         pswarpMapApply(&xIn, &yIn, map, xOut + 0.5, yOut + 0.5);
-        // XXX Why are we subtracting the *output* col0,row0 from the *input* coordinates?
-        // I expect these are zero, so that it makes no difference, but it's distracting.
-        xIn -= outCol0;
-        yIn -= outRow0;
         psKernel *kernel = psImageInterpolationKernel(xIn, yIn, args->interp->mode); // Interpolation kernel
         args->covariance = psImageCovarianceCalculate(kernel, args->input->covariance);
+        args->jacobian = sqrt(jacobian);
         psFree(kernel);
     }