MeshLib C++ Docs
Loading...
Searching...
No Matches
MRTriangleIntersection.h
Go to the documentation of this file.
1#pragma once
2
3#include "MRVector3.h"
5#include "MRTriPoint.h"
6
7#include <algorithm>
8#include <optional>
9
10namespace MR
11{
12
16
18{
19 // barycentric representation
20 TriPointf bary;
21 // distance from ray origin to p in dir length units
22 float t = 0;
23 TriIntersectResult(float U, float V, float dist)
24 {
25 bary.a = U; bary.b = V; t = dist;
26 }
27};
28
31template <typename T>
35)
36{
37 const auto abcd = mixed( a - d, b - d, c - d );
38 const auto abce = mixed( a - e, b - e, c - e );
39 const auto abcf = mixed( a - f, b - f, c - f );
40 const auto abc_de = abcd * abce >= 0; // segment DE is located at one side of the plane ABC
41 const auto abc_fd = abcf * abcd >= 0; // segment FD is located at one side of the plane ABC
42
43 if ( abc_de && abc_fd && abce * abcf >= 0 )
44 return false; // triangle DEF is located at one side of the plane ABC
45
46 const auto defa = mixed( d - a, e - a, f - a );
47 const auto defb = mixed( d - b, e - b, f - b );
48 const auto defc = mixed( d - c, e - c, f - c );
49 const auto def_ab = defa * defb >= 0; // segment AB is located at one side of the plane DEF
50 const auto def_ca = defc * defa >= 0; // segment CA is located at one side of the plane DEF
51
52 if ( def_ab && def_ca && defb * defc >= 0 )
53 return false; // triangle ABC is located at one side of the plane DEF
54
55 if ( abc_de )
56 std::swap( d, f );
57 else if( abc_fd )
58 std::swap( d, e );
59 // now segments DE and FD are crossed by the plane ABC: D at one side and EF at the other
60
61 if ( def_ab )
62 std::swap( a, c );
63 else if ( def_ca )
64 std::swap( a, b );
65 // now segments AB and CA are crossed by the plane DEF: A at one side and BC at the other
66
67 const auto abde = mixed( a - e, b - e, d - e );
68 const auto abdf = mixed( a - f, b - f, d - f );
69
70 if ( abde * abdf < 0 )
71 return true; // AB segment penetrates triangle DEF since points E and F are at distinct sides of ABD
72
73 const auto acde = mixed( a - e, c - e, d - e );
74
75 if ( abde * acde < 0 )
76 return true; // DE segment penetrates triangle ABC since points B and C are at distinct sides of ADE
77
78 if ( abdf == 0 && acde == 0 )
79 return true; // AB and DF segments are in the same plane, and AC and DE segments are in other same plane => triangles intersect, but no edge intersect the interior of other triangle
80
81 const auto acdf = mixed( a - f, c - f, d - f );
82
83 if ( acde * acdf < 0 )
84 return true; // AC segment penetrates triangle DEF since points E and F are at distinct sides of ACD
85
86 if ( abdf * acdf < 0 )
87 return true; // DF segment penetrates triangle ABC since points B and C are at distinct sides of ADF
88
89 if ( abde == 0 && acdf == 0 )
90 return true; // AB and DE segments are in the same plane, and AC and DF segments are in other same plane => triangles intersect, but no edge intersect the interior of other triangle
91
92 return false;
93}
94
96template<typename T>
97bool isPointInPlane( const Vector3<T>& p, const Vector3<T>& a, const Vector3<T>& b, const Vector3<T>& c )
98{
99 return mixed( p - a, p - b, p - c ) == T( 0 );
100}
101
103template<typename T>
104bool isPointInTriangle( const Vector3<T>& p, const Vector3<T>& a, const Vector3<T>& b, const Vector3<T>& c )
105{
106 if ( !isPointInPlane( p, a, b, c ) )
107 return false;
108 const auto normDir = cross( b - a, c - a );
109 if ( dot( normDir, cross( b - a, p - a ) ) < 0 )
110 return false;
111 if ( dot( normDir, cross( c - b, p - b ) ) < 0 )
112 return false;
113 if ( dot( normDir, cross( a - c, p - c ) ) < 0 )
114 return false;
115 return true;
116}
117
119template <typename T>
121 const Vector3<T> & x, const Vector3<T> & y, const Vector3<T> & z,
122 const Vector3<T> & u, const Vector3<T> & v, const Vector3<T> & w,
123 Vector3<T> d // approximate normal of the plane
124)
125{
126 const auto xy = ( y - x ).normalized();
127 d = ( d - xy * dot( xy, d ) ).normalized();
128 // now d is orthogonal to xy
129 const auto dz = dot( d, z - x );
130 return
131 dz * dot( d, u - x ) < 0 &&
132 dz * dot( d, v - x ) < 0 &&
133 dz * dot( d, w - x ) < 0;
134}
135
138template <typename T>
140 const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c,
141 const Vector3<T> & d, const Vector3<T> & e, const Vector3<T> & f
142)
143{
144 if ( !doTrianglesIntersect( a, b, c, d, e, f ) )
145 return false;
146
147 // direction from center to center
148 const auto dir = a + b + c - d - e - f;
149
150 return
151 !doesEdgeXySeparate( a, b, c, d, e, f, dir ) &&
152 !doesEdgeXySeparate( b, c, a, d, e, f, dir ) &&
153 !doesEdgeXySeparate( c, a, b, d, e, f, dir ) &&
154 !doesEdgeXySeparate( d, e, f, a, b, c, dir ) &&
155 !doesEdgeXySeparate( e, f, d, a, b, c, dir ) &&
156 !doesEdgeXySeparate( f, d, e, a, b, c, dir );
157}
158
160template <typename T>
162 const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c,
163 const Vector3<T> & d, const Vector3<T> & e
164)
165{
166 const auto dabe = mixed( d - e, a - e, b - e );
167 const auto dbce = mixed( d - e, b - e, c - e );
168 if ( dabe * dbce <= 0 )
169 return false; // segment AC is located at one side of the plane DEB
170
171 const auto dcae = mixed( d - e, c - e, a - e );
172 if ( dbce * dcae <= 0 )
173 return false; // segment AB is located at one side of the plane DEC
174
175 if ( dcae * dabe <= 0 )
176 return false; // segment BC is located at one side of the plane DEA
177
178 return true;
179}
180
182template <typename T>
184 const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c,
185 const Vector3<T> & d, const Vector3<T> & e
186)
187{
188 const auto abcd = mixed( a - d, b - d, c - d );
189 const auto abce = mixed( a - e, b - e, c - e );
190 if ( abcd * abce >= 0 )
191 return false; // segment DE is located at one side of the plane ABC
192 return doTriangleLineIntersect( a, b, c, d, e );
193}
194
196template <typename T>
198 const Vector3<T>& a, const Vector3<T>& b, const Vector3<T>& c,
199 const Vector3<T>& d, const Vector3<T>& e
200)
201{
202 const auto abcd = std::abs( mixed( a - d, b - d, c - d ) );
203 const auto abce = std::abs( mixed( a - e, b - e, c - e ) );
204 auto r = std::clamp( abcd / ( abcd + abce ), T( 0 ), T( 1 ) );
205 return r * e + ( 1 - r ) * d;
206}
207
208template <typename T>
209std::optional<TriIntersectResult> rayTriangleIntersect_( const Vector3<T>& oriA, const Vector3<T>& oriB, const Vector3<T>& oriC,
210 const IntersectionPrecomputes<T>& prec )
211{
212 const T& Sx = prec.Sx;
213 const T& Sy = prec.Sy;
214 const T& Sz = prec.Sz;
215
216 const T Ax = oriA[prec.idxX] - Sx * oriA[prec.maxDimIdxZ];
217 const T Ay = oriA[prec.idxY] - Sy * oriA[prec.maxDimIdxZ];
218 const T Bx = oriB[prec.idxX] - Sx * oriB[prec.maxDimIdxZ];
219 const T By = oriB[prec.idxY] - Sy * oriB[prec.maxDimIdxZ];
220 const T Cx = oriC[prec.idxX] - Sx * oriC[prec.maxDimIdxZ];
221 const T Cy = oriC[prec.idxY] - Sy * oriC[prec.maxDimIdxZ];
222
223 // due to fused multiply-add (FMA): (A*B-A*B) can be different from zero, so we need epsilon
224 const T eps = std::numeric_limits<T>::epsilon() * std::max( { Ax, Bx, Cx, Ay, By, Cy } );
225 T U = Cx * By - Cy * Bx;
226 T V = Ax * Cy - Ay * Cx;
227 T W = Bx * Ay - By * Ax;
228
229 if( U < -eps || V < -eps || W < -eps )
230 {
231 if( U > eps || V > eps || W > eps )
232 {
233 // U,V,W have clearly different signs, so the ray misses the triangle
234 return std::nullopt;
235 }
236 }
237
238 T det = U + V + W;
239 if( det == T( 0 ) )
240 return std::nullopt;
241 const T Az = Sz * oriA[prec.maxDimIdxZ];
242 const T Bz = Sz * oriB[prec.maxDimIdxZ];
243 const T Cz = Sz * oriC[prec.maxDimIdxZ];
244 const T t = U * Az + V * Bz + W * Cz;
245
246 auto invDet = T( 1 ) / det;
247 return TriIntersectResult( float( V * invDet ), float( W * invDet ), float( t * invDet ) );
248}
249
250inline std::optional<TriIntersectResult> rayTriangleIntersect( const Vector3f& oriA, const Vector3f& oriB, const Vector3f& oriC,
252{
253 return rayTriangleIntersect_( oriA, oriB, oriC, prec );
254}
255inline std::optional<TriIntersectResult> rayTriangleIntersect( const Vector3f& oriA, const Vector3f& oriB, const Vector3f& oriC,
256 const Vector3f& dir )
257{
258 const IntersectionPrecomputes<float> prec( dir );
259 return rayTriangleIntersect_( oriA, oriB, oriC, prec );
260}
261
262inline std::optional<TriIntersectResult> rayTriangleIntersect( const Vector3d& oriA, const Vector3d& oriB, const Vector3d& oriC,
264{
265 return rayTriangleIntersect_( oriA, oriB, oriC, prec );
266}
267
268inline std::optional<TriIntersectResult> rayTriangleIntersect( const Vector3d& oriA, const Vector3d& oriB, const Vector3d& oriC,
269 const Vector3d& dir )
270{
271 const IntersectionPrecomputes<double> prec( dir );
272 return rayTriangleIntersect_( oriA, oriB, oriC, prec );
273}
274
276
277} // namespace MR
int idxY
Definition MRIntersectionPrecomputes.h:124
T Sx
precomputed factors
Definition MRIntersectionPrecomputes.h:130
int idxX
Definition MRIntersectionPrecomputes.h:123
T Sz
Definition MRIntersectionPrecomputes.h:130
T Sy
Definition MRIntersectionPrecomputes.h:130
int maxDimIdxZ
Definition MRIntersectionPrecomputes.h:122
bool isPointInPlane(const Vector3< T > &p, const Vector3< T > &a, const Vector3< T > &b, const Vector3< T > &c)
returns true if ABC plane contains point P
Definition MRTriangleIntersection.h:97
bool doTrianglesIntersect(Vector3< T > a, Vector3< T > b, Vector3< T > c, Vector3< T > d, Vector3< T > e, Vector3< T > f)
Definition MRTriangleIntersection.h:32
bool doesEdgeXySeparate(const Vector3< T > &x, const Vector3< T > &y, const Vector3< T > &z, const Vector3< T > &u, const Vector3< T > &v, const Vector3< T > &w, Vector3< T > d)
returns true if a plane containing edge XY separates point Z from triangle UVW
Definition MRTriangleIntersection.h:120
std::optional< TriIntersectResult > rayTriangleIntersect(const Vector3f &oriA, const Vector3f &oriB, const Vector3f &oriC, const IntersectionPrecomputes< float > &prec)
Definition MRTriangleIntersection.h:250
bool isPointInTriangle(const Vector3< T > &p, const Vector3< T > &a, const Vector3< T > &b, const Vector3< T > &c)
returns true if ABC triangle contains point P
Definition MRTriangleIntersection.h:104
std::optional< TriIntersectResult > rayTriangleIntersect_(const Vector3< T > &oriA, const Vector3< T > &oriB, const Vector3< T > &oriC, const IntersectionPrecomputes< T > &prec)
Definition MRTriangleIntersection.h:209
bool doTriangleLineIntersect(const Vector3< T > &a, const Vector3< T > &b, const Vector3< T > &c, const Vector3< T > &d, const Vector3< T > &e)
checks whether triangle ABC and infinite line DE intersect
Definition MRTriangleIntersection.h:161
Vector3< T > findTriangleSegmentIntersection(const Vector3< T > &a, const Vector3< T > &b, const Vector3< T > &c, const Vector3< T > &d, const Vector3< T > &e)
this function input should have intersection
Definition MRTriangleIntersection.h:197
bool doTrianglesIntersectExt(const Vector3< T > &a, const Vector3< T > &b, const Vector3< T > &c, const Vector3< T > &d, const Vector3< T > &e, const Vector3< T > &f)
Definition MRTriangleIntersection.h:139
MRMESH_API TriangleSegmentIntersectResult doTriangleSegmentIntersect(const std::array< PreciseVertCoords, 5 > &vs)
Definition MRMesh/MRMeshFwd.h:429
Definition MRTriangleIntersection.h:18
TriIntersectResult(float U, float V, float dist)
Definition MRTriangleIntersection.h:23
float t
Definition MRTriangleIntersection.h:22
TriPointf bary
Definition MRTriangleIntersection.h:20
Definition MRMesh/MRVector3.h:26