source: molecuilder/src/tesselation.cpp@ d6b011

/*
 * tesselation.cpp
 *
 *  Created on: Aug 3, 2009
 *      Author: heber
 */

#include "tesselation.hpp"

// ======================================== 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)
{
  node = Walker;
  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()
{
  cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
  if (!lines.empty())
    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)
{
  cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."
      << endl;
  if (line->endpoints[0] == this)
    {
      lines.insert(LinePair(line->endpoints[1]->Nr, line));
    }
  else
    {
      lines.insert(LinePair(line->endpoints[0]->Nr, line));
    }
  LinesCount++;
};

/** output operator for BoundaryPointSet.
 * \param &ost output stream
 * \param &a boundary point
 */
ostream & operator <<(ostream &ost, BoundaryPointSet &a)
{
  ost << "[" << a.Nr << "|" << a.node->Name << "]";
  return ost;
}
;

// ======================================== Lines on Boundary =================================

/** Constructor of BoundaryLineSet.
 */
BoundaryLineSet::BoundaryLineSet()
{
  for (int i = 0; i < 2; i++)
    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)
{
  // set number
  Nr = number;
  // set endpoints in ascending order
  SetEndpointsOrdered(endpoints, Point[0], Point[1]);
  // add this line to the hash maps of both endpoints
  Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
  Point[1]->AddLine(this); //
  // 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];

  // 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++) {
    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;
        }
      }
    }
  }
  if (!triangles.empty())
    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;
};

/** 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.
 * \param *line other line to test
 * \return true - common endpoint present, false - not connected
 */
bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
{
  if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
    return true;
  else
    return false;
};

/** Checks whether the adjacent triangles of a baseline are convex or not.
 * We sum the two angles of each normal vector with a ficticious normnal vector from this baselinbe pointing outwards.
 * If greater/equal M_PI than we are convex.
 * \param *out output stream for debugging
 * \return true - triangles are convex, false - concave or less than two triangles connected
 */
bool BoundaryLineSet::CheckConvexityCriterion(ofstream *out)
{
  Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck;
  // get the two triangles
  if (triangles.size() != 2) {
    *out << Verbose(1) << "ERROR: Baseline " << *this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
    return false;
  }
  // check normal vectors
  // have a normal vector on the base line pointing outwards
  *out << Verbose(3) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl;
  BaseLineCenter.CopyVector(endpoints[0]->node->node);
  BaseLineCenter.AddVector(endpoints[1]->node->node);
  BaseLineCenter.Scale(1./2.);
  BaseLine.CopyVector(endpoints[0]->node->node);
  BaseLine.SubtractVector(endpoints[1]->node->node);
  *out << Verbose(3) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl;

  BaseLineNormal.Zero();
  NormalCheck.Zero();
  double sign = -1.;
  int i=0;
  class BoundaryPointSet *node = NULL;
  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
    *out << Verbose(3) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
    NormalCheck.AddVector(&runner->second->NormalVector);
    NormalCheck.Scale(sign);
    sign = -sign;
    BaseLineNormal.SubtractVector(&runner->second->NormalVector);   // we subtract as BaseLineNormal has to point inward in direction of [pi,2pi]
    node = runner->second->GetThirdEndpoint(this);
    if (node != NULL) {
      *out << Verbose(3) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;
      helper[i].CopyVector(node->node->node);
      helper[i].SubtractVector(&BaseLineCenter);
      helper[i].MakeNormalVector(&BaseLine);  // we want to compare the triangle's heights' angles!
      *out << Verbose(4) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
      i++;
    } else {
      *out << Verbose(2) << "WARNING: I cannot find third node in triangle, something's wrong." << endl;
      return true;
    }
  }
  *out << Verbose(3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
  if (NormalCheck.NormSquared() < MYEPSILON) {
    *out << Verbose(3) << "Normalvectors of both triangles are the same: convex." << endl;
    return true;
  }
  double angle = GetAngle(helper[0], helper[1], BaseLineNormal);
  if ((angle - M_PI) > -MYEPSILON)
    return true;
  else
    return false;
}

/** Checks whether point is any of the two endpoints this line contains.
 * \param *point point to test
 * \return true - point is of the line, false - is not
 */
bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
{
  for(int i=0;i<2;i++)
    if (point == endpoints[i])
      return true;
  return false;
};

/** Returns other endpoint of the line.
 * \param *point other endpoint
 * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
 */
class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point)
{
  if (endpoints[0] == point)
    return endpoints[1];
  else if (endpoints[1] == point)
    return endpoints[0];
  else
    return NULL;
};

/** 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()
{
  for (int i = 0; i < 3; i++)
    {
      endpoints[i] = NULL;
      lines[i] = NULL;
    }
  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)
{
  // set number
  Nr = number;
  // set lines
  cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
  for (int i = 0; i < 3; i++)
    {
      lines[i] = line[i];
      lines[i]->AddTriangle(this);
    }
  // get ascending order of endpoints
  map<int, class BoundaryPointSet *> OrderMap;
  for (int i = 0; i < 3; i++)
    // for all three lines
    for (int j = 0; j < 2; j++)
      { // for both endpoints
        OrderMap.insert(pair<int, class BoundaryPointSet *> (
            line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
        // and we don't care whether insertion fails
      }
  // set endpoints
  int Counter = 0;
  cout << Verbose(6) << " with end points ";
  for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner
      != OrderMap.end(); runner++)
    {
      endpoints[Counter] = runner->second;
      cout << " " << *endpoints[Counter];
      Counter++;
    }
  if (Counter < 3)
    {
      cerr << "ERROR! We have a triangle with only two distinct endpoints!"
          << endl;
      //exit(1);
    }
  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++) {
    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.
 * Is made unique by comparison with \a OtherVector to point in the other direction.
 * \param &OtherVector direction vector to make normal vector unique.
 */
void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
{
  // get normal vector
  NormalVector.MakeNormalVector(endpoints[0]->node->node, endpoints[1]->node->node, endpoints[2]->node->node);

  // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
  if (NormalVector.Projection(&OtherVector) > 0)
    NormalVector.Scale(-1.);
};

/** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
 * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
 * 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: if it's really outside, then for any endpoint of the triangle and it's opposite
 * base line, the intersection between the line from endpoint to intersection and the base line will have a Vector::NormSquared()
 * smaller than the first line.
 * \param *out output stream for debugging
 * \param *MolCenter offset vector of line
 * \param *x second endpoint of line, minus \a *MolCenter is directional vector of line
 * \param *Intersection intersection on plane on return
 * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
 */
bool BoundaryTriangleSet::GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection)
{
  Vector CrossPoint;
  Vector helper;

  if (!Intersection->GetIntersectionWithPlane(out, &NormalVector, endpoints[0]->node->node, MolCenter, x)) {
    *out << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl;
    return false;
  }

  // Calculate cross point between one baseline and the line from the third endpoint to intersection
  int i=0;
  do {
    if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(out, endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection, &NormalVector)) {
      helper.CopyVector(endpoints[(i+1)%3]->node->node);
      helper.SubtractVector(endpoints[i%3]->node->node);
    } else 
      i++;
    if (i>2)
      break;
  } while (CrossPoint.NormSquared() < MYEPSILON);
  if (i==3) {
    *out << Verbose(1) << "ERROR: Could not find any cross points, something's utterly wrong here!" << endl;
    exit(255);
  }
  CrossPoint.SubtractVector(endpoints[i%3]->node->node);

  // check whether intersection is inside or not by comparing length of intersection and length of cross point
  if ((CrossPoint.NormSquared() - helper.NormSquared()) > -MYEPSILON) { // inside
    return true;
  } else { // outside!
    Intersection->Zero();
    return false;
  }
};

/** 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(class BoundaryLineSet *line)
{
  for(int i=0;i<3;i++)
    if (line == lines[i])
      return true;
  return false;
};

/** Checks whether point is any of the three endpoints this triangle contains.
 * \param *point point to test
 * \return true - point is of the triangle, false - is not
 */
bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
{
  for(int i=0;i<3;i++)
    if (point == endpoints[i])
      return true;
  return false;
};

/** Checks whether three given \a *Points coincide with triangle's endpoints.
 * \param *Points[3] pointer to BoundaryPointSet
 * \return true - is the very triangle, false - is not
 */
bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3])
{
  return (((endpoints[0] == Points[0])
            || (endpoints[0] == Points[1])
            || (endpoints[0] == Points[2])
          ) && (
            (endpoints[1] == Points[0])
            || (endpoints[1] == Points[1])
            || (endpoints[1] == Points[2])
          ) && (
            (endpoints[2] == Points[0])
            || (endpoints[2] == Points[1])
            || (endpoints[2] == Points[2])

          ));
};

/** Returns the endpoint which is not contained in the given \a *line.
 * \param *line baseline defining two endpoints
 * \return pointer third endpoint or NULL if line does not belong to triangle.
 */
class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(class BoundaryLineSet *line)
{
  // sanity check
  if (!ContainsBoundaryLine(line))
    return NULL;
  for(int i=0;i<3;i++)
    if (!line->ContainsBoundaryPoint(endpoints[i]))
      return endpoints[i];
  // actually, that' impossible :)
  return NULL;
};

/** Calculates the center point of the triangle.
 * Is third of the sum of all endpoints.
 * \param *center central point on return.
 */
void BoundaryTriangleSet::GetCenter(Vector *center)
{
  center->Zero();
  for(int i=0;i<3;i++)
    center->AddVector(endpoints[i]->node->node);
  center->Scale(1./3.);
}

/** 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 ===========================================

/** Constructor of class TesselPoint.
 */
TesselPoint::TesselPoint()
{
  node = NULL;
  nr = -1;
  Name =  NULL;
};

/** Destructor for class TesselPoint.
 */
TesselPoint::~TesselPoint()
{
  Free((void **)&Name, "TesselPoint::~TesselPoint: *Name");
};

/** Prints LCNode to screen.
 */
ostream & operator << (ostream &ost, const TesselPoint &a)
{
  ost << "[" << (a.Name) << "|" << &a << "]";
  return ost;
};

/** Prints LCNode to screen.
 */
ostream & TesselPoint::operator << (ostream &ost)
{
  ost << "[" << (Name) << "|" << this << "]";
  return ost;
};


// =========================================================== class POINTCLOUD ============================================

/** Constructor of class PointCloud.
 */
PointCloud::PointCloud()
{

};

/** Destructor for class PointCloud.
 */
PointCloud::~PointCloud()
{

};

// ============================ CandidateForTesselation =============================

/** Constructor of class CandidateForTesselation.
 */
CandidateForTesselation::CandidateForTesselation(TesselPoint *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) {
  point = candidate;
  BaseLine = line;
  OptCenter.CopyVector(&OptCandidateCenter);
  OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
};

/** Destructor for class CandidateForTesselation.
 */
CandidateForTesselation::~CandidateForTesselation() {
  point = NULL;
  BaseLine = NULL;
};

// =========================================================== class TESSELATION ===========================================

/** Constructor of class Tesselation.
 */
Tesselation::Tesselation()
{
  PointsOnBoundaryCount = 0;
  LinesOnBoundaryCount = 0;
  TrianglesOnBoundaryCount = 0;
  InternalPointer = PointsOnBoundary.begin();
}
;

/** Destructor of class Tesselation.
 * We have to free all points, lines and triangles.
 */
Tesselation::~Tesselation()
{
  cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
    if (runner->second != NULL) {
      delete (runner->second);
      runner->second = NULL;
    } else
      cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
  }
}
;

/** PointCloud implementation of GetCenter
 * Uses PointsOnBoundary and STL stuff.
 */   
Vector * Tesselation::GetCenter(ofstream *out)
{
  Vector *Center = new Vector(0.,0.,0.);
  int num=0;
  for (GoToFirst(); (!IsEnd()); GoToNext()) {
    Center->AddVector(GetPoint()->node);
    num++;
  }
  Center->Scale(1./num);
  return Center;
};

/** PointCloud implementation of GoPoint
 * Uses PointsOnBoundary and STL stuff.
 */   
TesselPoint * Tesselation::GetPoint()
{
  return (InternalPointer->second->node);
};

/** PointCloud implementation of GetTerminalPoint.
 * Uses PointsOnBoundary and STL stuff.
 */   
TesselPoint * Tesselation::GetTerminalPoint()
{
  PointMap::iterator Runner = PointsOnBoundary.end();
  Runner--;
  return (Runner->second->node);
};

/** PointCloud implementation of GoToNext.
 * Uses PointsOnBoundary and STL stuff.
 */   
void Tesselation::GoToNext()
{
  if (InternalPointer != PointsOnBoundary.end())
    InternalPointer++;
};

/** PointCloud implementation of GoToPrevious.
 * Uses PointsOnBoundary and STL stuff.
 */   
void Tesselation::GoToPrevious()
{
  if (InternalPointer != PointsOnBoundary.begin())
    InternalPointer--;
};

/** PointCloud implementation of GoToFirst.
 * Uses PointsOnBoundary and STL stuff.
 */   
void Tesselation::GoToFirst()
{
  InternalPointer = PointsOnBoundary.begin();
};

/** PointCloud implementation of GoToLast.
 * Uses PointsOnBoundary and STL stuff.
 */   
void Tesselation::GoToLast()
{
  InternalPointer = PointsOnBoundary.end();
  InternalPointer--;
};

/** PointCloud implementation of IsEmpty.
 * Uses PointsOnBoundary and STL stuff.
 */   
bool Tesselation::IsEmpty()
{
  return (PointsOnBoundary.empty());
};

/** PointCloud implementation of IsLast.
 * Uses PointsOnBoundary and STL stuff.
 */   
bool Tesselation::IsEnd()
{
  return (InternalPointer == PointsOnBoundary.end());
};


/** Gueses first starting triangle of the convex envelope.
 * We guess the starting triangle by taking the smallest distance between two points and looking for a fitting third.
 * \param *out output stream for debugging
 * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
 */
void
Tesselation::GuessStartingTriangle(ofstream *out)
{
  // 4b. create a starting triangle
  // 4b1. create all distances
  DistanceMultiMap DistanceMMap;
  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->node->DistanceSquared(B->second->node->node);
              distance = tmp * tmp;
              tmp = A->second->node->node->DistanceSquared(C->second->node->node);
              distance += tmp * tmp;
              tmp = B->second->node->node->DistanceSquared(C->second->node->node);
              distance += tmp * tmp;
              DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
            }
        }
    }
  //    // listing distances
  //    *out << Verbose(1) << "Listing DistanceMMap:";
  //    for(DistanceMultiMap::iterator runner = DistanceMMap.begin(); runner != DistanceMMap.end(); runner++) {
  //      *out << " " << runner->first << "(" << *runner->second.first->second << ", " << *runner->second.second->second << ")";
  //    }
  //    *out << endl;
  // 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();
  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->node,
          baseline->second.first->second->node->node,
          baseline->second.second->second->node->node);
      *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->node);
          TrialVector.SubtractVector(A->second->node->node);
          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 << " leaves "
                  << 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->node->DistanceSquared(A->second->node->node);
          int innerpoint = 0;
          if ((tmp < A->second->node->node->DistanceSquared(
              baseline->second.first->second->node->node)) && (tmp
              < A->second->node->node->DistanceSquared(
                  baseline->second.second->second->node->node)))
            innerpoint++;
          tmp = checker->second->node->node->DistanceSquared(
              baseline->second.first->second->node->node);
          if ((tmp < baseline->second.first->second->node->node->DistanceSquared(
              A->second->node->node)) && (tmp
              < baseline->second.first->second->node->node->DistanceSquared(
                  baseline->second.second->second->node->node)))
            innerpoint++;
          tmp = checker->second->node->node->DistanceSquared(
              baseline->second.second->second->node->node);
          if ((tmp < baseline->second.second->second->node->node->DistanceSquared(
              baseline->second.first->second->node->node)) && (tmp
              < baseline->second.second->second->node->node->DistanceSquared(
                  A->second->node->node)))
            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);
    }
}
;

/** Tesselates the convex envelope of a cluster from a single starting triangle.
 * The starting triangle is made out of three baselines. Each line in the final tesselated cluster may belong to at most
 * 2 triangles. Hence, we go through all current lines:
 * -# if the lines contains to only one triangle
 * -# We search all points in the boundary
 *    -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
 *       baseline in triangle plane pointing out of the triangle and normal vector of new triangle)
 *    -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors)
 *    -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
 * \param *out output stream for debugging
 * \param *configuration for IsAngstroem
 * \param *cloud cluster of points
 */
void Tesselation::TesselateOnBoundary(ofstream *out, PointCloud *cloud)
{
  bool flag;
  PointMap::iterator winner;
  class BoundaryPointSet *peak = NULL;
  double SmallestAngle, TempAngle;
  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector, helper, PropagationVector, *Center = NULL;
  LineMap::iterator LineChecker[2];

  Center = cloud->GetCenter(out);
  // create a first tesselation with the given BoundaryPoints
  do {
    flag = false;
    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
      if (baseline->second->triangles.size() == 1) {
        // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
        SmallestAngle = M_PI;

        // get peak point with respect to this base line's only triangle
        BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
        *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl;
        for (int i = 0; i < 3; i++)
          if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1]))
            peak = BTS->endpoints[i];
        *out << Verbose(3) << " and has peak " << *peak << "." << endl;

        // prepare some auxiliary vectors
        Vector BaseLineCenter, BaseLine;
        BaseLineCenter.CopyVector(baseline->second->endpoints[0]->node->node);
        BaseLineCenter.AddVector(baseline->second->endpoints[1]->node->node);
        BaseLineCenter.Scale(1. / 2.); // points now to center of base line
        BaseLine.CopyVector(baseline->second->endpoints[0]->node->node);
        BaseLine.SubtractVector(baseline->second->endpoints[1]->node->node);

        // offset to center of triangle
        CenterVector.Zero();
        for (int i = 0; i < 3; i++)
          CenterVector.AddVector(BTS->endpoints[i]->node->node);
        CenterVector.Scale(1. / 3.);
        *out << Verbose(4) << "CenterVector of base triangle is " << CenterVector << endl;

        // normal vector of triangle
        NormalVector.CopyVector(Center);
        NormalVector.SubtractVector(&CenterVector);
        BTS->GetNormalVector(NormalVector);
        NormalVector.CopyVector(&BTS->NormalVector);
        *out << Verbose(4) << "NormalVector of base triangle is " << NormalVector << endl;

        // vector in propagation direction (out of triangle)
        // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
        PropagationVector.MakeNormalVector(&BaseLine, &NormalVector);
        TempVector.CopyVector(&CenterVector);
        TempVector.SubtractVector(baseline->second->endpoints[0]->node->node); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
        //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
        if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
          PropagationVector.Scale(-1.);
        *out << Verbose(4) << "PropagationVector of base triangle is " << PropagationVector << endl;
        winner = PointsOnBoundary.end();

        // loop over all points and calculate angle between normal vector of new and present triangle
        for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) {
          if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints
            *out << Verbose(3) << "Target point is " << *(target->second) << ":" << endl;

            // first check direction, so that triangles don't intersect
            VirtualNormalVector.CopyVector(target->second->node->node);
            VirtualNormalVector.SubtractVector(&BaseLineCenter); // points from center of base line to target
            VirtualNormalVector.ProjectOntoPlane(&NormalVector);
            TempAngle = VirtualNormalVector.Angle(&PropagationVector);
            *out << Verbose(4) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl;
            if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees)
              *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;
              continue;
            } else
              *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl;

            // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle)
            LineChecker[0] = baseline->second->endpoints[0]->lines.find(target->first);
            LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first);
            if (((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[0]->second->triangles.size() == 2))) {
              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->triangles.size() << " triangles." << endl;
              continue;
            }
            if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->triangles.size() == 2))) {
              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->triangles.size() << " triangles." << endl;
              continue;
            }

            // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
            if ((((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak)))) {
              *out << Verbose(4) << "Current target is peak!" << endl;
              continue;
            }

            // check for linear dependence
            TempVector.CopyVector(baseline->second->endpoints[0]->node->node);
            TempVector.SubtractVector(target->second->node->node);
            helper.CopyVector(baseline->second->endpoints[1]->node->node);
            helper.SubtractVector(target->second->node->node);
            helper.ProjectOntoPlane(&TempVector);
            if (fabs(helper.NormSquared()) < MYEPSILON) {
              *out << Verbose(4) << "Chosen set of vectors is linear dependent." << endl;
              continue;
            }

            // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle
            flag = true;
            VirtualNormalVector.MakeNormalVector(baseline->second->endpoints[0]->node->node, baseline->second->endpoints[1]->node->node, target->second->node->node);
            TempVector.CopyVector(baseline->second->endpoints[0]->node->node);
            TempVector.AddVector(baseline->second->endpoints[1]->node->node);
            TempVector.AddVector(target->second->node->node);
            TempVector.Scale(1./3.);
            TempVector.SubtractVector(Center);
            // make it always point outward
            if (VirtualNormalVector.Projection(&TempVector) < 0)
              VirtualNormalVector.Scale(-1.);
            // calculate angle
            TempAngle = NormalVector.Angle(&VirtualNormalVector);
            *out << Verbose(4) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl;
            if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner
              SmallestAngle = TempAngle;
              winner = target;
              *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
            } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
              // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
              helper.CopyVector(target->second->node->node);
              helper.SubtractVector(&BaseLineCenter);
              helper.ProjectOntoPlane(&BaseLine);
              // ...the one with the smaller angle is the better candidate
              TempVector.CopyVector(target->second->node->node);
              TempVector.SubtractVector(&BaseLineCenter);
              TempVector.ProjectOntoPlane(&VirtualNormalVector);
              TempAngle = TempVector.Angle(&helper);
              TempVector.CopyVector(winner->second->node->node);
              TempVector.SubtractVector(&BaseLineCenter);
              TempVector.ProjectOntoPlane(&VirtualNormalVector);
              if (TempAngle < TempVector.Angle(&helper)) {
                TempAngle = NormalVector.Angle(&VirtualNormalVector);
                SmallestAngle = TempAngle;
                winner = target;
                *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl;
              } else
                *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl;
            } else
              *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
          }
        } // end of loop over all boundary points

        // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
        if (winner != PointsOnBoundary.end()) {
          *out << Verbose(2) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl;
          // create the lins of not yet present
          BLS[0] = baseline->second;
          // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
          LineChecker[0] = baseline->second->endpoints[0]->lines.find(winner->first);
          LineChecker[1] = baseline->second->endpoints[1]->lines.find(winner->first);
          if (LineChecker[0] == baseline->second->endpoints[0]->lines.end()) { // create
            BPS[0] = baseline->second->endpoints[0];
            BPS[1] = winner->second;
            BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
            LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[1]));
            LinesOnBoundaryCount++;
          } else
            BLS[1] = LineChecker[0]->second;
          if (LineChecker[1] == baseline->second->endpoints[1]->lines.end()) { // create
            BPS[0] = baseline->second->endpoints[1];
            BPS[1] = winner->second;
            BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
            LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[2]));
            LinesOnBoundaryCount++;
          } else
            BLS[2] = LineChecker[1]->second;
          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
          BTS->GetCenter(&helper);
          helper.SubtractVector(Center);
          helper.Scale(-1);
          BTS->GetNormalVector(helper);
          TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
          TrianglesOnBoundaryCount++;
        } else {
          *out << Verbose(1) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl;
        }

        // 5d. If the set of lines is not yet empty, go to 5. and continue
      } else
        *out << Verbose(2) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->triangles.size() << "." << endl;
  } while (flag);

  // exit
  delete(Center);
};

/** Inserts all points outside of the tesselated surface into it by adding new triangles.
 * \param *out output stream for debugging
 * \param *cloud cluster of points
 * \param *LC LinkedCell structure to find nearest point quickly
 * \return true - all straddling points insert, false - something went wrong
 */
bool Tesselation::InsertStraddlingPoints(ofstream *out, PointCloud *cloud, LinkedCell *LC)
{
  Vector Intersection, Normal;
  TesselPoint *Walker = NULL;
  Vector *Center = cloud->GetCenter(out);
  list<BoundaryTriangleSet*> *triangles = NULL;

  *out << Verbose(1) << "Begin of InsertStraddlingPoints" << endl;

  cloud->GoToFirst();
  while (!cloud->IsEnd()) {  // we only have to go once through all points, as boundary can become only bigger
    LinkedCell BoundaryPoints(this, 5.);
    Walker = cloud->GetPoint();
    *out << Verbose(2) << "Current point is " << *Walker << "." << endl;
    // get the next triangle
    triangles = FindClosestTrianglesToPoint(out, Walker->node, &BoundaryPoints);
    if (triangles == NULL) {
      *out << Verbose(1) << "No triangles found, probably a tesselation point itself." << endl;
      cloud->GoToNext();
      continue;
    } else {
      BTS = triangles->front();
    }
    *out << Verbose(2) << "Closest triangle is " << *BTS << "." << endl;
    // get the intersection point
    if (BTS->GetIntersectionInsideTriangle(out, Center, Walker->node, &Intersection)) {
      *out << Verbose(2) << "We have an intersection at " << Intersection << "." << endl;
      // we have the intersection, check whether in- or outside of boundary
      if ((Center->DistanceSquared(Walker->node) - Center->DistanceSquared(&Intersection)) < -MYEPSILON) {
        // inside, next!
        *out << Verbose(2) << *Walker << " is inside wrt triangle " << *BTS << "." << endl;
      } else {
        // outside!
        *out << Verbose(2) << *Walker << " is outside wrt triangle " << *BTS << "." << endl;
        class BoundaryLineSet *OldLines[3], *NewLines[3];
        class BoundaryPointSet *OldPoints[3], *NewPoint;
        // store the three old lines and old points
        for (int i=0;i<3;i++) {
          OldLines[i] = BTS->lines[i];
          OldPoints[i] = BTS->endpoints[i];
        }
        Normal.CopyVector(&BTS->NormalVector);
        // add Walker to boundary points
        *out << Verbose(2) << "Adding " << *Walker << " to BoundaryPoints." << endl;
        if (AddBoundaryPoint(Walker,0))
          NewPoint = BPS[0];
        else
          continue;
        // remove triangle
        *out << Verbose(2) << "Erasing triangle " << *BTS << "." << endl;
        TrianglesOnBoundary.erase(BTS->Nr);
        delete(BTS);
        // create three new boundary lines
        for (int i=0;i<3;i++) {
          BPS[0] = NewPoint;
          BPS[1] = OldPoints[i];
          NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
          *out << Verbose(3) << "Creating new line " << *NewLines[i] << "." << endl;
          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one
          LinesOnBoundaryCount++;
        }
        // create three new triangle with new point
        for (int i=0;i<3;i++) { // find all baselines
          BLS[0] = OldLines[i];
          int n = 1;
          for (int j=0;j<3;j++) {
            if (NewLines[j]->IsConnectedTo(BLS[0])) {
              if (n>2) {
                *out << Verbose(1) << "ERROR: " << BLS[0] << " connects to all of the new lines?!" << endl;
                return false;
              } else
                BLS[n++] = NewLines[j];
            }
          }
          // create the triangle
          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
          Normal.Scale(-1.);
          BTS->GetNormalVector(Normal);
          Normal.Scale(-1.);
          *out << Verbose(2) << "Created new triangle " << *BTS << "." << endl;
          TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
          TrianglesOnBoundaryCount++;
        }
      }
    } else { // something is wrong with FindClosestTriangleToPoint!
      *out << Verbose(1) << "ERROR: The closest triangle did not produce an intersection!" << endl;
      return false;
    }
    cloud->GoToNext();
  }

  // exit
  delete(Center);
  *out << Verbose(1) << "End of InsertStraddlingPoints" << endl;
  return true;
};

/** Adds a point to the tesselation::PointsOnBoundary list.
 * \param *Walker point to add
 * \param n TesselStruct::BPS index to put pointer into
 * \return true - new point was added, false - point already present
 */
bool
Tesselation::AddBoundaryPoint(TesselPoint *Walker, int n)
{
  PointTestPair InsertUnique;
  BPS[n] = new class BoundaryPointSet(Walker);
  InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[n]));
  if (InsertUnique.second) { // if new point was not present before, increase counter
    PointsOnBoundaryCount++;
    return true;
  } else {
    delete(BPS[n]);
    BPS[n] = InsertUnique.first->second;
    return false;
  }
}
;

/** Adds point to Tesselation::PointsOnBoundary if not yet present.
 * Tesselation::TPS is set to either this new BoundaryPointSet or to the existing one of not unique.
 * @param Candidate point to add
 * @param n index for this point in Tesselation::TPS array
 */
void
Tesselation::AddTesselationPoint(TesselPoint* Candidate, int n)
{
  PointTestPair InsertUnique;
  TPS[n] = new class BoundaryPointSet(Candidate);
  InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
  if (InsertUnique.second) { // if new point was not present before, increase counter
    PointsOnBoundaryCount++;
  } else {
    delete TPS[n];
    cout << Verbose(3) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
    TPS[n] = (InsertUnique.first)->second;
  }
}
;

/** Function tries to add line from current Points in BPS to BoundaryLineSet.
 * If successful it raises the line count and inserts the new line into the BLS,
 * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one.
 * @param *a first endpoint
 * @param *b second endpoint
 * @param n index of Tesselation::BLS giving the line with both endpoints
 */
void Tesselation::AddTesselationLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) {
  bool insertNewLine = true;

  if (a->lines.find(b->node->nr) != a->lines.end()) {
    LineMap::iterator FindLine;
    pair<LineMap::iterator,LineMap::iterator> FindPair;
    FindPair = a->lines.equal_range(b->node->nr);

    for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
      // If there is a line with less than two attached triangles, we don't need a new line.
      if (FindLine->second->triangles.size() < 2) {
        insertNewLine = false;
        cout << Verbose(3) << "Using existing line " << *FindLine->second << endl;

        BPS[0] = FindLine->second->endpoints[0];
        BPS[1] = FindLine->second->endpoints[1];
        BLS[n] = FindLine->second;

        break;
      }
    }
  }

  if (insertNewLine) {
    AlwaysAddTesselationTriangleLine(a, b, n);
  }
}
;

/**
 * Adds lines from each of the current points in the BPS to BoundaryLineSet.
 * Raises the line count and inserts the new line into the BLS.
 *
 * @param *a first endpoint
 * @param *b second endpoint
 * @param n index of Tesselation::BLS giving the line with both endpoints
 */
void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
{
  cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
  BPS[0] = a;
  BPS[1] = b;
  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);  // this also adds the line to the local maps
  // add line to global map
  LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
  // increase counter
  LinesOnBoundaryCount++;
};

/** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm).
 * Furthermore it adds the triangle to all of its lines, in order to recognize those which are saturated later.
 */
void Tesselation::AddTesselationTriangle()
{
  cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;

  // add triangle to global map
  TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
  TrianglesOnBoundaryCount++;

  // 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.
 * Searches for the points in Tesselation::PointsOnBoundary and checks their
 * lines. If any of the three edges already has two triangles attached, false is
 * returned.
 * \param *out output stream for debugging
 * \param *Candidates endpoints of the triangle candidate
 * \return integer 0 if no triangle exists, 1 if one triangle exists, 2 if two
 *                 triangles exist which is the maximum for three points
 */
int Tesselation::CheckPresenceOfTriangle(ofstream *out, TesselPoint *Candidates[3]) {
  int adjacentTriangleCount = 0;
  class BoundaryPointSet *Points[3];

  *out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
  // builds a triangle point set (Points) of the end points
  for (int i = 0; i < 3; i++) {
    PointMap::iterator FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
    if (FindPoint != PointsOnBoundary.end()) {
      Points[i] = FindPoint->second;
    } else {
      Points[i] = NULL;
    }
  }

  // checks lines between the points in the Points for their adjacent triangles
  for (int i = 0; i < 3; i++) {
    if (Points[i] != NULL) {
      for (int j = i; j < 3; j++) {
        if (Points[j] != NULL) {
          LineMap::iterator FindLine = Points[i]->lines.find(Points[j]->node->nr);
          for (; (FindLine != Points[i]->lines.end()) && (FindLine->first == Points[j]->node->nr); FindLine++) {
            TriangleMap *triangles = &FindLine->second->triangles;
            *out << Verbose(3) << "Current line is " << FindLine->first << ": " << *(FindLine->second) << " with triangles " << triangles << "." << endl;
            for (TriangleMap::iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) {
              if (FindTriangle->second->IsPresentTupel(Points)) {
                adjacentTriangleCount++;
              }
            }
            *out << Verbose(3) << "end." << endl;
          }
          // Only one of the triangle lines must be considered for the triangle count.
          *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
          return adjacentTriangleCount;
        }
      }
    }
  }

  *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
  *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
  return adjacentTriangleCount;
};


/** Finds the starting triangle for FindNonConvexBorder().
 * Looks at the outermost point per axis, then FindSecondPointForTesselation()
 * for the second and FindNextSuitablePointViaAngleOfSphere() for the third
 * point are called.
 * \param *out output stream for debugging
 * \param RADIUS radius of virtual rolling sphere
 * \param *LC LinkedCell structure with neighbouring TesselPoint's
 */
void Tesselation::FindStartingTriangle(ofstream *out, const double RADIUS, LinkedCell *LC)
{
  cout << Verbose(1) << "Begin of FindStartingTriangle\n";
  int i = 0;
  LinkedNodes *List = NULL;
  TesselPoint* FirstPoint = NULL;
  TesselPoint* SecondPoint = NULL;
  TesselPoint* MaxPoint[NDIM];
  double maxCoordinate[NDIM];
  Vector Oben;
  Vector helper;
  Vector Chord;
  Vector SearchDirection;

  Oben.Zero();

  for (i = 0; i < 3; i++) {
    MaxPoint[i] = NULL;
    maxCoordinate[i] = -1;
  }

  // 1. searching topmost point with respect to each axis
  for (int i=0;i<NDIM;i++) { // each axis
    LC->n[i] = LC->N[i]-1; // current axis is topmost cell
    for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
      for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
        List = LC->GetCurrentCell();
        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
        if (List != NULL) {
          for (LinkedNodes::iterator Runner = List->begin();Runner != List->end();Runner++) {
            if ((*Runner)->node->x[i] > maxCoordinate[i]) {
              cout << Verbose(2) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl;
              maxCoordinate[i] = (*Runner)->node->x[i];
              MaxPoint[i] = (*Runner);
            }
          }
        } else {
          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
        }
      }
  }

  cout << Verbose(2) << "Found maximum coordinates: ";
  for (int i=0;i<NDIM;i++)
    cout << i << ": " << *MaxPoint[i] << "\t";
  cout << endl;

  BTS = NULL;
  CandidateList *OptCandidates = new CandidateList();
  for (int k=0;k<NDIM;k++) {
    Oben.x[k] = 1.;
    FirstPoint = MaxPoint[k];
    cout << Verbose(1) << "Coordinates of start node at " << *FirstPoint->node << "." << endl;

    double ShortestAngle;
    TesselPoint* OptCandidate = NULL;
    ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.

    FindSecondPointForTesselation(FirstPoint, NULL, Oben, OptCandidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
    SecondPoint = OptCandidate;
    if (SecondPoint == NULL)  // have we found a second point?
      continue;
    else
      cout << Verbose(1) << "Found second point is at " << *SecondPoint->node << ".\n";

    helper.CopyVector(FirstPoint->node);
    helper.SubtractVector(SecondPoint->node);
    helper.Normalize();
    Oben.ProjectOntoPlane(&helper);
    Oben.Normalize();
    helper.VectorProduct(&Oben);
    ShortestAngle = 2.*M_PI; // This will indicate the quadrant.

    Chord.CopyVector(FirstPoint->node); // bring into calling function
    Chord.SubtractVector(SecondPoint->node);
    double radius = Chord.ScalarProduct(&Chord);
    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
    helper.CopyVector(&Oben);
    helper.Scale(CircleRadius);
    // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)

    // look in one direction of baseline for initial candidate
    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...

    // adding point 1 and point 2 and add the line between them
    AddTesselationPoint(FirstPoint, 0);
    AddTesselationPoint(SecondPoint, 1);
    AddTesselationLine(TPS[0], TPS[1], 0);

    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n";
    FindThirdPointForTesselation(
      Oben, SearchDirection, helper, BLS[0], NULL, *&OptCandidates, &ShortestAngle, RADIUS, LC
    );
    cout << Verbose(1) << "List of third Points is ";
    for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
        cout << " " << *(*it)->point;
    }
    cout << endl;

    for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
      // add third triangle point
      AddTesselationPoint((*it)->point, 2);
      // add the second and third line
      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);
      AddTesselationTriangle();
      // ... and calculate its normal vector (with correct orientation)
      (*it)->OptCenter.Scale(-1.);
      cout << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl;
      BTS->GetNormalVector((*it)->OptCenter);  // vector to compare with should point inwards
      cout << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and "
      << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n";

      // if we do not reach the end with the next step of iteration, we need to setup a new first line
      if (it != OptCandidates->end()--) {
        FirstPoint = (*it)->BaseLine->endpoints[0]->node;
        SecondPoint = (*it)->point;
        // adding point 1 and point 2 and the line between them
        AddTesselationPoint(FirstPoint, 0);
        AddTesselationPoint(SecondPoint, 1);
        AddTesselationLine(TPS[0], TPS[1], 0);
      }
      cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
    }
    if (BTS != NULL) // we have created one starting triangle
      break;
    else {
      // remove all candidates from the list and then the list itself
      class CandidateForTesselation *remover = NULL;
      for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
        remover = *it;
        delete(remover);
      }
      OptCandidates->clear();
    }
  }

  // remove all candidates from the list and then the list itself
  class CandidateForTesselation *remover = NULL;
  for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
    remover = *it;
    delete(remover);
  }
  delete(OptCandidates);
  cout << Verbose(1) << "End of FindStartingTriangle\n";
};


/** This function finds a triangle to a line, adjacent to an existing one.
 * @param out output stream for debugging
 * @param Line current baseline to search from
 * @param T current triangle which \a Line is edge of
 * @param RADIUS radius of the rolling ball
 * @param N number of found triangles
 * @param *LC LinkedCell structure with neighbouring points
 */
bool Tesselation::FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, LinkedCell *LC)
{
  cout << Verbose(0) << "Begin of FindNextSuitableTriangle\n";
  bool result = true;
  CandidateList *OptCandidates = new CandidateList();

  Vector CircleCenter;
  Vector CirclePlaneNormal;
  Vector OldSphereCenter;
  Vector SearchDirection;
  Vector helper;
  TesselPoint *ThirdNode = NULL;
  LineMap::iterator testline;
  double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
  double radius, CircleRadius;

  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
  for (int i=0;i<3;i++)
    if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
      ThirdNode = T.endpoints[i]->node;

  // construct center of circle
  CircleCenter.CopyVector(Line.endpoints[0]->node->node);
  CircleCenter.AddVector(Line.endpoints[1]->node->node);
  CircleCenter.Scale(0.5);

  // construct normal vector of circle
  CirclePlaneNormal.CopyVector(Line.endpoints[0]->node->node);
  CirclePlaneNormal.SubtractVector(Line.endpoints[1]->node->node);

  // calculate squared radius of circle
  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
  if (radius/4. < RADIUS*RADIUS) {
    CircleRadius = RADIUS*RADIUS - radius/4.;
    CirclePlaneNormal.Normalize();
    cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;

    // construct old center
    GetCenterofCircumcircle(&OldSphereCenter, T.endpoints[0]->node->node, T.endpoints[1]->node->node, T.endpoints[2]->node->node);
    helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
    radius = Line.endpoints[0]->node->node->DistanceSquared(&OldSphereCenter);
    helper.Scale(sqrt(RADIUS*RADIUS - radius));
    OldSphereCenter.AddVector(&helper);
    OldSphereCenter.SubtractVector(&CircleCenter);
    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;

    // construct SearchDirection
    SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
    helper.CopyVector(Line.endpoints[0]->node->node);
    helper.SubtractVector(ThirdNode->node);
    if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
      SearchDirection.Scale(-1.);
    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
    SearchDirection.Normalize();
    cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
      // rotated the wrong way!
      cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
    }

    // add third point
    FindThirdPointForTesselation(
      T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, OptCandidates,
      &ShortestAngle, RADIUS, LC
    );

  } else {
    cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
  }

  if (OptCandidates->begin() == OptCandidates->end()) {
    cerr << "WARNING: Could not find a suitable candidate." << endl;
    return false;
  }
  cout << Verbose(1) << "Third Points are ";
  for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
    cout << " " << *(*it)->point;
  }
  cout << endl;

  BoundaryLineSet *BaseRay = &Line;
  for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
    cout << Verbose(1) << " Third point candidate is " << *(*it)->point
    << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
    cout << Verbose(1) << " Baseline is " << *BaseRay << endl;

    // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
    TesselPoint *PointCandidates[3];
    PointCandidates[0] = (*it)->point;
    PointCandidates[1] = BaseRay->endpoints[0]->node;
    PointCandidates[2] = BaseRay->endpoints[1]->node;
    int existentTrianglesCount = CheckPresenceOfTriangle(out, PointCandidates);

    BTS = NULL;
    // If there is no triangle, add it regularly.
    if (existentTrianglesCount == 0) {
      AddTesselationPoint((*it)->point, 0);
      AddTesselationPoint(BaseRay->endpoints[0]->node, 1);
      AddTesselationPoint(BaseRay->endpoints[1]->node, 2);

      if (CheckLineCriteriaForDegeneratedTriangle(TPS)) {
        AddTesselationLine(TPS[0], TPS[1], 0);
        AddTesselationLine(TPS[0], TPS[2], 1);
        AddTesselationLine(TPS[1], TPS[2], 2);

        BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
        AddTesselationTriangle();
        (*it)->OptCenter.Scale(-1.);
        BTS->GetNormalVector((*it)->OptCenter);
        (*it)->OptCenter.Scale(-1.);

        cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector
          << " for this triangle ... " << endl;
      //cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl;
      } else {
        cout << Verbose(1) << "WARNING: This triangle consisting of ";
        cout << *(*it)->point << ", ";
        cout << *BaseRay->endpoints[0]->node << " and ";
        cout << *BaseRay->endpoints[1]->node << " ";
        cout << "exists and is not added, as it does not seem helpful!" << endl;
        result = false;
      }
    } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
        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)) {
          AddTesselationLine(TPS[0], TPS[1], 0);
          AddTesselationLine(TPS[0], TPS[2], 1);
          AddTesselationLine(TPS[1], TPS[2], 2);

          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
          AddTesselationTriangle();  // add to global map

          (*it)->OtherOptCenter.Scale(-1.);
          BTS->GetNormalVector((*it)->OtherOptCenter);
          (*it)->OtherOptCenter.Scale(-1.);

          cout << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector
          << " for this triangle ... " << endl;
          cout << Verbose(1) << "We have "<< BaseRay->triangles.size() << " for line " << BaseRay << "." << endl;
        } else {
          cout << Verbose(1) << "WARNING: This triangle consisting of ";
          cout << *(*it)->point << ", ";
          cout << *BaseRay->endpoints[0]->node << " and ";
          cout << *BaseRay->endpoints[1]->node << " ";
          cout << "exists and is not added, as it does not seem helpful!" << endl;
          result = false;
        }
    } else {
      cout << Verbose(1) << "This triangle consisting of ";
      cout << *(*it)->point << ", ";
      cout << *BaseRay->endpoints[0]->node << " and ";
      cout << *BaseRay->endpoints[1]->node << " ";
      cout << "is invalid!" << endl;
      result = false;
    }

    // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point))
    BaseRay = BLS[0];
  }

  // remove all candidates from the list and then the list itself
  class CandidateForTesselation *remover = NULL;
  for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
    remover = *it;
    delete(remover);
  }
  delete(OptCandidates);
  cout << Verbose(0) << "End of FindNextSuitableTriangle\n";
  return result;
};

/** 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
 * \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) {
      *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()

    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";
      }
  *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;

  // 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.
 * \param *a first node
 * \param *Candidate pointer to candidate node on return
 * \param Oben vector indicating the outside
 * \param OptCandidate reference to recommended candidate on return
 * \param Storage[3] array storing angles and other candidate information
 * \param RADIUS radius of virtual sphere
 * \param *LC LinkedCell structure with neighbouring points
 */
void Tesselation::FindSecondPointForTesselation(TesselPoint* a, TesselPoint* Candidate, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, LinkedCell *LC)
{
  cout << Verbose(2) << "Begin of FindSecondPointForTesselation" << endl;
  Vector AngleCheck;
  double norm = -1., angle;
  LinkedNodes *List = NULL;
  int N[NDIM], Nlower[NDIM], Nupper[NDIM];

  if (LC->SetIndexToNode(a)) {  // get cell for the starting point
    for(int i=0;i<NDIM;i++) // store indices of this cell
      N[i] = LC->n[i];
  } else {
    cerr << "ERROR: Point " << *a << " is not found in cell " << LC->index << "." << endl;
    return;
  }
  // then go through the current and all neighbouring cells and check the contained points for possible candidates
  cout << Verbose(3) << "LC Intervals from [";
  for (int i=0;i<NDIM;i++) {
  cout << " " << N[i] << "<->" << LC->N[i];
  }
  cout << "] :";
  for (int i=0;i<NDIM;i++) {
    Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
    Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
    cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
  }
  cout << endl;


  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
        List = LC->GetCurrentCell();
        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
        if (List != NULL) {
          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
            Candidate = (*Runner);
            // check if we only have one unique point yet ...
            if (a != Candidate) {
              // Calculate center of the circle with radius RADIUS through points a and Candidate
              Vector OrthogonalizedOben, aCandidate, Center;
              double distance, scaleFactor;

              OrthogonalizedOben.CopyVector(&Oben);
              aCandidate.CopyVector(a->node);
              aCandidate.SubtractVector(Candidate->node);
              OrthogonalizedOben.ProjectOntoPlane(&aCandidate);
              OrthogonalizedOben.Normalize();
              distance = 0.5 * aCandidate.Norm();
              scaleFactor = sqrt(((RADIUS * RADIUS) - (distance * distance)));
              OrthogonalizedOben.Scale(scaleFactor);

              Center.CopyVector(Candidate->node);
              Center.AddVector(a->node);
              Center.Scale(0.5);
              Center.AddVector(&OrthogonalizedOben);

              AngleCheck.CopyVector(&Center);
              AngleCheck.SubtractVector(a->node);
              norm = aCandidate.Norm();
              // second point shall have smallest angle with respect to Oben vector
              if (norm < RADIUS*2.) {
                angle = AngleCheck.Angle(&Oben);
                if (angle < Storage[0]) {
                  //cout << Verbose(3) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
                  cout << Verbose(3) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";
                  OptCandidate = Candidate;
                  Storage[0] = angle;
                  //cout << Verbose(3) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
                } else {
                  //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *OptCandidate << endl;
                }
              } else {
                //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
              }
            } else {
              //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
            }
          }
        } else {
          cout << Verbose(3) << "Linked cell list is empty." << endl;
        }
      }
  cout << Verbose(2) << "End of FindSecondPointForTesselation" << endl;
};


/** This recursive function finds a third point, to form a triangle with two given ones.
 * Note that this function is for the starting triangle.
 * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
 * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
 * the center of the sphere is still fixed up to a single parameter. The band of possible values
 * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
 * us the "null" on this circle, the new center of the candidate point will be some way along this
 * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
 * by the normal vector of the base triangle that always points outwards by construction.
 * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
 * We construct the normal vector that defines the plane this circle lies in, it is just in the
 * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
 * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
 * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
 * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
 * both.
 * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
 * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
 * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
 * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
 * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
 * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
 * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa FindStartingTriangle())
 * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
 * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
 * @param BaseLine BoundaryLineSet with the current base line
 * @param ThirdNode third point to avoid in search
 * @param candidates list of equally good candidates to return
 * @param ShortestAngle the current path length on this circle band for the current OptCandidate
 * @param RADIUS radius of sphere
 * @param *LC LinkedCell structure with neighbouring points
 */
void Tesselation::FindThirdPointForTesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, class TesselPoint  *ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
{
  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
  Vector SphereCenter;
  Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
  Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
  Vector NewNormalVector;   // normal vector of the Candidate's triangle
  Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
  LinkedNodes *List = NULL;
  double CircleRadius; // radius of this circle
  double radius;
  double alpha, Otheralpha; // angles (i.e. parameter for the circle).
  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
  TesselPoint *Candidate = NULL;
  CandidateForTesselation *optCandidate = NULL;

  cout << Verbose(1) << "Begin of FindThirdPointForTesselation" << endl;

  //cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;

  // construct center of circle
  CircleCenter.CopyVector(BaseLine->endpoints[0]->node->node);
  CircleCenter.AddVector(BaseLine->endpoints[1]->node->node);
  CircleCenter.Scale(0.5);

  // construct normal vector of circle
  CirclePlaneNormal.CopyVector(BaseLine->endpoints[0]->node->node);
  CirclePlaneNormal.SubtractVector(BaseLine->endpoints[1]->node->node);

  // calculate squared radius TesselPoint *ThirdNode,f circle
  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
  if (radius/4. < RADIUS*RADIUS) {
    CircleRadius = RADIUS*RADIUS - radius/4.;
    CirclePlaneNormal.Normalize();
    //cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;

    // test whether old center is on the band's plane
    if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
      cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
      OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
    }
    radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
    if (fabs(radius - CircleRadius) < HULLEPSILON) {

      // check SearchDirection
      //cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
        cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
      }

      // get cell for the starting point
      if (LC->SetIndexToVector(&CircleCenter)) {
        for(int i=0;i<NDIM;i++) // store indices of this cell
        N[i] = LC->n[i];
        //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
      } else {
        cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
        return;
      }
      // then go through the current and all neighbouring cells and check the contained points for possible candidates
      //cout << Verbose(2) << "LC Intervals:";
      for (int i=0;i<NDIM;i++) {
        Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
        Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
        //cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
      }
      //cout << endl;
      for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
        for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
          for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
            List = LC->GetCurrentCell();
            //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
            if (List != NULL) {
              for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
                Candidate = (*Runner);

                // check for three unique points
                //cout << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->node << "." << endl;
                if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) ){

                  // construct both new centers
                  GetCenterofCircumcircle(&NewSphereCenter, BaseLine->endpoints[0]->node->node, BaseLine->endpoints[1]->node->node, Candidate->node);
                  OtherNewSphereCenter.CopyVector(&NewSphereCenter);

                  if ((NewNormalVector.MakeNormalVector(BaseLine->endpoints[0]->node->node, BaseLine->endpoints[1]->node->node, Candidate->node))
                  && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
                  ) {
                    helper.CopyVector(&NewNormalVector);
                    //cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
                    radius = BaseLine->endpoints[0]->node->node->DistanceSquared(&NewSphereCenter);
                    if (radius < RADIUS*RADIUS) {
                      helper.Scale(sqrt(RADIUS*RADIUS - radius));
                      //cout << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;
                      NewSphereCenter.AddVector(&helper);
                      NewSphereCenter.SubtractVector(&CircleCenter);
                      //cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;

                      // OtherNewSphereCenter is created by the same vector just in the other direction
                      helper.Scale(-1.);
                      OtherNewSphereCenter.AddVector(&helper);
                      OtherNewSphereCenter.SubtractVector(&CircleCenter);
                      //cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;

                      alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
                      Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
                      alpha = min(alpha, Otheralpha);
                      // if there is a better candidate, drop the current list and add the new candidate
                      // otherwise ignore the new candidate and keep the list
                      if (*ShortestAngle > (alpha - HULLEPSILON)) {
                        optCandidate = new CandidateForTesselation(Candidate, BaseLine, OptCandidateCenter, OtherOptCandidateCenter);
                        if (fabs(alpha - Otheralpha) > MYEPSILON) {
                          optCandidate->OptCenter.CopyVector(&NewSphereCenter);
                          optCandidate->OtherOptCenter.CopyVector(&OtherNewSphereCenter);
                        } else {
                          optCandidate->OptCenter.CopyVector(&OtherNewSphereCenter);
                          optCandidate->OtherOptCenter.CopyVector(&NewSphereCenter);
                        }
                        // if there is an equal candidate, add it to the list without clearing the list
                        if ((*ShortestAngle - HULLEPSILON) < alpha) {
                          candidates->push_back(optCandidate);
                          cout << Verbose(2) << "ACCEPT: We have found an equally good candidate: " << *(optCandidate->point) << " with "
                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
                        } else {
                          // remove all candidates from the list and then the list itself
                          class CandidateForTesselation *remover = NULL;
                          for (CandidateList::iterator it = candidates->begin(); it != candidates->end(); ++it) {
                            remover = *it;
                            delete(remover);
                          }
                          candidates->clear();
                          candidates->push_back(optCandidate);
                          cout << Verbose(2) << "ACCEPT: We have found a better candidate: " << *(optCandidate->point) << " with "
                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
                        }
                        *ShortestAngle = alpha;
                        //cout << Verbose(2) << "INFO: There are " << candidates->size() << " candidates in the list now." << endl;
                      } else {
                        if ((optCandidate != NULL) && (optCandidate->point != NULL)) {
                          //cout << Verbose(2) << "REJECT: Old candidate " << *(optCandidate->point) << " with " << *ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl;
                        } else {
                          //cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
                        }
                      }

                    } else {
                      //cout << Verbose(2) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl;
                    }
                  } else {
                    //cout << Verbose(2) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
                  }
                } else {
                  if (ThirdNode != NULL) {
                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
                  } else {
                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
                  }
                }
              }
            }
          }
    } else {
      cerr << Verbose(2) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
    }
  } else {
    if (ThirdNode != NULL)
      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
    else
      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
  }

  //cout << Verbose(2) << "INFO: Sorting candidate list ..." << endl;
  if (candidates->size() > 1) {
    candidates->unique();
    candidates->sort(SortCandidates);
  }

  cout << Verbose(1) << "End of FindThirdPointForTesselation" << endl;
};

/** Finds the endpoint two lines are sharing.
 * \param *line1 first line
 * \param *line2 second line
 * \return point which is shared or NULL if none
 */
class BoundaryPointSet *Tesselation::GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
{
  class BoundaryLineSet * lines[2] =
    { line1, line2 };
  class BoundaryPointSet *node = NULL;
  map<int, class BoundaryPointSet *> OrderMap;
  pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
  for (int i = 0; i < 2; i++)
    // for both lines
    for (int j = 0; j < 2; j++)
      { // for both endpoints
        OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (
            lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
        if (!OrderTest.second)
          { // if insertion fails, we have common endpoint
            node = OrderTest.first->second;
            cout << Verbose(5) << "Common endpoint of lines " << *line1
                << " and " << *line2 << " is: " << *node << "." << endl;
            j = 2;
            i = 2;
            break;
          }
      }
  return node;
};

/** Finds the triangle that is closest to a given Vector \a *x.
 * \param *out output stream for debugging
 * \param *x Vector to look from
 * \return list of BoundaryTriangleSet of nearest triangles or NULL in degenerate case.
 */
list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(ofstream *out, Vector *x, LinkedCell* LC)
{
  TesselPoint *trianglePoints[3];
  TesselPoint *SecondPoint = NULL;

  if (LinesOnBoundary.empty()) {
    *out << Verbose(0) << "Error: There is no tesselation structure to compare the point with, please create one first.";
    return NULL;
  }

  trianglePoints[0] = FindClosestPoint(x, SecondPoint, LC);
  
  // check whether closest point is "too close" :), then it's inside
  if (trianglePoints[0] == NULL) {
    *out << Verbose(1) << "Is the only point, no one else is closeby." << endl;
    return NULL;
  }
  if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) {
    *out << Verbose(1) << "Point is right on a tesselation point, no nearest triangle." << endl;
    return NULL;
  }
  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();
  for (int i=0;i<3;i++) {
    if (trianglePoints[i] == NULL) {
      *out << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl;
    }
    //*out << Verbose(1) << "List of possible points:" << endl;
    //*out << Verbose(2) << *trianglePoints[i] << endl;
  }

  list<BoundaryTriangleSet*> *triangles = FindTriangles(trianglePoints);

  delete(connectedClosestPoints);
  
  if (triangles->empty()) {
    *out << Verbose(0) << "Error: There is no nearest triangle. Please check the tesselation structure.";
    return NULL;
  } else
    return triangles;
};

/** Finds closest triangle to a point.
 * This basically just takes care of the degenerate case, which is not handled in FindClosestTrianglesToPoint().
 * \param *out output stream for debugging
 * \param *x Vector to look from
 * \return list of BoundaryTriangleSet of nearest triangles or NULL.
 */
class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(ofstream *out, Vector *x, LinkedCell* LC)
{
  class BoundaryTriangleSet *result = NULL;
  list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(out, x, LC);

  if (triangles == NULL)
    return NULL;

  if (x->ScalarProduct(&triangles->front()->NormalVector) < 0)
    result = triangles->back();
  else
    result = triangles->front();

  delete(triangles);
  return result;
};

/** Checks whether the provided Vector is within the tesselation structure.
 *
 * @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(ofstream *out, Vector Point, LinkedCell* LC)
{
  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(out, &Point, LC);
  if (result == NULL)
    return true;
  if (Point.ScalarProduct(&result->NormalVector) < 0)
    return true;
  else
    return false;
}

/** Checks whether the provided TesselPoint is within the tesselation structure.
 *
 * @param *Point of which to check the position
 * @param *LC Linked Cell structure
 *
 * @return true if the point is inside the tesselation structure, false otherwise
 */
bool Tesselation::IsInnerPoint(ofstream *out, TesselPoint *Point, LinkedCell* LC)
{
  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(out, Point->node, LC);
  if (result == NULL)
    return true;
  if (Point->node->ScalarProduct(&result->NormalVector) < 0)
    return true;
  else
    return false;
}

/** Gets all points connected to the provided point by triangulation lines.
 *
 * @param *Point of which get all connected points
 *
 * @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;
  TesselPoint* current;
  bool takePoint = false;

  // find the respective boundary point
  PointMap::iterator PointRunner = PointsOnBoundary.find(Point->nr);
  if (PointRunner != PointsOnBoundary.end()) {
    ReferencePoint = PointRunner->second;
  } else {
    *out << Verbose(2) << "getCircleOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
    ReferencePoint = NULL;
  }

  // little trick so that we look just through lines connect to the BoundaryPoint 
  // OR fall-back to look through all lines if there is no such BoundaryPoint
  LineMap *Lines = &LinesOnBoundary;
  if (ReferencePoint != NULL)
    Lines = &(ReferencePoint->lines);
  LineMap::iterator findLines = Lines->begin();
  while (findLines != Lines->end()) {
    takePoint = false;

    if (findLines->second->endpoints[0]->Nr == Point->nr) {
      takePoint = true;
      current = findLines->second->endpoints[1]->node;
    } else if (findLines->second->endpoints[1]->Nr == Point->nr) {
      takePoint = true;
      current = findLines->second->endpoints[0]->node;
    }
    
    if (takePoint) {
      *out << Verbose(3) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is taken into the circle." << endl;
      connectedPoints->push_back(current);
    }

    findLines++;
  }

  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();
  center.Scale(1.0/connectedPoints->size());
  *out << Verbose(4) << "INFO: Calculated center of all circle points is " << center << "." << endl;

  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
  PlaneNormal.CopyVector(Point->node);
  PlaneNormal.SubtractVector(&center);
  PlaneNormal.Normalize();
  *out << Verbose(4) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;

  // construct one orthogonal vector
  AngleZero.CopyVector(Reference);
  AngleZero.SubtractVector(Point->node);
  AngleZero.ProjectOntoPlane(&PlaneNormal);
  *out << Verbose(4) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl;
  OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
  *out << Verbose(4) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;

  // go through all connected points and calculate angle
  for (list<TesselPoint*>::iterator listRunner = connectedPoints->begin(); listRunner != connectedPoints->end(); listRunner++) {
    helper.CopyVector((*listRunner)->node);
    helper.SubtractVector(Point->node);
    helper.ProjectOntoPlane(&PlaneNormal);
    double angle = GetAngle(helper, AngleZero, OrthogonalVector);
    *out << Verbose(2) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
    anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
  }

  list<TesselPoint*> *result = new list<TesselPoint*>;
  runner = anglesOfPoints.begin();
  result->push_back(runner->second);
  runner = anglesOfPoints.end();
  runner--;
  result->push_back(runner->second);

  *out << Verbose(2) << "List of closest points has " << result->size() << " elements, which are "
    << *(result->front()) << " and " << *(result->back()) << endl;

  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.
 * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
 * make it bigger (i.e. closing one (the baseline) and opening two new ones).
 * \param TPS[3] nodes of the triangle
 * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
 */
bool CheckLineCriteriaForDegeneratedTriangle(class BoundaryPointSet *nodes[3])
{
  bool result = false;
  int counter = 0;

  // check all three points
  for (int i=0;i<3;i++)
    for (int j=i+1; j<3; j++) {
      if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
        LineMap::iterator FindLine;
        pair<LineMap::iterator,LineMap::iterator> FindPair;
        FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
        for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
          // If there is a line with less than two attached triangles, we don't need a new line.
          if (FindLine->second->triangles.size() < 2) {
            counter++;
            break;  // increase counter only once per edge
          }
        }
      } else { // no line
        cout << Verbose(1) << "The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle." << endl;
        result = true;
      }
    }
  if (counter > 1) {
    cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
    result = true;
  }
  return result;
};


/** Sort function for the candidate list.
 */
bool SortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2)
{
  Vector BaseLineVector, OrthogonalVector, helper;
  if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
    cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
    //return false;
    exit(1);
  }
  // create baseline vector
  BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node);
  BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
  BaseLineVector.Normalize();

  // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!)
  helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node);
  helper.SubtractVector(candidate1->point->node);
  OrthogonalVector.CopyVector(&helper);
  helper.VectorProduct(&BaseLineVector);
  OrthogonalVector.SubtractVector(&helper);
  OrthogonalVector.Normalize();

  // calculate both angles and correct with in-plane vector
  helper.CopyVector(candidate1->point->node);
  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
  double phi = BaseLineVector.Angle(&helper);
  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    phi = 2.*M_PI - phi;
  }
  helper.CopyVector(candidate2->point->node);
  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
  double psi = BaseLineVector.Angle(&helper);
  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    psi = 2.*M_PI - psi;
  }

  cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
  cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl;

  // return comparison
  return phi < psi;
};

/**
 * Finds the point which is second closest to the provided one.
 *
 * @param Point to which to find the second closest other point
 * @param linked cell structure
 *
 * @return point which is second closest to the provided one
 */
TesselPoint* FindSecondClosestPoint(const Vector* Point, LinkedCell* LC)
{
  LinkedNodes *List = NULL;
  TesselPoint* closestPoint = NULL;
  TesselPoint* secondClosestPoint = NULL;
  double distance = 1e16;
  double secondDistance = 1e16;
  Vector helper;
  int N[NDIM], Nlower[NDIM], Nupper[NDIM];

  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
  for(int i=0;i<NDIM;i++) // store indices of this cell
    N[i] = LC->n[i];
  cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;

  LC->GetNeighbourBounds(Nlower, Nupper);
  //cout << endl;
  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
        List = LC->GetCurrentCell();
        cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
        if (List != NULL) {
          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
            helper.CopyVector(Point);
            helper.SubtractVector((*Runner)->node);
            double currentNorm = helper. Norm();
            if (currentNorm < distance) {
              // remember second point
              secondDistance = distance;
              secondClosestPoint = closestPoint;
              // mark down new closest point
              distance = currentNorm;
              closestPoint = (*Runner);
              cout << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
            }
          }
        } else {
          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
            << LC->n[2] << " is invalid!" << endl;
        }
      }

  return secondClosestPoint;
};

/**
 * Finds the point which is closest to the provided one.
 *
 * @param Point to which to find the closest other point
 * @param SecondPoint the second closest other point on return, NULL if none found
 * @param linked cell structure
 *
 * @return point which is closest to the provided one, NULL if none found
 */
TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC)
{
  LinkedNodes *List = NULL;
  TesselPoint* closestPoint = NULL;
  SecondPoint = NULL;
  double distance = 1e16;
  double secondDistance = 1e16;
  Vector helper;
  int N[NDIM], Nlower[NDIM], Nupper[NDIM];

  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
  for(int i=0;i<NDIM;i++) // store indices of this cell
    N[i] = LC->n[i];
  cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;

  LC->GetNeighbourBounds(Nlower, Nupper);
  //cout << endl;
  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
        List = LC->GetCurrentCell();
        cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
        if (List != NULL) {
          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
            helper.CopyVector(Point);
            helper.SubtractVector((*Runner)->node);
            double currentNorm = helper. Norm();
            if (currentNorm < distance) {
              secondDistance = distance;
              SecondPoint = closestPoint;
              distance = currentNorm;
              closestPoint = (*Runner);
              cout << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
            } else if (currentNorm < secondDistance) {
              secondDistance = currentNorm;
              SecondPoint = (*Runner);
              cout << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
            }
          }
        } else {
          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
            << LC->n[2] << " is invalid!" << endl;
        }
      }

  return closestPoint;
};

/**
 * Finds triangles belonging to the three provided points.
 *
 * @param *Points[3] list, is expected to contain three points
 *
 * @return triangles which belong to the provided points, will be empty if there are none,
 *         will usually be one, in case of degeneration, there will be two
 */
list<BoundaryTriangleSet*> *Tesselation::FindTriangles(TesselPoint* Points[3]) 
{
  list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
  LineMap::iterator FindLine;
  PointMap::iterator FindPoint;
  TriangleMap::iterator FindTriangle;
  class BoundaryPointSet *TrianglePoints[3];

  for (int i = 0; i < 3; i++) {
    FindPoint = PointsOnBoundary.find(Points[i]->nr);
    if (FindPoint != PointsOnBoundary.end()) {
      TrianglePoints[i] = FindPoint->second;
    } else {
      TrianglePoints[i] = NULL;
    }
  }

  // checks lines between the points in the Points for their adjacent triangles
  for (int i = 0; i < 3; i++) {
    if (TrianglePoints[i] != NULL) {
      for (int j = i; j < 3; j++) {
        if (TrianglePoints[j] != NULL) {
          FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr);
          if (FindLine != TrianglePoints[i]->lines.end()) {
            for (; FindLine->first == TrianglePoints[j]->node->nr; FindLine++) {
              FindTriangle = FindLine->second->triangles.begin();
              for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
                if ((
                  (FindTriangle->second->endpoints[0] == TrianglePoints[0])
                    || (FindTriangle->second->endpoints[0] == TrianglePoints[1])
                    || (FindTriangle->second->endpoints[0] == TrianglePoints[2])
                  ) && (
                    (FindTriangle->second->endpoints[1] == TrianglePoints[0])
                    || (FindTriangle->second->endpoints[1] == TrianglePoints[1])
                    || (FindTriangle->second->endpoints[1] == TrianglePoints[2])
                  ) && (
                    (FindTriangle->second->endpoints[2] == TrianglePoints[0])
                    || (FindTriangle->second->endpoints[2] == TrianglePoints[1])
                    || (FindTriangle->second->endpoints[2] == TrianglePoints[2])
                  )
                ) {
                  result->push_back(FindTriangle->second);
                }
              }
            }
            // Is it sufficient to consider one of the triangle lines for this.
            return result;

          }
        }
      }
    }
  }

  return result;
}

/**
 * Finds all degenerated triangles within the tesselation structure.
 *
 * @return map of keys of degenerated triangle pairs, each triangle occurs twice
 *         in the list, once as key and once as value
 */
map<int, int> Tesselation::FindAllDegeneratedTriangles()
{
  map<int, int> DegeneratedTriangles;

  // sanity check
  if (LinesOnBoundary.empty()) {
    cout << Verbose(1) << "Warning: FindAllDegeneratedTriangles() was called without any tesselation structure.";
    return DegeneratedTriangles;
  }

  LineMap::iterator LineRunner1, LineRunner2;

  for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
    for (LineRunner2 = LinesOnBoundary.begin(); LineRunner2 != LinesOnBoundary.end(); ++LineRunner2) {
      if ((LineRunner1->second != LineRunner2->second)
        && (LineRunner1->second->endpoints[0] == LineRunner2->second->endpoints[0])
        && (LineRunner1->second->endpoints[1] == LineRunner2->second->endpoints[1])
      ) {
        TriangleMap::iterator TriangleRunner1 = LineRunner1->second->triangles.begin(),
          TriangleRunner2 = LineRunner2->second->triangles.begin();

        for (; TriangleRunner1 != LineRunner1->second->triangles.end(); ++TriangleRunner1) {
          for (; TriangleRunner2 != LineRunner2->second->triangles.end(); ++TriangleRunner2) {
            if ((TriangleRunner1->second != TriangleRunner2->second)
              && (TriangleRunner1->second->endpoints[0] == TriangleRunner2->second->endpoints[0])
              && (TriangleRunner1->second->endpoints[1] == TriangleRunner2->second->endpoints[1])
              && (TriangleRunner1->second->endpoints[2] == TriangleRunner2->second->endpoints[2])
            ) {
              DegeneratedTriangles[TriangleRunner1->second->Nr] = TriangleRunner2->second->Nr;
              DegeneratedTriangles[TriangleRunner2->second->Nr] = TriangleRunner1->second->Nr;
            }
          }
        }
      }
    }
  }

  cout << Verbose(1) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles.size() << " triangles." << endl;
  map<int,int>::iterator it;
  for (it = DegeneratedTriangles.begin(); it != DegeneratedTriangles.end(); it++)
      cout << Verbose(2) << (*it).first << " => " << (*it).second << endl;

  return DegeneratedTriangles;
}

/**
 * Purges degenerated triangles from the tesselation structure if they are not
 * necessary to keep a single point within the structure.
 */
void Tesselation::RemoveDegeneratedTriangles()
{
  map<int, int> DegeneratedTriangles = FindAllDegeneratedTriangles();

  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles.begin();
    TriangleKeyRunner != DegeneratedTriangles.end(); ++TriangleKeyRunner
  ) {
    BoundaryTriangleSet *triangle = TrianglesOnBoundary.find(TriangleKeyRunner->first)->second,
      *partnerTriangle = TrianglesOnBoundary.find(TriangleKeyRunner->second)->second;

    bool trianglesShareLine = false;
    for (int i = 0; i < 3; ++i)
      for (int j = 0; j < 3; ++j)
        trianglesShareLine = trianglesShareLine || triangle->lines[i] == partnerTriangle->lines[j];

    if (trianglesShareLine
      && (triangle->endpoints[1]->LinesCount > 2)
      && (triangle->endpoints[2]->LinesCount > 2)
      && (triangle->endpoints[0]->LinesCount > 2)
    ) {
      cout << Verbose(1) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl;
      RemoveTesselationTriangle(triangle);
      cout << Verbose(1) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl;
      RemoveTesselationTriangle(partnerTriangle);
      DegeneratedTriangles.erase(DegeneratedTriangles.find(partnerTriangle->Nr));
    } else {
      cout << Verbose(1) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle
        << " and its partner " << *partnerTriangle << " because it is essential for at"
        << " least one of the endpoints to be kept in the tesselation structure." << endl;
    }
  }
}

/** Gets the angle between a point and a reference relative to the provided center.
 * We have two shanks point and reference between which the angle is calculated
 * and by scalar product with OrthogonalVector we decide the interval.
 * @param point to calculate the angle for
 * @param reference to which to calculate the angle
 * @param OrthogonalVector points in direction of [pi,2pi] interval
 *
 * @return angle between point and reference
 */
double GetAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector)
{
  if (reference.IsNull())
    return M_PI;

  // calculate both angles and correct with in-plane vector
  if (point.IsNull())
    return M_PI;
  double phi = point.Angle(&reference);
  if (OrthogonalVector.ScalarProduct(&point) > 0) {
    phi = 2.*M_PI - phi;
  }

  cout << Verbose(3) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;

  return phi;
}