[cig-commits] r7949 - cs/spatialdata-0.1/trunk/libsrc/spatialdb

[brad at geodynamics.org]

Author: brad
Date: 2007-09-12 11:14:36 -0700 (Wed, 12 Sep 2007)
New Revision: 7949

Modified:
   cs/spatialdata-0.1/trunk/libsrc/spatialdb/SimpleDBQuery.cc
Log:
Fixed bug in colinear test in query with area interpolation.

Modified: cs/spatialdata-0.1/trunk/libsrc/spatialdb/SimpleDBQuery.cc
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/SimpleDBQuery.cc	2007-09-12 04:16:26 UTC (rev 7948)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/SimpleDBQuery.cc	2007-09-12 18:14:36 UTC (rev 7949)
@@ -429,6 +429,7 @@
     double dirABC[numCoords];
     _area(&areaABC, dirABC, ptA, ptB, ptC);
 
+#if 0
     // Alternate method of determining collinearity.
     // Compute unit vectors AB and AC, then compute the dot product.
     // If the absolute value of the dot product is somewhat less than 1,
@@ -448,21 +449,23 @@
     double abdotac = 0.0;
     for (int iDir=0; iDir < numCoords; ++iDir) {
       abdotac += (vecAB[iDir]/magAB) * (vecAC[iDir]/magAC);
-    }
-
+    } // for
+    
     const double tolerance = 0.98;
-
-    /* *** comment this out for now
+    if (fabs(abdotac) < tolerance) {
+#else
     // make sure A,B,C are not collinear by checking if area of
     // triangle ABC is not a tiny fraction of the distance AB
-
+    
     // length(ab)^2
-    const double ab2 = ptA[0]*ptB[0] + ptA[1]*ptB[1] + ptA[2]*ptB[2];
+    const double ab2 = 
+      pow(ptA[0]-ptB[0], 2) +
+      pow(ptA[1]-ptB[1], 2) +
+      pow(ptA[2]-ptB[2], 2);
 
     const double tolerance = 1.0e-06;
     if (areaABC > tolerance*ab2) {
-    */
-    if (fabs(abdotac) < tolerance) {
+#endif
       // project P onto abc plane
       double qProj[numCoords];
       const double qmod =