From 727a39823c4a6ed8e8b2686fa69f3a867d78b818 Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Thu, 28 Jul 2016 16:38:15 +0200
Subject: [PATCH] [MaL] Impl. of isPointInTriangleXY().

---
 MathLib/GeometricBasics.cpp | 25 +++++++++++++++++++++++++
 MathLib/GeometricBasics.h   |  7 +++++++
 2 files changed, 32 insertions(+)

diff --git a/MathLib/GeometricBasics.cpp b/MathLib/GeometricBasics.cpp
index a90c5bf86f0..255564779e7 100644
--- a/MathLib/GeometricBasics.cpp
+++ b/MathLib/GeometricBasics.cpp
@@ -169,6 +169,31 @@ bool barycentricPointInTriangle(MathLib::Point3d const& p,
     return true;
 }
 
+bool isPointInTriangleXY(MathLib::Point3d const& p,
+                         MathLib::Point3d const& a,
+                         MathLib::Point3d const& b,
+                         MathLib::Point3d const& c)
+{
+    // criterion: p-a = u0 * (b-a) + u1 * (c-a); 0 <= u0, u1 <= 1, u0+u1 <= 1
+    MathLib::DenseMatrix<double> mat(2, 2);
+    mat(0, 0) = b[0] - a[0];
+    mat(0, 1) = c[0] - a[0];
+    mat(1, 0) = b[1] - a[1];
+    mat(1, 1) = c[1] - a[1];
+    double y[2] = {p[0] - a[0], p[1] - a[1]};
+
+    MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss;
+    gauss.solve(mat, y, true);
+
+    // check if u0 and u1 fulfills the condition
+    if (0 <= y[0] && y[0] <= 1 && 0 <= y[1] && y[1] <= 1 && y[0] + y[1] <= 1)
+    {
+        return true;
+    }
+    return false;
+
+}
+
 bool dividedByPlane(const MathLib::Point3d& a, const MathLib::Point3d& b,
                     const MathLib::Point3d& c, const MathLib::Point3d& d)
 {
diff --git a/MathLib/GeometricBasics.h b/MathLib/GeometricBasics.h
index 80a5783e559..e87d8b553aa 100644
--- a/MathLib/GeometricBasics.h
+++ b/MathLib/GeometricBasics.h
@@ -154,6 +154,13 @@ bool barycentricPointInTriangle(
     double eps_pnt_out_of_plane = std::numeric_limits<float>::epsilon(),
     double eps_pnt_out_of_tri = std::numeric_limits<float>::epsilon());
 
+/// Checks if the point \f$p'\f$ is in the triangle defined by the points
+/// \f$a', b', c'\f$, where the \f$p', a', b', c' \f$ are the orthogonal
+/// projections to the \f$x\f$-\f$y\f$ plane of the points \f$p, a, b, c\f$,
+/// respectively.
+bool isPointInTriangleXY(MathLib::Point3d const& p, MathLib::Point3d const& a,
+MathLib::Point3d const& b, MathLib::Point3d const& c);
+
 /**
  * Checks if a and b can be placed on a plane such that c and d lie on different
  * sides of said plane. (In 2D space this checks if c and d are on different
-- 
GitLab