Changeset 28225

Timestamp:

Jun 4, 2010, 2:23:58 PM (16 years ago)

Author:

Paul Price

Message:

Fix multiple problems:

1. Some functions were not pruning, hence were searching EVERY point (exactly what we're trying to avoid).
2. Distance to branch was not considering the ra=0 boundary.
3. Was passing psVectorAppend a long that it was reading as psS64, which is wrong when sizeof(long) != sizeof(psS64), e.g., on 32-bit machines.

File:

• trunk/psLib/src/types/psTree.c (modified) (10 diffs)

trunk/psLib/src/types/psTree.c

@@ -342 +342 @@
         fprintf(fptr, " ");
     }
-    fprintf(fptr, "(");
     for (int i = 0; i < node->num; i++) {
+        fprintf(fptr, "%ld=(",node->contents[i]);
         psF64 *coords = tree->data->F64[node->contents[i]]; // Coordinates
         for (int j = 0; j < dim; j++) {
@@ -351 +351 @@
             }
         }
-        fprintf(fptr, ")");
-        if (i < node->num - 1) {
-            fprintf(fptr, " (");
-        }
+        fprintf(fptr, ") ");
     }
     fprintf(fptr, "\n");
@@ -451 +448 @@
 }
 
+static inline double sin2diff(double coord1, double coord2)
+{
+    double sinDiff = sin(coord1 - coord2);
+    return PS_SQR(sinDiff);
+}
+
 // Returns the square of the distance between a leaf and an arbitrary point
 static inline double treeBranchDistance(const psTree *tree, long index, const psVector *coords)
@@ -456 +459 @@
     int dim = tree->dim;                // Dimensionality
     double distance = 0.0;              // Distance to box
-    for (int i = 0; i < dim; i++) {
-        double minDiff = tree->min->F64[index][i] - coords->data.F64[i];
-        if (minDiff > 0) {
-            switch (tree->type) {
-              case PS_TREE_EUCLIDEAN:
+    switch (tree->type) {
+      case PS_TREE_EUCLIDEAN:
+        for (int i = 0; i < dim; i++) {
+            double minDiff = tree->min->F64[index][i] - coords->data.F64[i];
+            if (minDiff > 0) {
                 distance += PS_SQR(minDiff);
-                break;
-              case PS_TREE_SPHERICAL: {
-                  double sinDiff = sin(minDiff / 2.0);
-                  distance += PS_SQR(sinDiff);
-                  break;
+                continue;
+            }
+            double maxDiff = coords->data.F64[i] - tree->max->F64[index][i];
+            if (maxDiff > 0) {
+                distance += PS_SQR(maxDiff);
+                continue;
+            }
+        }
+        break;
+      case PS_TREE_SPHERICAL: {
+          double ra = coords->data.F64[0], dec = coords->data.F64[1];                 // Coords of interest
+          double raMin = tree->min->F64[index][0], raMax = tree->max->F64[index][0]; // RA bounds
+          double decMin = tree->min->F64[index][1], decMax = tree->max->F64[index][1]; // Dec bounds
+
+          double raDist = 0.0;
+          if (ra < raMin || ra > raMax) {
+              double ra1 = sin2diff(ra, raMin);
+              double ra2 = sin2diff(ra + 2.0 * M_PI, raMin);
+              raDist = PS_MIN(ra1, ra2);
+              double ra3 = sin2diff(ra, raMax);
+              if (ra3 < raDist) {
+                  raDist = ra3;
               }
-              default:
-                psAbort("Unrecognised type: %x", tree->type);
-            }
-            continue;
-        }
-        double maxDiff = coords->data.F64[i] - tree->max->F64[index][i];
-        if (maxDiff > 0) {
-            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;
+              double ra4 = sin2diff(ra + 2.0 * M_PI, raMax);
+              if (ra4 < raDist) {
+                  raDist = ra4;
               }
-              default:
-                psAbort("Unrecognised type: %x", tree->type);
-            }
-            continue;
-        }
+          }
+
+          double decDist = 0.0;
+          if (dec < decMin || dec > decMax) {
+              double dec1 = sin2diff(dec, decMin);
+              double dec2 = sin2diff(dec, decMax);
+              decDist = PS_MIN(dec1, dec2);
+          }
+          break;
+      }
+      default:
+        psAbort("Unrecognised type: %x", tree->type);
     }
     return distance;
@@ -653 +668 @@
     long num = 0;                       // Number of points in circle
 
-    // This is essentially the same as psTreeNearest, except we're not allowed to prune
+    // This is essentially the same as psTreeNearest
 
     // Find the closest point in the leaf that contains the point of interest
@@ -667 +682 @@
         while (workIndex >= 0) {
             psTreeNode *node = work->data[workIndex];
+            if (treeBranchDistance(tree, node->index, coords) > distance) {
+                // No need to investigate
+                work->data[workIndex] = NULL;
+                workIndex--;
+                continue;
+            }
             if (node->left) {
                 // Branch node
@@ -696 +717 @@
 {
     for (int i = 0; i < leaf->num; i++) {
-        long index = leaf->contents[i]; // Index of point
+        psS64 index = leaf->contents[i]; // Index of point
         if (treeContentDistance(tree, index, coords) < distance) {
             psVectorAppend(result, index);
@@ -720 +741 @@
     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
+    // This is essentially the same as psTreeNearest, except pruning based on the search radius
 
     // Find the closest point in the leaf that contains the point of interest
@@ -734 +755 @@
         while (workIndex >= 0) {
             psTreeNode *node = work->data[workIndex];
+            if (treeBranchDistance(tree, node->index, coords) > distance) {
+                // No need to investigate
+                work->data[workIndex] = NULL;
+                workIndex--;
+                continue;
+            }
             if (node->left) {
                 // Branch node
@@ -801 +828 @@
         while (workIndex >= 0) {
             psTreeNode *node = work->data[workIndex];
+            if (treeBranchDistance(tree, node->index, coords) > distance) {
+                // No need to investigate
+                work->data[workIndex] = NULL;
+                workIndex--;
+                continue;
+            }
             if (node->left) {
                 // Branch node