Changes in src/tesselation.cpp [9473f6:241485]

File:

• src/tesselation.cpp (modified) (29 diffs)

src/tesselation.cpp

@@ -35 +35 @@
  * \param *Walker TesselPoint this boundary point represents
  */
-BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
+BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker) :
   LinesCount(0),
   node(Walker),
@@ -61 +61 @@
  * \param *line line to add
  */
-void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
+void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -105 +105 @@
  * \param number number of the list
  */
-BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number)
+BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], const int number)
 {
         Info FunctionInfo(__func__);
@@ -115 +115 @@
   Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
   Point[1]->AddLine(this); //
-  // set skipped to false
-  skipped = false;
-  // clear triangles list
-  Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl;
-};
-
-/** Constructor of BoundaryLineSet with two endpoints.
- * Adds line automatically to each endpoints' LineMap
- * \param *Point1 first boundary point
- * \param *Point2 second boundary point
- * \param number number of the list
- */
-BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point1, BoundaryPointSet * const Point2, const int number)
-{
-  Info FunctionInfo(__func__);
-  // set number
-  Nr = number;
-  // set endpoints in ascending order
-  SetEndpointsOrdered(endpoints, Point1, Point2);
-  // add this line to the hash maps of both endpoints
-  Point1->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
-  Point2->AddLine(this); //
   // set skipped to false
   skipped = false;
@@ -193 +171 @@
  * \param *triangle to add
  */
-void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
+void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
 {
         Info FunctionInfo(__func__);
@@ -204 +182 @@
  * \return true - common endpoint present, false - not connected
  */
-bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
+bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -219 +197 @@
  * \return true - triangles are convex, false - concave or less than two triangles connected
  */
-bool BoundaryLineSet::CheckConvexityCriterion() const
+bool BoundaryLineSet::CheckConvexityCriterion()
 {
         Info FunctionInfo(__func__);
@@ -243 +221 @@
   int i=0;
   class BoundaryPointSet *node = NULL;
-  for(TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
+  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
     //Log() << Verbose(0) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
     NormalCheck.AddVector(&runner->second->NormalVector);
@@ -286 +264 @@
  * \return true - point is of the line, false - is not
  */
-bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
 {
         Info FunctionInfo(__func__);
@@ -299 +277 @@
  * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
  */
-class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
+class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point)
 {
         Info FunctionInfo(__func__);
@@ -339 +317 @@
  * \param number number of triangle
  */
-BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number) :
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number) :
   Nr(number)
 {
@@ -398 +376 @@
  * \param &OtherVector direction vector to make normal vector unique.
  */
-void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
+void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
 {
         Info FunctionInfo(__func__);
@@ -412 +390 @@
 /** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
  * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
- * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
  * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
  * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
@@ -422 +400 @@
  * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
  */
-bool BoundaryTriangleSet::GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const
+bool BoundaryTriangleSet::GetIntersectionInsideTriangle(Vector *MolCenter, Vector *x, Vector *Intersection)
 {
   Info FunctionInfo(__func__);
@@ -474 +452 @@
 };
 
-/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
- * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
- * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
- * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
- * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
- * the first two basepoints) or not.
- * \param *x point
- * \param *ClosestPoint desired closest point inside triangle to \a *x, is absolute vector
- * \return Distance squared between \a *x and closest point inside triangle
- */
-double BoundaryTriangleSet::GetClosestPointInsideTriangle(const Vector * const x, Vector * const ClosestPoint) const
-{
-  Info FunctionInfo(__func__);
-  Vector Direction;
-
-  // 1. get intersection with plane
-  Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl;
-  GetCenter(&Direction);
-  if (!ClosestPoint->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, x, &Direction)) {
-    ClosestPoint->CopyVector(x);
-  }
-
-  // 2. Calculate in plane part of line (x, intersection)
-  Vector InPlane;
-  InPlane.CopyVector(x);
-  InPlane.SubtractVector(ClosestPoint);  // points from plane intersection to straight-down point
-  InPlane.ProjectOntoPlane(&NormalVector);
-  InPlane.AddVector(ClosestPoint);
-
-  Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl;
-  Log() << Verbose(2) << "INFO: Line is from " << Direction << " to " << *x << "." << endl;
-  Log() << Verbose(2) << "INFO: In-plane part is " << InPlane << "." << endl;
-
-  // Calculate cross point between one baseline and the desired point such that distance is shortest
-  double ShortestDistance = -1.;
-  bool InsideFlag = false;
-  Vector CrossDirection[3];
-  Vector CrossPoint[3];
-  Vector helper;
-  for (int i=0;i<3;i++) {
-    // treat direction of line as normal of a (cut)plane and the desired point x as the plane offset, the intersect line with point
-    Direction.CopyVector(endpoints[(i+1)%3]->node->node);
-    Direction.SubtractVector(endpoints[i%3]->node->node);
-    // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
-    CrossPoint[i].GetIntersectionWithPlane(&Direction, &InPlane, endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node);
-    CrossDirection[i].CopyVector(&CrossPoint[i]);
-    CrossDirection[i].SubtractVector(&InPlane);
-    CrossPoint[i].SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
-    const double s = CrossPoint[i].ScalarProduct(&Direction)/Direction.NormSquared();
-    Log() << Verbose(2) << "INFO: Factor s is " << s << "." << endl;
-    if ((s >= -MYEPSILON) && ((s-1.) <= MYEPSILON)) {
-      CrossPoint[i].AddVector(endpoints[i%3]->node->node);  // make cross point absolute again
-      Log() << Verbose(2) << "INFO: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between " << *endpoints[i%3]->node->node << " and " << *endpoints[(i+1)%3]->node->node << "." << endl;
-      const double distance = CrossPoint[i].DistanceSquared(x);
-      if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
-        ShortestDistance = distance;
-        ClosestPoint->CopyVector(&CrossPoint[i]);
-      }
-    } else
-      CrossPoint[i].Zero();
-  }
-  InsideFlag = true;
-  for (int i=0;i<3;i++) {
-    const double sign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+1)%3]);
-    const double othersign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+2)%3]);;
-    if ((sign > -MYEPSILON) && (othersign > -MYEPSILON))  // have different sign
-      InsideFlag = false;
-  }
-  if (InsideFlag) {
-    ClosestPoint->CopyVector(&InPlane);
-    ShortestDistance = InPlane.DistanceSquared(x);
-  } else {  // also check endnodes
-    for (int i=0;i<3;i++) {
-      const double distance = x->DistanceSquared(endpoints[i]->node->node);
-      if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
-        ShortestDistance = distance;
-        ClosestPoint->CopyVector(endpoints[i]->node->node);
-      }
-    }
-  }
-  Log() << Verbose(1) << "INFO: Closest Point is " << *ClosestPoint << " with shortest squared distance is " << ShortestDistance << "." << endl;
-  return ShortestDistance;
-};
-
 /** Checks whether lines is any of the three boundary lines this triangle contains.
  * \param *line line to test
  * \return true - line is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
+bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -575 +469 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
 {
         Info FunctionInfo(__func__);
@@ -588 +482 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(const TesselPoint * const point) const
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class TesselPoint *point)
 {
         Info FunctionInfo(__func__);
@@ -601 +495 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(const BoundaryPointSet * const Points[3]) const
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3])
 {
         Info FunctionInfo(__func__);
@@ -624 +518 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(const BoundaryTriangleSet * const T) const
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T)
 {
         Info FunctionInfo(__func__);
@@ -646 +540 @@
  * \return pointer third endpoint or NULL if line does not belong to triangle.
  */
-class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const
+class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -663 +557 @@
  * \param *center central point on return.
  */
-void BoundaryTriangleSet::GetCenter(Vector * const center) const
+void BoundaryTriangleSet::GetCenter(Vector *center)
 {
         Info FunctionInfo(__func__);
@@ -670 +564 @@
     center->AddVector(endpoints[i]->node->node);
   center->Scale(1./3.);
-  Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl;
 }
 
@@ -3242 +3135 @@
   // for each point, check its lines, remember closest
   Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl;
-  LineSet ClosestLines;
+  BoundaryLineSet *ClosestLine = NULL;
   double MinDistance = 1e+16;
   Vector BaseLineIntersection;
@@ -3264 +3157 @@
 
       if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) {  // intersection would be outside, take closer endpoint
-        const double lengthEnd = Min(lengthEndA, lengthEndB);
-        if (lengthEnd - MinDistance < -MYEPSILON) { // new best line
-          ClosestLines.clear();
-          ClosestLines.insert(LineRunner->second);
-          MinDistance = lengthEnd;
-          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEnd << "." << endl;
-        } else if  (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate
-          ClosestLines.insert(LineRunner->second);
-          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "." << endl;
-        } else { // line is worse
+        if ((lengthEndA < MinDistance) && (lengthEndA < lengthEndB)) {
+          ClosestLine = LineRunner->second;
+          MinDistance = lengthEndA;
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEndA << "." << endl;
+        } else if (lengthEndB < MinDistance) {
+          ClosestLine = LineRunner->second;
+          MinDistance = lengthEndB;
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is closer with " << lengthEndB << "." << endl;
+        } else {
           Log() << Verbose(1) << "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "." << endl;
         }
@@ -3287 +3179 @@
           eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl;
         }
-        if ((ClosestLines.empty()) || (distance < MinDistance)) {
-          ClosestLines.insert(LineRunner->second);
+        if ((ClosestLine == NULL) || (distance < MinDistance)) {
+          ClosestLine = LineRunner->second;
           MinDistance = distance;
-          Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLines.begin() << " with projected distance " << MinDistance << "." << endl;
+          Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
         } else {
           Log() << Verbose(2) << "REJECT: Point is further away from line " << *LineRunner->second << " than present closest line: " << distance << " >> " << MinDistance << "." << endl;
@@ -3300 +3192 @@
 
   // check whether closest line is "too close" :), then it's inside
-  if (ClosestLines.empty()) {
+  if (ClosestLine == NULL) {
     Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
     return NULL;
   }
   TriangleList * candidates = new TriangleList;
-  for (LineSet::iterator LineRunner = ClosestLines.begin(); LineRunner != ClosestLines.end(); LineRunner++)
-    for (TriangleMap::iterator Runner = (*LineRunner)->triangles.begin(); Runner != (*LineRunner)->triangles.end(); Runner++) {
+  for (TriangleMap::iterator Runner = ClosestLine->triangles.begin(); Runner != ClosestLine->triangles.end(); Runner++) {
     candidates->push_back(Runner->second);
   }
@@ -3350 +3241 @@
 };
 
-/** Checks whether the provided Vector is within the Tesselation structure.
- * Basically calls Tesselation::GetDistanceToSurface() and checks the sign of the return value.
- * @param point of which to check the position
- * @param *LC LinkedCell structure
- *
- * @return true if the point is inside the Tesselation structure, false otherwise
- */
-bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
-{
-  return (GetDistanceSquaredToSurface(Point, LC) > MYEPSILON);
-}
-
-/** Returns the distance to the surface given by the tesselation.
+/** Checks whether the provided Vector is within the tesselation structure.
  * Calls FindClosestTriangleToVector() and checks whether the resulting triangle's BoundaryTriangleSet#NormalVector points
- * towards or away from the given \a &Point. Additionally, we check whether it's normal to the normal vector, i.e. on the
- * closest triangle's plane. Then, we have to check whether \a Point is inside the triangle or not to determine whether it's
- * an inside or outside point. This is done by calling BoundaryTriangleSet::GetIntersectionInsideTriangle().
- * In the end we additionally find the point on the triangle who was smallest distance to \a Point:
- *  -# Separate distance from point to center in vector in NormalDirection and on the triangle plane.
- *  -# Check whether vector on triangle plane points inside the triangle or crosses triangle bounds.
- *  -# If inside, take it to calculate closest distance
- *  -# If not, take intersection with BoundaryLine as distance
- *
- * @note distance is squared despite it still contains a sign to determine in-/outside!
+ * towards or away from the given \a &Point.
  *
  * @param point of which to check the position
  * @param *LC LinkedCell structure
  *
- * @return >0 if outside, ==0 if on surface, <0 if inside (Note that distance can be at most LinkedCell::RADIUS.)
- */
-double Tesselation::GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const
+ * @return true if the point is inside the tesselation structure, false otherwise
+ */
+bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
 {
   Info FunctionInfo(__func__);
@@ -3388 +3258 @@
   Vector DistanceToCenter;
   Vector Intersection;
-  double distance = 0.;
 
   if (result == NULL) {// is boundary point or only point in point cloud?
     Log() << Verbose(1) << Point << " is the only point in vicinity." << endl;
-    return LC->RADIUS;
+    return false;
   } else {
     Log() << Verbose(1) << "INFO: Closest triangle found is " << *result << " with normal vector " << result->NormalVector << "." << endl;
@@ -3413 +3282 @@
     if (result->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
       Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl;
-      return 0.;
+      return true;
     } else {
       Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl;
       return false;
     }
+  }
+
+  // then check direction to boundary
+  if (DistanceToCenter.ScalarProduct(&result->NormalVector) > MYEPSILON) {
+    Log() << Verbose(1) << Point << " is an inner point." << endl;
+    return true;
   } else {
-    // calculate smallest distance
-    distance = result->GetClosestPointInsideTriangle(&Point, &Intersection);
-    Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl;
-    distance = Min(distance, (LC->RADIUS*LC->RADIUS));
-
-    // then check direction to boundary
-    if (DistanceToCenter.ScalarProduct(&result->NormalVector) > MYEPSILON) {
-      Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl;
-      return -distance;
-    } else {
-      Log() << Verbose(1) << Point << " is NOT an inner point, " << distance << " above surface." << endl;
-      return +distance;
-    }
-  }
-};
+    Log() << Verbose(1) << Point << " is NOT an inner point." << endl;
+    return false;
+  }
+}
 
 /** Gets all points connected to the provided point by triangulation lines.