Index: src/Line.cpp
===================================================================
--- src/Line.cpp	(revision 45ef76d8bee7dfcd142bef7a7cea8b1bc62fbafd)
+++ src/Line.cpp	(revision 643e7632e285a017dfb1b785b9f930ba145de897)
@@ -68,4 +68,13 @@
 }
 
+/** Calculates the intersection of the two lines that are both on the same plane.
+ * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html
+ * \param *out output stream for debugging
+ * \param *Line1a first vector of first line
+ * \param *Line1b second vector of first line
+ * \param *Line2a first vector of second line
+ * \param *Line2b second vector of second line
+ * \return true - \a this will contain the intersection on return, false - lines are parallel
+ */
 Vector Line::getIntersection(const Line& otherLine) const{
   Info FunctionInfo(__func__);
Index: src/tesselation.cpp
===================================================================
--- src/tesselation.cpp	(revision 45ef76d8bee7dfcd142bef7a7cea8b1bc62fbafd)
+++ src/tesselation.cpp	(revision 643e7632e285a017dfb1b785b9f930ba145de897)
@@ -17,4 +17,5 @@
 #include "triangleintersectionlist.hpp"
 #include "vector.hpp"
+#include "Line.hpp"
 #include "vector_ops.hpp"
 #include "verbose.hpp"
@@ -460,8 +461,7 @@
     int i = 0;
     do {
-      CrossPoint = GetIntersectionOfTwoLinesOnPlane(*(endpoints[i%3]->node->node),
-                                                    *(endpoints[(i+1)%3]->node->node),
-                                                    *(endpoints[(i+2)%3]->node->node),
-                                                    *Intersection);
+      Line line1 = makeLineThrough(*(endpoints[i%3]->node->node),*(endpoints[(i+1)%3]->node->node));
+      Line line2 = makeLineThrough(*(endpoints[(i+2)%3]->node->node),*Intersection);
+      CrossPoint = line1.getIntersection(line2);
       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
Index: src/tesselationhelpers.cpp
===================================================================
--- src/tesselationhelpers.cpp	(revision 45ef76d8bee7dfcd142bef7a7cea8b1bc62fbafd)
+++ src/tesselationhelpers.cpp	(revision 643e7632e285a017dfb1b785b9f930ba145de897)
@@ -14,4 +14,5 @@
 #include "tesselationhelpers.hpp"
 #include "vector.hpp"
+#include "Line.hpp"
 #include "vector_ops.hpp"
 #include "verbose.hpp"
@@ -666,7 +667,7 @@
   // calculate the intersection between this projected baseline and Base
   Vector *Intersection = new Vector;
-  *Intersection = GetIntersectionOfTwoLinesOnPlane(*(Base->endpoints[0]->node->node),
-                                                   *(Base->endpoints[1]->node->node),
-                                                     NewOffset, NewDirection);
+  Line line1 = makeLineThrough(*(Base->endpoints[0]->node->node),*(Base->endpoints[1]->node->node));
+  Line line2 = makeLineThrough(NewOffset, NewDirection);
+  *Intersection = line1.getIntersection(line2);
   Normal = (*Intersection) - (*Base->endpoints[0]->node->node);
   DoLog(1) && (Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(Baseline)/Baseline.NormSquared()) << "." << endl);
Index: src/vector_ops.cpp
===================================================================
--- src/vector_ops.cpp	(revision 45ef76d8bee7dfcd142bef7a7cea8b1bc62fbafd)
+++ src/vector_ops.cpp	(revision 643e7632e285a017dfb1b785b9f930ba145de897)
@@ -143,100 +143,2 @@
   return res;
 };
-
-/** Calculates the intersection of the two lines that are both on the same plane.
- * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html
- * \param *out output stream for debugging
- * \param *Line1a first vector of first line
- * \param *Line1b second vector of first line
- * \param *Line2a first vector of second line
- * \param *Line2b second vector of second line
- * \return true - \a this will contain the intersection on return, false - lines are parallel
- */
-Vector GetIntersectionOfTwoLinesOnPlane(const Vector &Line1a, const Vector &Line1b, const Vector &Line2a, const Vector &Line2b)
-{
-  Info FunctionInfo(__func__);
-
-  Vector res;
-
-  auto_ptr<GSLMatrix> M = auto_ptr<GSLMatrix>(new GSLMatrix(4,4));
-
-  M->SetAll(1.);
-  for (int i=0;i<3;i++) {
-    M->Set(0, i, Line1a[i]);
-    M->Set(1, i, Line1b[i]);
-    M->Set(2, i, Line2a[i]);
-    M->Set(3, i, Line2b[i]);
-  }
-
-  //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
-  //for (int i=0;i<4;i++) {
-  //  for (int j=0;j<4;j++)
-  //    cout << "\t" << M->Get(i,j);
-  //  cout << endl;
-  //}
-  if (fabs(M->Determinant()) > MYEPSILON) {
-    Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
-    throw SkewException(__FILE__,__LINE__);
-  }
-
-  Log() << Verbose(1) << "INFO: Line1a = " << Line1a << ", Line1b = " << Line1b << ", Line2a = " << Line2a << ", Line2b = " << Line2b << "." << endl;
-
-
-  // constuct a,b,c
-  Vector a = Line1b - Line1a;
-  Vector b = Line2b - Line2a;
-  Vector c = Line2a - Line1a;
-  Vector d = Line2b - Line1b;
-  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
-  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
-   res.Zero();
-   Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
-   throw LinearDependenceException(__FILE__,__LINE__);
-  }
-
-  // check for parallelity
-  Vector parallel;
-  double factor = 0.;
-  if (fabs(a.ScalarProduct(b)*a.ScalarProduct(b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
-    parallel = Line1a - Line2a;
-    factor = parallel.ScalarProduct(a)/a.Norm();
-    if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
-      res = Line2a;
-      Log() << Verbose(1) << "Lines conincide." << endl;
-      return res;
-    } else {
-      parallel = Line1a - Line2b;
-      factor = parallel.ScalarProduct(a)/a.Norm();
-      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
-        res = Line2b;
-        Log() << Verbose(1) << "Lines conincide." << endl;
-        return res;
-      }
-    }
-    Log() << Verbose(1) << "Lines are parallel." << endl;
-    res.Zero();
-    throw LinearDependenceException(__FILE__,__LINE__);
-  }
-
-  // obtain s
-  double s;
-  Vector temp1, temp2;
-  temp1 = c;
-  temp1.VectorProduct(b);
-  temp2 = a;
-  temp2.VectorProduct(b);
-  Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
-  if (fabs(temp2.NormSquared()) > MYEPSILON)
-    s = temp1.ScalarProduct(temp2)/temp2.NormSquared();
-  else
-    s = 0.;
-  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
-
-  // construct intersection
-  res = a;
-  res.Scale(s);
-  res += Line1a;
-  Log() << Verbose(1) << "Intersection is at " << res << "." << endl;
-
-  return res;
-};
Index: src/vector_ops.hpp
===================================================================
--- src/vector_ops.hpp	(revision 45ef76d8bee7dfcd142bef7a7cea8b1bc62fbafd)
+++ src/vector_ops.hpp	(revision 643e7632e285a017dfb1b785b9f930ba145de897)
@@ -11,5 +11,4 @@
 bool LSQdistance(Vector &res,const Vector **vectors, int num);
 Vector RotateVector(const Vector &vec,const Vector &axis, const double alpha);
-Vector GetIntersectionOfTwoLinesOnPlane(const Vector &Line1a, const Vector &Line1b, const Vector &Line2a, const Vector &Line2b);
 
 #endif /* VECTOR_OPS_HPP_ */
