Index: /branches/eam_branches/20090715/psLib/src/db/psDB.c
===================================================================
--- /branches/eam_branches/20090715/psLib/src/db/psDB.c	(revision 25404)
+++ /branches/eam_branches/20090715/psLib/src/db/psDB.c	(revision 25405)
@@ -2261,4 +2261,5 @@
     PS_DB_OP_LE,
     PS_DB_OP_GE,
+    PS_DB_OP_NE,
 } psDBOpValue;
 
@@ -2308,6 +2309,10 @@
             opStr = "<";
             op = PS_DB_OP_LT;
-        }
-    }
+        } else if (strstr(item->comment, "!=")) {
+            opStr = "!=";
+            op = PS_DB_OP_NE;
+        }
+    }
+
 
     // XXX why are >, < searches not supported here????
@@ -2341,4 +2346,7 @@
         case PS_DB_OP_EQ:
             psStringAppend(&query, "(ABS(%s - %.8f) < %.8f)", itemName, (float)(item->data.F32), PS_DB_FLT_PAD);
+            break;
+        case PS_DB_OP_NE:
+            psStringAppend(&query, "(ABS(%s - %.8f) >= %.8f)", itemName, (float)(item->data.F32), PS_DB_FLT_PAD);
             break;
         case PS_DB_OP_LE:
@@ -2360,4 +2368,7 @@
             psStringAppend(&query, "(ABS(%s - %.17f) < %.17f)", itemName, (float)(item->data.F64), PS_DB_DBL_PAD);
             break;
+        case PS_DB_OP_NE:
+            psStringAppend(&query, "(ABS(%s - %.17f) >= %.17f)", itemName, (float)(item->data.F64), PS_DB_DBL_PAD);
+            break;
         case PS_DB_OP_LE:
         case PS_DB_OP_LT:
@@ -2376,4 +2387,7 @@
         case PS_DB_OP_EQ:
             psStringAppend(&query, "%s = %d", itemName, (int)(item->data.B));
+            break;
+        case PS_DB_OP_NE:
+            psStringAppend(&query, "%s != %d", itemName, (int)(item->data.B));
             break;
         default:
@@ -2397,5 +2411,5 @@
                 psStringAppend(&query, "%s LIKE '%s'", itemName, item->data.str);
             } else {
-                psStringAppend(&query, "%s = '%s'", itemName, item->data.str);
+                psStringAppend(&query, "%s %s '%s'", itemName, opStr, item->data.str);
             }
         }
Index: /branches/eam_branches/20090715/psLib/src/fits/psFits.c
===================================================================
--- /branches/eam_branches/20090715/psLib/src/fits/psFits.c	(revision 25404)
+++ /branches/eam_branches/20090715/psLib/src/fits/psFits.c	(revision 25405)
@@ -150,8 +150,8 @@
             char errorBuf[MAX_STRING_LENGTH], *errorMsg;
 #if ((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE)
-            errorMsg = strerror_r(thisErrno, errorBuf, MAX_STRING_LENGTH);
-#else
             strerror_r(thisErrno, errorBuf, MAX_STRING_LENGTH);
             errorMsg = errorBuf;
+#else
+            errorMsg = strerror_r(thisErrno, errorBuf, MAX_STRING_LENGTH);
 #endif
             psError(PS_ERR_IO, true, "Directory (%s) for requested file is not accessible: %s",
@@ -209,8 +209,8 @@
                 char errorBuf[64], *errorMsg;
 # if ((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE)
-                errorMsg = strerror_r (errno, errorBuf, 64);
-# else
                 strerror_r (errno, errorBuf, 64);
                 errorMsg = errorBuf;
+# else
+                errorMsg = strerror_r (errno, errorBuf, 64);
 # endif
                 psError(PS_ERR_IO, true, "Failed to delete a previously-existing file (%s), error %d: %s",
@@ -344,6 +344,8 @@
 // Therefore, we implement our own version of moving to an extension specified by name.  The pure cfitsio
 // version is used if "conventions.compression" handling is turned off in the psFits structure.
-bool psFitsMoveExtName(const psFits* fits,
-                       const char* extname)
+static bool fitsMoveExtName(const psFits* fits, // FITS file
+                            const char* extname, // Extension name
+                            bool errors // Generate errors?
+    )
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
@@ -356,5 +358,7 @@
         // User wants to use cfitsio.  Good luck to them!
         if (fits_movnam_hdu(fits->fd, ANY_HDU, (char*)extname, 0, &status) != 0) {
-            psFitsError(status, true, _("Could not find HDU '%s'"), extname);
+            if (errors) {
+                psFitsError(status, true, _("Could not find HDU '%s'"), extname);
+            }
             return false;
         }
@@ -378,5 +382,7 @@
         if (fits_movabs_hdu(fits->fd, i, &hdutype, &status)) {
             // We've run off the end
-            psFitsError(status, true, _("Could not find HDU with %s = '%s'"), extword, extname);
+            if (errors) {
+                psFitsError(status, true, _("Could not find HDU with %s = '%s'"), extword, extname);
+            }
             return false;
         }
@@ -406,4 +412,15 @@
     }
     psAbort("Should never reach here.");
+}
+
+
+bool psFitsMoveExtName(const psFits* fits, const char* extname)
+{
+    return fitsMoveExtName(fits, extname, true);
+}
+
+bool psFitsMoveExtNameClean(const psFits* fits, const char* extname)
+{
+    return fitsMoveExtName(fits, extname, false);
 }
 
Index: /branches/eam_branches/20090715/psLib/src/fits/psFits.h
===================================================================
--- /branches/eam_branches/20090715/psLib/src/fits/psFits.h	(revision 25404)
+++ /branches/eam_branches/20090715/psLib/src/fits/psFits.h	(revision 25405)
@@ -245,4 +245,14 @@
 );
 
+/** Moves the FITS HDU to the specified extension name without generating errors.
+ *
+ *  @return bool        TRUE if the extension name was found and move was
+ *                      successful, otherwise FALSE
+ */
+bool psFitsMoveExtNameClean(
+    const psFits* fits,                ///< the psFits object to move
+    const char* extname                ///< the extension name
+);
+
 /** Moves the FITS HDU to the specified extension number
  *
Index: /branches/eam_branches/20090715/psLib/src/fits/psFitsHeader.c
===================================================================
--- /branches/eam_branches/20090715/psLib/src/fits/psFitsHeader.c	(revision 25404)
+++ /branches/eam_branches/20090715/psLib/src/fits/psFitsHeader.c	(revision 25405)
@@ -50,5 +50,5 @@
 
 // List of the start of FITS header keys not to write (handled by cfitsio); NULL-terminated
-static const char *noWriteFitsKeyStarts[] = { "NAXIS", "TTYPE", "TFORM", NULL };
+static const char *noWriteFitsKeyStarts[] = { "NAXIS", "TTYPE", "TFORM", "TZERO", "TSCAL", NULL };
 
 // List of compressed FITS header keys not to write (handled by cfitsio); NULL-terminated
@@ -305,4 +305,5 @@
                 keyType = 'C';
             }
+            psTrace("psLib.fits", 3, "Reading keyword %s, type %c\n", keyName, keyType);
             if (status != 0) {
                 break;
@@ -327,6 +328,12 @@
                                                -INFINITY);
                 } else {
-                    success = psMetadataAddS32(header, PS_LIST_TAIL, keyNameTrans, dupFlag, keyComment,
-                                               atoi(keyValue));
+                    long long value = atoll(keyValue); // Value
+                    if (value > PS_MIN_S32 && value < PS_MAX_S32) {
+                        success = psMetadataAddS32(header, PS_LIST_TAIL, keyNameTrans, dupFlag, keyComment,
+                                                   value);
+                    } else {
+                        success = psMetadataAddS64(header, PS_LIST_TAIL, keyNameTrans, dupFlag, keyComment,
+                                                   value);
+                    }
                 }
                 break;
@@ -509,16 +516,10 @@
         psMetadataItem *simpleItem = psMetadataLookup(output, "SIMPLE"); // SIMPLE in the header
         if (simpleItem) {
-            if (simpleItem->type != PS_DATA_BOOL) {
-                psError(PS_ERR_BAD_PARAMETER_TYPE, true, "SIMPLE in a FITS header must be of boolean type: "
-                        "not %x --- assuming FALSE.\n", simpleItem->type);
+            if (simpleItem->type != PS_DATA_BOOL || !simpleItem->data.B) {
                 int value = false;          // Temporary holder for boolean
+                psWarning("Writing SIMPLE=F to FITS header by request");
                 fits_update_key(fits->fd, TLOGICAL, "SIMPLE", &value,
                                 "File does not conform to FITS standard", &status);
                 simple = false;
-            } else if (!simpleItem->data.B) {
-                simple = false;
-                int value = false;      // Temporary holder for boolean
-                fits_update_key(fits->fd, TLOGICAL, "SIMPLE", &value,
-                                "File does not conform to FITS standard", &status);
             }
             // Uncompressed SIMPLE = T is taken care of by cfitsio.
@@ -546,4 +547,5 @@
             char comment[FLEN_CARD];    // Comment for ZIMAGE; unused
             int value;                  // Value for ZIMAGE; unused
+            psTrace("psLib.fits", 3, "Writing header ZSIMPLE to preserve PHU");
             fits_read_key(fits->fd, TLOGICAL, "ZIMAGE", &value, comment, &status);
             fits_insert_key_log(fits->fd, "ZSIMPLE", simple, "Uncompressed file's conforms to FITS", &status);
@@ -575,8 +577,10 @@
                 if (keywordInList(name, noWriteCompressedKeys) ||
                     (keyStarts && keywordStartsWith(name, noWriteCompressedKeyStarts))) {
+                    psTrace("psLib.fits", 3, "Not writing FITS keyword %s", name);
                     continue;
                 }
             } else if (keywordInList(name, noWriteCompressedKeys) ||
                        (keyStarts && keywordStartsWith(name, noWriteCompressedKeyStarts))) {
+                psTrace("psLib.fits", 3, "Not writing FITS keyword %s", name);
                 continue;
             }
@@ -584,4 +588,5 @@
             if (keywordInList(name, noWriteFitsKeys) ||
                 (keyStarts && keywordStartsWith(name, noWriteFitsKeyStarts))) {
+                psTrace("psLib.fits", 3, "Not writing FITS keyword %s", name);
                 continue;
             }
@@ -592,4 +597,5 @@
                 psWarning("COMMENT header is not of type STRING (%x) --- ignored.", item->type);
             } else {
+                psTrace("psLib.fits", 5, "Writing header COMMENT: %s", item->data.str);
                 fits_write_comment(fits->fd, item->data.str, &status);
             }
@@ -598,4 +604,5 @@
                 psWarning("COMMENT header is not of type STRING (%x) --- ignored.", item->type);
             } else {
+                psTrace("psLib.fits", 5, "Writing header HISTORY: %s", item->data.str);
                 fits_write_history(fits->fd, item->data.str, &status);
             }
@@ -605,28 +612,36 @@
               case PS_DATA_BOOL: {
                   int value = item->data.B;
+                  psTrace("psLib.fits", 5, "Writing BOOL header %s: %d", name, value);
                   fits_update_key(fits->fd, TLOGICAL, name, &value, item->comment, &status);
                   break;
               }
               case PS_DATA_S8:
+                psTrace("psLib.fits", 5, "Writing S8 header %s: %d", name, (int)item->data.S8);
                 fits_update_key(fits->fd, TBYTE, name, &item->data.S8, item->comment, &status);
                 break;
               case PS_DATA_S16:
+                psTrace("psLib.fits", 5, "Writing S16 header %s: %d", name, (int)item->data.S16);
                 fits_update_key(fits->fd, TSHORT, name, &item->data.S16, item->comment, &status);
                 break;
               case PS_DATA_S32:
+                psTrace("psLib.fits", 5, "Writing S32 header %s: %d", name, (int)item->data.S32);
                 fits_update_key(fits->fd, TINT, name, &item->data.S32, item->comment, &status);
                 break;
               case PS_DATA_S64:
+                psTrace("psLib.fits", 5, "Writing S64 header %s: %" PRId64, name, item->data.S64);
                 fits_update_key(fits->fd, TLONGLONG, name, &item->data.S64, item->comment, &status);
                 break;
               case PS_DATA_U8: {
                   unsigned short int temp = item->data.U8;
+                psTrace("psLib.fits", 5, "Writing U8 header %s: %d", name, (int)item->data.U8);
                   fits_update_key(fits->fd, TUSHORT, name, &temp, item->comment, &status);
                   break;
               }
               case PS_DATA_U16:
+                psTrace("psLib.fits", 5, "Writing U16 header %s: %d", name, (int)item->data.U16);
                 fits_update_key(fits->fd, TUSHORT, name, &item->data.U16, item->comment, &status);
                 break;
               case PS_DATA_U32:
+                psTrace("psLib.fits", 5, "Writing U32 header %s: %d", name, (unsigned int)item->data.U32);
                 fits_update_key(fits->fd, TUINT, name, &item->data.U32, item->comment, &status);
                 break;
@@ -639,8 +654,10 @@
                 }
                 psS64 temp = item->data.U64; // Signed version
+                psTrace("psLib.fits", 5, "Writing U64 header %s: %" PRIu64, name, item->data.U64);
                 fits_update_key(fits->fd, TLONGLONG, name, &temp, item->comment, &status);
                 break;
               case PS_DATA_F32: {
                   int infCheck = 0;         // Result of isinf()
+                  psTrace("psLib.fits", 5, "Writing F32 header %s: %f", name, item->data.F32);
                   if (isnan(item->data.F32)) {
                       fits_update_key(fits->fd, TSTRING, name, "NaN", item->comment, &status);
@@ -659,4 +676,5 @@
               case PS_DATA_F64: {
                   int infCheck = 0;         // Result of isinf()
+                  psTrace("psLib.fits", 5, "Writing F32 header %s: %lf", name, item->data.F64);
                   if (isnan(item->data.F64)) {
                       fits_update_key(fits->fd, TSTRING, name, "NaN", item->comment, &status);
@@ -674,4 +692,5 @@
               }
               case PS_DATA_STRING:
+                psTrace("psLib.fits", 5, "Writing STR header %s: %s", name, item->data.str);
                 fits_update_key(fits->fd, TSTRING, name, item->data.V, item->comment, &status);
                 break;
Index: /branches/eam_branches/20090715/psLib/src/imageops/psImageConvolve.c
===================================================================
--- /branches/eam_branches/20090715/psLib/src/imageops/psImageConvolve.c	(revision 25404)
+++ /branches/eam_branches/20090715/psLib/src/imageops/psImageConvolve.c	(revision 25405)
@@ -678,4 +678,137 @@
 }
 
+static bool imageSmoothMaskPixels(psVector *out, const psImage *image, const psImage *mask,
+                                  psImageMaskType maskVal, const psVector *x, const psVector *y,
+                                  const psVector *gaussNorm, float minGauss, int size, int start, int stop)
+{
+    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    int xLast = numCols - 1, yLast = numRows - 1; // Last index
+    for (int i = start; i < stop; i++) {
+        int xPix = x->data.S32[i], yPix = y->data.S32[i]; // Pixel coordinates for smoothing
+
+        int yMin = PS_MAX(yPix - size, 0);
+        int yMax = PS_MIN(yPix + size, yLast);
+        int xMin = PS_MAX(xPix - size, 0);
+        int xMax = PS_MIN(xPix + size, xLast);
+
+        const float *yGauss = &gauss[yMin - yPix];
+
+        double ySumIG = 0.0, ySumG = 0.0;
+        for (int v = yMin; v <= yMax; v++, yGauss++) {
+            const float *xGauss = &gauss[xMin - xPix];
+            double xSumIG = 0.0, xSumG = 0.0;
+            const psImageMaskType *maskData = &mask->data.PS_TYPE_IMAGE_MASK_DATA[v][xMin];
+            const psF32 *imageData = &image->data.F32[v][xMin];
+            for (int u = xMin; u <= xMax; u++, xGauss++, imageData++, maskData++) {
+                if (*maskData & maskVal) {
+                    continue;
+                }
+                xSumIG += *imageData * *xGauss;
+                xSumG += *xGauss;
+            }
+            if (xSumG > minGauss) {
+                ySumIG += xSumIG * *yGauss;
+                ySumG += xSumG * *yGauss;
+            }
+        }
+
+        out->data.F32[i] = ySumG > minGauss ? ySumIG / ySumG : NAN;
+    }
+
+    return true;
+}
+
+static bool psImageSmoothMaskPixelsThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert(job->args, "programming error: no job arguments");
+    psAssert(job->args->n == 11, "programming error: wrong number of job arguments");
+
+    psVector *out = job->args->data[0]; // Output vector
+    const psImage *image  = job->args->data[1]; // Input image
+    const psImage *mask   = job->args->data[2]; // Input mask
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[3], PS_TYPE_IMAGE_MASK_DATA);
+    const psVector *x = job->args->data[4];
+    const psVector *y = job->args->data[5];
+    const psVector *gaussNorm = job->args->data[6];
+    float minGauss = PS_SCALAR_VALUE(job->args->data[7], F32);
+    int size = PS_SCALAR_VALUE(job->args->data[8], S32);
+    int start = PS_SCALAR_VALUE(job->args->data[9], S32);
+    int stop = PS_SCALAR_VALUE(job->args->data[10], S32);
+    return imageSmoothMaskPixels(out, image, mask, maskVal, x, y, gaussNorm,
+                                 minGauss, size, start, stop);
+}
+
+
+psVector *psImageSmoothMaskPixels(const psImage *image, const psImage *mask, psImageMaskType maskVal,
+                                  const psVector *x, const psVector *y,
+                                  float sigma, float numSigma, float minGauss)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+    PS_ASSERT_IMAGE_NON_NULL(mask, NULL);
+    PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_IMAGE_MASK, NULL);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(image, mask, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
+
+    int size = sigma * numSigma + 0.5;  // Half-size of kernel
+
+    int num = x->n;                     // Number of pixels to smooth
+    psVector *out = psVectorAlloc(num, PS_TYPE_F32); // Output results
+
+    // Generate normalized gaussian
+    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
+
+    // Columns
+    if (threaded) {
+        int numThreads = psThreadPoolSize(); // Number of threads
+        int delta = (numThreads) ? num / numThreads + 1 : num; // Block of cols to do at once
+        for (int start = 0; start < num; start += delta) {
+            int stop = PS_MIN(start + delta, num);  // End of block
+
+            psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS");
+            psArrayAdd(job->args, 1, out);
+            psArrayAdd(job->args, 1, (psImage*)image);
+            psArrayAdd(job->args, 1, (psImage*)mask);
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
+            psArrayAdd(job->args, 1, (psVector*)x);
+            psArrayAdd(job->args, 1, (psVector*)y);
+            psArrayAdd(job->args, 1, gaussNorm);
+            PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, size, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, start, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, stop, PS_TYPE_S32);
+            if (!psThreadJobAddPending(job)) {
+                psFree(job);
+                psFree(gaussNorm);
+                psFree(out);
+                return NULL;
+            }
+            psFree(job);
+        }
+    } else if (!imageSmoothMaskPixels(out, image, mask, maskVal, x, y,
+                                      gaussNorm, minGauss, size, 0, num)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to smooth pixels.");
+        psFree(gaussNorm);
+        psFree(out);
+        return NULL;
+    }
+
+    if (threaded && !psThreadPoolWait(true)) {
+        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+        psFree(gaussNorm);
+        psFree(out);
+        return NULL;
+    }
+
+    psFree(gaussNorm);
+    return out;
+}
+
+
 // Smooth an image with masked pixels
 // The calculation and calcMask images are *deliberately* backwards (row,col instead of col,row or
@@ -1325,7 +1458,7 @@
     if (threaded) {
         int numThreads = psThreadPoolSize(); // Number of threads
-        int deltaRows = (numThreads) ? numRows / numThreads : numRows; // Block of rows to do at once
+        int deltaRows = (numThreads) ? numRows / numThreads + 1 : numRows; // Block of rows to do at once
         for (int start = 0; start < numRows; start += deltaRows) {
-            int stop = PS_MIN (start + deltaRows, numRows);  // end of row block
+            int stop = PS_MIN(start + deltaRows, numRows);  // end of row block
 
             psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
@@ -1363,7 +1496,7 @@
     if (threaded) {
         int numThreads = psThreadPoolSize(); // Number of threads
-        int deltaCols = (numThreads) ? numCols / numThreads : numCols; // Block of cols to do at once
+        int deltaCols = (numThreads) ? numCols / numThreads + 1 : numCols; // Block of cols to do at once
         for (int start = 0; start < numCols; start += deltaCols) {
-            int stop = PS_MIN (start + deltaCols, numCols);  // end of col block
+            int stop = PS_MIN(start + deltaCols, numCols);  // end of col block
 
             psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
@@ -1486,4 +1619,10 @@
             psFree(task);
         }
+        {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS", 9);
+            task->function = &psImageSmoothMaskPixelsThread;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
     } else if (!set && threaded) {
         psThreadTaskRemove("PSLIB_IMAGE_CONVOLVE_MASK");
Index: /branches/eam_branches/20090715/psLib/src/imageops/psImageConvolve.h
===================================================================
--- /branches/eam_branches/20090715/psLib/src/imageops/psImageConvolve.h	(revision 25404)
+++ /branches/eam_branches/20090715/psLib/src/imageops/psImageConvolve.h	(revision 25405)
@@ -222,4 +222,18 @@
     );
 
+/// Smooth particular pixels on an image, allowing for masked pixels
+///
+/// Applies a circularly symmetric Gaussian smoothing first in x and then in y
+/// directions with just a vector.  This process is 2N faster than 2D convolutions (in general).
+psVector *psImageSmoothMaskPixels(
+    const psImage *image,               ///< Input image (F32)
+    const psImage *mask,                ///< Mask image
+    psImageMaskType maskVal,            ///< Value to mask
+    const psVector *x,                  ///< x coordinates
+    const psVector *y,                  ///< y coordinates
+    float sigma,                        ///< Width of the smoothing kernel (pixels)
+    float numSigma,                     ///< Size of the smoothing box (sigma)
+    float minGauss                      ///< Minimum fraction of Gaussian to accept
+    );
 
 psImage *psImageSmoothMask_Threaded(psImage *output,
Index: /branches/eam_branches/20090715/psLib/src/sys/psType.h
===================================================================
--- /branches/eam_branches/20090715/psLib/src/sys/psType.h	(revision 25404)
+++ /branches/eam_branches/20090715/psLib/src/sys/psType.h	(revision 25405)
@@ -141,5 +141,5 @@
 } psDataType;
 
-// macros to abstract the generic mask type : these values must be consistent 
+// macros to abstract the generic mask type : these values must be consistent
 #define PS_TYPE_MASK PS_TYPE_U8        /**< the psElemType to use for mask image */
 #define PS_TYPE_MASK_DATA U8           /**< the data member to use for mask image */
@@ -152,8 +152,8 @@
 // alternate versions if needed
 // #define PS_NOT_MASK(A)(UINT16_MAX-(A))
-// #define PS_NOT_MASK(A)(UINT32_MAX-(A)) 
+// #define PS_NOT_MASK(A)(UINT32_MAX-(A))
 // #define PS_NOT_MASK(A)(UINT64_MAX-(A))
 
-// macros to abstract the vector mask type : these values must be consistent 
+// macros to abstract the vector mask type : these values must be consistent
 #define PS_TYPE_VECTOR_MASK PS_TYPE_U8        /**< the psElemType to use for mask image */
 #define PS_TYPE_VECTOR_MASK_DATA U8           /**< the data member to use for mask image */
@@ -161,11 +161,11 @@
 #define PS_MIN_VECTOR_MASK_TYPE 0             /**< minimum valid Vector Mask value */
 #define PS_MAX_VECTOR_MASK_TYPE UINT8_MAX     /**< maximum valid Vector Mask value */
-typedef psU8 psVectorMaskType;			  ///< the C datatype for a mask image
+typedef psU8 psVectorMaskType;                    ///< the C datatype for a mask image
 #define PS_NOT_VECTOR_MASK(A)(UINT8_MAX-(A))
 
-// macros to abstract the image mask type : these values must be consistent 
-#define PS_TYPE_IMAGE_MASK PS_TYPE_U16	     /**< the psElemType to use for mask image */
-#define PS_TYPE_IMAGE_MASK_DATA U16	     /**< the data member to use for mask image */
-#define PS_TYPE_IMAGE_MASK_NAME "psU16"	     /**< the data type for mask as a string */
+// macros to abstract the image mask type : these values must be consistent
+#define PS_TYPE_IMAGE_MASK PS_TYPE_U16       /**< the psElemType to use for mask image */
+#define PS_TYPE_IMAGE_MASK_DATA U16          /**< the data member to use for mask image */
+#define PS_TYPE_IMAGE_MASK_NAME "psU16"      /**< the data type for mask as a string */
 #define PS_MIN_IMAGE_MASK_TYPE 0             /**< minimum valid Image Mask value */
 #define PS_MAX_IMAGE_MASK_TYPE UINT16_MAX    /**< maximum valid Image Mask value */
@@ -246,16 +246,4 @@
 };
 
-/// Macro to get the bad pixel reason code (stored as part of mask value)
-#define PS_BADPIXEL_BITMASK 0x0f
-#define PS_GET_BADPIXEL(maskValue) (maskValue & PS_BADPIXEL_BITMASK)
-
-#define PS_IS_BADPIXEL(maskValue) (PS_GET_BADPIXEL(maskValue) != 0)
-
-/// Macro to apply a bad pixel reason code to mask image
-#define PS_SET_BADPIXEL(maskValue, reasonCode) \
-{ \
-    maskValue = (psMaskType)((reasonCode & PS_BADPIXEL_BITMASK) | (maskValue & ~PS_BADPIXEL_BITMASK)); \
-}
-
 /// Macro to determine if the psElemType is an integer.
 #define PS_IS_PSELEMTYPE_INT(x) ((x & 0x100) == 0x100)
Index: /branches/eam_branches/20090715/psLib/src/types/psMetadata.c
===================================================================
--- /branches/eam_branches/20090715/psLib/src/types/psMetadata.c	(revision 25404)
+++ /branches/eam_branches/20090715/psLib/src/types/psMetadata.c	(revision 25405)
@@ -996,5 +996,4 @@
     if (metadataItem->type == PS_DATA_METADATA_MULTI) {
         // if multiple keys found, use the first.
-        //        metadataItem = (psMetadataItem*)((metadataItem->data.list)->head);
         metadataItem = (psMetadataItem*)(metadataItem->data.list->head->data);
         if (status) {
@@ -1082,5 +1081,4 @@
     } \
     \
-    /* psFree(metadataItem); currently, the lookup doesn't increment the ref count */ \
     return value; \
 }
@@ -1101,4 +1099,44 @@
 psMetadataLookupNumTYPE(VectorMaskType,VectorMask)
 psMetadataLookupNumTYPE(ImageMaskType,ImageMask)
+
+#define psMetadataLookupPtrTYPE(TYPENAME,NAME,TYPE,VAL) \
+TYPENAME psMetadataLookup##NAME(bool *status, const psMetadata *md, const char *key) \
+{ \
+    PS_ASSERT_METADATA_NON_NULL(md, NULL); \
+    PS_ASSERT_STRING_NON_EMPTY(key, NULL); \
+    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); /* The item of interest */ \
+    if (!item) { \
+        if (status) { \
+            *status = false; \
+        } else { \
+            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key); \
+        } \
+        return NULL; \
+    } \
+    \
+    if (item->type == PS_DATA_METADATA_MULTI) { \
+        /* if multiple keys found, use the first. */ \
+        item = item->data.list->head->data; \
+    } \
+    if (item->type != TYPE) { \
+        if (status) { \
+            *status = false; \
+        } else { \
+            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_META, as expected.\n", key); \
+        } \
+        return NULL; \
+    } \
+    \
+    if (status) { \
+        *status = true; \
+    } \
+    return item->data.VAL; \
+}
+
+psMetadataLookupPtrTYPE(psMetadata*, Metadata, PS_DATA_METADATA, md)
+psMetadataLookupPtrTYPE(psString, Str, PS_DATA_STRING, str)
+psMetadataLookupPtrTYPE(psTime*, Time, PS_DATA_TIME, V)
+psMetadataLookupPtrTYPE(psVector*, Vector, PS_DATA_VECTOR, V)
+
 
 psMetadataItem* psMetadataGet(const psMetadata *md,
@@ -1257,105 +1295,7 @@
 }
 
-psMetadata *psMetadataLookupMetadata(bool *status,
-                                     const psMetadata *md,
-                                     const char *key)
-{
-    PS_ASSERT_METADATA_NON_NULL(md,NULL);
-    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); // The metadata with instruments
-    psMetadata *value = NULL;  // The value to return
-    if (!item) {
-        // The given key isn't in the metadata
-        if (status) {
-            *status = false;
-        } else {
-            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key);
-        }
-    } else if (item->type != PS_DATA_METADATA) {
-        // The value at the key isn't metadata
-        if (status) {
-            *status = false;
-        } else {
-            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_META, as expected.\n", key);
-        }
-        //        value = NULL;
-    } else {
-        // We have the requested metadata
-        if (status) {
-            *status = true;
-        }
-        value = item->data.md; // The requested metadata
-    }
-    return value;
-}
-
-psTime *psMetadataLookupTime(bool *status,
-                             const psMetadata *md,
-                             const char *key)
-{
-    PS_ASSERT_METADATA_NON_NULL(md,NULL);
-
-    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key);
-    if (!item) {
-        // The given key isn't in the metadata
-        if (status) {
-            *status = false;
-        } else {
-            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key);
-        }
-        return NULL;
-    }
-
-    if (item->type != PS_DATA_TIME) {
-        // The value at the key isn't metadata
-        if (status) {
-            *status = false;
-        } else {
-            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_TIME, as expected.\n", key);
-        }
-        return NULL;
-    }
-
-    // We have the requested metadata
-    if (status) {
-        *status = true;
-    }
-
-    return item->data.V;
-}
-
-
-psString psMetadataLookupStr(bool *status,
-                             const psMetadata *md,
-                             const char *key)
-{
-    PS_ASSERT_METADATA_NON_NULL(md,NULL);
-
-    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); // The metadata with instruments
-    //    char *value = NULL;   // The value to return
-    psString value = NULL;
-    if (!item) {
-        // The given key isn't in the metadata
-        if (status) {
-            *status = false;
-        } else {
-            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key);
-        }
-    } else if (item->type != PS_DATA_STRING) {
-        // The value at the key isn't of the desired type
-        if (status) {
-            *status = false;
-        } else {
-            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_STRING, as expected.\n", key);
-        }
-        //        value = NULL;
-    } else {
-        // We have the requested metadata
-        if (status) {
-            *status = true;
-        }
-        value = item->data.V;
-    }
-    return value;
-}
+
+
+
 
 psList *psMetadataKeys(psMetadata *md)
Index: /branches/eam_branches/20090715/psLib/src/types/psMetadata.h
===================================================================
--- /branches/eam_branches/20090715/psLib/src/types/psMetadata.h	(revision 25404)
+++ /branches/eam_branches/20090715/psLib/src/types/psMetadata.h	(revision 25405)
@@ -717,4 +717,5 @@
 PS_METADATA_LOOKUP_TYPE_DECL(Ptr, psPtr);
 PS_METADATA_LOOKUP_TYPE_DECL(Str, psString);
+PS_METADATA_LOOKUP_TYPE_DECL(Vector, psVector*);
 PS_METADATA_LOOKUP_TYPE_DECL(Metadata, psMetadata*);
 PS_METADATA_LOOKUP_TYPE_DECL(Time, psTime*);
Index: /branches/eam_branches/20090715/psLib/src/types/psTree.c
===================================================================
--- /branches/eam_branches/20090715/psLib/src/types/psTree.c	(revision 25404)
+++ /branches/eam_branches/20090715/psLib/src/types/psTree.c	(revision 25405)
@@ -14,4 +14,6 @@
 #include "psTree.h"
 
+//#define INPUT_CHECK                   // Check inputs for functions that may be in a tight loop?
+
 
 // XXX Upgrades:
@@ -84,5 +86,5 @@
 }
 
-psTree *psTreeAlloc(int dim, int maxLeafContents, long numNodes)
+psTree *psTreeAlloc(int dim, int maxLeafContents, psTreeType type, long numNodes)
 {
     psAssert(dim > 0, "Dimensionality (%d) must be positive", dim);
@@ -95,4 +97,5 @@
     tree->dim = dim;
     tree->maxLeafContents = maxLeafContents;
+    tree->type = type;
 
     tree->numNodes = numNodes;
@@ -115,5 +118,5 @@
 
 
-psTree *psTreePlant(int dim, int maxLeafContents, ...)
+psTree *psTreePlant(int dim, int maxLeafContents, psTreeType type, ...)
 {
     PS_ASSERT_INT_POSITIVE(dim, NULL);
@@ -122,5 +125,5 @@
     // Parse coordinate list
     va_list args;                       // Variable argument list
-    va_start(args, maxLeafContents);
+    va_start(args, type);
     psArray *coords = psArrayAlloc(dim); // Array of coordinates
     long numData = 0;                   // Number of data points
@@ -145,8 +148,8 @@
         long pow2;                      // Smallest power of two >= numData
         for (pow2 = 1; pow2 < numData; pow2 <<= 1);
-        numNodes = PS_MIN(pow2 - 1, 2 * numData - pow2 / 2 - 1);
-    }
-
-    psTree *tree = psTreeAlloc(dim, maxLeafContents, numNodes);
+        numNodes = PS_MAX(PS_MIN(pow2 - 1, 2 * numData - pow2 / 2 - 1), 1);
+    }
+
+    psTree *tree = psTreeAlloc(dim, maxLeafContents, type, numNodes);
     tree->data = psTreeCoordArrayAlloc(numData, dim);
     tree->numData = numData;
@@ -199,4 +202,9 @@
     }
 
+    if (numData <= maxLeafContents) {
+        // Don't need to do any more work
+        return tree;
+    }
+
     psArray *work = psArrayAlloc(numNodes); // Work queue
     work->data[0] = root;
@@ -365,5 +373,5 @@
 psTreeNode *psTreeLeaf(const psTree *tree, const psVector *coords)
 {
-#if 1 // Might be in a tight loop
+#ifdef INPUT_CHECK // Might be in a tight loop
     PS_ASSERT_TREE_NON_NULL(tree, NULL);
     PS_ASSERT_VECTOR_NON_NULL(coords, NULL);
@@ -403,20 +411,41 @@
 {
     int dim = tree->dim;                // Dimensionality
-    switch (dim) {
-      case 2:
-        return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) +
-            PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]);
-      case 3:
-        return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) +
-            PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]) +
-            PS_SQR(coords->data.F64[2] - tree->data->F64[index][2]);
-      default: {
-          double distance2 = 0.0;             // Distance of interest
-          for (int i = 0; i < dim; i++) {
-              distance2 += PS_SQR(coords->data.F64[i] - tree->data->F64[index][i]);
+    switch (tree->type) {
+      case PS_TREE_EUCLIDEAN:
+        switch (dim) {
+          case 2:
+            return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) +
+                PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]);
+          case 3:
+            return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) +
+                PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]) +
+                PS_SQR(coords->data.F64[2] - tree->data->F64[index][2]);
+          default: {
+              double distance2 = 0.0;             // Distance of interest
+              for (int i = 0; i < dim; i++) {
+                  distance2 += PS_SQR(coords->data.F64[i] - tree->data->F64[index][i]);
+              }
+              return distance2;
           }
-          return distance2;
-      }
-    }
+        }
+        break;
+      case PS_TREE_SPHERICAL:
+        switch (dim) {
+          case 2: {
+              // Haversine formula
+              double dphi = coords->data.F64[1] - tree->data->F64[index][1];
+              double sindphi = sin(dphi / 2.0);
+              double dlambda = coords->data.F64[0] - tree->data->F64[index][0];
+              double sindlambda = sin(dlambda / 2.0);
+              return PS_SQR(sindphi) +
+                  cos(coords->data.F64[1]) * cos(tree->data->F64[index][1]) * PS_SQR(sindlambda);
+          }
+          default:
+            psAbort("Spherical distances not supported for more than 2 dimensions");
+        }
+      default:
+        psAbort("Unrecognised type: %x", tree->type);
+    }
+
     return NAN;
 }
@@ -430,10 +459,32 @@
         double minDiff = tree->min->F64[index][i] - coords->data.F64[i];
         if (minDiff > 0) {
-            distance += PS_SQR(minDiff);
+            switch (tree->type) {
+              case PS_TREE_EUCLIDEAN:
+                distance += PS_SQR(minDiff);
+                break;
+              case PS_TREE_SPHERICAL: {
+                  double sinDiff = sin(minDiff / 2.0);
+                  distance += PS_SQR(sinDiff);
+                  break;
+              }
+              default:
+                psAbort("Unrecognised type: %x", tree->type);
+            }
             continue;
         }
         double maxDiff = coords->data.F64[i] - tree->max->F64[index][i];
         if (maxDiff > 0) {
-            distance += PS_SQR(maxDiff);
+            switch (tree->type) {
+              case PS_TREE_EUCLIDEAN:
+                distance += PS_SQR(maxDiff);
+                break;
+              case PS_TREE_SPHERICAL: {
+                  double sinDiff = sin(maxDiff / 2.0);
+                  distance += PS_SQR(sinDiff);
+                  break;
+              }
+              default:
+                psAbort("Unrecognised type: %x", tree->type);
+            }
             continue;
         }
@@ -479,12 +530,12 @@
 }
 
-// Return the index of the nearest neighbour to given coordinates, within some radius
+// Return the index of the nearest neighbour to given coordinates, within some distance measure
 // This is the engine for psTreeNearest() and psTreeNearestWithin()
 static inline long treeNearestWithin(const psTree *tree, // Tree
                                      const psVector *coordinates, // Coordinates of interest
-                                     double bestDistance // Distance (radius-squared) to best point
+                                     double bestDistance // Distance measure to best point
     )
 {
-#if 1 // Might be in a tight loop
+#ifdef INPUT_CHECK // Might be in a tight loop
     PS_ASSERT_TREE_NON_NULL(tree, -1);
     PS_ASSERT_VECTOR_NON_NULL(coordinates, -1);
@@ -545,12 +596,31 @@
 
 
+// Convert a radius to our internal "distance measure"
+// Often, we're given a search radius, but for efficiency reasons, we don't use that internally.
+static double treeRadiusToDistance(const psTree *tree, double radius)
+{
+    switch (tree->type) {
+      case PS_TREE_EUCLIDEAN:
+        // Using the square of the distance as the distance measure
+        return PS_SQR(radius);
+      case PS_TREE_SPHERICAL: {
+          // Using a rearrangement of the Haversine formula
+          double sindist = sin(radius / 2.0);
+          return PS_SQR(sindist);
+      }
+      default:
+        psAbort("Unrecognised type: %x", tree->type);
+    }
+}
+
+
 long psTreeNearestWithin(const psTree *tree, const psVector *coords, double radius)
 {
-    return treeNearestWithin(tree, coords, PS_SQR(radius));
-}
-
-
-// Search a leaf node for points within radius squared
-static inline long treeLeafSearchWithin(double radius2, // Radius squared to search
+    return treeNearestWithin(tree, coords, treeRadiusToDistance(tree, radius));
+}
+
+
+// Search a leaf node for points within distance
+static inline long treeLeafSearchWithin(double distance, // Distance to search
                                         const psTree *tree, // Tree of interest
                                         const psTreeNode *leaf, // Leaf to search
@@ -561,6 +631,5 @@
     for (int i = 0; i < leaf->num; i++) {
         long index = leaf->contents[i]; // Index of point
-        double distance = treeContentDistance(tree, index, coords); // Distance to point
-        if (distance < radius2) {
+        if (treeContentDistance(tree, index, coords) < distance) {
             num++;
         }
@@ -572,5 +641,5 @@
 long psTreeWithin(const psTree *tree, const psVector *coordinates, double radius)
 {
-#if 1 // Might be in a tight loop
+#ifdef INPUT_CHECK // Might be in a tight loop
     PS_ASSERT_TREE_NON_NULL(tree, -1);
     PS_ASSERT_VECTOR_NON_NULL(coordinates, -1);
@@ -581,5 +650,5 @@
                         psVectorCopy(NULL, coordinates, PS_TYPE_F64)); // F64 version of coordinates
 
-    radius *= radius;                   // We work with the radius-squared
+    double distance = treeRadiusToDistance(tree, radius); // Distance measure
     long num = 0;                       // Number of points in circle
 
@@ -588,5 +657,5 @@
     // Find the closest point in the leaf that contains the point of interest
     psTreeNode *leaf = psTreeLeaf(tree, coords); // Leaf containing the point of interest
-    num += treeLeafSearchWithin(radius, tree, leaf, coords);
+    num += treeLeafSearchWithin(distance, tree, leaf, coords);
 
     psArray *work = psArrayAlloc(tree->numNodes); // Work queue
@@ -605,5 +674,5 @@
             }
             // Leaf node
-            num += treeLeafSearchWithin(radius, tree, node, coords);
+            num += treeLeafSearchWithin(distance, tree, node, coords);
             work->data[workIndex] = NULL;
             workIndex--;
@@ -618,28 +687,28 @@
 }
 
-// Search a leaf node for any points within radius squared
-static inline bool treeLeafSearchWithinAny(double radius2, // Radius squared to search
-                                           const psTree *tree, // Tree of interest
-                                           const psTreeNode *leaf, // Leaf to search
-                                           const psVector *coords // Coordinates of interest
+// Search a leaf node for points within distance
+static inline void treeLeafSearchAllWithin(psVector *result,       // Result vector
+                                          double distance, // Distance to search
+                                          const psTree *tree, // Tree of interest
+                                          const psTreeNode *leaf, // Leaf to search
+                                          const psVector *coords // Coordinates of interest
     )
 {
     for (int i = 0; i < leaf->num; i++) {
         long index = leaf->contents[i]; // Index of point
-        double distance = treeContentDistance(tree, index, coords); // Distance to point
-        if (distance < radius2) {
-            return true;
-        }
-    }
-    return false;
-}
-
-// Given an arbitrary point and a matching radius, return whether there are any points within that radius
-bool psTreeWithinAny(const psTree *tree, const psVector *coordinates, double radius)
-{
-#if 1 // Might be in a tight loop
-    PS_ASSERT_TREE_NON_NULL(tree, false);
-    PS_ASSERT_VECTOR_NON_NULL(coordinates, false);
-    PS_ASSERT_VECTOR_SIZE(coordinates, (long)tree->dim, false);
+        if (treeContentDistance(tree, index, coords) < distance) {
+            psVectorAppend(result, index);
+        }
+    }
+    return;
+}
+
+// Given an arbitrary point and a matching radius, return the index of all points within that radius
+psVector *psTreeAllWithin(const psTree *tree, const psVector *coordinates, double radius)
+{
+#ifdef INPUT_CHECK // Might be in a tight loop
+    PS_ASSERT_TREE_NON_NULL(tree, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(coordinates, NULL);
+    PS_ASSERT_VECTOR_SIZE(coordinates, (long)tree->dim, NULL);
 #endif
 
@@ -647,14 +716,13 @@
                         psVectorCopy(NULL, coordinates, PS_TYPE_F64)); // F64 version of coordinates
 
-    radius *= radius;                   // We work with the radius-squared
-
-    // This is essentially the same as psTreeWithin, except we can bail as soon as we find something
+    double distance = treeRadiusToDistance(tree, radius); // Distance measure
+
+    psVector *result = psVectorAllocEmpty(4, PS_TYPE_S64); // Indices of points within match radius
+
+    // This is essentially the same as psTreeNearest, except we're not allowed to prune
 
     // Find the closest point in the leaf that contains the point of interest
     psTreeNode *leaf = psTreeLeaf(tree, coords); // Leaf containing the point of interest
-    if (treeLeafSearchWithinAny(radius, tree, leaf, coords)) {
-        psFree(coords);
-        return true;
-    }
+    treeLeafSearchAllWithin(result, distance, tree, leaf, coords);
 
     psArray *work = psArrayAlloc(tree->numNodes); // Work queue
@@ -673,5 +741,72 @@
             }
             // Leaf node
-            if (treeLeafSearchWithinAny(radius, tree, node, coords)) {
+            treeLeafSearchAllWithin(result, distance, tree, node, coords);
+            work->data[workIndex] = NULL;
+            workIndex--;
+        }
+
+        leaf = up;
+    }
+    psFree(work);
+    psFree(coords);
+
+    return result;
+}
+
+// Search a leaf node for any points within distance measure
+static inline bool treeLeafSearchWithinAny(double distance, // Distance to search
+                                           const psTree *tree, // Tree of interest
+                                           const psTreeNode *leaf, // Leaf to search
+                                           const psVector *coords // Coordinates of interest
+    )
+{
+    for (int i = 0; i < leaf->num; i++) {
+        long index = leaf->contents[i]; // Index of point
+        if (treeContentDistance(tree, index, coords) < distance) {
+            return true;
+        }
+    }
+    return false;
+}
+
+// Given an arbitrary point and a matching radius, return whether there are any points within that radius
+bool psTreeWithinAny(const psTree *tree, const psVector *coordinates, double radius)
+{
+#ifdef INPUT_CHECK // Might be in a tight loop
+    PS_ASSERT_TREE_NON_NULL(tree, false);
+    PS_ASSERT_VECTOR_NON_NULL(coordinates, false);
+    PS_ASSERT_VECTOR_SIZE(coordinates, (long)tree->dim, false);
+#endif
+
+    psVector *coords = (coordinates->type.type == PS_TYPE_F64 ? psMemIncrRefCounter((psVector*)coordinates) :
+                        psVectorCopy(NULL, coordinates, PS_TYPE_F64)); // F64 version of coordinates
+
+    double distance = treeRadiusToDistance(tree, radius); // Distance measure
+
+    // This is essentially the same as psTreeWithin, except we can bail as soon as we find something
+
+    // Find the closest point in the leaf that contains the point of interest
+    psTreeNode *leaf = psTreeLeaf(tree, coords); // Leaf containing the point of interest
+    if (treeLeafSearchWithinAny(distance, tree, leaf, coords)) {
+        psFree(coords);
+        return true;
+    }
+
+    psArray *work = psArrayAlloc(tree->numNodes); // Work queue
+    while (leaf->up) {
+        psTreeNode *up = leaf->up;      // Parent node
+
+        long workIndex = 0;
+        work->data[workIndex] = (leaf == up->left ? up->right : up->left); // The other node
+        while (workIndex >= 0) {
+            psTreeNode *node = work->data[workIndex];
+            if (node->left) {
+                // Branch node
+                work->data[workIndex] = node->right;
+                work->data[++workIndex] = node->left;
+                continue;
+            }
+            // Leaf node
+            if (treeLeafSearchWithinAny(distance, tree, node, coords)) {
                 // Clear out the work queue
                 memset(work->data, 0, (workIndex + 1) * sizeof(void));
@@ -695,5 +830,5 @@
 psVector *psTreeCoords(psVector *out, const psTree *tree, long index)
 {
-#if 1 // Might be in a tight loop
+#ifdef INPUT_CHECK // Might be in a tight loop
     PS_ASSERT_TREE_NON_NULL(tree, NULL);
     PS_ASSERT_INT_NONNEGATIVE(index, NULL);
Index: /branches/eam_branches/20090715/psLib/src/types/psTree.h
===================================================================
--- /branches/eam_branches/20090715/psLib/src/types/psTree.h	(revision 25404)
+++ /branches/eam_branches/20090715/psLib/src/types/psTree.h	(revision 25405)
@@ -16,4 +16,12 @@
 } psTreeCoordArray;
 
+/// Type of tree
+///
+/// Specifies how distances are measured
+typedef enum {
+    PS_TREE_EUCLIDEAN,                  // d^2 = dx^2 + dy^2 + ...
+    PS_TREE_SPHERICAL,                  // sin(dist/2)^2 = sin(dphi/2)^2 + cos(phi1)cos(phi2)sin(dlambda/2)^2
+} psTreeType;
+
 /// A simple kd-tree implementation
 ///
@@ -23,4 +31,5 @@
     int dim;                            ///< Dimensionality
     int maxLeafContents;                ///< Maximum number of points on a leaf
+    psTreeType type;                    ///< Type of tree
     long numNodes;                      ///< Number of nodes
     long numData;                       ///< Number of data points
@@ -67,4 +76,5 @@
 psTree *psTreeAlloc(int dim,            ///< Dimensionality
                     int maxLeafContents,///< Maximum number of points on a leaf
+                    psTreeType type,    ///< Type of tree
                     long numNodes       ///< Number of nodes in tree
                     );
@@ -75,4 +85,5 @@
 psTree *psTreePlant(int dim,            ///< Dimensionality
                     int maxLeafContents,///< Maximum number of points on a leaf
+                    psTreeType type,    ///< Type of tree
                     ...                 ///< psVector for each coordinate
                     );
@@ -108,4 +119,10 @@
                      );
 
+/// Return the index of all points within some radius of given coordinates
+psVector *psTreeAllWithin(const psTree *tree,          ///< Tree
+                          const psVector *coordinates, ///< Coordinates of interest
+                          double radius                ///< Radius of interest
+                          );
+
 /// Return the coordinates of a point in the tree, specified by its index
 psVector *psTreeCoords(psVector *out,   ///< Output vector, or NULL
Index: /branches/eam_branches/20090715/psLib/test/types/tap_psTree.c
===================================================================
--- /branches/eam_branches/20090715/psLib/test/types/tap_psTree.c	(revision 25404)
+++ /branches/eam_branches/20090715/psLib/test/types/tap_psTree.c	(revision 25405)
@@ -3,11 +3,13 @@
 #include "pstap.h"
 
-#define NUM 10000                       // Number of points
+#define NUM 1000000                      // Number of points
+#define SPHERICAL_DISTANCE 3.0          // Distance of interest for spherical test
 
 int main(int argc, char *argv[])
 {
     psLibInit(NULL);
-    plan_tests(6);
+    plan_tests(13);
 
+    // Euclidean geometry: 6 tests
     {
         psMemId id = psMemGetId();
@@ -23,5 +25,5 @@
         psFree(rng);
 
-        psTree *tree = psTreePlant(2, 2, x, y);
+        psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y);
 
         ok(tree, "Tree planted");
@@ -35,5 +37,5 @@
             long closeIndex = psTreeNearest(tree, coords);
             psFree(coords);
-            ok(closeIndex >= 0 && closeIndex < tree->numNodes, "found point: %ld", closeIndex);
+            ok(closeIndex >= 0 && closeIndex < tree->numNodes, "found closest point: %ld", closeIndex);
 
             long bestIndex = -1;
@@ -68,4 +70,89 @@
     }
 
+    // Spherical geometry: 7 tests
+    {
+        psMemId id = psMemGetId();
+
+        psVector *ra = psVectorAlloc(NUM, PS_TYPE_F64);
+        psVector *dec = psVectorAlloc(NUM, PS_TYPE_F64);
+
+        psRandom *rng = psRandomAllocSpecific(PS_RANDOM_TAUS, 0);
+        for (int i = 0; i < NUM; i++) {
+            // Using http://mathworld.wolfram.com/SpherePointPicking.html
+            ra->data.F64[i] = psRandomUniform(rng);
+            dec->data.F64[i] = acos(2.0 * psRandomUniform(rng) - 1.0) - M_PI_2;
+        }
+
+        psTree *tree = psTreePlant(2, 2, PS_TREE_SPHERICAL, ra, dec);
+
+        ok(tree, "Tree planted");
+        skip_start(!tree, 4, "tree died");
+        {
+            //            psTreePrint(stderr, tree);
+
+            psVector *coords = psVectorAlloc(2, PS_TYPE_F64);
+            coords->data.F64[0] = psRandomUniform(rng);
+            coords->data.F64[1] = acos(2.0 * psRandomUniform(rng) - 1.0) - M_PI_2;
+
+            psVector *indices = psTreeAllWithin(tree, coords, DEG_TO_RAD(SPHERICAL_DISTANCE));
+            ok(indices && indices->type.type == PS_TYPE_S64, "got list of indices (%ld points)", indices->n);
+            long closeIndex = psTreeNearest(tree, coords);
+            ok(closeIndex >= 0 && closeIndex < tree->numNodes, "found closest point: %ld", closeIndex);
+
+            ok(psVectorSortInPlace(indices), "sorted indices");
+
+            bool allgood = true;        // All points in the appropriate place?
+            double bestDistance = INFINITY; // Distance to best point
+            long bestIndex = -1;        // Index of best point
+            for (long i = 0, j = 0; i < NUM; i++) {
+#if 1
+                // Traditional formula
+                double dist = acos(sin(coords->data.F64[1]) * sin(dec->data.F64[i]) +
+                                   cos(coords->data.F64[1]) * cos(dec->data.F64[i]) *
+                                   cos(coords->data.F64[0] - ra->data.F64[i]));
+#else
+                // Haversine formula: used in psTree
+                double dphi = coords->data.F64[1] - dec->data.F64[i];
+                double sindphi = sin(dphi/2.0);
+                double dlambda = coords->data.F64[0] - ra->data.F64[i];
+                double sindlambda = sin(dlambda/2.0);
+                double dist = 2.0 * asin(sqrt(PS_SQR(sindphi) +
+                                              cos(coords->data.F64[1]) * cos(dec->data.F64[i]) *
+                                              PS_SQR(sindlambda)));
+#endif
+                                              if (dist < bestDistance) {
+                    bestDistance = dist;
+                    bestIndex = i;
+                }
+                if (i == indices->data.S64[j]) {
+                    j++;
+                    if (dist > DEG_TO_RAD(SPHERICAL_DISTANCE)) {
+                        diag("%ld is in the list, but shouldn't be (%lf vs %lf)",
+                             i, RAD_TO_DEG(dist), SPHERICAL_DISTANCE);
+                        allgood = false;
+                    }
+                } else if (dist <= DEG_TO_RAD(SPHERICAL_DISTANCE)) {
+                    diag("%ld is not in the list, but should be (%lf vs %lf)",
+                         i, RAD_TO_DEG(dist), SPHERICAL_DISTANCE);
+                    allgood = false;
+                }
+            }
+            ok(allgood, "list is acurate");
+            ok(bestIndex == closeIndex, "correct point: %ld vs %ld", closeIndex, bestIndex);
+
+            psFree(coords);
+            psFree(indices);
+
+        }
+        skip_end();
+
+        psFree(rng);
+        psFree(tree);
+        psFree(ra);
+        psFree(dec);
+
+        ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks");
+    }
+
     psLibFinalize();
 }