diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index c9d9a01740cf7e77f5cd071da670f006e0c01803..a5d9f28d1cc7cbb35989c35feb558602e28d6fb7 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -21,6 +21,8 @@
#include <limits>
#include <list>
+#include "logog/include/logog.hpp"
+
// BaseLib
#include "quicksort.h"
@@ -169,8 +171,8 @@ bool lineSegmentIntersect(
}
bool lineSegmentsIntersect(const GeoLib::Polyline* ply,
- size_t &idx0,
- size_t &idx1,
+ size_t &idx0,
+ size_t &idx1,
GeoLib::Point& intersection_pnt)
{
size_t n_segs(ply->getNumberOfPoints() - 1);
@@ -197,17 +199,28 @@ bool lineSegmentsIntersect(const GeoLib::Polyline* ply,
return false;
}
-bool isPointInTriangle(const GeoLib::Point* p, const GeoLib::Point* a, const GeoLib::Point* b,
- const GeoLib::Point* c)
+bool isPointInTriangle(GeoLib::Point const& p,
+ GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c,
+ double eps_pnt_out_of_plane,
+ double eps_pnt_out_of_tri,
+ GeoLib::TriangleTest algorithm)
{
- return isPointInTriangle(*p, *a, *b, *c);
+ switch (algorithm)
+ {
+ case GeoLib::GAUSS:
+ return gaussPointInTriangle(p, a, b, c, eps_pnt_out_of_plane, eps_pnt_out_of_tri);
+ case GeoLib::BARYCENTRIC:
+ return barycentricPointInTriangle(p, a, b, c, eps_pnt_out_of_plane, eps_pnt_out_of_tri);
+ default:
+ ERR ("Selected algorithm for point in triangle testing not found, falling back on default.");
+ }
+ return gaussPointInTriangle(p, a, b, c, eps_pnt_out_of_plane, eps_pnt_out_of_tri);
}
-bool isPointInTriangle(GeoLib::Point const& q,
- GeoLib::Point const& a,
- GeoLib::Point const& b,
- GeoLib::Point const& c,
- double eps)
+bool gaussPointInTriangle(GeoLib::Point const& q,
+ GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c,
+ double eps_pnt_out_of_plane,
+ double eps_pnt_out_of_tri)
{
MathLib::Vector3 const v(a, b);
MathLib::Vector3 const w(a, c);
@@ -225,7 +238,7 @@ bool isPointInTriangle(GeoLib::Point const& q,
MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss(mat);
gauss.solve(y);
- const double lower (std::numeric_limits<float>::epsilon());
+ const double lower (eps_pnt_out_of_tri);
const double upper (1 + lower);
if (-lower <= y[0] && y[0] <= upper && -lower <= y[1] && y[1] <= upper && y[0] + y[1] <=
@@ -235,13 +248,38 @@ bool isPointInTriangle(GeoLib::Point const& q,
a[1] + y[0] * v[1] + y[1] * w[1],
a[2] + y[0] * v[2] + y[1] * w[2]
);
- if (MathLib::sqrDist(q, q_projected) < eps)
+ if (MathLib::sqrDist(q, q_projected) < eps_pnt_out_of_plane)
return true;
}
return false;
}
+bool barycentricPointInTriangle(GeoLib::Point const& p,
+ GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c,
+ double eps_pnt_out_of_plane,
+ double eps_pnt_out_of_tri)
+{
+ if (std::abs(orientation3d(p, a, b, c)) > eps_pnt_out_of_plane)
+ return false;
+
+ MathLib::Vector3 const pa (p,a);
+ MathLib::Vector3 const pb (p,b);
+ MathLib::Vector3 const pc (p,c);
+ double const area_x_2 (calcTriangleArea(a,b,c) * 2);
+
+ double const alpha ((MathLib::crossProduct(pb,pc).getLength()) / area_x_2);
+ if (alpha < -eps_pnt_out_of_tri || alpha > 1+eps_pnt_out_of_tri)
+ return false;
+ double const beta ((MathLib::crossProduct(pc,pa).getLength()) / area_x_2);
+ if (beta < -eps_pnt_out_of_tri || beta > 1+eps_pnt_out_of_tri)
+ return false;
+ double const gamma (1 - alpha - beta);
+ if (gamma < -eps_pnt_out_of_tri || gamma > 1+eps_pnt_out_of_tri)
+ return false;
+ return true;
+}
+
bool isPointInTetrahedron(GeoLib::Point const& p, GeoLib::Point const& a, GeoLib::Point const& b,
GeoLib::Point const& c, GeoLib::Point const& d, double eps)
{
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index 2c89448f72a7d97805c22f6baad7e77a49392418..892c333553144aed0b66ac7cf584cdcb1b08acf1 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -26,6 +26,11 @@ namespace GeoLib
{
class Polyline;
+enum TriangleTest
+{
+ GAUSS, BARYCENTRIC
+};
+
enum Orientation
{
CW = 1, CCW = 2, COLLINEAR = 3
@@ -111,22 +116,60 @@ double calcTriangleArea(GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::
*/
double calcTetrahedronVolume(const double* x1, const double* x2, const double* x3, const double* x4);
-bool isPointInTriangle (const GeoLib::Point* p,
- const GeoLib::Point* a, const GeoLib::Point* b, const GeoLib::Point* c);
-
/**
* Tests if the given point p is within the triangle, defined by its edge nodes a, b and c.
- * Using the eps it is possible to test a 'epsilon' neighbourhood around the triangle.
+ * Using two eps-values it is possible to test an 'epsilon' neighbourhood around the triangle
+ * as well as an 'epsilon' outside the triangles plane.
* @param p test point
* @param a edge node of triangle
* @param b edge node of triangle
* @param c edge node of triangle
- * @param eps size of neighbourhood (orthogonal distance to the plane spaned by triangle)
+ * @param eps_pnt_out_of_plane eps allowing for p to be slightly off the plane spanned by abc
+ * @param eps_pnt_out_of_tri eps allowing for p to be slightly off outside of abc
+ * @param algorithm defines the method to use
* @return true if the test point p is within the 'epsilon'-neighbourhood of the triangle
*/
bool isPointInTriangle(GeoLib::Point const& p,
- GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c,
- double eps = std::numeric_limits<float>::epsilon());
+ GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c,
+ double eps_pnt_out_of_plane = std::numeric_limits<float>::epsilon(),
+ double eps_pnt_out_of_tri = std::numeric_limits<float>::epsilon(),
+ GeoLib::TriangleTest algorithm = GeoLib::GAUSS);
+
+/**
+ * Tests if the given point p is within the triangle, defined by its edge nodes a, b and c.
+ * Using two eps-values it is possible to test an 'epsilon' neighbourhood around the triangle
+ * as well as an 'epsilon' outside the triangles plane.
+ * @param p test point
+ * @param a edge node of triangle
+ * @param b edge node of triangle
+ * @param c edge node of triangle
+ * @param eps_pnt_out_of_plane eps allowing for p to be slightly off the plane spanned by abc
+ * ((orthogonal distance to the plane spaned by triangle)
+ * @param eps_pnt_out_of_tri eps allowing for p to be slightly off outside of abc
+ * @return true if the test point p is within the 'epsilon'-neighbourhood of the triangle
+ */
+bool gaussPointInTriangle(GeoLib::Point const& p,
+ GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c,
+ double eps_pnt_out_of_plane = std::numeric_limits<float>::epsilon(),
+ double eps_pnt_out_of_tri = std::numeric_limits<float>::epsilon());
+
+/**
+ * Tests if the given point p is within the triangle, defined by its edge nodes a, b and c.
+ * Using two eps-values it is possible to test an 'epsilon' neighbourhood around the triangle
+ * as well as an 'epsilon' outside the triangles plane.
+ * Algorithm based on "Fundamentals of Computer Graphics" by Peter Shirley.
+ * @param p test point
+ * @param a edge node of triangle
+ * @param b edge node of triangle
+ * @param c edge node of triangle
+ * @param eps_pnt_out_of_plane eps allowing for p to be slightly off the plane spanned by abc
+ * @param eps_pnt_out_of_tri eps allowing for p to be slightly off outside of abc
+ * @return true if the test point p is within the 'epsilon'-neighbourhood of the triangle
+ */
+bool barycentricPointInTriangle(GeoLib::Point const& p,
+ GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c,
+ double eps_pnt_out_of_plane = std::numeric_limits<float>::epsilon(),
+ double eps_pnt_out_of_tri = std::numeric_limits<float>::epsilon());
/**
* Tests if the given point p is located within a tetrahedron spanned by points a, b, c, d.
diff --git a/GeoLib/EarClippingTriangulation.cpp b/GeoLib/EarClippingTriangulation.cpp
index eec05c6ba7960b7be1a6dc3c535f74003f636cef..b14fc0de6ed8e1bf1722203ce9b8e34288088105 100644
--- a/GeoLib/EarClippingTriangulation.cpp
+++ b/GeoLib/EarClippingTriangulation.cpp
@@ -119,7 +119,7 @@ bool EarClippingTriangulation::isEar(std::size_t v0, std::size_t v1, std::size_t
for (std::list<std::size_t>::const_iterator it (_vertex_list.begin ());
it != _vertex_list.end(); ++it) {
if (*it != v0 && *it != v1 && *it != v2) {
- if (isPointInTriangle (_pnts[*it], _pnts[v0], _pnts[v1], _pnts[v2])) {
+ if (isPointInTriangle (*_pnts[*it], *_pnts[v0], *_pnts[v1], *_pnts[v2])) {
return false;
}
}
diff --git a/Tests/GeoLib/TestIsPointInTriangle.cpp b/Tests/GeoLib/TestIsPointInTriangle.cpp
index 3171a674df4872e3569303a7d20c9be00981891d..ab1e9378216c2631154aeb0e9944dc7ebdd3905e 100644
--- a/Tests/GeoLib/TestIsPointInTriangle.cpp
+++ b/Tests/GeoLib/TestIsPointInTriangle.cpp
@@ -19,42 +19,60 @@ TEST(GeoLib, IsPointInTriangle)
GeoLib::Point const b(100.0, 0.0, 0.0);
GeoLib::Point const c(0.0, 100.0, 0.0);
+ double const default_eps (std::numeric_limits<double>::epsilon());
+
// check point on corner points of triangle
GeoLib::Point q(a);
- EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::gaussPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::barycentricPointInTriangle(q, a, b, c));
q = b;
- EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::gaussPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::barycentricPointInTriangle(q, a, b, c));
q = c;
- EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::gaussPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::barycentricPointInTriangle(q, a, b, c));
// check points on edges of triangle
q = GeoLib::Point(0.5*(b[0]+a[0]), 0.5*(b[1]+a[1]), 0.5*(b[2]+a[2]));
- EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::gaussPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::barycentricPointInTriangle(q, a, b, c));
q = GeoLib::Point(0.5*(c[0]+a[0]), 0.5*(c[1]+a[1]), 0.5*(c[2]+a[2]));
- EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::gaussPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::barycentricPointInTriangle(q, a, b, c));
q = GeoLib::Point(0.5*(c[0]+b[0]), 0.5*(c[1]+b[1]), 0.5*(c[2]+b[2]));
- EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::gaussPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::barycentricPointInTriangle(q, a, b, c));
// check points inside
q = GeoLib::Point (0.1, 0.1, 0.0);
- EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::gaussPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::barycentricPointInTriangle(q, a, b, c));
q = GeoLib::Point (0.1, 0.1, 1e-10);
- EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ EXPECT_TRUE(GeoLib::gaussPointInTriangle(q, a, b, c));
+ // here is a higher eps value needed for the second algorithm
+ EXPECT_TRUE(GeoLib::barycentricPointInTriangle(q, a, b, c, 1e-5));
// check points outside
q = GeoLib::Point (-0.1, 0.1, 0.0);
- EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c));
+ EXPECT_FALSE(GeoLib::gaussPointInTriangle(q, a, b, c));
+ EXPECT_FALSE(GeoLib::barycentricPointInTriangle(q, a, b, c));
q = GeoLib::Point (0.1, 0.1, 0.0005);
- EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c)); // default eps from float
+ EXPECT_FALSE(GeoLib::gaussPointInTriangle(q, a, b, c));
+ EXPECT_FALSE(GeoLib::barycentricPointInTriangle(q, a, b, c));
q = GeoLib::Point (0.1, 0.1, 0.0001);
- EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c, std::numeric_limits<double>::epsilon()));
+ EXPECT_FALSE(GeoLib::gaussPointInTriangle(q, a, b, c, std::numeric_limits<double>::epsilon()));
+ EXPECT_FALSE(GeoLib::barycentricPointInTriangle(q, a, b, c, std::numeric_limits<double>::epsilon()));
q = GeoLib::Point (0.1, 0.1, 0.000001);
- EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c, std::numeric_limits<double>::epsilon()));
+ EXPECT_FALSE(GeoLib::gaussPointInTriangle(q, a, b, c, std::numeric_limits<double>::epsilon()));
+ EXPECT_FALSE(GeoLib::barycentricPointInTriangle(q, a, b, c, std::numeric_limits<double>::epsilon()));
q = GeoLib::Point (0.1, 0.1, 1e-7);
- EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c, std::numeric_limits<double>::epsilon()));
+ EXPECT_FALSE(GeoLib::gaussPointInTriangle(q, a, b, c, std::numeric_limits<double>::epsilon()));
+ EXPECT_FALSE(GeoLib::barycentricPointInTriangle(q, a, b, c, std::numeric_limits<double>::epsilon()));
q = GeoLib::Point (0.1, 0.1, 0.001);
- EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c));
+ EXPECT_FALSE(GeoLib::gaussPointInTriangle(q, a, b, c));
+ EXPECT_FALSE(GeoLib::barycentricPointInTriangle(q, a, b, c));
q = GeoLib::Point (0.1, 0.1, 0.1);
- EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c));
+ EXPECT_FALSE(GeoLib::gaussPointInTriangle(q, a, b, c));
+ EXPECT_FALSE(GeoLib::barycentricPointInTriangle(q, a, b, c));
}