diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index 37771bc670bbc9f4d9a0b1dc7d929204b01cdd22..b8fa5b673470556943795b1fb9d4b272e3ad10c1 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -184,70 +184,59 @@ bool lineSegmentsIntersect(const GeoLib::Polyline* ply,
return false;
}
-static
-bool isPointInTriangle(const double p[3], const double a[3], const double b[3], const double c[3])
-{
- // criterion: p-b = u0 * (b - a) + u1 * (b - c); 0 <= u0, u1 <= 1, u0+u1 <= 1
- MathLib::DenseMatrix<double> mat(2, 2);
- mat(0, 0) = a[0] - b[0];
- mat(0, 1) = c[0] - b[0];
- mat(1, 0) = a[1] - b[1];
- mat(1, 1) = c[1] - b[1];
- double rhs[2] = { p[0] - b[0], p[1] - b[1] };
-
- MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss(mat);
- gauss.solve(rhs);
-
- if (0 <= rhs[0] && rhs[0] <= 1 && 0 <= rhs[1] && rhs[1] <= 1 && rhs[0] + rhs[1] <= 1)
- return true;
- return false;
-}
-
bool isPointInTriangle(const GeoLib::Point* p, const GeoLib::Point* a, const GeoLib::Point* b,
const GeoLib::Point* c)
{
- return isPointInTriangle(p->getCoords(), a->getCoords(), b->getCoords(), c->getCoords());
+ return isPointInTriangle(*p, *a, *b, *c);
}
-static
-double getOrientedTriArea(GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c)
+bool isPointInTriangle(GeoLib::Point const& q,
+ GeoLib::Point const& a,
+ GeoLib::Point const& b,
+ GeoLib::Point const& c,
+ double eps)
{
- const MathLib::Vector3 u(a,c);
- const MathLib::Vector3 v(a,b);
- const MathLib::Vector3 w(MathLib::crossProduct(u, v));
- return 0.5 * sqrt(MathLib::scalarProduct(w, w));
-}
+ MathLib::Vector3 const v(a, b);
+ MathLib::Vector3 const w(a, c);
-bool isPointInTriangle(GeoLib::Point const& p, GeoLib::Point const& a, GeoLib::Point const& b,
- GeoLib::Point const& c, double eps)
-{
- const unsigned dim(3);
- MathLib::DenseMatrix<double> m(dim, dim);
- for (unsigned i(0); i < dim; i++)
- m(i, 0) = b[i] - a[i];
- for (unsigned i(0); i < dim; i++)
- m(i, 1) = c[i] - a[i];
- for (unsigned i(0); i < dim; i++)
- m(i, 2) = p[i] - a[i];
-
- // point p is in the same plane as the triangle if and only if
- // the following determinate of the 3x3 matrix equals zero (up to an eps)
- double det3x3(m(0, 0) * (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2))
- - m(1, 0) * (m(2, 1) * m(0, 2) - m(0, 1) * m(2, 2))
- + m(2, 0) * (m(0, 1) * m(1, 2) - m(1, 1) * m(0, 2)));
- if (fabs(det3x3) > eps)
- return false;
+ MathLib::DenseMatrix<double> mat (2,2);
+ mat(0,0) = v.getSqrLength();
+ mat(0,1) = v[0]*w[0] + v[1]*w[1] + v[2]*w[2];
+ mat(1,0) = mat(0,1);
+ mat(1,1) = w.getSqrLength();
+ double y[2] = {
+ v[0] * (q[0] - a[0]) + v[1] * (q[1] - a[1]) + v[2] * (q[2] - a[2]),
+ w[0] * (q[0] - a[0]) + w[1] * (q[1] - a[1]) + w[2] * (q[2] - a[2])
+ };
- double total_area(getOrientedTriArea(a, b, c));
- double abp_area(getOrientedTriArea(a, b, p));
- double bcp_area(getOrientedTriArea(b, c, p));
- double cap_area(getOrientedTriArea(c, a, p));
+ MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss(mat);
+ gauss.solve(y);
+
+ const double lower (std::numeric_limits<float>::epsilon());
+ const double upper (1 + lower);
+
+ if (-lower <= y[0] && y[0] <= upper && -lower <= y[1] && y[1] <= upper && y[0] + y[1] <=
+ upper) {
+ GeoLib::Point const q_projected(
+ a[0] + y[0] * v[0] + y[1] * w[0],
+ 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)
+ return true;
+ }
- if (fabs(abp_area + bcp_area + cap_area - total_area) < eps)
- return true;
return false;
}
+double calcTriangleArea(GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c)
+{
+ MathLib::Vector3 const u(a,c);
+ MathLib::Vector3 const v(a,b);
+ MathLib::Vector3 const w(MathLib::crossProduct(u, v));
+ return 0.5 * w.getLength();
+}
+
// NewellPlane from book Real-Time Collision detection p. 494
void getNewellPlane(const std::vector<GeoLib::Point*>& pnts, MathLib::Vector3 &plane_normal, double& d)
{
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index 9aed682d823bf508076880a638387cc83735f726..1a6b1a48d9dd33275de281d3b8fc77d547cf23ae 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -98,12 +98,29 @@ void rotatePointsToXY(std::vector<GeoLib::Point*> &pnts);
*/
void rotatePointsToXZ(std::vector<GeoLib::Point*> &pnts);
+/**
+ * Calculates the area of the triangle defined by its edge nodes a, b and c..
+ * The formula is \f$A= \frac{1}{2} \cdot |u \times v|\f$, i.e. half of the area of the
+ * parallelogram specified by the vectors\f$u=b-a\f$ and \f$v=c-a\f.
+ */
+double calcTriangleArea(GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c);
+
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.
+ * @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)
+ * @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<double>::epsilon());
+ double eps = std::numeric_limits<float>::epsilon());
/**
* test for intersections of the line segments of the Polyline
diff --git a/GeoLib/Triangle.cpp b/GeoLib/Triangle.cpp
index 14ea90d525007cc60466d73dc450c49eba242c91..f915d3da8b4f05e5d2f125f0441c287b87cc5682 100644
--- a/GeoLib/Triangle.cpp
+++ b/GeoLib/Triangle.cpp
@@ -14,6 +14,9 @@
#include "Triangle.h"
+// GeoLib
+#include "AnalyticalGeometry.h"
+
// MathLib
#include "LinAlg/Solvers/GaussAlgorithm.h"
#include "MathTools.h"
@@ -59,106 +62,12 @@ void Triangle::setTriangle (size_t pnt_a, size_t pnt_b, size_t pnt_c)
_longest_edge = sqrt (_longest_edge);
}
-bool Triangle::containsPoint (Point const& pnt) const
+bool Triangle::containsPoint(Point const& q, double eps) const
{
- GeoLib::Point const& a_tmp (*(_pnts[_pnt_ids[0]]));
- GeoLib::Point const& b_tmp (*(_pnts[_pnt_ids[1]]));
- GeoLib::Point const& c_tmp (*(_pnts[_pnt_ids[2]]));
-
- GeoLib::Point s(a_tmp);
- for (size_t k(0); k<3; k++) {
- s[k] += b_tmp[k] + c_tmp[k];
- s[k] /= 3.0;
- }
-
- double eps (1e-2);
- GeoLib::Point const a (a_tmp[0] + eps *(a_tmp[0]-s[0]),
- a_tmp[1] + eps *(a_tmp[1]-s[1]),
- a_tmp[2] + eps *(a_tmp[2]-s[2]));
- GeoLib::Point const b (b_tmp[0] + eps *(b_tmp[0]-s[0]),
- b_tmp[1] + eps *(b_tmp[1]-s[1]),
- b_tmp[2] + eps *(b_tmp[2]-s[2]));
- GeoLib::Point const c (c_tmp[0] + eps *(c_tmp[0]-s[0]),
- c_tmp[1] + eps *(c_tmp[1]-s[1]),
- c_tmp[2] + eps *(c_tmp[2]-s[2]));
-
- const double delta (std::numeric_limits<double>::epsilon());
- const double upper (1+delta);
-
- // check special case where points of triangle have the same x-coordinate
- if (fabs(b[0]-a[0]) <= std::numeric_limits<double>::epsilon() &&
- fabs(c[0]-a[0]) <= std::numeric_limits<double>::epsilon()) {
- // all points of triangle have same x-coordinate
- if (fabs(pnt[0]-a[0]) / _longest_edge <= 1e-3) {
- // 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[1] - a[1];
- mat(0,1) = c[1] - a[1];
- mat(1,0) = b[2] - a[2];
- mat(1,1) = c[2] - a[2];
- double y[2] = {pnt[1]-a[1], pnt[2]-a[2]};
-
- MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss (mat);
- gauss.solve (y);
-
- if (-delta <= y[0] && y[0] <= upper && -delta <= y[1] && y[1] <= upper
- && y[0] + y[1] <= upper) {
- return true;
- } else {
- return false;
- }
- } else {
- return false;
- }
- }
-
- // check special case where points of triangle have the same y-coordinate
- if (fabs(b[1]-a[1]) <= std::numeric_limits<double>::epsilon() &&
- fabs(c[1]-a[1]) <= std::numeric_limits<double>::epsilon()) {
- // all points of triangle have same y-coordinate
- if (fabs(pnt[1]-a[1]) / _longest_edge <= 1e-3) {
- // 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[2] - a[2];
- mat(1,1) = c[2] - a[2];
- double y[2] = {pnt[0]-a[0], pnt[2]-a[2]};
-
- MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss (mat);
- gauss.solve (y);
-
- if (-delta <= y[0] && y[0] <= upper && -delta <= y[1] && y[1] <= upper && y[0] + y[1] <= upper) {
- return true;
- } else {
- return false;
- }
- } else {
- return false;
- }
- }
-
- // 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] = {pnt[0]-a[0], pnt[1]-a[1]};
-
- MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss (mat);
- gauss.solve (y);
-
- // check if the solution fulfills the third equation
- if (fabs((b[2]-a[2]) * y[0] + (c[2]-a[2]) * y[1] - (pnt[2] - a[2])) < 1e-3) {
- if (-delta <= y[0] && y[0] <= upper && -delta <= y[1] && y[1] <= upper &&
- y[0] + y[1] <= upper) {
- return true;
- }
- return false;
- } else {
- return false;
- }
+ GeoLib::Point const& a(*(_pnts[_pnt_ids[0]]));
+ GeoLib::Point const& b(*(_pnts[_pnt_ids[1]]));
+ GeoLib::Point const& c(*(_pnts[_pnt_ids[2]]));
+ return GeoLib::isPointInTriangle(q, a, b, c, eps);
}
bool Triangle::containsPoint2D (Point const& pnt) const
diff --git a/GeoLib/Triangle.h b/GeoLib/Triangle.h
index 3b0dcb07dfd741cd3256315b885b618519af8d89..caa14f01d0a956856e336e2c17d5ea55fa0adb8c 100644
--- a/GeoLib/Triangle.h
+++ b/GeoLib/Triangle.h
@@ -71,11 +71,12 @@ public:
}
/**
- * checks if point is in triangle
- * @param pnt The input point
+ * Checks if point q is within the triangle, uses GeoLib::isPointInTriangle().
+ * @param q The input point.
+ * @param eps Checks the 'epsilon'-neighbourhood
* @return true, if point is in triangle, else false
*/
- bool containsPoint (Point const& pnt) const;
+ bool containsPoint(Point const& q, double eps = std::numeric_limits<float>::epsilon()) const;
/**
* projects the triangle points to the x-y-plane and
diff --git a/MathLib/MathTools.cpp b/MathLib/MathTools.cpp
index 651fa75369b10a4218753378150339db1694d346..a086660e996d0f8a01256e1c90eb0ba8a49489c7 100644
--- a/MathLib/MathTools.cpp
+++ b/MathLib/MathTools.cpp
@@ -62,23 +62,6 @@ double getAngle (const double p0[3], const double p1[3], const double p2[3])
return acos (scalarProduct<double,3> (v0,v1) / (sqrt(scalarProduct<double,3>(v0,v0)) * sqrt(scalarProduct<double,3>(v1,v1))));
}
-double calcTriangleArea(const double* p0, const double* p1, const double* p2)
-{
- const double u0 (p2[0] - p0[0]);
- const double u1 (p2[1] - p0[1]);
- const double u2 (p2[2] - p0[2]);
-
- const double v0 (p1[0] - p0[0]);
- const double v1 (p1[1] - p0[1]);
- const double v2 (p1[2] - p0[2]);
-
- const double z0 (u1 * v2 - u2 * v1);
- const double z1 (u2 * v0 - u0 * v2);
- const double z2 (u0 * v1 - u1 * v0);
-
- return 0.5 * sqrt(z0 * z0 + z1 * z1 + z2 * z2);
-}
-
double calcTetrahedronVolume(const double* x1, const double* x2, const double* x3, const double* x4)
{
return fabs((x1[0] - x4[0]) * ((x2[1] - x4[1]) * (x3[2] - x4[2]) - (x2[2] - x4[2]) * (x3[1] - x4[1]))
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index 585d141b14c405b6143eeb70d2ac50a59fc559ab..8c78efca61690d9c1f4cc4dca30118855648d0af 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -148,12 +148,6 @@ float normalize(float min, float max, float val);
*/
double getAngle (const double p0[3], const double p1[3], const double p2[3]);
-/**
- * Calculates the area of a triangle.
- * The formula is A=.5*|u x v|, i.e. half of the area of the parallelogram specified by u=p0->p1 and v=p0->p2.
- */
-double calcTriangleArea(const double* p0, const double* p1, const double* p2);
-
/**
* Calculates the volume of a tetrahedron.
* The formula is V=1/6*|a(b x c)| with a=x1->x2, b=x1->x3 and c=x1->x4.
diff --git a/MeshLib/Elements/TemplateQuad-impl.h b/MeshLib/Elements/TemplateQuad-impl.h
index afd624856b42bf684cc3ca4911314694970c72c8..fff572636170247c3932311c97c6387247f765fb 100644
--- a/MeshLib/Elements/TemplateQuad-impl.h
+++ b/MeshLib/Elements/TemplateQuad-impl.h
@@ -90,8 +90,8 @@ TemplateQuad<NNODES,CELLQUADTYPE>::~TemplateQuad()
template <unsigned NNODES, CellType CELLQUADTYPE>
double TemplateQuad<NNODES,CELLQUADTYPE>::computeVolume()
{
- return MathLib::calcTriangleArea(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords())
- + MathLib::calcTriangleArea(_nodes[2]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords());
+ return GeoLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2])
+ + GeoLib::calcTriangleArea(*_nodes[2], *_nodes[3], *_nodes[0]);
}
template <unsigned NNODES, CellType CELLQUADTYPE>
diff --git a/MeshLib/Elements/TemplateTri.h b/MeshLib/Elements/TemplateTri.h
index a0c05f27d38dd4ccb0cf4a630e3eb3d951b56cb4..c91290193e21420ee7fa232d2bc5aab822f78fa7 100644
--- a/MeshLib/Elements/TemplateTri.h
+++ b/MeshLib/Elements/TemplateTri.h
@@ -22,6 +22,7 @@
#include "MeshEnums.h"
#include "MathTools.h"
+#include "AnalyticalGeometry.h"
namespace MeshLib {
@@ -121,7 +122,7 @@ protected:
/// Calculates the area of the triangle by returning half of the area of the corresponding parallelogram.
double computeVolume()
{
- return MathLib::calcTriangleArea(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords());
+ return GeoLib::calcTriangleArea(*_nodes[0], *_nodes[1], *_nodes[2]);
}
protected:
diff --git a/NumLib/Function/LinearInterpolationOnSurface.cpp b/NumLib/Function/LinearInterpolationOnSurface.cpp
index d71716a7b856ac5400e1ea1ccacd5f5b105e328d..18aa67f09ae1f3a19cf025b85722bd54756d1429 100644
--- a/NumLib/Function/LinearInterpolationOnSurface.cpp
+++ b/NumLib/Function/LinearInterpolationOnSurface.cpp
@@ -70,10 +70,10 @@ double LinearInterpolationOnSurface::interpolateInTri(const GeoLib::Triangle &tr
std::vector<GeoLib::Point*> p_pnts = {{&pnts[0], &pnts[1], &pnts[2]}};
GeoLib::rotatePointsToXY(p_pnts);
- const double* v1 = pnts[0].getCoords();
- const double* v2 = pnts[1].getCoords();
- const double* v3 = pnts[2].getCoords();
- const double area = MathLib::calcTriangleArea(v1, v2, v3);
+ GeoLib::Point const& v1(pnts[0]);
+ GeoLib::Point const& v2(pnts[1]);
+ GeoLib::Point const& v3(pnts[2]);
+ const double area = GeoLib::calcTriangleArea(v1, v2, v3);
if (area==.0) {
// take average if all points have the same coordinates
diff --git a/Tests/GeoLib/TestIsPointInTriangle.cpp b/Tests/GeoLib/TestIsPointInTriangle.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ce90557fa3a5b82c894d91f512649f7ed03dc1cf
--- /dev/null
+++ b/Tests/GeoLib/TestIsPointInTriangle.cpp
@@ -0,0 +1,60 @@
+/**
+ * @file TestIsPointInTriangle.cpp
+ * @date 2014-05-27
+ *
+ * @copyright
+ * Copyright (c) 2014, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include "gtest/gtest.h"
+
+#include "AnalyticalGeometry.h"
+
+TEST(GeoLib, IsPointInTriangle)
+{
+ GeoLib::Point const a(0.0, 0.0, 0.0);
+ GeoLib::Point const b(100.0, 0.0, 0.0);
+ GeoLib::Point const c(0.0, 100.0, 0.0);
+
+ // check point on corner points of triangle
+ GeoLib::Point q(a);
+ EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ q = b;
+ EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ q = c;
+ EXPECT_TRUE(GeoLib::isPointInTriangle(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));
+ 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));
+ 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));
+
+ // check points inside
+ q = GeoLib::Point (0.1, 0.1, 0.0);
+ EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+ q = GeoLib::Point (0.1, 0.1, 1e-10);
+ EXPECT_TRUE(GeoLib::isPointInTriangle(q, a, b, c));
+
+ // check points outside
+ q = GeoLib::Point (-0.1, 0.1, 0.0);
+ EXPECT_FALSE(GeoLib::isPointInTriangle(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
+ q = GeoLib::Point (0.1, 0.1, 0.0001);
+ EXPECT_FALSE(GeoLib::isPointInTriangle(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()));
+ q = GeoLib::Point (0.1, 0.1, 1e-7);
+ EXPECT_FALSE(GeoLib::isPointInTriangle(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));
+ q = GeoLib::Point (0.1, 0.1, 0.1);
+ EXPECT_FALSE(GeoLib::isPointInTriangle(q, a, b, c));
+}
+