Index: /branches/eam_branch_20080719/pswarp/src/pswarpThreads.c
===================================================================
--- /branches/eam_branch_20080719/pswarp/src/pswarpThreads.c	(revision 18666)
+++ /branches/eam_branch_20080719/pswarp/src/pswarpThreads.c	(revision 18666)
@@ -0,0 +1,52 @@
+# include "pswarp.h"
+
+typedef enum {
+  PSWARP_TRANSFORM_TILE
+} pswarpJobType;
+
+typedef struct {
+  pswarpJobType type;
+  void *opts;
+} pswarpJob;
+
+static psArray *jobs = NULL;
+static pthread_t *threads = NULL;
+
+bool pswarpThreadsAddJob (pwarpJob *job) {
+
+  if (jobs == NULL) {
+    jobs = psArrayAlloc (16);
+  }
+
+  psArrayAdd (jobs, job);
+  return true;
+}
+
+bool pswarpThreadsLaunchJobs () {
+
+  while () {
+    while ((job = psListGetAndRemove (jobs, PS_LIST_HEAD)) == NULL) {
+      usleep (50000);
+    }
+
+    switch (job->type) {
+      case PSWARP_TRANSFORM_TILE:
+	status = pswarpTransformTile ((pswarpTransformTileOpts *)job->opts);
+	// send a message somewhere if the job fails
+	break;
+    }
+  }  
+}
+
+// create a pool of Nthreads
+bool pswarpThreadsCreate (int nThreads) {
+
+  threads = (pthread_t *) psAlloc (sizeof(pthread_t));
+
+  for (int i = 0; i < nThreads; i++) {
+    pthread_create (&threads[i], NULL, &function, NULL);
+  }
+
+  return true;
+}
+
Index: /branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout_Threaded.c
===================================================================
--- /branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout_Threaded.c	(revision 18666)
+++ /branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout_Threaded.c	(revision 18666)
@@ -0,0 +1,287 @@
+# 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
+bool pswarpTransformReadout_Threaded(pmReadout *output, pmReadout *input, pmConfig *config)
+{
+    // XXX this implementation currently ignores the use of the region
+    psImage *region = NULL;
+
+    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");
+    // load the recipe
+    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
+    if (!recipe) {
+        psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE);
+        return false;
+    }
+
+    // output mask bits
+    psMaskType maskIn   = psMetadataLookupU8(&mdok, recipe, "MASK.INPUT"); 
+    psMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config); 
+    psMaskType maskBad  = pmConfigMaskGet("BAD.WARP", config); 
+    psAssert (mdok, "MASK.INPUT was not defined");
+
+    float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); // Flux fraction for "poor"
+
+    // we need to apply the offset to convert parent coordinates to child coordinates for
+    // psImagePixelInterpolate below
+    int inCol0 = input->image->col0;
+    int inRow0 = input->image->row0;
+    int outCol0 = output->image->col0;
+    int outRow0 = output->image->row0;
+
+    // we might want to do the rectangular regions outside of the selection independently
+    // psImageInit (output->image, NAN);
+
+    // find the output pixel range
+    int minX, minY, maxX, maxY;
+    pswarpMatchRange (&minX, &minY, &maxX, &maxY, input, output);
+
+    // pswarpMapGridFromImage builds a set of locally-linear maps which convert the
+    // output coordinates to input coordinates
+    pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, nGridX, nGridY);
+
+    // XXX need to 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);
+
+    int gridX, gridY, nextGridX, nextGridY;
+    pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY);
+    nextGridY = pswarpMapGridNextGrid_Y (grid, gridY);
+    pswarpMap *map = grid->maps[gridX][gridY];
+
+    assert ((int)(minX - outCol0) >= 0);
+    assert ((int)(maxX - outCol0) <= output->image->numCols);
+    assert ((int)(minY - outRow0) >= 0);
+    assert ((int)(maxY - outRow0) <= output->image->numRows);
+
+    int gridXo = gridX;
+    int nextGridXo = pswarpMapGridNextGrid_X (grid, gridX);
+
+    psImage *inImage = input->image;    // Input image
+    psImage *inVar   = input->weight;   // Input weight map
+    psImage *inMask  = input->mask;     // Input mask
+
+    // Interpolation options
+    psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(interpolationMode, inImage, 
+                                                                       inVar, inMask, maskIn, NAN, NAN,
+                                                                       maskBad, maskPoor, poorFrac);
+
+    psPlane *inPix = psPlaneAlloc();    // Coordinates on the input detector
+    psF32 **outImageData = output->image->data.F32; // Output image pixels
+    psF32 **outVarData = NULL;          // Output variance pixels
+    psMaskType **outMaskData = NULL;    // Output mask pixels
+    if (inVar) {
+        if (!output->weight) {
+            output->weight = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_F32);
+            psImageInit(output->weight, NAN);
+        }
+        outVarData = output->weight->data.F32;
+    }
+    if (inMask || maskPoor || maskBad) {
+        if (!output->mask) {
+            output->mask = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_MASK);
+            psImageInit(output->mask, maskBad);
+        }
+        outMaskData = output->mask->data.PS_TYPE_MASK_DATA;
+    }
+
+    // create N threads
+    for (gridY = gridYmin; gridY < gridYmax; gridY++) {
+	for (gridX = gridXmin; gridX < gridXmax; gridX++) {
+	    pswarpTransformTile (gridX, gridY);
+	    // send function and arguments to thread manager
+	    // threadFunction (pswarpTransformTile)
+	    // while (threadFull()) {
+	    // usleep (50000);
+	    // }
+	}
+    }
+
+
+    // Iterate over the output image pixels (parent frame)
+    long goodPixels = 0;                // Number of input pixels landing on the output image
+    for (int y = minY; y < maxY; y++) {
+        if (y >= nextGridY) {
+            gridY ++;
+            nextGridY = pswarpMapGridNextGrid_Y (grid, gridY);
+            map = grid->maps[gridX][gridY];
+        }
+
+        gridX = gridXo;
+        nextGridX = nextGridXo;
+        map = grid->maps[gridX][gridY];
+        int yOut = y - outRow0;         // Position on image
+        for (int x = minX; x < maxX; x++) {
+            if (x >= nextGridX) {
+                gridX ++;
+                nextGridX = pswarpMapGridNextGrid_X (grid, gridX);
+                map = grid->maps[gridX][gridY];
+            }
+
+            // Only transform those pixels requested
+            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 (&inPix->x, &inPix->y, map, x + 0.5, y + 0.5);
+
+            if (inPix->x - inCol0 < 0) continue;
+            if (inPix->x - inCol0 >= inImage->numCols) continue;
+            if (inPix->y - inRow0 < 0) continue;
+            if (inPix->y - inRow0 >= inImage->numRows) continue;
+
+            goodPixels++;
+
+            // XXX include mask
+            // XXX apply scale and offset?
+            // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates
+            double imageValue, varValue;
+            psMaskType maskValue;
+            if (!psImageInterpolate(&imageValue, &varValue, &maskValue,
+                                    inPix->x - inCol0, inPix->y - inRow0, interp)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
+                psFree(interp);
+                psFree(inPix);
+                psFree(grid);
+                return false;
+            }
+            int xOut = x - outCol0;     // Position on image
+            outImageData[yOut][xOut] = imageValue;
+            if (inVar) {
+                outVarData[yOut][xOut] = varValue;
+            }
+            if (outMaskData) {
+                outMaskData[yOut][xOut] = maskValue;
+            }
+        }
+    }
+    psFree(interp);
+    psFree(inPix);
+    psFree(grid);
+
+    // Store the variance factor and number of good pixels
+    if (goodPixels > 0) {
+        float varFactor = psImageInterpolateVarianceFactor(inImage->numCols / 2.0 + inCol0,
+                                                           inImage->numRows / 2.0 + inRow0,
+                                                           interp); // Variance factor: large --> small scale
+        psMetadataItem *vfItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_VARFACTOR);
+        if (vfItem) {
+            psMetadataItem *goodpixItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_GOODPIX);
+            psAssert(goodpixItem, "It should be where we left it!");
+            psAssert(vfItem->type == PS_TYPE_F32 && goodpixItem->type == PS_TYPE_S64,
+                     "Should be the type we said.");
+
+            vfItem->data.F32 += varFactor * goodPixels;
+            goodpixItem->data.S64 += goodPixels;
+        } else {
+            psMetadataAddF32(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_VARFACTOR, 0,
+                             "Variance factor weighted by the good pixels", varFactor * goodPixels);
+            psMetadataAddS64(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_GOODPIX, 0,
+                             "Number of good pixels", goodPixels);
+        }
+    }
+
+    // Transform sources
+    psArray *inSources = psMetadataLookupPtr(&mdok, input->analysis, "PSPHOT.SOURCES"); // Sources in source
+    if (goodPixels > 0 && 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 + 0.5, yIn + 0.5)) {
+                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 - 0.5;
+            yOut += outRow0 - 0.5;
+            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));
+            dyOut = sqrt(PS_SQR(map->Yx * dxIn) + PS_SQR(map->Yy * 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;
+    output->parent->data_exists = true;
+    output->parent->parent->data_exists = true;
+
+    psLogMsg ("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp"));
+
+    return true;
+}
Index: /branches/eam_branch_20080719/pswarp/src/pswarpTransformTile.c
===================================================================
--- /branches/eam_branch_20080719/pswarp/src/pswarpTransformTile.c	(revision 18666)
+++ /branches/eam_branch_20080719/pswarp/src/pswarpTransformTile.c	(revision 18666)
@@ -0,0 +1,66 @@
+
+// XXX need to determine all of the needed inputs and put them in a structure
+
+typedef struct {
+    int gridX;
+    int gridY;
+} pswarpTransformTileOpts;
+
+pswarpTransformTile (int gridX, int gridY) {
+
+    minX = pswarpMapGridNextGrid_Y (grid, gridY);
+    minY = pswarpMapGridNextGrid_X (grid, gridX);
+
+    maxX = pswarpMapGridNextGrid_Y (grid, gridY);
+    maxY = pswarpMapGridNextGrid_X (grid, gridX);
+
+    map = grid->maps[gridX][gridY];
+
+    psPlane *inPix = psPlaneAlloc();    // Coordinates on the input detector
+
+    // Iterate over the output image pixels (parent frame)
+    long goodPixels = 0;                // Number of input pixels landing on the output image
+    for (int y = minY; y < maxY; y++) {
+        int yOut = y - outRow0;         // Position on image
+        for (int x = minX; x < maxX; x++) {
+
+            // Only transform those pixels requested
+            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 (&inPix->x, &inPix->y, map, x + 0.5, y + 0.5);
+
+            if (inPix->x - inCol0 < 0) continue;
+            if (inPix->x - inCol0 >= inImage->numCols) continue;
+            if (inPix->y - inRow0 < 0) continue;
+            if (inPix->y - inRow0 >= inImage->numRows) continue;
+
+            goodPixels++;
+
+            // XXX include mask
+            // XXX apply scale and offset?
+            // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates
+            double imageValue, varValue;
+            psMaskType maskValue;
+            if (!psImageInterpolate(&imageValue, &varValue, &maskValue,
+                                    inPix->x - inCol0, inPix->y - inRow0, interp)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
+                psFree(interp);
+                psFree(inPix);
+                psFree(grid);
+                return false;
+            }
+            int xOut = x - outCol0;     // Position on image
+            outImageData[yOut][xOut] = imageValue;
+            if (inVar) {
+                outVarData[yOut][xOut] = varValue;
+            }
+            if (outMaskData) {
+                outMaskData[yOut][xOut] = maskValue;
+            }
+        }
+    }
+
+    psFree(inPix);
+}