Index: src/tesselation.cpp
===================================================================
--- src/tesselation.cpp	(revision 86381b38132f1af7b1e83f6a7c271dc5decbac5c)
+++ src/tesselation.cpp	(revision 16d86622dc2c759142b589c64308c884f9339eb6)
@@ -10,11 +10,15 @@
 // ======================================== Points on Boundary =================================
 
+/** Constructor of BoundaryPointSet.
+ */
 BoundaryPointSet::BoundaryPointSet()
 {
   LinesCount = 0;
   Nr = -1;
-}
-;
-
+};
+
+/** Constructor of BoundaryPointSet with Tesselpoint.
+ * \param *Walker TesselPoint this boundary point represents
+ */
 BoundaryPointSet::BoundaryPointSet(TesselPoint *Walker)
 {
@@ -22,7 +26,10 @@
   LinesCount = 0;
   Nr = Walker->nr;
-}
-;
-
+};
+
+/** Destructor of BoundaryPointSet.
+ * Sets node to NULL to avoid removing the original, represented TesselPoint.
+ * \note When removing point from a class Tesselation, use RemoveTesselationPoint()
+ */
 BoundaryPointSet::~BoundaryPointSet()
 {
@@ -31,7 +38,9 @@
     cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
   node = NULL;
-}
-;
-
+};
+
+/** Add a line to the LineMap of this point.
+ * \param *line line to add
+ */
 void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
 {
@@ -47,9 +56,11 @@
     }
   LinesCount++;
-}
-;
-
-ostream &
-operator <<(ostream &ost, BoundaryPointSet &a)
+};
+
+/** output operator for BoundaryPointSet.
+ * \param &ost output stream
+ * \param &a boundary point
+ */
+ostream & operator <<(ostream &ost, BoundaryPointSet &a)
 {
   ost << "[" << a.Nr << "|" << a.node->Name << "]";
@@ -60,4 +71,6 @@
 // ======================================== Lines on Boundary =================================
 
+/** Constructor of BoundaryLineSet.
+ */
 BoundaryLineSet::BoundaryLineSet()
 {
@@ -65,7 +78,11 @@
     endpoints[i] = NULL;
   Nr = -1;
-}
-;
-
+};
+
+/** Constructor of BoundaryLineSet with two endpoints.
+ * Adds line automatically to each endpoints' LineMap
+ * \param *Point[2] array of two boundary points
+ * \param number number of the list
+ */
 BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
 {
@@ -79,44 +96,60 @@
   // clear triangles list
   cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
-}
-;
-
+};
+
+/** Destructor for BoundaryLineSet.
+ * Removes itself from each endpoints' LineMap, calling RemoveTrianglePoint() when point not connected anymore.
+ * \note When removing lines from a class Tesselation, use RemoveTesselationLine()
+ */
 BoundaryLineSet::~BoundaryLineSet()
 {
   int Numbers[2];
-  Numbers[0] = endpoints[1]->Nr;
-  Numbers[1] = endpoints[0]->Nr;
+
+  // get other endpoint number of finding copies of same line
+  if (endpoints[1] != NULL)
+    Numbers[0] = endpoints[1]->Nr;
+  else
+    Numbers[0] = -1;
+  if (endpoints[0] != NULL)
+    Numbers[1] = endpoints[0]->Nr;
+  else
+    Numbers[1] = -1;
+
   for (int i = 0; i < 2; i++) {
-    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
-    // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
-    pair<LineMap::iterator, LineMap::iterator> erasor = endpoints[i]->lines.equal_range(Numbers[i]);
-    for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
-      if ((*Runner).second == this) {
-        endpoints[i]->lines.erase(Runner);
-        break;
+    if (endpoints[i] != NULL) {
+      if (Numbers[i] != -1) { // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
+        pair<LineMap::iterator, LineMap::iterator> erasor = endpoints[i]->lines.equal_range(Numbers[i]);
+        for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
+          if ((*Runner).second == this) {
+            cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+            endpoints[i]->lines.erase(Runner);
+            break;
+          }
+      } else { // there's just a single line left
+        if (endpoints[i]->lines.erase(Nr))
+          cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
       }
-    if (endpoints[i]->lines.empty()) {
-      cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
-      if (endpoints[i] != NULL) {
-        delete(endpoints[i]);
-        endpoints[i] = NULL;
-      } else
-        cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
-    } else
-      cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
+      if (endpoints[i]->lines.empty()) {
+        cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+        if (endpoints[i] != NULL) {
+          delete(endpoints[i]);
+          endpoints[i] = NULL;
+        }
+      }
+    }
   }
   if (!triangles.empty())
     cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;
-}
-;
-
-void
-BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
+};
+
+/** Add triangle to TriangleMap of this boundary line.
+ * \param *triangle to add
+ */
+void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
 {
   cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
       << endl;
   triangles.insert(TrianglePair(triangle->Nr, triangle));
-}
-;
+};
 
 /** Checks whether we have a common endpoint with given \a *line.
@@ -218,16 +251,18 @@
 };
 
-
-ostream &
-operator <<(ostream &ost, BoundaryLineSet &a)
+/** output operator for BoundaryLineSet.
+ * \param &ost output stream
+ * \param &a boundary line
+ */
+ostream & operator <<(ostream &ost, BoundaryLineSet &a)
 {
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "]";
   return ost;
-}
-;
+};
 
 // ======================================== Triangles on Boundary =================================
 
-
+/** Constructor for BoundaryTriangleSet.
+ */
 BoundaryTriangleSet::BoundaryTriangleSet()
 {
@@ -238,7 +273,10 @@
     }
   Nr = -1;
-}
-;
-
+};
+
+/** Constructor for BoundaryTriangleSet with three lines.
+ * \param *line[3] lines that make up the triangle
+ * \param number number of triangle
+ */
 BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number)
 {
@@ -279,24 +317,25 @@
     }
   cout << "." << endl;
-}
-;
-
+};
+
+/** Destructor of BoundaryTriangleSet.
+ * Removes itself from each of its lines' LineMap and removes them if necessary.
+ * \note When removing triangles from a class Tesselation, use RemoveTesselationTriangle()
+ */
 BoundaryTriangleSet::~BoundaryTriangleSet()
 {
   for (int i = 0; i < 3; i++) {
-    cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
-    lines[i]->triangles.erase(Nr);
-    if (lines[i]->triangles.empty()) {
-      if (lines[i] != NULL) {
-        cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
-        delete (lines[i]);
-        lines[i] = NULL;
-      } else
-        cerr << "ERROR: This line " << i << " has already been free'd." << endl;
-    } else
-      cout << Verbose(5) << *lines[i] << " is still attached to another triangle." << endl;
-  }
-}
-;
+    if (lines[i] != NULL) {
+      if (lines[i]->triangles.erase(Nr))
+        cout << Verbose(5) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;
+      if (lines[i]->triangles.empty()) {
+          cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+          delete (lines[i]);
+          lines[i] = NULL;
+      }
+    }
+  }
+  cout << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
+};
 
 /** Calculates the normal vector for this triangle.
@@ -435,12 +474,14 @@
 }
 
-ostream &
-operator <<(ostream &ost, BoundaryTriangleSet &a)
+/** output operator for BoundaryTriangleSet.
+ * \param &ost output stream
+ * \param &a boundary triangle
+ */
+ostream &operator <<(ostream &ost, BoundaryTriangleSet &a)
 {
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
       << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
   return ost;
-}
-;
+};
 
 // =========================================================== class TESSELPOINT ===========================================
@@ -1044,5 +1085,5 @@
         // add Walker to boundary points
         *out << Verbose(2) << "Adding " << *Walker << " to BoundaryPoints." << endl;
-        if (AddPoint(Walker,0))
+        if (AddBoundaryPoint(Walker,0))
           NewPoint = BPS[0];
         else
@@ -1097,5 +1138,5 @@
 };
 
-/** Adds an point to the tesselation::PointsOnBoundary list.
+/** Adds a point to the tesselation::PointsOnBoundary list.
  * \param *Walker point to add
  * \param n TesselStruct::BPS index to put pointer into
@@ -1103,5 +1144,5 @@
  */
 bool
-Tesselation::AddPoint(TesselPoint *Walker, int n)
+Tesselation::AddBoundaryPoint(TesselPoint *Walker, int n)
 {
   PointTestPair InsertUnique;
@@ -1125,5 +1166,5 @@
  */
 void
-Tesselation::AddTrianglePoint(TesselPoint* Candidate, int n)
+Tesselation::AddTesselationPoint(TesselPoint* Candidate, int n)
 {
   PointTestPair InsertUnique;
@@ -1147,5 +1188,5 @@
  * @param n index of Tesselation::BLS giving the line with both endpoints
  */
-void Tesselation::AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) {
+void Tesselation::AddTesselationLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) {
   bool insertNewLine = true;
 
@@ -1171,5 +1212,5 @@
 
   if (insertNewLine) {
-    AlwaysAddTriangleLine(a, b, n);
+    AlwaysAddTesselationTriangleLine(a, b, n);
   }
 }
@@ -1184,5 +1225,5 @@
  * @param n index of Tesselation::BLS giving the line with both endpoints
  */
-void Tesselation::AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
+void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
 {
   cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
@@ -1199,6 +1240,5 @@
  * Furthermore it adds the triangle to all of its lines, in order to recognize those which are saturated later.
  */
-void
-Tesselation::AddTriangle()
+void Tesselation::AddTesselationTriangle()
 {
   cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
@@ -1209,6 +1249,98 @@
 
   // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
-}
-;
+};
+
+/** Removes a triangle from the tesselation.
+ * Removes itself from the TriangleMap's of its lines, calls for them RemoveTriangleLine() if they are no more connected.
+ * Removes itself from memory.
+ * \param *triangle to remove
+ */
+void Tesselation::RemoveTesselationTriangle(class BoundaryTriangleSet *triangle)
+{
+  if (triangle == NULL)
+    return;
+  for (int i = 0; i < 3; i++) {
+    if (triangle->lines[i] != NULL) {
+      cout << Verbose(5) << "Removing triangle Nr." << triangle->Nr << " in line " << *triangle->lines[i] << "." << endl;
+      triangle->lines[i]->triangles.erase(triangle->Nr);
+      if (triangle->lines[i]->triangles.empty()) {
+          cout << Verbose(5) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl;
+          RemoveTesselationLine(triangle->lines[i]);
+          triangle->lines[i] = NULL;
+      } else
+        cout << Verbose(5) << *triangle->lines[i] << " is still attached to another triangle." << endl;
+    } else
+      cerr << "ERROR: This line " << i << " has already been free'd." << endl;
+  }
+
+  if (TrianglesOnBoundary.erase(triangle->Nr))
+    cout << Verbose(5) << "Removing triangle Nr. " << triangle->Nr << "." << endl;
+  delete(triangle);
+};
+
+/** Removes a line from the tesselation.
+ * Removes itself from each endpoints' LineMap, then removes itself from global LinesOnBoundary list and free's the line.
+ * \param *line line to remove
+ */
+void Tesselation::RemoveTesselationLine(class BoundaryLineSet *line)
+{
+  int Numbers[2];
+
+  if (line == NULL)
+    return;
+  // get other endpoint number of finding copies of same line
+  if (line->endpoints[1] != NULL)
+    Numbers[0] = line->endpoints[1]->Nr;
+  else
+    Numbers[0] = -1;
+  if (line->endpoints[0] != NULL)
+    Numbers[1] = line->endpoints[0]->Nr;
+  else
+    Numbers[1] = -1;
+
+  for (int i = 0; i < 2; i++) {
+    if (line->endpoints[i] != NULL) {
+      if (Numbers[i] != -1) { // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
+        pair<LineMap::iterator, LineMap::iterator> erasor = line->endpoints[i]->lines.equal_range(Numbers[i]);
+        for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
+          if ((*Runner).second == line) {
+            cout << Verbose(5) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl;
+            line->endpoints[i]->lines.erase(Runner);
+            break;
+          }
+      } else { // there's just a single line left
+        if (line->endpoints[i]->lines.erase(line->Nr))
+          cout << Verbose(5) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl;
+      }
+      if (line->endpoints[i]->lines.empty()) {
+        cout << Verbose(5) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+        RemoveTesselationPoint(line->endpoints[i]);
+        line->endpoints[i] = NULL;
+      } else
+        cout << Verbose(5) << *line->endpoints[i] << " has still lines it's attached to." << endl;
+    } else
+      cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
+  }
+  if (!line->triangles.empty())
+    cerr << "WARNING: Memory Leak! I " << *line << " am still connected to some triangles." << endl;
+
+  if (LinesOnBoundary.erase(line->Nr))
+    cout << Verbose(5) << "Removing line Nr. " << line->Nr << "." << endl;
+  delete(line);
+};
+
+/** Removes a point from the tesselation.
+ * Checks whether there are still lines connected, removes from global PointsOnBoundary list, then free's the point.
+ * \note If a point should be removed, while keep the tesselated surface intact (i.e. closed), use RemovePointFromTesselatedSurface()
+ * \param *point point to remove
+ */
+void Tesselation::RemoveTesselationPoint(class BoundaryPointSet *point)
+{
+  if (point == NULL)
+    return;
+  if (PointsOnBoundary.erase(point->Nr))
+    cout << Verbose(5) << "Removing point Nr. " << point->Nr << "." << endl;
+  delete(point);
+};
 
 /** Checks whether the triangle consisting of the three points is already present.
@@ -1359,7 +1491,7 @@
 
     // adding point 1 and point 2 and add the line between them
-    AddTrianglePoint(FirstPoint, 0);
-    AddTrianglePoint(SecondPoint, 1);
-    AddTriangleLine(TPS[0], TPS[1], 0);
+    AddTesselationPoint(FirstPoint, 0);
+    AddTesselationPoint(SecondPoint, 1);
+    AddTesselationLine(TPS[0], TPS[1], 0);
 
     //cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n";
@@ -1375,11 +1507,11 @@
     for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
       // add third triangle point
-      AddTrianglePoint((*it)->point, 2);
+      AddTesselationPoint((*it)->point, 2);
       // add the second and third line
-      AddTriangleLine(TPS[1], TPS[2], 1);
-      AddTriangleLine(TPS[0], TPS[2], 2);
+      AddTesselationLine(TPS[1], TPS[2], 1);
+      AddTesselationLine(TPS[0], TPS[2], 2);
       // ... and triangles to the Maps of the Tesselation class
       BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-      AddTriangle();
+      AddTesselationTriangle();
       // ... and calculate its normal vector (with correct orientation)
       (*it)->OptCenter.Scale(-1.);
@@ -1394,7 +1526,7 @@
         SecondPoint = (*it)->point;
         // adding point 1 and point 2 and the line between them
-        AddTrianglePoint(FirstPoint, 0);
-        AddTrianglePoint(SecondPoint, 1);
-        AddTriangleLine(TPS[0], TPS[1], 0);
+        AddTesselationPoint(FirstPoint, 0);
+        AddTesselationPoint(SecondPoint, 1);
+        AddTesselationLine(TPS[0], TPS[1], 0);
       }
       cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
@@ -1528,15 +1660,15 @@
     // If there is no triangle, add it regularly.
     if (existentTrianglesCount == 0) {
-      AddTrianglePoint((*it)->point, 0);
-      AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
-      AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
+      AddTesselationPoint((*it)->point, 0);
+      AddTesselationPoint(BaseRay->endpoints[0]->node, 1);
+      AddTesselationPoint(BaseRay->endpoints[1]->node, 2);
 
       if (CheckLineCriteriaforDegeneratedTriangle(TPS)) {
-        AddTriangleLine(TPS[0], TPS[1], 0);
-        AddTriangleLine(TPS[0], TPS[2], 1);
-        AddTriangleLine(TPS[1], TPS[2], 2);
+        AddTesselationLine(TPS[0], TPS[1], 0);
+        AddTesselationLine(TPS[0], TPS[2], 1);
+        AddTesselationLine(TPS[1], TPS[2], 2);
 
         BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-        AddTriangle();
+        AddTesselationTriangle();
         (*it)->OptCenter.Scale(-1.);
         BTS->GetNormalVector((*it)->OptCenter);
@@ -1555,17 +1687,17 @@
       }
     } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
-        AddTrianglePoint((*it)->point, 0);
-        AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
-        AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
+        AddTesselationPoint((*it)->point, 0);
+        AddTesselationPoint(BaseRay->endpoints[0]->node, 1);
+        AddTesselationPoint(BaseRay->endpoints[1]->node, 2);
 
         // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
         // i.e. at least one of the three lines must be present with TriangleCount <= 1
         if (CheckLineCriteriaforDegeneratedTriangle(TPS)) {
-          AddTriangleLine(TPS[0], TPS[1], 0);
-          AddTriangleLine(TPS[0], TPS[2], 1);
-          AddTriangleLine(TPS[1], TPS[2], 2);
+          AddTesselationLine(TPS[0], TPS[1], 0);
+          AddTesselationLine(TPS[0], TPS[2], 1);
+          AddTesselationLine(TPS[1], TPS[2], 2);
 
           BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          AddTriangle();  // add to global map
+          AddTesselationTriangle();  // add to global map
 
           (*it)->OtherOptCenter.Scale(-1.);
@@ -1608,27 +1740,98 @@
 };
 
-
-/** Goes over all baselines and checks whether adjacent triangles and convex to each other.
+/** Checks whether the quadragon of the two triangles connect to \a *Base is convex.
+ * We look whether the closest point on \a *Base with respect to the other baseline is outside
+ * of the segment formed by both endpoints (concave) or not (convex).
  * \param *out output stream for debugging
- * \return true - all baselines were corrected, false - there are still concave pieces
- */
-bool Tesselation::CorrectConcaveBaselines(ofstream *out)
-{
-  class BoundaryLineSet *OldLines[4], *NewLine;
-  class BoundaryPointSet *OldPoints[2];
-  Vector BaseLineNormal;
-  class BoundaryLineSet *Base = NULL;
-  int OldTriangles[2], OldBaseLine;
-  int i,m;
-
-  *out << Verbose(1) << "Begin of CorrectConcaveBaselines" << endl;
-
-  // go through every baseline
-  LineMap::iterator Advancer = LinesOnBoundary.begin();
-  for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline = Advancer) {
-    Advancer++;  // as the base line might be removed, we rather have the next one already present
-    Base = baseline->second;
-    *out << Verbose(2) << "Current baseline is " << *Base << " ... " << endl;
-    // calculate NormalVector for later use
+ * \param *Base line to be flipped
+ * \return NULL - concave, otherwise endpoint that makes it concave
+ */
+class BoundaryPointSet *Tesselation::IsConvexRectangle(ofstream *out, class BoundaryLineSet *Base)
+{
+  class BoundaryPointSet *Spot = NULL;
+  class BoundaryLineSet *OtherBase;
+  Vector *ClosestPoint[2];
+
+  int m=0;
+  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+    for (int j=0;j<3;j++) // all of their endpoints and baselines
+      if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
+        BPS[m++] = runner->second->endpoints[j];
+  OtherBase = new class BoundaryLineSet(BPS,-1);
+
+  *out << Verbose(3) << "INFO: Current base line is " << *Base << "." << endl;
+  *out << Verbose(3) << "INFO: Other base line is " << *OtherBase << "." << endl;
+
+  // get the closest point on each line to the other line
+  ClosestPoint[0] = GetClosestPointBetweenLine(out, Base, OtherBase);
+  ClosestPoint[1] = GetClosestPointBetweenLine(out, OtherBase, Base);
+
+  // delete the temporary other base line
+  delete(OtherBase);
+
+  // get the distance vector from Base line to OtherBase line
+  Vector DistanceToIntersection, BaseLine;
+  BaseLine.CopyVector(Base->endpoints[1]->node->node);
+  BaseLine.SubtractVector(Base->endpoints[0]->node->node);
+  DistanceToIntersection.CopyVector(ClosestPoint[0]);
+  DistanceToIntersection.SubtractVector(Base->endpoints[0]->node->node);
+  Spot = Base->endpoints[1];
+  if (BaseLine.ScalarProduct(&DistanceToIntersection) < -MYEPSILON) {
+    DistanceToIntersection.CopyVector(ClosestPoint[0]);
+    DistanceToIntersection.SubtractVector(Base->endpoints[1]->node->node);
+    Spot = Base->endpoints[0];
+    if (BaseLine.ScalarProduct(&DistanceToIntersection) < -MYEPSILON) {
+      *out << Verbose(3) << "ERROR: Something's very wrong here, both SKPs return negative?!" << endl;
+      return NULL;
+    }
+  }
+  if ((BaseLine.NormSquared() - DistanceToIntersection.NormSquared()) < -MYEPSILON) {  // distance is smaller: concave
+    *out << Verbose(3) << "INFO: Rectangle of triangles of base line " << *Base << " is concave: Base line squared length " << BaseLine.NormSquared() << " against Distance to intersection squared " << DistanceToIntersection.NormSquared() << "." << endl;
+    return Spot;
+  } else {  // base line is longer: convex
+    *out << Verbose(3) << "INFO: Rectangle of triangles of base line " << *Base << " is concave." << endl;
+    return NULL;
+  }
+
+};
+
+
+/** For a given boundary line \a *Base and its two triangles, picks the central baseline that is "higher".
+ * \param *out output stream for debugging
+ * \param *Base line to be flipped
+ * \return true - line was changed, false - same line as before
+ */
+bool Tesselation::PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base)
+{
+  class BoundaryLineSet *OtherBase;
+  Vector *ClosestPoint[2];
+
+  int m=0;
+  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+    for (int j=0;j<3;j++) // all of their endpoints and baselines
+      if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
+        BPS[m++] = runner->second->endpoints[j];
+  OtherBase = new class BoundaryLineSet(BPS,-1);
+
+  *out << Verbose(3) << "INFO: Current base line is " << *Base << "." << endl;
+  *out << Verbose(3) << "INFO: Other base line is " << *OtherBase << "." << endl;
+
+  // get the closest point on each line to the other line
+  ClosestPoint[0] = GetClosestPointBetweenLine(out, Base, OtherBase);
+  ClosestPoint[1] = GetClosestPointBetweenLine(out, OtherBase, Base);
+
+  // get the distance vector from Base line to OtherBase line
+  Vector Distance;
+  Distance.CopyVector(ClosestPoint[1]);
+  Distance.SubtractVector(ClosestPoint[0]);
+
+  // delete the temporary other base line
+  delete(OtherBase);
+
+  if (Distance.NormSquared() < MYEPSILON) { // check for intersection
+    *out << Verbose(3) << "REJECT: Both lines have an intersection: Nothing to do." << endl;
+    return false;
+  } else { // check for sign against BaseLineNormal
+    Vector BaseLineNormal;
     BaseLineNormal.Zero();
     if (Base->triangles.size() < 2) {
@@ -1641,97 +1844,167 @@
     }
     BaseLineNormal.Scale(-1./2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
-    // check convexity
-    if (Base->CheckConvexityCriterion(out)) { // triangles are convex
-      *out << Verbose(2) << "... has two convex triangles." << endl;
-    } else { // not convex!
-      *out << Verbose(2) << "... has two concave triangles!" << endl;
-      // get the two triangles
-      // gather four endpoints and four lines
-      for (int j=0;j<4;j++)
-        OldLines[j] = NULL;
-      for (int j=0;j<2;j++)
-        OldPoints[j] = NULL;
-      i=0;
-      m=0;
-      *out << Verbose(3) << "The four old lines are: ";
-      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-        for (int j=0;j<3;j++) // all of their endpoints and baselines
-          if (runner->second->lines[j] != Base) { // pick not the central baseline
-            OldLines[i++] = runner->second->lines[j];
-            *out << *runner->second->lines[j] << "\t";
-          }
-      *out << endl;
-      *out << Verbose(3) << "The two old points are: ";
-      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-        for (int j=0;j<3;j++) // all of their endpoints and baselines
-          if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) { // and neither of its endpoints
-            OldPoints[m++] = runner->second->endpoints[j];
-            *out << *runner->second->endpoints[j] << "\t";
-          }
-      *out << endl;
-      if (i<4) {
-        *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-        return false;
+
+    if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
+      *out << Verbose(2) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl;
+      FlipBaseline(out, Base);
+      return true;
+    } else {  // Base higher than OtherBase -> do nothing
+      *out << Verbose(2) << "REJECT: Base line is higher: Nothing to do." << endl;
+      return false;
+    }
+  }
+};
+
+/** Returns the closest point on \a *Base with respect to \a *OtherBase.
+ * \param *out output stream for debugging
+ * \param *Base reference line
+ * \param *OtherBase other base line
+ * \return Vector on reference line that has closest distance
+ */
+Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase)
+{
+  // construct the plane of the two baselines (i.e. take both their directional vectors)
+  Vector Normal;
+  Vector OtherBaseline;
+  Normal.CopyVector(Base->endpoints[1]->node->node);
+  Normal.SubtractVector(Base->endpoints[0]->node->node);
+  OtherBaseline.CopyVector(OtherBase->endpoints[1]->node->node);
+  OtherBaseline.SubtractVector(OtherBase->endpoints[0]->node->node);
+  Normal.VectorProduct(&OtherBaseline);
+  Normal.Normalize();
+
+  // project one offset point of OtherBase onto this plane
+  Vector NewOffset;
+  NewOffset.CopyVector(OtherBase->endpoints[0]->node->node);
+  NewOffset.ProjectOntoPlane(&Normal);
+  Vector NewDirection;
+  NewDirection.CopyVector(OtherBase->endpoints[1]->node->node);
+  NewDirection.ProjectOntoPlane(&Normal);
+
+  // calculate the intersection between this projected baseline and Base
+  Vector *Intersection = new Vector;
+  Intersection->GetIntersectionOfTwoLinesOnPlane(out, Base->endpoints[0]->node->node, Base->endpoints[1]->node->node, &NewOffset, &NewDirection, &Normal);
+
+  return Intersection;
+};
+
+/** For a given baseline and its two connected triangles, flips the baseline.
+ * I.e. we create the new baseline between the other two endpoints of these four
+ * endpoints and reconstruct the two triangles accordingly.
+ * \param *out output stream for debugging
+ * \param *Base line to be flipped
+ * \return true - flipping successful, false - something went awry
+ */
+bool Tesselation::FlipBaseline(ofstream *out, class BoundaryLineSet *Base)
+{
+  class BoundaryLineSet *OldLines[4], *NewLine;
+  class BoundaryPointSet *OldPoints[2];
+  Vector BaseLineNormal;
+  int OldTriangleNrs[2], OldBaseLineNr;
+  int i,m;
+
+  *out << Verbose(1) << "Begin of FlipBaseline" << endl;
+
+  // calculate NormalVector for later use
+  BaseLineNormal.Zero();
+  if (Base->triangles.size() < 2) {
+    *out << Verbose(2) << "ERROR: Less than two triangles are attached to this baseline!" << endl;
+    return false;
+  }
+  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
+    *out << Verbose(4) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl;
+    BaseLineNormal.AddVector(&(runner->second->NormalVector));
+  }
+  BaseLineNormal.Scale(-1./2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
+
+  // get the two triangles
+  // gather four endpoints and four lines
+  for (int j=0;j<4;j++)
+    OldLines[j] = NULL;
+  for (int j=0;j<2;j++)
+    OldPoints[j] = NULL;
+  i=0;
+  m=0;
+  *out << Verbose(3) << "The four old lines are: ";
+  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+    for (int j=0;j<3;j++) // all of their endpoints and baselines
+      if (runner->second->lines[j] != Base) { // pick not the central baseline
+        OldLines[i++] = runner->second->lines[j];
+        *out << *runner->second->lines[j] << "\t";
       }
-      for (int j=0;j<4;j++)
-        if (OldLines[j] == NULL) {
-          *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-          return false;
-        }
-      for (int j=0;j<2;j++)
-        if (OldPoints[j] == NULL) {
-          *out << Verbose(1) << "ERROR: We have not gathered enough endpoints!" << endl;
-          return false;
-        }
-
-      // remove triangles and baseline removes itself
-      *out << Verbose(3) << "INFO: Deleting baseline " << *Base << "." << endl;
-      OldBaseLine = Base->Nr;
-      LinesOnBoundary.erase(OldBaseLine);
-      m=0;
-      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
-        *out << Verbose(3) << "INFO: Deleting triangle " << *(runner->second) << "." << endl;
-        TrianglesOnBoundary.erase(OldTriangles[m++] = runner->second->Nr);
-        delete(runner->second);
+  *out << endl;
+  *out << Verbose(3) << "The two old points are: ";
+  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+    for (int j=0;j<3;j++) // all of their endpoints and baselines
+      if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) { // and neither of its endpoints
+        OldPoints[m++] = runner->second->endpoints[j];
+        *out << *runner->second->endpoints[j] << "\t";
       }
-
-      // construct new baseline (with same number as old one)
-      BPS[0] = OldPoints[0];
-      BPS[1] = OldPoints[1];
-      NewLine = new class BoundaryLineSet(BPS, OldBaseLine);
-      LinesOnBoundary.insert(LinePair(OldBaseLine, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
-      *out << Verbose(3) << "INFO: Created new baseline " << *NewLine << "." << endl;
-
-      // construct new triangles with flipped baseline
-      i=-1;
-      if (OldLines[0]->IsConnectedTo(OldLines[2]))
-        i=2;
-      if (OldLines[0]->IsConnectedTo(OldLines[3]))
-        i=3;
-      if (i!=-1) {
-        BLS[0] = OldLines[0];
-        BLS[1] = OldLines[i];
-        BLS[2] = NewLine;
-        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[0]);
-        BTS->GetNormalVector(BaseLineNormal);
-        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[0], BTS));
-        *out << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl;
-
-        BLS[0] = (i==2 ? OldLines[3] : OldLines[2]);
-        BLS[1] = OldLines[1];
-        BLS[2] = NewLine;
-        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[1]);
-        BTS->GetNormalVector(BaseLineNormal);
-        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[1], BTS));
-        *out << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl;
-      } else {
-        *out << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;
-        return false;
-      }
-    }
-  }
-  *out << Verbose(1) << "End of CorrectConcaveBaselines" << endl;
+  *out << endl;
+
+  // check whether everything is in place to create new lines and triangles
+  if (i<4) {
+    *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
+    return false;
+  }
+  for (int j=0;j<4;j++)
+    if (OldLines[j] == NULL) {
+      *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
+      return false;
+    }
+  for (int j=0;j<2;j++)
+    if (OldPoints[j] == NULL) {
+      *out << Verbose(1) << "ERROR: We have not gathered enough endpoints!" << endl;
+      return false;
+    }
+
+  // remove triangles and baseline removes itself
+  *out << Verbose(3) << "INFO: Deleting baseline " << *Base << " from global list." << endl;
+  OldBaseLineNr = Base->Nr;
+  m=0;
+  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
+    *out << Verbose(3) << "INFO: Deleting triangle " << *(runner->second) << "." << endl;
+    OldTriangleNrs[m++] = runner->second->Nr;
+    RemoveTesselationTriangle(runner->second);
+  }
+
+  // construct new baseline (with same number as old one)
+  BPS[0] = OldPoints[0];
+  BPS[1] = OldPoints[1];
+  NewLine = new class BoundaryLineSet(BPS, OldBaseLineNr);
+  LinesOnBoundary.insert(LinePair(OldBaseLineNr, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
+  *out << Verbose(3) << "INFO: Created new baseline " << *NewLine << "." << endl;
+
+  // construct new triangles with flipped baseline
+  i=-1;
+  if (OldLines[0]->IsConnectedTo(OldLines[2]))
+    i=2;
+  if (OldLines[0]->IsConnectedTo(OldLines[3]))
+    i=3;
+  if (i!=-1) {
+    BLS[0] = OldLines[0];
+    BLS[1] = OldLines[i];
+    BLS[2] = NewLine;
+    BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[0]);
+    BTS->GetNormalVector(BaseLineNormal);
+    TrianglesOnBoundary.insert(TrianglePair(OldTriangleNrs[0], BTS));
+    *out << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl;
+
+    BLS[0] = (i==2 ? OldLines[3] : OldLines[2]);
+    BLS[1] = OldLines[1];
+    BLS[2] = NewLine;
+    BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[1]);
+    BTS->GetNormalVector(BaseLineNormal);
+    TrianglesOnBoundary.insert(TrianglePair(OldTriangleNrs[1], BTS));
+    *out << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl;
+  } else {
+    *out << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;
+    return false;
+  }
+
+  *out << Verbose(1) << "End of FlipBaseline" << endl;
   return true;
 };
+
 
 /** Finds the second point of starting triangle.
@@ -2100,5 +2373,7 @@
     return NULL;
   }
-  list<TesselPoint*> *connectedClosestPoints = getCircleOfConnectedPoints(out, trianglePoints[0], x);
+  list<TesselPoint*> *connectedPoints = getCircleOfConnectedPoints(out, trianglePoints[0]);
+  list<TesselPoint*> *connectedClosestPoints = getNeighboursonCircleofConnectedPoints(out, connectedPoints, trianglePoints[0], x);
+  delete(connectedPoints);
   trianglePoints[1] = connectedClosestPoints->front();
   trianglePoints[2] = connectedClosestPoints->back();
@@ -2182,22 +2457,13 @@
 
 /** Gets all points connected to the provided point by triangulation lines.
- * Maps them down onto the plane designated by the axis \a *Point and \a *Reference. The center of all points
- * connected in the tesselation to \a *Point is mapped to spherical coordinates with the zero angle being given
- * by the mapped down \a *Reference. Hence, the biggest and the smallest angles are those of the two shanks of the
- * triangle we are looking for.
  *
  * @param *Point of which get all connected points
- * @param *Reference Vector to be checked whether it is an inner point
  *
- * @return list of the two points linked to the provided one and closest to the point to be checked,
- */
-list<TesselPoint*> * Tesselation::getCircleOfConnectedPoints(ofstream *out, TesselPoint* Point, Vector* Reference)
-{
-  list<TesselPoint*> connectedPoints;
-  map<double, TesselPoint*> anglesOfPoints;
-  map<double, TesselPoint*>::iterator runner;
-  list<TesselPoint*>::iterator listRunner;
+ * @return list of the all points linked to the provided one
+ */
+list<TesselPoint*> * Tesselation::getCircleOfConnectedPoints(ofstream *out, TesselPoint* Point)
+{
+  list<TesselPoint*> *connectedPoints = new list<TesselPoint*>;
   class BoundaryPointSet *ReferencePoint = NULL;
-  Vector center, PlaneNormal, currentPoint, OrthogonalVector, helper, AngleZero;
   TesselPoint* current;
   bool takePoint = false;
@@ -2231,8 +2497,5 @@
     if (takePoint) {
       *out << Verbose(3) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is taken into the circle." << endl;
-      connectedPoints.push_back(current);
-      currentPoint.CopyVector(current->node);
-      currentPoint.ProjectOntoPlane(&PlaneNormal);
-      center.AddVector(&currentPoint);
+      connectedPoints->push_back(current);
     }
 
@@ -2240,13 +2503,42 @@
   }
 
+  if (connectedPoints->size() == 0) { // if have not found any points
+    *out << Verbose(1) << "ERROR: We have not found any connected points to " << *Point<< "." << endl;
+    return NULL;
+  }
+  return connectedPoints;
+}
+
+/** Gets the two neighbouring points with respect to a reference line to the provided point.
+ * Maps them down onto the plane designated by the axis \a *Point and \a *Reference. The center of all points
+ * connected in the tesselation to \a *Point is mapped to spherical coordinates with the zero angle being given
+ * by the mapped down \a *Reference. Hence, the biggest and the smallest angles are those of the two shanks of the
+ * triangle we are looking for.
+ *
+ * @param *out output stream for debugging
+ * @param *connectedPoints list of connected points to the central \a *Point
+ * @param *Point of which get all connected points
+ * @param *Reference Vector to be checked whether it is an inner point
+ *
+ * @return list of the two points linked to the provided one and closest to the point to be checked,
+ */
+list<TesselPoint*> * Tesselation::getNeighboursonCircleofConnectedPoints(ofstream *out, list<TesselPoint*> *connectedPoints, TesselPoint* Point, Vector* Reference)
+{
+  map<double, TesselPoint*> anglesOfPoints;
+  map<double, TesselPoint*>::iterator runner;
+  ;
+  Vector center, PlaneNormal, OrthogonalVector, helper, AngleZero;
+
+  if (connectedPoints->size() == 0) { // if have not found any points
+    *out << Verbose(1) << "ERROR: We have not found any connected points to " << *Point<< "." << endl;
+    return NULL;
+  }
+
+  // calculate central point
+  for (list<TesselPoint*>::iterator TesselRunner = connectedPoints->begin(); TesselRunner != connectedPoints->end(); TesselRunner++)
+    center.AddVector((*TesselRunner)->node);
   //*out << "Summed vectors " << center << "; number of points " << connectedPoints.size()
   //  << "; scale factor " << 1.0/connectedPoints.size();
-
-  if (connectedPoints.size() == 0) { // if have not found any points
-    *out << Verbose(1) << "ERROR: We have not found any connected points to " << *Point<< "." << endl;
-    return NULL;
-  }
-  center.Scale(1.0/connectedPoints.size());
-
+  center.Scale(1.0/connectedPoints->size());
   *out << Verbose(4) << "INFO: Calculated center of all circle points is " << center << "." << endl;
 
@@ -2264,8 +2556,7 @@
   OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
   *out << Verbose(4) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;
-  
+
   // go through all connected points and calculate angle
-  listRunner = connectedPoints.begin();
-  while (listRunner != connectedPoints.end()) {
+  for (list<TesselPoint*>::iterator listRunner = connectedPoints->begin(); listRunner != connectedPoints->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
     helper.SubtractVector(Point->node);
@@ -2274,5 +2565,4 @@
     *out << Verbose(2) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
     anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
-    listRunner++;
   }
 
@@ -2289,4 +2579,90 @@
   return result;
 }
+
+/** Removes a boundary point from the envelope while keeping it closed.
+ * We create new triangles and remove the old ones connected to the point.
+ * \param *out output stream for debugging
+ * \param *point point to be removed
+ * \return volume added to the volume inside the tesselated surface by the removal
+ */
+double Tesselation::RemovePointFromTesselatedSurface(ofstream *out, class BoundaryPointSet *point) {
+  class BoundaryLineSet *line = NULL;
+  class BoundaryTriangleSet *triangle = NULL;
+  Vector OldPoint, TetraederVector[3];
+  double volume = 0;
+  int *numbers = NULL;
+  int count = 0;
+  int i;
+
+  // copy old location for the volume
+  OldPoint.CopyVector(point->node->node);
+
+  // get list of connected points
+  if (point->lines.empty()) {
+    *out << Verbose(1) << "ERROR: Cannot remove the point " << *point << ", it's connected to no lines!" << endl;
+    return 0.;
+  }
+  list<TesselPoint*> *CircleofPoints = getCircleOfConnectedPoints(out, point->node);
+
+  // remove all triangles
+  for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
+    count+=LineRunner->second->triangles.size();
+  numbers = new int[count];
+  i=0;
+  for (LineMap::iterator LineRunner = point->lines.begin(); (point != NULL) && (LineRunner != point->lines.end()); LineRunner++) {
+    line = LineRunner->second;
+    for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
+      triangle = TriangleRunner->second;
+      *out << Verbose(2) << "Erasing triangle " << *triangle << "." << endl;
+      numbers[i++] = triangle->Nr;
+      RemoveTesselationTriangle(triangle);
+      triangle = NULL;
+    }
+  }
+  *out << Verbose(1) << i << " triangles were removed." << endl;
+
+  // re-create all triangles by going through connected points list
+  list<TesselPoint*>::iterator CircleRunner = CircleofPoints->begin();
+  list<TesselPoint*>::iterator OtherCircleRunner = CircleofPoints->begin();
+  class TesselPoint *CentralNode = *CircleRunner;
+  // advance two with CircleRunner and one with OtherCircleRunner
+  CircleRunner++;
+  CircleRunner++;
+  OtherCircleRunner++;
+  i=0;
+  cout << Verbose(2) << "INFO: CentralNode is " << *CentralNode << "." << endl;
+  for (; (OtherCircleRunner != CircleofPoints->end()) && (CircleRunner != CircleofPoints->end()); (CircleRunner++), (OtherCircleRunner++)) {
+    cout << Verbose(3) << "INFO: CircleRunner's node is " << **CircleRunner << "." << endl;
+    cout << Verbose(3) << "INFO: OtherCircleRunner's node is " << **OtherCircleRunner << "." << endl;
+    *out << Verbose(4) << "Adding new triangle points."<< endl;
+    AddTesselationPoint(CentralNode, 0);
+    AddTesselationPoint(*OtherCircleRunner, 1);
+    AddTesselationPoint(*CircleRunner, 2);
+    *out << Verbose(4) << "Adding new triangle lines."<< endl;
+    AddTesselationLine(TPS[0], TPS[1], 0);
+    AddTesselationLine(TPS[0], TPS[2], 1);
+    AddTesselationLine(TPS[1], TPS[2], 2);
+    BTS = new class BoundaryTriangleSet(BLS, numbers[i]);
+    TrianglesOnBoundary.insert(TrianglePair(numbers[i], BTS));
+    *out << Verbose(4) << "Created triangle " << *BTS << "." << endl;
+    // calculate volume summand as a general tetraeder
+    for (int j=0;j<3;j++) {
+      TetraederVector[j].CopyVector(TPS[j]->node->node);
+      TetraederVector[j].SubtractVector(&OldPoint);
+    }
+    OldPoint.CopyVector(&TetraederVector[0]);
+    OldPoint.VectorProduct(&TetraederVector[1]);
+    volume += 1./6. * fabs(OldPoint.ScalarProduct(&TetraederVector[2]));
+    // advance number
+    i++;
+    if (i >= count)
+      *out << Verbose(2) << "WARNING: Maximum of numbers reached!" << endl;
+  }
+  *out << Verbose(1) << i << " triangles were created." << endl;
+
+  delete[](numbers);
+
+  return volume;
+};
 
 /** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
Index: src/tesselation.hpp
===================================================================
--- src/tesselation.hpp	(revision 86381b38132f1af7b1e83f6a7c271dc5decbac5c)
+++ src/tesselation.hpp	(revision 16d86622dc2c759142b589c64308c884f9339eb6)
@@ -78,4 +78,5 @@
     int LinesCount;
     TesselPoint *node;
+    double value;
     int Nr;
 };
@@ -190,9 +191,12 @@
     virtual ~Tesselation();
 
-    bool AddPoint(TesselPoint *Walker, int n);
-    void AddTrianglePoint(TesselPoint* Candidate, int n);
-    void AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
-    void AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
-    void AddTriangle();
+    void AddTesselationPoint(TesselPoint* Candidate, int n);
+    void AddTesselationLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
+    void AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
+    void AddTesselationTriangle();
+    void RemoveTesselationTriangle(class BoundaryTriangleSet *triangle);
+    void RemoveTesselationLine(class BoundaryLineSet *line);
+    void RemoveTesselationPoint(class BoundaryPointSet *point);
+
     bool IsInside(Vector *pointer);
     class BoundaryPointSet *GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2);
@@ -209,7 +213,12 @@
     void GuessStartingTriangle(ofstream *out);
     bool InsertStraddlingPoints(ofstream *out, PointCloud *cloud, LinkedCell *LC);
-    bool CorrectConcaveBaselines(ofstream *out);
-
-    list<TesselPoint*> * getCircleOfConnectedPoints(ofstream *out, TesselPoint* Point, Vector* Reference);
+    double RemovePointFromTesselatedSurface(ofstream *out, class BoundaryPointSet *point);
+    bool FlipBaseline(ofstream *out, class BoundaryLineSet *Base);
+    bool PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base);
+    class BoundaryPointSet *IsConvexRectangle(ofstream *out, class BoundaryLineSet *Base);
+
+
+    list<TesselPoint*> * getCircleOfConnectedPoints(ofstream *out, TesselPoint* Point);
+    list<TesselPoint*> * getNeighboursonCircleofConnectedPoints(ofstream *out, list<TesselPoint*> *connectedPoints, TesselPoint* Point, Vector* Reference);
     list<BoundaryTriangleSet*> *FindTriangles(TesselPoint* Points[3]);
     list<BoundaryTriangleSet*> * FindClosestTrianglesToPoint(ofstream *out, Vector *x, LinkedCell* LC);
@@ -217,4 +226,5 @@
     bool IsInnerPoint(ofstream *out, Vector Point, LinkedCell* LC);
     bool IsInnerPoint(ofstream *out, TesselPoint *Point, LinkedCell* LC);
+    bool AddBoundaryPoint(TesselPoint *Walker, int n);
 
     PointMap PointsOnBoundary;
@@ -251,4 +261,5 @@
 TesselPoint* findSecondClosestPoint(const Vector*, LinkedCell*);
 double getAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector);
+Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase);
 
 #endif /* TESSELATION_HPP_ */
