Changeset d403a1 for molecuilder/src

Timestamp:

Jan 11, 2010, 3:06:30 AM (16 years ago)

Author:

Frederik Heber <heber@…>

Children:

91b1e79

Parents:

23f6ec

Message:

InsideOutside unit test of tesselation is working correctly.

• FIX: BoundaryTriangleSet::GetIntersectionInsideTriangle() - don't need helper, just check whether CrossPoint is returned (and true) for all of the three sides.
• FIX: Tesselation::IsInnerPoint() - projection onto plane and stuff was nonsense, just take the Point ans Intersection which is on the plane anyway.
• FIX: Vector::GetIntersectionOfTwoLinesOnPlane() - coefficient MUST be zero (then vectors are coplanar), but parallelity check was missing. Also, we have to check whether s is in [0,1] in order to see whether we are inside the triangle side or outside.

Signed-off-by: Frederik Heber <heber@tabletINS.(none)>

Location:

molecuilder/src

Files:

• tesselation.cpp (modified) (4 diffs)
• vector.cpp (modified) (4 diffs)

molecuilder/src/tesselation.cpp

@@ -420 +420 @@
     if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection, &NormalVector)) {
       CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
-      helper.CopyVector(endpoints[(i+1)%3]->node->node);
-      helper.SubtractVector(endpoints[i%3]->node->node);
+      i++;
+    } else 
       break;
-    } else 
-      i++;
   } while (i<3);
   if (i==3) {
-    eLog() << Verbose(0) << "Could not find any cross points, something's utterly wrong here!" << endl;
-  }
-  Log() << Verbose(1) << "INFO: Crosspoint is " << CrossPoint << "." << endl;
-
-  // check whether intersection is inside or not by comparing length of intersection and length of cross point
-  if ((CrossPoint.NormSquared() - helper.NormSquared()) < MYEPSILON) { // inside
+    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl;
     return true;
-  } else { // outside!
-    Intersection->Zero();
+  } else {
+    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl;
     return false;
   }
@@ -3241 +3234 @@
 bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC, const double epsilon) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   class BoundaryTriangleSet *result = FindClosestTriangleToVector(&Point, LC);
   Vector Center;
@@ -3252 +3245 @@
     return false;
   } else {
-    Log() << Verbose(1) << "INFO: Closest triangle found is " << *result << "." << endl;
+    Log() << Verbose(1) << "INFO: Closest triangle found is " << *result << " with normal vector " << result->NormalVector << "." << endl;
   }
 
@@ -3259 +3252 @@
   DistanceToCenter.CopyVector(&Center);
   DistanceToCenter.SubtractVector(&Point);
-  Log() << Verbose(2) << "INFO: Vector from point to test to center is " << Center << "." << endl;
+  Log() << Verbose(2) << "INFO: Vector from point to test to center is " << DistanceToCenter << "." << endl;
 
   // check whether we are on boundary
   if (fabs(DistanceToCenter.ScalarProduct(&result->NormalVector)) < MYEPSILON) {
     // calculate whether inside of triangle
-    DistanceToCenter.ProjectOntoPlane(&result->NormalVector);
-    DistanceToCenter.AddVector(&Center);
-    Center.CopyVector(&DistanceToCenter);
+    DistanceToCenter.CopyVector(&Point);
+    Center.CopyVector(&Point);
     Center.SubtractVector(&result->NormalVector); // points towards MolCenter
     DistanceToCenter.AddVector(&result->NormalVector); // points outside

molecuilder/src/vector.cpp

@@ -319 +319 @@
     M->Set(3, i, Line2b->x[i]);
   }
-  if (fabs(M->Determinant()) < MYEPSILON) 
+  Log() << Verbose(1) << "Coefficent matrix is:" << endl;
+  for (int i=0;i<4;i++) {
+    for (int j=0;j<4;j++)
+      cout << "\t" << M->Get(i,j);
+    cout << endl;
+  }
+  if (fabs(M->Determinant()) > MYEPSILON) {
+    Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
     return false;
+  }
+  Log() << Verbose(1) << "INFO: Line1a = " << *Line1a << ", Line1b = " << *Line1b << ", Line2a = " << *Line2a << ", Line2b = " << *Line2b << "." << endl;
+
 
   // constuct a,b,c
@@ -329 +339 @@
   c.CopyVector(Line2a);
   c.SubtractVector(Line1a);
+  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
+
+  // check for parallelity
+  Vector parallel;
+  parallel.CopyVector(&a);
+  parallel.SubtractVector(&b);
+  if (parallel.NormSquared() < MYEPSILON) {
+    Log() << Verbose(1) << "Lines are parallel." << endl;
+    return false;
+  }
 
   // obtain s
@@ -337 +357 @@
   temp2.CopyVector(&a);
   temp2.VectorProduct(&b);
-  s = temp1.ScalarProduct(&temp2)/temp2.NormSquared();
-  Log() << Verbose(1) << "Factor s is " << s << "." << endl;
+  Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
+  if (fabs(temp2.NormSquared()) > MYEPSILON)
+    s = temp1.ScalarProduct(&temp2)/temp2.NormSquared();
+  else
+    s = 0.;
+  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
 
   // construct intersection
@@ -346 +370 @@
   Log() << Verbose(1) << "Intersection is at " << *this << "." << endl;
 
-  return true;
+  if ((s >=0 ) && (s<=1))
+    return true;
+  else
+    return false;
 };