Changeset 0d111b for molecuilder/src/tesselation.cpp

Timestamp:

Apr 29, 2010, 1:55:21 PM (16 years ago)

Author:

Tillmann Crueger <crueger@…>

Children:

070651, 5d1a94

Parents:

90c4460 (diff), 32842d8 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'VectorRefactoring' into StructureRefactoring

Conflicts:

molecuilder/src/Legacy/oldmenu.cpp
molecuilder/src/Makefile.am
molecuilder/src/analysis_correlation.cpp
molecuilder/src/boundary.cpp
molecuilder/src/builder.cpp
molecuilder/src/config.cpp
molecuilder/src/ellipsoid.cpp
molecuilder/src/linkedcell.cpp
molecuilder/src/molecule.cpp
molecuilder/src/molecule_fragmentation.cpp
molecuilder/src/molecule_geometry.cpp
molecuilder/src/molecule_graph.cpp
molecuilder/src/moleculelist.cpp
molecuilder/src/tesselation.cpp
molecuilder/src/tesselationhelpers.cpp
molecuilder/src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp
molecuilder/src/unittests/bondgraphunittest.cpp
molecuilder/src/vector.cpp
molecuilder/src/vector.hpp

File:

• molecuilder/src/tesselation.cpp (modified) (81 diffs)

molecuilder/src/tesselation.cpp

@@ -17 +17 @@
 #include "triangleintersectionlist.hpp"
 #include "vector.hpp"
+#include "vector_ops.hpp"
 #include "verbose.hpp"
+#include "Plane.hpp"
+#include "Exceptions/LinearDependenceException.hpp"
 
 class molecule;
@@ -234 +237 @@
   // have a normal vector on the base line pointing outwards
   //Log() << Verbose(0) << "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);
+  BaseLineCenter = (1./2.)*((*endpoints[0]->node->node) + (*endpoints[1]->node->node));
+  BaseLine = (*endpoints[0]->node->node) - (*endpoints[1]->node->node);
+
   //Log() << Verbose(0) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl;
 
@@ -248 +249 @@
   for (TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
     //Log() << Verbose(0) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
-    NormalCheck.AddVector(&runner->second->NormalVector);
-    NormalCheck.Scale(sign);
+    NormalCheck += runner->second->NormalVector;
+    NormalCheck *= sign;
     sign = -sign;
     if (runner->second->NormalVector.NormSquared() > MYEPSILON)
-      BaseLineNormal.CopyVector(&runner->second->NormalVector); // yes, copy second on top of first
+      BaseLineNormal = runner->second->NormalVector;   // yes, copy second on top of first
     else {
       DoeLog(0) && (eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl);
@@ -259 +260 @@
     if (node != NULL) {
       //Log() << Verbose(0) << "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!
+      helper[i] = (*node->node->node) - BaseLineCenter;
+      helper[i].MakeNormalTo(BaseLine);  // we want to compare the triangle's heights' angles!
       //Log() << Verbose(0) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
       i++;
@@ -409 +409 @@
   Info FunctionInfo(__func__);
   // get normal vector
-  NormalVector.MakeNormalVector(endpoints[0]->node->node, endpoints[1]->node->node, endpoints[2]->node->node);
+  NormalVector = Plane(*(endpoints[0]->node->node),
+                       *(endpoints[1]->node->node),
+                       *(endpoints[2]->node->node)).getNormal();
 
   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
-  if (NormalVector.ScalarProduct(&OtherVector) > 0.)
+  if (NormalVector.ScalarProduct(OtherVector) > 0.)
     NormalVector.Scale(-1.);
   DoLog(1) && (Log() << Verbose(1) << "Normal Vector is " << NormalVector << "." << endl);
@@ -430 +432 @@
  * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
  */
+
 bool BoundaryTriangleSet::GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const
 {
@@ -436 +439 @@
   Vector helper;
 
-  if (!Intersection->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, MolCenter, x)) {
+  try {
+    *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*MolCenter, *x);
+  }
+  catch (LinearDependenceException &excp) {
+    Log() << Verbose(1) << excp;
     DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
     return false;
@@ -445 +452 @@
   DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);
 
-  if (Intersection->DistanceSquared(endpoints[0]->node->node) < MYEPSILON) {
+  if (Intersection->DistanceSquared(*endpoints[0]->node->node) < MYEPSILON) {
     DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);
     return true;
-  } else if (Intersection->DistanceSquared(endpoints[1]->node->node) < MYEPSILON) {
+  }   else if (Intersection->DistanceSquared(*endpoints[1]->node->node) < MYEPSILON) {
     DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);
     return true;
-  } else if (Intersection->DistanceSquared(endpoints[2]->node->node) < MYEPSILON) {
+  }   else if (Intersection->DistanceSquared(*endpoints[2]->node->node) < MYEPSILON) {
     DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl);
     return true;
@@ -458 +465 @@
   int i = 0;
   do {
-    if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(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);
-      CrossPoint.SubtractVector(endpoints[i % 3]->node->node); // cross point was returned as absolute vector
-      const double s = CrossPoint.ScalarProduct(&helper) / helper.NormSquared();
+    try {
+      CrossPoint = GetIntersectionOfTwoLinesOnPlane(*(endpoints[i%3]->node->node),
+                                                    *(endpoints[(i+1)%3]->node->node),
+                                                    *(endpoints[(i+2)%3]->node->node),
+                                                    *Intersection);
+      helper = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
+      CrossPoint -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
+      const double s = CrossPoint.ScalarProduct(helper)/helper.NormSquared();
       DoLog(1) && (Log() << Verbose(1) << "INFO: Factor s is " << s << "." << endl);
-      if ((s < -MYEPSILON) || ((s - 1.) > MYEPSILON)) {
+      if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
         DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl);
-        i = 4;
+        i=4;
         break;
       }
       i++;
-    } else
+    } catch (LinearDependenceException &excp){
       break;
+    }
   } while (i < 3);
   if (i == 3) {
@@ -503 +514 @@
   DoLog(1) && (Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl);
   GetCenter(&Direction);
-  if (!ClosestPoint->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, x, &Direction)) {
-    ClosestPoint->CopyVector(x);
+  try {
+    *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*x, Direction);
+  }
+  catch (LinearDependenceException &excp) {
+    (*ClosestPoint) = (*x);
   }
 
   // 2. Calculate in plane part of line (x, intersection)
-  Vector InPlane;
-  InPlane.CopyVector(x);
-  InPlane.SubtractVector(ClosestPoint); // points from plane intersection to straight-down point
-  InPlane.ProjectOntoPlane(&NormalVector);
-  InPlane.AddVector(ClosestPoint);
+  Vector InPlane = (*x) - (*ClosestPoint); // points from plane intersection to straight-down point
+  InPlane.ProjectOntoPlane(NormalVector);
+  InPlane += *ClosestPoint;
 
   DoLog(2) && (Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl);
@@ -526 +538 @@
   for (int i = 0; i < 3; i++) {
     // treat direction of line as normal of a (cut)plane and the desired point x as the plane offset, the intersect line with point
-    Direction.CopyVector(endpoints[(i + 1) % 3]->node->node);
-    Direction.SubtractVector(endpoints[i % 3]->node->node);
+    Direction = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
     // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
-    CrossPoint[i].GetIntersectionWithPlane(&Direction, &InPlane, endpoints[i % 3]->node->node, endpoints[(i + 1) % 3]->node->node);
-    CrossDirection[i].CopyVector(&CrossPoint[i]);
-    CrossDirection[i].SubtractVector(&InPlane);
-    CrossPoint[i].SubtractVector(endpoints[i % 3]->node->node); // cross point was returned as absolute vector
-    const double s = CrossPoint[i].ScalarProduct(&Direction) / Direction.NormSquared();
+    CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(*(endpoints[i%3]->node->node), *(endpoints[(i+1)%3]->node->node));
+    CrossDirection[i] = CrossPoint[i] - InPlane;
+    CrossPoint[i] -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
+    const double s = CrossPoint[i].ScalarProduct(Direction)/Direction.NormSquared();
     DoLog(2) && (Log() << Verbose(2) << "INFO: Factor s is " << s << "." << endl);
-    if ((s >= -MYEPSILON) && ((s - 1.) <= MYEPSILON)) {
-      CrossPoint[i].AddVector(endpoints[i % 3]->node->node); // make cross point absolute again
+    if ((s >= -MYEPSILON) && ((s-1.) <= MYEPSILON)) {
+          CrossPoint[i] += (*endpoints[i%3]->node->node);  // make cross point absolute again
       DoLog(2) && (Log() << Verbose(2) << "INFO: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between " << *endpoints[i % 3]->node->node << " and " << *endpoints[(i + 1) % 3]->node->node << "." << endl);
-      const double distance = CrossPoint[i].DistanceSquared(x);
+      const double distance = CrossPoint[i].DistanceSquared(*x);
       if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
         ShortestDistance = distance;
-        ClosestPoint->CopyVector(&CrossPoint[i]);
+        (*ClosestPoint) = CrossPoint[i];
       }
     } else
@@ -548 +558 @@
   InsideFlag = true;
   for (int i = 0; i < 3; i++) {
-    const double sign = CrossDirection[i].ScalarProduct(&CrossDirection[(i + 1) % 3]);
-    const double othersign = CrossDirection[i].ScalarProduct(&CrossDirection[(i + 2) % 3]);
-    ;
+    const double sign = CrossDirection[i].ScalarProduct(CrossDirection[(i + 1) % 3]);
+    const double othersign = CrossDirection[i].ScalarProduct(CrossDirection[(i + 2) % 3]);
+
     if ((sign > -MYEPSILON) && (othersign > -MYEPSILON)) // have different sign
       InsideFlag = false;
   }
   if (InsideFlag) {
-    ClosestPoint->CopyVector(&InPlane);
-    ShortestDistance = InPlane.DistanceSquared(x);
+    (*ClosestPoint) = InPlane;
+    ShortestDistance = InPlane.DistanceSquared(*x);
   } else { // also check endnodes
     for (int i = 0; i < 3; i++) {
-      const double distance = x->DistanceSquared(endpoints[i]->node->node);
+      const double distance = x->DistanceSquared(*endpoints[i]->node->node);
       if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
         ShortestDistance = distance;
-        ClosestPoint->CopyVector(endpoints[i]->node->node);
+        (*ClosestPoint) = (*endpoints[i]->node->node);
       }
     }
@@ -667 +677 @@
   center->Zero();
   for (int i = 0; i < 3; i++)
-    center->AddVector(endpoints[i]->node->node);
+    (*center) += (*endpoints[i]->node->node);
   center->Scale(1. / 3.);
   DoLog(1) && (Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl);
@@ -733 +743 @@
   int counter = 0;
   for (; Runner[2] != endpoints.end();) {
-    TemporaryNormal.MakeNormalVector((*Runner[0])->node->node, (*Runner[1])->node->node, (*Runner[2])->node->node);
+    TemporaryNormal = Plane(*((*Runner[0])->node->node),
+                            *((*Runner[1])->node->node),
+                            *((*Runner[2])->node->node)).getNormal();
     for (int i = 0; i < 3; i++) // increase each of them
       Runner[i]++;
-    TotalNormal->AddVector(&TemporaryNormal);
+    (*TotalNormal) += TemporaryNormal;
   }
   TotalNormal->Scale(1. / (double) counter);
 
   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
-  if (TotalNormal->ScalarProduct(&OtherVector) > 0.)
+  if (TotalNormal->ScalarProduct(OtherVector) > 0.)
     TotalNormal->Scale(-1.);
   DoLog(1) && (Log() << Verbose(1) << "Normal Vector is " << *TotalNormal << "." << endl);
@@ -758 +770 @@
   center->Zero();
   int counter = 0;
-  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
-    center->AddVector((*Runner)->node->node);
+  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
+    (*center) += (*(*Runner)->node->node);
     counter++;
   }
@@ -1020 +1032 @@
   BaseLine(line), ThirdPoint(point), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI)
 {
-  Info FunctionInfo(__func__);
-  OptCenter.CopyVector(&OptCandidateCenter);
-  OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
-}
-;
+        Info FunctionInfo(__func__);
+  OptCenter = OptCandidateCenter;
+  OtherOptCenter = OtherOptCandidateCenter;
+};
+
 
 /** Destructor for class CandidateForTesselation.
@@ -1055 +1067 @@
   for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
     for (int i = 0; i < 2; i++) {
-      const double distance = fabs((*VRunner)->DistanceSquared(BaseLine->endpoints[i]->node->node) - radiusSquared);
+      const double distance = fabs((*VRunner)->DistanceSquared(*BaseLine->endpoints[i]->node->node) - radiusSquared);
       if (distance > HULLEPSILON) {
         DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << *BaseLine->endpoints[i] << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl);
@@ -1067 +1079 @@
     const TesselPoint *Walker = *Runner;
     for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
-      const double distance = fabs((*VRunner)->DistanceSquared(Walker->node) - radiusSquared);
+      const double distance = fabs((*VRunner)->DistanceSquared(*Walker->node) - radiusSquared);
       if (distance > HULLEPSILON) {
         DoeLog(1) && (eLog() << Verbose(1) << "Candidate " << *Walker << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl);
@@ -1085 +1097 @@
     DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << OtherOptCenter << ":" << endl);
     for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(&OtherOptCenter) << "." << endl);
+      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(OtherOptCenter) << "." << endl);
 
     // remove baseline's endpoints and candidates
@@ -1107 +1119 @@
     // check with animate_sphere.tcl VMD script
     if (ThirdPoint != NULL) {
-      DoLog(1) && (Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " " << ThirdPoint->Nr + 1 << " " << RADIUS << " " << OldCenter.x[0] << " " << OldCenter.x[1] << " " << OldCenter.x[2] << " " << (*VRunner)->x[0] << " " << (*VRunner)->x[1] << " " << (*VRunner)->x[2] << endl);
+      DoLog(1) && (Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " " << ThirdPoint->Nr + 1 << " " << RADIUS << " " << OldCenter[0] << " " << OldCenter[1] << " " << OldCenter[2] << " " << (*VRunner)->at(0) << " " << (*VRunner)->at(1) << " " << (*VRunner)->at(2) << endl);
     } else {
       DoLog(1) && (Log() << Verbose(1) << "Check by: ... missing third point ..." << endl);
-      DoLog(1) && (Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " ??? " << RADIUS << " " << OldCenter.x[0] << " " << OldCenter.x[1] << " " << OldCenter.x[2] << " " << (*VRunner)->x[0] << " " << (*VRunner)->x[1] << " " << (*VRunner)->x[2] << endl);
+      DoLog(1) && (Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " ??? " << RADIUS << " " << OldCenter[0] << " " << OldCenter[1] << " " << OldCenter[2] << " " << (*VRunner)->at(0) << " " << (*VRunner)->at(1) << " " << (*VRunner)->at(2) << endl);
     }
   }
@@ -1179 +1191 @@
   int num = 0;
   for (GoToFirst(); (!IsEnd()); GoToNext()) {
-    Center->AddVector(GetPoint()->node);
+    (*Center) += (*GetPoint()->node);
     num++;
   }
@@ -1296 +1308 @@
       C++;
       for (; C != PointsOnBoundary.end(); C++) {
-        tmp = A->second->node->node->DistanceSquared(B->second->node->node);
+        tmp = A->second->node->node->DistanceSquared(*B->second->node->node);
         distance = tmp * tmp;
-        tmp = A->second->node->node->DistanceSquared(C->second->node->node);
+        tmp = A->second->node->node->DistanceSquared(*C->second->node->node);
         distance += tmp * tmp;
-        tmp = B->second->node->node->DistanceSquared(C->second->node->node);
+        tmp = B->second->node->node->DistanceSquared(*C->second->node->node);
         distance += tmp * tmp;
         DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
@@ -1319 +1331 @@
     // 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);
+    PlaneVector = Plane(*A->second->node->node,
+                        *baseline->second.first->second->node->node,
+                        *baseline->second.second->second->node->node).getNormal();
     DoLog(2) && (Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl);
     // 4. loop over all points
@@ -1329 +1343 @@
         continue;
       // 4a. project onto plane vector
-      TrialVector.CopyVector(checker->second->node->node);
-      TrialVector.SubtractVector(A->second->node->node);
-      distance = TrialVector.ScalarProduct(&PlaneVector);
+      TrialVector = (*checker->second->node->node);
+      TrialVector.SubtractVector(*A->second->node->node);
+      distance = TrialVector.ScalarProduct(PlaneVector);
       if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
         continue;
@@ -1346 +1360 @@
       }
       // 4d. Check whether the point is inside the triangle (check distance to each node
-      tmp = checker->second->node->node->DistanceSquared(A->second->node->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)))
+      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)))
+      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)))
+      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
@@ -1435 +1449 @@
         // 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);
+        BaseLineCenter = 0.5 * ((*baseline->second->endpoints[0]->node->node) +
+                                (*baseline->second->endpoints[1]->node->node));
+        BaseLine = (*baseline->second->endpoints[0]->node->node) - (*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 += (*BTS->endpoints[i]->node->node);
         CenterVector.Scale(1. / 3.);
         DoLog(2) && (Log() << Verbose(2) << "CenterVector of base triangle is " << CenterVector << endl);
 
         // normal vector of triangle
-        NormalVector.CopyVector(Center);
-        NormalVector.SubtractVector(&CenterVector);
+        NormalVector = (*Center) - CenterVector;
         BTS->GetNormalVector(NormalVector);
-        NormalVector.CopyVector(&BTS->NormalVector);
+        NormalVector = BTS->NormalVector;
         DoLog(2) && (Log() << Verbose(2) << "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!
+        PropagationVector = Plane(BaseLine, NormalVector,0).getNormal();
+        TempVector = CenterVector - (*baseline->second->endpoints[0]->node->node); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
         //Log() << Verbose(0) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
-        if (PropagationVector.ScalarProduct(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
+        if (PropagationVector.ScalarProduct(TempVector) > 0) // make sure normal propagation vector points outward from baseline
           PropagationVector.Scale(-1.);
         DoLog(2) && (Log() << Verbose(2) << "PropagationVector of base triangle is " << PropagationVector << endl);
@@ -1472 +1482 @@
 
             // 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);
+            VirtualNormalVector = (*target->second->node->node) - BaseLineCenter;
+            VirtualNormalVector.ProjectOntoPlane(NormalVector);
+            TempAngle = VirtualNormalVector.Angle(PropagationVector);
             DoLog(2) && (Log() << Verbose(2) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl);
             if (TempAngle > (M_PI / 2.)) { // no bends bigger than Pi/2 (90 degrees)
@@ -1502 +1511 @@
 
             // 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);
+            TempVector = (*baseline->second->endpoints[0]->node->node) - (*target->second->node->node);
+            helper = (*baseline->second->endpoints[1]->node->node) - (*target->second->node->node);
+            helper.ProjectOntoPlane(TempVector);
             if (fabs(helper.NormSquared()) < MYEPSILON) {
               DoLog(2) && (Log() << Verbose(2) << "Chosen set of vectors is linear dependent." << endl);
@@ -1514 +1521 @@
             // 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);
+            VirtualNormalVector = Plane(*(baseline->second->endpoints[0]->node->node),
+                                        *(baseline->second->endpoints[1]->node->node),
+                                        *(target->second->node->node)).getNormal();
+            TempVector = (1./3.) * ((*baseline->second->endpoints[0]->node->node) +
+                                    (*baseline->second->endpoints[1]->node->node) +
+                                    (*target->second->node->node));
+            TempVector -= (*Center);
             // make it always point outward
-            if (VirtualNormalVector.ScalarProduct(&TempVector) < 0)
+            if (VirtualNormalVector.ScalarProduct(TempVector) < 0)
               VirtualNormalVector.Scale(-1.);
             // calculate angle
-            TempAngle = NormalVector.Angle(&VirtualNormalVector);
+            TempAngle = NormalVector.Angle(VirtualNormalVector);
             DoLog(2) && (Log() << Verbose(2) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl);
             if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner
@@ -1532 +1540 @@
             } 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);
+              helper = (*target->second->node->node) - 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);
+              TempVector = (*target->second->node->node) - BaseLineCenter;
+              TempVector.ProjectOntoPlane(VirtualNormalVector);
+              TempAngle = TempVector.Angle(helper);
+              TempVector = (*winner->second->node->node) - BaseLineCenter;
+              TempVector.ProjectOntoPlane(VirtualNormalVector);
+              if (TempAngle < TempVector.Angle(helper)) {
+                TempAngle = NormalVector.Angle(VirtualNormalVector);
                 SmallestAngle = TempAngle;
                 winner = target;
@@ -1581 +1586 @@
           BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
           BTS->GetCenter(&helper);
-          helper.SubtractVector(Center);
-          helper.Scale(-1);
+          helper -= (*Center);
+          helper *= -1;
           BTS->GetNormalVector(helper);
           TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
@@ -1640 +1645 @@
       DoLog(0) && (Log() << Verbose(0) << "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) {
+      if ((Center->DistanceSquared(*Walker->node) - Center->DistanceSquared(Intersection)) < -MYEPSILON) {
         // inside, next!
         DoLog(0) && (Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl);
@@ -1653 +1658 @@
           OldPoints[i] = BTS->endpoints[i];
         }
-        Normal.CopyVector(&BTS->NormalVector);
+        Normal = BTS->NormalVector;
         // add Walker to boundary points
         DoLog(0) && (Log() << Verbose(0) << "Adding " << *Walker << " to BoundaryPoints." << endl);
@@ -1800 +1805 @@
         // get open line
         for (TesselPointList::const_iterator CandidateChecker = Finder->second->pointlist.begin(); CandidateChecker != Finder->second->pointlist.end(); ++CandidateChecker) {
-          if ((*(CandidateChecker) == candidate->node) && (OptCenter == NULL || OptCenter->DistanceSquared(&Finder->second->OptCenter) < MYEPSILON )) { // stop searching if candidate matches
+          if ((*(CandidateChecker) == candidate->node) && (OptCenter == NULL || OptCenter->DistanceSquared(Finder->second->OptCenter) < MYEPSILON )) { // stop searching if candidate matches
             DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Candidate " << *(*CandidateChecker) << " has the right center " << Finder->second->OptCenter << "." << endl);
             insertNewLine = false;
@@ -2039 +2044 @@
   bool flag = true;
 
-  DoLog(1) && (Log() << Verbose(1) << "Check by: draw sphere {" << CandidateLine.OtherOptCenter.x[0] << " " << CandidateLine.OtherOptCenter.x[1] << " " << CandidateLine.OtherOptCenter.x[2] << "} radius " << RADIUS << " resolution 30" << endl);
+  DoLog(1) && (Log() << Verbose(1) << "Check by: draw sphere {" << CandidateLine.OtherOptCenter[0] << " " << CandidateLine.OtherOptCenter[1] << " " << CandidateLine.OtherOptCenter[2] << "} radius " << RADIUS << " resolution 30" << endl);
   // get all points inside the sphere
   TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, &CandidateLine.OtherOptCenter);
@@ -2045 +2050 @@
   DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":" << endl);
   for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-    DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(&CandidateLine.OtherOptCenter) << "." << endl);
+    DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(CandidateLine.OtherOptCenter) << "." << endl);
 
   // remove triangles's endpoints
@@ -2061 +2066 @@
     DoLog(1) && (Log() << Verbose(1) << "External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":" << endl);
     for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(&CandidateLine.OtherOptCenter) << "." << endl);
+      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(CandidateLine.OtherOptCenter) << "." << endl);
   }
   delete (ListofPoints);
@@ -2215 +2220 @@
         if (List != NULL) {
           for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            if ((*Runner)->node->x[i] > maxCoordinate[i]) {
+            if ((*Runner)->node->at(i) > maxCoordinate[i]) {
               DoLog(1) && (Log() << Verbose(1) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl);
-              maxCoordinate[i] = (*Runner)->node->x[i];
+              maxCoordinate[i] = (*Runner)->node->at(i);
               MaxPoint[i] = (*Runner);
             }
@@ -2235 +2240 @@
   for (int k = 0; k < NDIM; k++) {
     NormalVector.Zero();
-    NormalVector.x[k] = 1.;
+    NormalVector[k] = 1.;
     BaseLine = new BoundaryLineSet();
     BaseLine->endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
@@ -2253 +2258 @@
 
     // construct center of circle
-    CircleCenter.CopyVector(BaseLine->endpoints[0]->node->node);
-    CircleCenter.AddVector(BaseLine->endpoints[1]->node->node);
-    CircleCenter.Scale(0.5);
+    CircleCenter = 0.5 * ((*BaseLine->endpoints[0]->node->node) + (*BaseLine->endpoints[1]->node->node));
 
     // construct normal vector of circle
-    CirclePlaneNormal.CopyVector(BaseLine->endpoints[0]->node->node);
-    CirclePlaneNormal.SubtractVector(BaseLine->endpoints[1]->node->node);
+    CirclePlaneNormal = (*BaseLine->endpoints[0]->node->node) - (*BaseLine->endpoints[1]->node->node);
 
     double radius = CirclePlaneNormal.NormSquared();
     double CircleRadius = sqrt(RADIUS * RADIUS - radius / 4.);
 
-    NormalVector.ProjectOntoPlane(&CirclePlaneNormal);
+    NormalVector.ProjectOntoPlane(CirclePlaneNormal);
     NormalVector.Normalize();
     ShortestAngle = 2. * M_PI; // This will indicate the quadrant.
 
-    SphereCenter.CopyVector(&NormalVector);
-    SphereCenter.Scale(CircleRadius);
-    SphereCenter.AddVector(&CircleCenter);
+    SphereCenter = (CircleRadius * NormalVector) +  CircleCenter;
     // Now, NormalVector and SphereCenter are two orthonormalized vectors in the plane defined by CirclePlaneNormal (not normalized)
 
     // look in one direction of baseline for initial candidate
-    SearchDirection.MakeNormalVector(&CirclePlaneNormal, &NormalVector); // whether we look "left" first or "right" first is not important ...
+    SearchDirection = Plane(CirclePlaneNormal, NormalVector,0).getNormal();  // whether we look "left" first or "right" first is not important ...
 
     // adding point 1 and point 2 and add the line between them
@@ -2478 +2478 @@
 
   // construct center of circle
-  CircleCenter.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
-  CircleCenter.AddVector(CandidateLine.BaseLine->endpoints[1]->node->node);
-  CircleCenter.Scale(0.5);
+  CircleCenter = 0.5 * ((*CandidateLine.BaseLine->endpoints[0]->node->node) +
+                        (*CandidateLine.BaseLine->endpoints[1]->node->node));
 
   // construct normal vector of circle
-  CirclePlaneNormal.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
-  CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node);
+  CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
+                      (*CandidateLine.BaseLine->endpoints[1]->node->node);
 
   // calculate squared radius of circle
-  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+  radius = CirclePlaneNormal.ScalarProduct(CirclePlaneNormal);
   if (radius / 4. < RADIUS * RADIUS) {
     // construct relative sphere center with now known CircleCenter
-    RelativeSphereCenter.CopyVector(&T.SphereCenter);
-    RelativeSphereCenter.SubtractVector(&CircleCenter);
+    RelativeSphereCenter = T.SphereCenter - CircleCenter;
 
     CircleRadius = RADIUS * RADIUS - radius / 4.;
@@ -2500 +2498 @@
 
     // construct SearchDirection and an "outward pointer"
-    SearchDirection.MakeNormalVector(&RelativeSphereCenter, &CirclePlaneNormal);
-    helper.CopyVector(&CircleCenter);
-    helper.SubtractVector(ThirdPoint->node->node);
-    if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
+    SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal,0).getNormal();
+    helper = CircleCenter - (*ThirdPoint->node->node);
+    if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
       SearchDirection.Scale(-1.);
     DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
-    if (fabs(RelativeSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
+    if (fabs(RelativeSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
       DoeLog(1) && (eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl);
@@ -2643 +2640 @@
         Finder->second->pointlist.push_back(Sprinter);
         Finder->second->ShortestAngle = 0.;
-        Finder->second->OptCenter.CopyVector(OptCenter);
+        Finder->second->OptCenter = *OptCenter;
       }
     }
@@ -2678 +2675 @@
   // create normal vector
   BTS->GetCenter(&Center);
-  Center.SubtractVector(&CandidateLine.OptCenter);
-  BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
+  Center -= CandidateLine.OptCenter;
+  BTS->SphereCenter = CandidateLine.OptCenter;
   BTS->GetNormalVector(Center);
   // give some verbose output about the whole procedure
@@ -2712 +2709 @@
   AddTesselationTriangle();
 
-  BTS->SphereCenter.CopyVector(&CandidateLine.OtherOptCenter);
+  BTS->SphereCenter = CandidateLine.OtherOptCenter;
   // create normal vector in other direction
-  BTS->GetNormalVector(&triangle->NormalVector);
+  BTS->GetNormalVector(triangle->NormalVector);
   BTS->NormalVector.Scale(-1.);
   // give some verbose output about the whole procedure
@@ -2757 +2754 @@
   // create normal vector
   BTS->GetCenter(&Center);
-  Center.SubtractVector(OptCenter);
-  BTS->SphereCenter.CopyVector(OptCenter);
+  Center.SubtractVector(*OptCenter);
+  BTS->SphereCenter = *OptCenter;
   BTS->GetNormalVector(Center);
 
@@ -2802 +2799 @@
   Vector DistanceToIntersection[2], BaseLine;
   double distance[2];
-  BaseLine.CopyVector(Base->endpoints[1]->node->node);
-  BaseLine.SubtractVector(Base->endpoints[0]->node->node);
+  BaseLine = (*Base->endpoints[1]->node->node) - (*Base->endpoints[0]->node->node);
   for (int i = 0; i < 2; i++) {
-    DistanceToIntersection[i].CopyVector(ClosestPoint);
-    DistanceToIntersection[i].SubtractVector(Base->endpoints[i]->node->node);
-    distance[i] = BaseLine.ScalarProduct(&DistanceToIntersection[i]);
+    DistanceToIntersection[i] = (*ClosestPoint) - (*Base->endpoints[i]->node->node);
+    distance[i] = BaseLine.ScalarProduct(DistanceToIntersection[i]);
   }
   delete (ClosestPoint);
@@ -2883 +2878 @@
 
   // get the distance vector from Base line to OtherBase line
-  Vector Distance;
-  Distance.CopyVector(ClosestPoint[1]);
-  Distance.SubtractVector(ClosestPoint[0]);
+  Vector Distance = (*ClosestPoint[1]) - (*ClosestPoint[0]);
 
   // calculate volume
@@ -2906 +2899 @@
     }
     for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
-      DoLog(1) && (Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl);
-      BaseLineNormal.AddVector(&(runner->second->NormalVector));
+    DoLog(1) && (Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl);
+      BaseLineNormal += (runner->second->NormalVector);
     }
     BaseLineNormal.Scale(1. / 2.);
 
-    if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
+    if (Distance.ScalarProduct(BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
       DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl);
       // calculate volume summand as a general tetraeder
@@ -2947 +2940 @@
   for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
     DoLog(1) && (Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl);
-    BaseLineNormal.AddVector(&(runner->second->NormalVector));
+    BaseLineNormal += (runner->second->NormalVector);
   }
   BaseLineNormal.Scale(-1. / 2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
@@ -3087 +3080 @@
               double distance, scaleFactor;
 
-              OrthogonalizedOben.CopyVector(&Oben);
-              aCandidate.CopyVector(a->node);
-              aCandidate.SubtractVector(Candidate->node);
-              OrthogonalizedOben.ProjectOntoPlane(&aCandidate);
+              OrthogonalizedOben = Oben;
+              aCandidate = (*a->node) - (*Candidate->node);
+              OrthogonalizedOben.ProjectOntoPlane(aCandidate);
               OrthogonalizedOben.Normalize();
               distance = 0.5 * aCandidate.Norm();
@@ -3096 +3088 @@
               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);
+              Center = 0.5 * ((*Candidate->node) + (*a->node));
+              Center += OrthogonalizedOben;
+
+              AngleCheck = Center - (*a->node);
               norm = aCandidate.Norm();
               // second point shall have smallest angle with respect to Oben vector
               if (norm < RADIUS * 2.) {
-                angle = AngleCheck.Angle(&Oben);
+                angle = AngleCheck.Angle(Oben);
                 if (angle < Storage[0]) {
                   //Log() << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
@@ -3182 +3171 @@
 
   // copy old center
-  CandidateLine.OldCenter.CopyVector(&OldSphereCenter);
+  CandidateLine.OldCenter = OldSphereCenter;
   CandidateLine.ThirdPoint = ThirdPoint;
   CandidateLine.pointlist.clear();
 
   // construct center of circle
-  CircleCenter.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
-  CircleCenter.AddVector(CandidateLine.BaseLine->endpoints[1]->node->node);
-  CircleCenter.Scale(0.5);
+  CircleCenter = 0.5 * ((*CandidateLine.BaseLine->endpoints[0]->node->node) +
+                        (*CandidateLine.BaseLine->endpoints[1]->node->node));
 
   // construct normal vector of circle
-  CirclePlaneNormal.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
-  CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node);
-
-  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
-  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
+  CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
+                      (*CandidateLine.BaseLine->endpoints[1]->node->node);
+
+  RelativeOldSphereCenter = OldSphereCenter - CircleCenter;
 
   // calculate squared radius TesselPoint *ThirdPoint,f circle
@@ -3206 +3193 @@
 
     // test whether old center is on the band's plane
-    if (fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-      DoeLog(1) && (eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl);
-      RelativeOldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+    if (fabs(RelativeOldSphereCenter.ScalarProduct(CirclePlaneNormal)) > HULLEPSILON) {
+      DoeLog(1) && (eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(CirclePlaneNormal)) << "!" << endl);
+      RelativeOldSphereCenter.ProjectOntoPlane(CirclePlaneNormal);
     }
     radius = RelativeOldSphereCenter.NormSquared();
@@ -3216 +3203 @@
       // check SearchDirection
       DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
-      if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way!
+      if (fabs(RelativeOldSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) { // rotated the wrong way!
         DoeLog(1) && (eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl);
       }
@@ -3254 +3241 @@
                   DoLog(1) && (Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl);
 
-                  if (NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node) && (fabs(NewNormalVector.NormSquared()) > HULLEPSILON)) {
+                  try {
+                    NewNormalVector = Plane(*(CandidateLine.BaseLine->endpoints[0]->node->node),
+                                            *(CandidateLine.BaseLine->endpoints[1]->node->node),
+                                            *(Candidate->node)).getNormal();
                     DoLog(1) && (Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl);
-                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewPlaneCenter);
+                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(NewPlaneCenter);
                     DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);
                     DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
                     DoLog(1) && (Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl);
                     if (radius < RADIUS * RADIUS) {
-                      otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(&NewPlaneCenter);
+                      otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(NewPlaneCenter);
                       if (fabs(radius - otherradius) < HULLEPSILON) {
                         // construct both new centers
-                        NewSphereCenter.CopyVector(&NewPlaneCenter);
-                        OtherNewSphereCenter.CopyVector(&NewPlaneCenter);
-                        helper.CopyVector(&NewNormalVector);
+                        NewSphereCenter = NewPlaneCenter;
+                        OtherNewSphereCenter= NewPlaneCenter;
+                        helper = NewNormalVector;
                         helper.Scale(sqrt(RADIUS * RADIUS - radius));
                         DoLog(2) && (Log() << Verbose(2) << "INFO: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << "." << endl);
-                        NewSphereCenter.AddVector(&helper);
+                        NewSphereCenter += helper;
                         DoLog(2) && (Log() << 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 += helper;
                         DoLog(2) && (Log() << 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);
                         if ((ThirdPoint != NULL) && (Candidate == ThirdPoint->node)) { // in that case only the other circlecenter is valid
-                          if (OldSphereCenter.DistanceSquared(&NewSphereCenter) < OldSphereCenter.DistanceSquared(&OtherNewSphereCenter))
+                          if (OldSphereCenter.DistanceSquared(NewSphereCenter) < OldSphereCenter.DistanceSquared(OtherNewSphereCenter))
                             alpha = Otheralpha;
                         } else
                           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 (CandidateLine.ShortestAngle > (alpha - HULLEPSILON)) {
                           if (fabs(alpha - Otheralpha) > MYEPSILON) {
-                            CandidateLine.OptCenter.CopyVector(&NewSphereCenter);
-                            CandidateLine.OtherOptCenter.CopyVector(&OtherNewSphereCenter);
+                            CandidateLine.OptCenter = NewSphereCenter;
+                            CandidateLine.OtherOptCenter = OtherNewSphereCenter;
                           } else {
-                            CandidateLine.OptCenter.CopyVector(&OtherNewSphereCenter);
-                            CandidateLine.OtherOptCenter.CopyVector(&NewSphereCenter);
+                            CandidateLine.OptCenter = OtherNewSphereCenter;
+                            CandidateLine.OtherOptCenter = NewSphereCenter;
                           }
                           // if there is an equal candidate, add it to the list without clearing the list
@@ -3319 +3307 @@
                       DoLog(1) && (Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl);
                     }
-                  } else {
-                    DoLog(1) && (Log() << Verbose(1) << "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl);
+                  }
+                  catch (LinearDependenceException &excp){
+                    Log() << Verbose(1) << excp;
+                    Log() << Verbose(1) << "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
                   }
                 } else {
@@ -3416 +3406 @@
             FindPoint = PointsOnBoundary.find((*Runner)->nr);
             if (FindPoint != PointsOnBoundary.end()) {
-              points->insert(DistanceToPointPair(FindPoint->second->node->node->DistanceSquared(x), FindPoint->second));
+              points->insert(DistanceToPointPair(FindPoint->second->node->node->DistanceSquared(*x), FindPoint->second));
               DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl);
             }
@@ -3460 +3450 @@
     for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
       // calculate closest point on line to desired point
-      helper.CopyVector((LineRunner->second)->endpoints[0]->node->node);
-      helper.AddVector((LineRunner->second)->endpoints[1]->node->node);
-      helper.Scale(0.5);
-      Center.CopyVector(x);
-      Center.SubtractVector(&helper);
-      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
-      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
-      Center.ProjectOntoPlane(&BaseLine);
+      helper = 0.5 * ((*(LineRunner->second)->endpoints[0]->node->node) +
+                      (*(LineRunner->second)->endpoints[1]->node->node));
+      Center = (*x) - helper;
+      BaseLine = (*(LineRunner->second)->endpoints[0]->node->node) -
+                 (*(LineRunner->second)->endpoints[1]->node->node);
+      Center.ProjectOntoPlane(BaseLine);
       const double distance = Center.NormSquared();
       if ((ClosestLine == NULL) || (distance < MinDistance)) {
         // additionally calculate intersection on line (whether it's on the line section or not)
-        helper.CopyVector(x);
-        helper.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
-        helper.SubtractVector(&Center);
-        const double lengthA = helper.ScalarProduct(&BaseLine);
-        helper.CopyVector(x);
-        helper.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
-        helper.SubtractVector(&Center);
-        const double lengthB = helper.ScalarProduct(&BaseLine);
+        helper = (*x) - (*(LineRunner->second)->endpoints[0]->node->node) - Center;
+        const double lengthA = helper.ScalarProduct(BaseLine);
+        helper = (*x) - (*(LineRunner->second)->endpoints[1]->node->node) - Center;
+        const double lengthB = helper.ScalarProduct(BaseLine);
         if (lengthB * lengthA < 0) { // if have different sign
           ClosestLine = LineRunner->second;
@@ -3527 +3511 @@
     for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
 
-      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
-      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      BaseLine = (*(LineRunner->second)->endpoints[0]->node->node) -
+                 (*(LineRunner->second)->endpoints[1]->node->node);
       const double lengthBase = BaseLine.NormSquared();
 
-      BaseLineIntersection.CopyVector(x);
-      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
+      BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[0]->node->node);
       const double lengthEndA = BaseLineIntersection.NormSquared();
 
-      BaseLineIntersection.CopyVector(x);
-      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[1]->node->node);
       const double lengthEndB = BaseLineIntersection.NormSquared();
 
@@ -3554 +3536 @@
       } else { // intersection is closer, calculate
         // calculate closest point on line to desired point
-        BaseLineIntersection.CopyVector(x);
-        BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
-        Center.CopyVector(&BaseLineIntersection);
-        Center.ProjectOntoPlane(&BaseLine);
-        BaseLineIntersection.SubtractVector(&Center);
+        BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[1]->node->node);
+        Center = BaseLineIntersection;
+        Center.ProjectOntoPlane(BaseLine);
+        BaseLineIntersection -= Center;
         const double distance = BaseLineIntersection.NormSquared();
         if (Center.NormSquared() > BaseLine.NormSquared()) {
@@ -3612 +3593 @@
   for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
     (*Runner)->GetCenter(&Center);
-    helper.CopyVector(x);
-    helper.SubtractVector(&Center);
-    const double Alignment = helper.Angle(&(*Runner)->NormalVector);
+    helper = (*x) - Center;
+    const double Alignment = helper.Angle((*Runner)->NormalVector);
     if (Alignment < MinAlignment) {
       result = *Runner;
@@ -3681 +3661 @@
   triangle->GetCenter(&Center);
   DoLog(2) && (Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl);
-  DistanceToCenter.CopyVector(&Center);
-  DistanceToCenter.SubtractVector(&Point);
+  DistanceToCenter = Center - Point;
   DoLog(2) && (Log() << Verbose(2) << "INFO: Vector from point to test to center is " << DistanceToCenter << "." << endl);
 
   // check whether we are on boundary
-  if (fabs(DistanceToCenter.ScalarProduct(&triangle->NormalVector)) < MYEPSILON) {
+  if (fabs(DistanceToCenter.ScalarProduct(triangle->NormalVector)) < MYEPSILON) {
     // calculate whether inside of triangle
-    DistanceToCenter.CopyVector(&Point);
-    Center.CopyVector(&Point);
-    Center.SubtractVector(&triangle->NormalVector); // points towards MolCenter
-    DistanceToCenter.AddVector(&triangle->NormalVector); // points outside
+    DistanceToCenter = Point + triangle->NormalVector; // points outside
+    Center = Point - triangle->NormalVector; // points towards MolCenter
     DoLog(1) && (Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl);
     if (triangle->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
@@ -3706 +3683 @@
 
     // then check direction to boundary
-    if (DistanceToCenter.ScalarProduct(&triangle->NormalVector) > MYEPSILON) {
+    if (DistanceToCenter.ScalarProduct(triangle->NormalVector) > MYEPSILON) {
       DoLog(1) && (Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl);
       return -distance;
@@ -3840 +3817 @@
   if ((triangles != NULL) && (!triangles->empty())) {
     for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
-      PlaneNormal.AddVector(&(*Runner)->NormalVector);
+      PlaneNormal += (*Runner)->NormalVector;
   } else {
     DoeLog(0) && (eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl);
@@ -3851 +3828 @@
   // construct one orthogonal vector
   if (Reference != NULL) {
-    AngleZero.CopyVector(Reference);
-    AngleZero.SubtractVector(Point->node);
-    AngleZero.ProjectOntoPlane(&PlaneNormal);
+    AngleZero = (*Reference) - (*Point->node);
+    AngleZero.ProjectOntoPlane(PlaneNormal);
   }
   if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON)) {
     DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl);
-    AngleZero.CopyVector((*SetOfNeighbours->begin())->node);
-    AngleZero.SubtractVector(Point->node);
-    AngleZero.ProjectOntoPlane(&PlaneNormal);
+    AngleZero = (*(*SetOfNeighbours->begin())->node) - (*Point->node);
+    AngleZero.ProjectOntoPlane(PlaneNormal);
     if (AngleZero.NormSquared() < MYEPSILON) {
       DoeLog(0) && (eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
@@ -3867 +3842 @@
   DoLog(1) && (Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl);
   if (AngleZero.NormSquared() > MYEPSILON)
-    OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
+    OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
   else
-    OrthogonalVector.MakeNormalVector(&PlaneNormal);
+    OrthogonalVector.MakeNormalTo(PlaneNormal);
   DoLog(1) && (Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl);
 
   // go through all connected points and calculate angle
   for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
-    helper.CopyVector((*listRunner)->node);
-    helper.SubtractVector(Point->node);
-    helper.ProjectOntoPlane(&PlaneNormal);
+    helper = (*(*listRunner)->node) - (*Point->node);
+    helper.ProjectOntoPlane(PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
     DoLog(0) && (Log() << Verbose(0) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl);
@@ -3934 +3908 @@
   int counter = 0;
   while (TesselC != SetOfNeighbours->end()) {
-    helper.MakeNormalVector((*TesselA)->node, (*TesselB)->node, (*TesselC)->node);
+    helper = Plane(*((*TesselA)->node),
+                   *((*TesselB)->node),
+                   *((*TesselC)->node)).getNormal();
     DoLog(0) && (Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl);
     counter++;
@@ -3940 +3916 @@
     TesselB++;
     TesselC++;
-    PlaneNormal.AddVector(&helper);
+    PlaneNormal += helper;
   }
   //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
@@ -3955 +3931 @@
   // construct one orthogonal vector
   if (Reference != NULL) {
-    AngleZero.CopyVector(Reference);
-    AngleZero.SubtractVector(Point->node);
-    AngleZero.ProjectOntoPlane(&PlaneNormal);
-  }
-  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON)) {
+    AngleZero = (*Reference) - (*Point->node);
+    AngleZero.ProjectOntoPlane(PlaneNormal);
+  }
+  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {
     DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl);
-    AngleZero.CopyVector((*SetOfNeighbours->begin())->node);
-    AngleZero.SubtractVector(Point->node);
-    AngleZero.ProjectOntoPlane(&PlaneNormal);
+    AngleZero = (*(*SetOfNeighbours->begin())->node) - (*Point->node);
+    AngleZero.ProjectOntoPlane(PlaneNormal);
     if (AngleZero.NormSquared() < MYEPSILON) {
       DoeLog(0) && (eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
@@ -3971 +3945 @@
   DoLog(1) && (Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl);
   if (AngleZero.NormSquared() > MYEPSILON)
-    OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
+    OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
   else
-    OrthogonalVector.MakeNormalVector(&PlaneNormal);
+    OrthogonalVector.MakeNormalTo(PlaneNormal);
   DoLog(1) && (Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl);
 
@@ -3979 +3953 @@
   pair<map<double, TesselPoint*>::iterator, bool> InserterTest;
   for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
-    helper.CopyVector((*listRunner)->node);
-    helper.SubtractVector(Point->node);
-    helper.ProjectOntoPlane(&PlaneNormal);
+    helper = (*(*listRunner)->node) - (*Point->node);
+    helper.ProjectOntoPlane(PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
     if (angle > M_PI) // the correction is of no use here (and not desired)
@@ -4225 +4198 @@
 
   // copy old location for the volume
-  OldPoint.CopyVector(point->node->node);
+  OldPoint = (*point->node->node);
 
   // get list of connected points
@@ -4253 +4226 @@
   for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
     DoLog(1) && (Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->second) << "." << endl);
-    NormalVector.SubtractVector(&Runner->second->NormalVector); // has to point inward
+    NormalVector -= Runner->second->NormalVector; // has to point inward
     RemoveTesselationTriangle(Runner->second);
     count++;
@@ -4288 +4261 @@
           StartNode--;
           //Log() << Verbose(3) << "INFO: StartNode is " << **StartNode << "." << endl;
-          Point.CopyVector((*StartNode)->node);
-          Point.SubtractVector((*MiddleNode)->node);
+          Point = (*(*StartNode)->node) - (*(*MiddleNode)->node);
           StartNode = MiddleNode;
           StartNode++;
@@ -4295 +4267 @@
             StartNode = connectedPath->begin();
           //Log() << Verbose(3) << "INFO: EndNode is " << **StartNode << "." << endl;
-          Reference.CopyVector((*StartNode)->node);
-          Reference.SubtractVector((*MiddleNode)->node);
-          OrthogonalVector.CopyVector((*MiddleNode)->node);
-          OrthogonalVector.SubtractVector(&OldPoint);
-          OrthogonalVector.MakeNormalVector(&Reference);
+          Reference = (*(*StartNode)->node) - (*(*MiddleNode)->node);
+          OrthogonalVector = (*(*MiddleNode)->node) - OldPoint;
+          OrthogonalVector.MakeNormalTo(Reference);
           angle = GetAngle(Point, Reference, OrthogonalVector);
           //if (angle < M_PI)  // no wrong-sided triangles, please?
@@ -4752 +4722 @@
   class BoundaryLineSet *BestLine = NULL;
   for (LineMap::iterator Runner = NearestBoundaryPoint->lines.begin(); Runner != NearestBoundaryPoint->lines.end(); Runner++) {
-    BaseLine.CopyVector(Runner->second->endpoints[0]->node->node);
-    BaseLine.SubtractVector(Runner->second->endpoints[1]->node->node);
-    CenterToPoint.CopyVector(Runner->second->endpoints[0]->node->node);
-    CenterToPoint.AddVector(Runner->second->endpoints[1]->node->node);
-    CenterToPoint.Scale(0.5);
-    CenterToPoint.SubtractVector(point->node);
-    angle = CenterToPoint.Angle(&BaseLine);
-    if (fabs(angle - M_PI / 2.) < fabs(BestAngle - M_PI / 2.)) {
+    BaseLine = (*Runner->second->endpoints[0]->node->node) -
+               (*Runner->second->endpoints[1]->node->node);
+    CenterToPoint = 0.5 * ((*Runner->second->endpoints[0]->node->node) +
+                           (*Runner->second->endpoints[1]->node->node));
+    CenterToPoint -= (*point->node);
+    angle = CenterToPoint.Angle(BaseLine);
+    if (fabs(angle - M_PI/2.) < fabs(BestAngle - M_PI/2.)) {
       BestAngle = angle;
       BestLine = Runner->second;
@@ -4925 +4894 @@
       for (map<int, Vector *>::iterator VectorRunner = VectorWalker; VectorRunner != TriangleVectors.end(); VectorRunner++)
         if (VectorWalker != VectorRunner) { // skip equals
-          const double SCP = VectorWalker->second->ScalarProduct(VectorRunner->second); // ScalarProduct should result in -1. for degenerated triangles
+          const double SCP = VectorWalker->second->ScalarProduct(*VectorRunner->second); // ScalarProduct should result in -1. for degenerated triangles
           DoLog(1) && (Log() << Verbose(1) << "Checking " << *VectorWalker->second << " against " << *VectorRunner->second << ": " << SCP << endl);
           if (fabs(SCP + 1.) < ParallelEpsilon) {
@@ -5010 +4979 @@
     TriangleSet::iterator TriangleWalker = T->begin(); // is the inner iterator
     /// 4a. Get NormalVector for one side (this is "front")
-    NormalVector.CopyVector(&(*TriangleWalker)->NormalVector);
+    NormalVector = (*TriangleWalker)->NormalVector;
     DoLog(1) && (Log() << Verbose(1) << "\"front\" defining triangle is " << **TriangleWalker << " and Normal vector of \"front\" side is " << NormalVector << endl);
     TriangleWalker++;
@@ -5021 +4990 @@
       TriangleSprinter++;
       DoLog(1) && (Log() << Verbose(1) << "Current triangle to test for removal: " << *triangle << endl);
-      if (triangle->NormalVector.ScalarProduct(&NormalVector) < 0) { // if from other side, then delete and remove from list
+      if (triangle->NormalVector.ScalarProduct(NormalVector) < 0) { // if from other side, then delete and remove from list
         DoLog(1) && (Log() << Verbose(1) << " Removing ... " << endl);
         TriangleNrs.push(triangle->Nr);
@@ -5043 +5012 @@
       AddTesselationTriangle(); // ... and add
       TriangleNrs.pop();
-      BTS->NormalVector.CopyVector(&(*TriangleWalker)->NormalVector);
-      BTS->NormalVector.Scale(-1.);
+      BTS->NormalVector = -1 * (*TriangleWalker)->NormalVector;
       TriangleWalker++;
     }