Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision ca8073901157e1302ad5c69a06c379e27ec66590)
+++ src/boundary.cpp	(revision 2b4a4087bfef04889a36e1b91688f2a56bea8e3a)
@@ -673,14 +673,27 @@
   // 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,47 +703,88 @@
 //    *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);
+  }
 };
 
