Index: trunk/pswarp/src/pswarpTransformReadout_Opt.c
===================================================================
--- trunk/pswarp/src/pswarpTransformReadout_Opt.c	(revision 13903)
+++ trunk/pswarp/src/pswarpTransformReadout_Opt.c	(revision 14768)
@@ -1,3 +1,5 @@
 # include "pswarp.h"
+
+#define SOURCE_ARRAY_BUFFER 100         // Size to grow the array of sources at a time
 
 // NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT
@@ -83,4 +85,5 @@
 
     // Iterate over the output image pixels (parent frame)
+    bool goodPixels = false;            // Any input pixels landing on the output image?
     for (int y = minY; y < maxY; y++) {
         if (y >= nextGridY) {
@@ -111,4 +114,6 @@
             if (inPix->y - inRow0 < 0) continue;
             if (inPix->y - inRow0 >= inImage->numRows) continue;
+
+            goodPixels = true;
 
             // XXX include mask
@@ -135,9 +140,98 @@
     }
 
-    output->data_exists = true;
-
     psFree(interp);
     psFree (inPix);
     psFree (grid);
+
+    // Transform sources
+    bool mdok;
+    psArray *inSources = psMetadataLookupPtr(&mdok, input->analysis, "PSPHOT.SOURCES"); // Sources in source
+    if (goodPixels && mdok && inSources) {
+        pswarpMapGrid *sourceGrid = pswarpMapGridFromImage(output, input, nGridX, nGridY); // Grid for sources
+
+        psArray *outSources = psMemIncrRefCounter(psMetadataLookupPtr(&mdok, output->analysis,
+                                                                      "PSPHOT.SOURCES")); // Target sources
+        if (!outSources) {
+            outSources = psArrayAllocEmpty(SOURCE_ARRAY_BUFFER);
+            psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY,
+                             "Warped sources", outSources);
+        }
+
+        for (int i = 0; i < inSources->n; i++) {
+            pmSource *source = inSources->data[i]; // Source of interest
+            pmModel *model = source->modelPSF; // Model for this source
+            float xIn, yIn;             // Coordinates of source
+            xIn = model->params->data.F32[PM_PAR_XPOS] - inCol0;
+            yIn = model->params->data.F32[PM_PAR_YPOS] - inRow0;
+            int xGrid, yGrid;           // Grid coordinates for local map
+            if (!pswarpMapGridSetGrid(sourceGrid, xIn + 0.5, yIn + 0.5, &xGrid, &yGrid)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to get grid coordinates for source at %f,%f\n",
+                        xIn, yIn);
+                psFree(outSources);
+                psFree(sourceGrid);
+                return false;
+            }
+            if (xGrid < 0 || xGrid >= sourceGrid->nXpts || yGrid < 0 || yGrid >= sourceGrid->nYpts) {
+                // It's not even on the grid
+                continue;
+            }
+
+            pswarpMap *map = sourceGrid->maps[xGrid][yGrid]; // Locally linear transformation
+            double xOut, yOut;          // Output coordinates
+            if (!pswarpMapApply(&xOut, &yOut, map, xIn, yIn)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to transform coordinates for source at %f,%f\n",
+                        xIn, yIn);
+                psFree(outSources);
+                psFree(sourceGrid);
+                return false;
+            }
+            xOut += outCol0;
+            yOut += outRow0;
+            if (xOut < minX || xOut > maxX || yOut < minY || yOut > maxY) {
+                // It's not in the output image
+                continue;
+            }
+
+            // Generate the new source in the output frame
+            pmSource *new = pmSourceAlloc(); // New source
+            new->peak = pmPeakAlloc(xOut, yOut, source->peak->flux, PM_PEAK_LONE);
+            new->peak->flux = source->peak->flux;
+            new->type = PM_SOURCE_TYPE_STAR;
+            new->psfMag = source->psfMag;
+            new->errMag = source->errMag;
+            new->sky = source->sky;
+            new->skyErr = source->skyErr;
+            new->pixWeight = source->pixWeight;
+            new->modelPSF = pmModelAlloc(source->modelPSF->type);
+
+#if 0
+            // XXX Note that this will not set the correct axes
+            pmPSF_AxesToModel(new->modelPSF->params->data.F32,
+                              pmPSF_ModelToAxes(source->modelPSF->params->data.F32, 20.0));
+#endif
+
+            // Propagate the position erorrs
+            float dxIn, dyIn;           // Errors in input coordinates
+            dxIn = model->dparams->data.F32[PM_PAR_XPOS];
+            dyIn = model->dparams->data.F32[PM_PAR_YPOS];
+
+            float dxOut, dyOut;         // Errors in output coordinates
+            dxOut = sqrt(PS_SQR(map->Xx * dxIn) + PS_SQR(map->Xy * dyIn));
+            dxOut = sqrt(PS_SQR(map->Xx * dxIn) + PS_SQR(map->Xy * dyIn));
+
+            new->modelPSF->params->data.F32[PM_PAR_XPOS] = xOut;
+            new->modelPSF->params->data.F32[PM_PAR_YPOS] = yOut;
+            new->modelPSF->dparams->data.F32[PM_PAR_XPOS] = dxOut;
+            new->modelPSF->dparams->data.F32[PM_PAR_YPOS] = dyOut;
+
+            psArrayAdd(outSources, SOURCE_ARRAY_BUFFER, new);
+            psFree(new);                // Drop reference
+        }
+        psFree(sourceGrid);
+        psFree(outSources);             // Drop reference
+    }
+
+    output->data_exists = true;
+
     psLogMsg ("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp"));