Index: /branches/pap_branch_20090128/pswarp/src/pswarp.h
===================================================================
--- /branches/pap_branch_20090128/pswarp/src/pswarp.h	(revision 21281)
+++ /branches/pap_branch_20090128/pswarp/src/pswarp.h	(revision 21282)
@@ -19,4 +19,5 @@
 #define PSWARP_ANALYSIS_VARFACTOR "PSWARP.VARFACTOR" // Name for variance factor in analysis metadata
 #define PSWARP_ANALYSIS_GOODPIX   "PSWARP.GOODPIX" // Name for number of good pixels in analysis metadata
+#define PSWARP_ANALYSIS_COVARIANCES "PSWARP.COVARIANCES" // Name for covariance matrices on analysis MD
 
 // a single pswarpMap converts coordinates from one image to a second image
@@ -53,4 +54,5 @@
     long goodPixels;                    // Number of good pixels
     int xMin, xMax, yMin, yMax;         // Bounds of tile
+    psKernel *covariance;               // Covariance matrix
 } pswarpTransformTileArgs;
 
Index: /branches/pap_branch_20090128/pswarp/src/pswarpLoop.c
===================================================================
--- /branches/pap_branch_20090128/pswarp/src/pswarpLoop.c	(revision 21281)
+++ /branches/pap_branch_20090128/pswarp/src/pswarpLoop.c	(revision 21282)
@@ -273,4 +273,15 @@
     }
 
+    // Set covariance matrix for output
+    {
+        psList *covariances = psMetadataLookupPtr(&mdok, output->analysis,
+                                                  PSWARP_ANALYSIS_COVARIANCES); // Covariance matrices
+        psAssert(covariances, "Should be there");
+        psArray *covars = psListToArray(covariances); // Array of covariance matrices
+        output->covariance = psImageCovarianceAverage(covars);
+        psFree(covars);
+        psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES);
+    }
+
     if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, false)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts.");
Index: /branches/pap_branch_20090128/pswarp/src/pswarpTransformReadout.c
===================================================================
--- /branches/pap_branch_20090128/pswarp/src/pswarpTransformReadout.c	(revision 21281)
+++ /branches/pap_branch_20090128/pswarp/src/pswarpTransformReadout.c	(revision 21282)
@@ -26,4 +26,8 @@
     psImageMaskType maskBad  = pmConfigMaskGet("BAD.WARP", config);
     psAssert(mdok, "MASK.INPUT was not defined");
+
+    if (!input->covariance) {
+        input->covariance = psImageCovarianceNone();
+    }
 
     int nThreads = psMetadataLookupS32(&mdok, config->arguments, "NTHREADS"); // Number of threads
@@ -122,4 +126,11 @@
     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
+    psList *covariances = psMetadataLookupPtr(&mdok, output->analysis,
+                                              PSWARP_ANALYSIS_COVARIANCES); // Collection of covar. matrices
+    if (!covariances) {
+        covariances = psListAlloc(NULL);
+        psMetadataAddList(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_COVARIANCES, 0,
+                          "Collection of covariance matrices", covariances);
+    }
     while ((job = psThreadJobGetDone()) != NULL) {
         if (job->args->n < 1) {
@@ -133,4 +144,7 @@
             yMin = PS_MIN(args->yMin, yMin);
             yMax = PS_MAX(args->yMax, yMax);
+            if (args->covariance) {
+                psListAdd(covariances, PS_LIST_TAIL, args->covariance);
+            }
         }
         psFree(job);
Index: /branches/pap_branch_20090128/pswarp/src/pswarpTransformTile.c
===================================================================
--- /branches/pap_branch_20090128/pswarp/src/pswarpTransformTile.c	(revision 21281)
+++ /branches/pap_branch_20090128/pswarp/src/pswarpTransformTile.c	(revision 21282)
@@ -8,4 +8,5 @@
     psFree(args->interp);
     psFree(args->region);
+    psFree(args->covariance);
     return;
 }
@@ -30,4 +31,5 @@
     args->yMin = PS_MAX_S32;
     args->yMax = PS_MIN_S32;
+    args->covariance = NULL;
 
     return args;
@@ -104,4 +106,17 @@
     }
 
+    if (goodPixels > 0) {
+        float xOut = 0.5 * (xMin + xMax), yOut = 0.5 * (yMin + yMax); // Position of interest on output
+        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->output->covariance);
+        psFree(kernel);
+    }
+
     args->goodPixels = goodPixels;
     args->xMin = xMin;