Changeset 7a3724 for molecuilder/src

Timestamp:

Oct 14, 2008, 7:14:10 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

46cd68

Parents:

51d33c5

Message:

new version of Tesselation::GuessStartingTriangle() seems to fix problems

+ now, we find the starting triangle in the following manner (point = boundary point!):

1. pick a point A
2. for all unequal pairs B,C (not A), calculate AB^2+BC2+AC^{2 and sort into ascending list}
3. Go through list, picking thus baseline B,C along with peak A of a candidate triangle
4. check whether all other boundary points lie only on one side
5. check whether no boundary point is inside the triangle
6. If all successful, starting triangle is found With this procedure, the convex hull of all SiOHCa3 clusters are found perfectly.

File:

• molecuilder/src/boundary.cpp (modified) (2 diffs)

molecuilder/src/boundary.cpp

@@ -673 +673 @@
   // 4b1. create all distances
   DistanceMultiMap DistanceMMap;
-  double distance;
-  for (PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
-    for(PointMap::iterator sprinter = PointsOnBoundary.begin(); sprinter != PointsOnBoundary.end(); sprinter++) {
-      if (runner->first < sprinter->first) {
-        distance = runner->second->node->x.Distance(&sprinter->second->node->x);
-        DistanceMMap.insert( DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator>(runner,sprinter) ) );
+  double distance, tmp;
+  Vector PlaneVector, TrialVector;
+  PointMap::iterator A, B, C; // three nodes of the first triangle
+  A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
+
+  // with A chosen, take each pair B,C and sort 
+  if (A != PointsOnBoundary.end()) {
+    B = A;
+    B++;
+    for (; B != PointsOnBoundary.end(); B++) {
+      C = B;
+      C++;
+      for (; C != PointsOnBoundary.end(); C++) {
+        tmp = A->second->node->x.Distance(&B->second->node->x);
+        distance = tmp*tmp;
+        tmp = A->second->node->x.Distance(&C->second->node->x);
+        distance += tmp*tmp;
+        tmp = B->second->node->x.Distance(&C->second->node->x);
+        distance += tmp*tmp;
+        DistanceMMap.insert( DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator>(B,C) ) );
       }
     }
   }
-
 //    // listing distances
 //    *out << Verbose(1) << "Listing DistanceMMap:";
@@ -690 +703 @@
 //    *out << endl;
   
-  // 4b2. take three smallest distance that form a triangle
-  // we take the smallest distance as the base line
+  // 4b2. pick three baselines forming a triangle
+  // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
   DistanceMultiMap::iterator baseline = DistanceMMap.begin();
-  BPS[0] = baseline->second.first->second;
-  BPS[1] = baseline->second.second->second;
-  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-
-  // take the second smallest as the second base line
-  DistanceMultiMap::iterator secondline = DistanceMMap.begin();
-  do {
-    secondline++;
-  } while (!(
-      ((BPS[0] == secondline->second.first->second) && (BPS[1] != secondline->second.second->second)) ||
-      ((BPS[0] == secondline->second.second->second) && (BPS[1] != secondline->second.first->second)) ||
-      ((BPS[1] == secondline->second.first->second) && (BPS[0] != secondline->second.second->second)) ||
-      ((BPS[1] == secondline->second.second->second) && (BPS[0] != secondline->second.first->second))
-    ));
-  BPS[0] = secondline->second.first->second;
-  BPS[1] = secondline->second.second->second;
-  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-
-  // connection yields the third line (note: first and second endpoint are sorted!)
-  if (baseline->second.first->second == secondline->second.first->second) {
-    SetEndpointsOrdered(BPS, baseline->second.second->second, secondline->second.second->second);
-  } else if (baseline->second.first->second == secondline->second.second->second) {
-    SetEndpointsOrdered(BPS, baseline->second.second->second, secondline->second.first->second);
-  } else if (baseline->second.second->second == secondline->second.first->second) {
-    SetEndpointsOrdered(BPS, baseline->second.first->second, baseline->second.second->second);
-  } else if (baseline->second.second->second == secondline->second.second->second) {
-    SetEndpointsOrdered(BPS, baseline->second.first->second, baseline->second.first->second);
-  }
-  BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-  
-  // 4b3. insert created triangle
-  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-  TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
-  TrianglesOnBoundaryCount++;
-  for(int i=0;i<NDIM;i++) {
-    LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
-    LinesOnBoundaryCount++;
-  }
-
-  *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
+  for (; baseline != DistanceMMap.end(); baseline++) {
+    // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
+    // 2. next, we have to check whether all points reside on only one side of the triangle
+    // 3. construct plane vector
+    PlaneVector.MakeNormalVector(&A->second->node->x, &baseline->second.first->second->node->x, &baseline->second.second->second->node->x);
+    *out << Verbose(2) << "Plane vector of candidate triangle is ";
+    PlaneVector.Output(out);
+    *out << endl;
+    // 4. loop over all points 
+    double sign = 0.;
+    PointMap::iterator checker = PointsOnBoundary.begin();
+    for (; checker != PointsOnBoundary.end(); checker++) {
+      // (neglecting A,B,C)
+      if ((checker == A) || (checker == baseline->second.first) || (checker == baseline->second.second))
+        continue;
+      // 4a. project onto plane vector
+      TrialVector.CopyVector(&checker->second->node->x);
+      TrialVector.SubtractVector(&A->second->node->x);
+      distance = TrialVector.Projection(&PlaneVector);
+      if (fabs(distance) < 1e-4)  // we need to have a small epsilon around 0 which is still ok
+        continue;
+      *out << Verbose(3) << "Projection of " << checker->second->node->Name << " yields distance of " << distance << "." << endl;
+      tmp = distance/fabs(distance);
+      // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
+      if ((sign != 0) && (tmp != sign)) {
+        // 4c. If so, break 4. loop and continue with next candidate in 1. loop
+        *out << Verbose(2) << "Current candidates: " << A->second->node->Name << "," << baseline->second.first->second->node->Name  << "," << baseline->second.second->second->node->Name << " leave " << checker->second->node->Name << " outside the convex hull." << endl;
+        break;
+      } else { // note the sign for later
+        *out << Verbose(2) << "Current candidates: " << A->second->node->Name << "," << baseline->second.first->second->node->Name  << "," << baseline->second.second->second->node->Name << " leave " << checker->second->node->Name << " inside the convex hull." << endl;
+        sign = tmp;
+      }
+      // 4d. Check whether the point is inside the triangle (check distance to each node
+      tmp = checker->second->node->x.Distance(&A->second->node->x);
+      int innerpoint = 0;
+      if ((tmp < A->second->node->x.Distance(&baseline->second.first->second->node->x)) 
+          && (tmp < A->second->node->x.Distance(&baseline->second.second->second->node->x)))
+        innerpoint++;
+      tmp = checker->second->node->x.Distance(&baseline->second.first->second->node->x);
+      if ((tmp < baseline->second.first->second->node->x.Distance(&A->second->node->x)) 
+          && (tmp < baseline->second.first->second->node->x.Distance(&baseline->second.second->second->node->x)))
+        innerpoint++;
+      tmp = checker->second->node->x.Distance(&baseline->second.second->second->node->x);
+      if ((tmp < baseline->second.second->second->node->x.Distance(&baseline->second.first->second->node->x)) 
+          && (tmp < baseline->second.second->second->node->x.Distance(&A->second->node->x)))
+        innerpoint++;
+      // 4e. If so, break 4. loop and continue with next candidate in 1. loop
+      if (innerpoint == 3)
+        break;
+    }
+    // 5. come this far, all on same side? Then break 1. loop and construct triangle
+    if (checker == PointsOnBoundary.end()) {
+      *out << "Looks like we have a candidate!" << endl;
+      break;
+    }
+  }
+  if (baseline != DistanceMMap.end()) {
+    BPS[0] = baseline->second.first->second;
+    BPS[1] = baseline->second.second->second;
+    BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+    BPS[0] = A->second;
+    BPS[1] = baseline->second.second->second;
+    BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+    BPS[0] = baseline->second.first->second;
+    BPS[1] = A->second;
+    BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);  
+   
+    // 4b3. insert created triangle
+    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+    TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
+    TrianglesOnBoundaryCount++;
+    for(int i=0;i<NDIM;i++) {
+      LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
+      LinesOnBoundaryCount++;
+    }
+  
+    *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
+  } else {
+    *out << Verbose(1) << "No starting triangle found." << endl;
+    exit(255);
+  }
 };