MeshLib C++ Docs
Loading...
Searching...
No Matches
MRTriMath.h
Go to the documentation of this file.
1#pragma once
2// triangle-related mathematical functions are here
3
4#include "MRVector2.h"
5#include "MRVector3.h"
6#include <algorithm>
7#include <cassert>
8#include <cmath>
9#include <limits>
10#include <optional>
11
12namespace MR
13{
14
17template <typename T>
18[[nodiscard]] T circumcircleDiameterSq( const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c )
19{
20 const auto ab = ( b - a ).lengthSq();
21 const auto ca = ( a - c ).lengthSq();
22 const auto bc = ( c - b ).lengthSq();
23 if ( ab <= 0 )
24 return ca;
25 if ( ca <= 0 )
26 return bc;
27 if ( bc <= 0 )
28 return ab;
29 const auto f = cross( b - a, c - a ).lengthSq();
30 if ( f <= 0 )
31 return std::numeric_limits<T>::infinity();
32 return ab * ca * bc / f;
33}
34
37template <typename T>
38[[nodiscard]] inline T circumcircleDiameter( const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c )
39{
40 return std::sqrt( circumcircleDiameterSq( a, b, c ) );
41}
42
44template <typename T>
45[[nodiscard]] Vector3<T> circumcircleCenter( const Vector3<T> & a, const Vector3<T> & b )
46{
47 const auto xabSq = cross( a, b ).lengthSq();
48 const auto aa = a.lengthSq();
49 const auto bb = b.lengthSq();
50 if ( xabSq <= 0 )
51 {
52 if ( aa <= 0 )
53 return b / T(2);
54 // else b == 0 || a == b
55 return a / T(2);
56 }
57 const auto ab = dot( a, b );
58 return ( bb * ( aa - ab ) * a + aa * ( bb - ab ) * b ) / ( 2 * xabSq );
59}
60
62template <typename T>
63[[nodiscard]] inline Vector3<T> circumcircleCenter( const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c )
64{
65 return circumcircleCenter( a - c, b - c ) + c;
66}
67
70template <typename T>
71[[nodiscard]] bool circumballCenters( const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c, T radius,
72 Vector3<T> & centerPos, // ball's center from the positive side of triangle
73 Vector3<T> & centerNeg )// ball's center from the negative side of triangle
74{
75 const auto rr = sqr( radius );
76 const auto circRadSq = circumcircleDiameterSq( a, b, c ) / T( 4 );
77 if ( rr < circRadSq )
78 return false;
79
80 const auto x = std::sqrt( rr - circRadSq );
81 const auto xn = x * normal( a, b, c );
82 const auto circCenter = circumcircleCenter( a, b, c );
83 centerPos = circCenter + xn;
84 centerNeg = circCenter - xn;
85
86 return true;
87}
88
91template <typename T>
92[[nodiscard]] T minTriangleAngleSin( const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c )
93{
94 const auto ab = ( b - a ).length();
95 const auto ca = ( a - c ).length();
96 const auto bc = ( c - b ).length();
97 if ( ab <= 0 || ca <= 0 || bc <= 0 )
98 return 0;
99 const auto f = cross( b - a, c - a ).length();
100 return f * std::min( { ab, ca, bc } ) / ( ab * ca * bc );
101}
102
103template <typename T>
104[[nodiscard]] T minTriangleAngle( const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c )
105{
106 return std::asin( minTriangleAngleSin( a, b, c ) );
107}
108
111template<typename T>
112[[nodiscard]] T triangleAspectRatio( const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c )
113{
114 const auto bc = ( c - b ).length();
115 const auto ca = ( a - c ).length();
116 const auto ab = ( b - a ).length();
117 auto halfPerimeter = ( bc + ca + ab ) / 2;
118 auto den = 8 * ( halfPerimeter - bc ) * ( halfPerimeter - ca ) * ( halfPerimeter - ab );
119 if ( den <= 0 )
120 return std::numeric_limits<T>::max();
121
122 return bc * ca * ab / den;
123}
124
126template<typename T>
127[[nodiscard]] inline Vector3<T> dirDblArea( const Triangle3<T> & t )
128{
129 return cross( t[1] - t[0], t[2] - t[0] );
130}
131
133template<typename T>
134[[nodiscard]] inline Vector3<T> dirDblArea( const Vector3<T> & q, const Vector3<T> & r )
135{
136 return cross( q, r );
137}
138
140template<typename T>
141[[nodiscard]] inline Vector3<T> dirDblArea( const Vector3<T> & p, const Vector3<T> & q, const Vector3<T> & r )
142{
143 return cross( q - p, r - p );
144}
145
147template<typename T>
148[[nodiscard]] inline Vector3<T> normal( const Vector3<T> & q, const Vector3<T> & r )
149{
150 return dirDblArea( q, r ).normalized();
151}
152
154template<typename T>
155[[nodiscard]] inline Vector3<T> normal( const Vector3<T> & p, const Vector3<T> & q, const Vector3<T> & r )
156{
157 return normal( q - p, r - p );
158}
159
161template<typename T>
162[[nodiscard]] inline T dblAreaSq( const Vector3<T> & p, const Vector3<T> & q, const Vector3<T> & r )
163{
164 return dirDblArea( p, q, r ).lengthSq();
165}
166
168template<typename T>
169[[nodiscard]] inline T dblArea( const Triangle3<T> & t )
170{
171 return dirDblArea( t ).length();
172}
173
175template<typename T>
176[[nodiscard]] inline T dblArea( const Vector3<T> & p, const Vector3<T> & q, const Vector3<T> & r )
177{
178 return dirDblArea( p, q, r ).length();
179}
180
182template<typename T>
183[[nodiscard]] inline T area( const Vector3<T> & p, const Vector3<T> & q, const Vector3<T> & r )
184{
185 return dblArea( p, q, r ) / 2;
186}
187
189template<typename T>
190[[nodiscard]] inline T dblArea( const Vector2<T> & p, const Vector2<T> & q, const Vector2<T> & r )
191{
192 return std::abs( cross( q - p, r - p ) );
193}
194
196template<typename T>
197[[nodiscard]] inline T area( const Vector2<T> & p, const Vector2<T> & q, const Vector2<T> & r )
198{
199 return dblArea( p, q, r ) / 2;
200}
201
203template <typename T>
204[[nodiscard]] Triangle3<T> makeDegenerate( const Triangle3<T> & t )
205{
206 const auto c = ( t[0] + t[1] + t[2] ) / T(3);
207 int longest = 0;
208 T longestSq = 0;
209 for ( int i = 0; i < 3; ++i )
210 {
211 const auto sq = ( t[i] - c ).lengthSq();
212 if ( longestSq >= sq )
213 continue;
214 longest = i;
215 longestSq = sq;
216 }
217 const auto d = ( t[longest] - c ).normalized();
218
219 // project triangle on the line (c, d)
220 Triangle3<T> res;
221 for ( int i = 0; i < 3; ++i )
222 res[i] = c + d * dot( d, t[i] - c );
223 return res;
224}
225
228template <typename T>
229[[nodiscard]] Triangle3<T> triangleWithNormal( const Triangle3<T> & t, const Vector3<T> & n )
230{
231 const auto c = ( t[0] + t[1] + t[2] ) / T(3);
232 Triangle3<T> res;
233 for ( int i = 0; i < 3; ++i )
234 res[i] = t[i] - n * dot( n, t[i] - c );
235
236 if ( dot( n, dirDblArea( res ) ) < 0 ) // projected triangle has inversed normal
237 res = makeDegenerate( res );
238 return res;
239}
240
246template <typename T>
247[[nodiscard]] T dihedralAngleSin( const Vector3<T>& leftNorm, const Vector3<T>& rightNorm, const Vector3<T>& edgeVec )
248{
249 auto edgeDir = edgeVec.normalized();
250 return dot( edgeDir, cross( leftNorm, rightNorm ) );
251}
252
258template <typename T>
259[[nodiscard]] T dihedralAngleCos( const Vector3<T>& leftNorm, const Vector3<T>& rightNorm )
260{
261 return dot( leftNorm, rightNorm );
262}
263
270template <typename T>
271[[nodiscard]] T dihedralAngle( const Vector3<T>& leftNorm, const Vector3<T>& rightNorm, const Vector3<T>& edgeVec )
272{
273 auto sin = dihedralAngleSin( leftNorm, rightNorm, edgeVec );
274 auto cos = dihedralAngleCos( leftNorm, rightNorm );
275 return std::atan2( sin, cos );
276}
277
281template <typename T>
282[[nodiscard]] std::optional<Vector2<T>> posFromTriEdgeLengths( T a, T b, T c )
283{
284 if ( c == 0 )
285 {
286 if ( a == b )
287 return Vector2<T>{ a, 0 };
288 else
289 return {};
290 }
291 const auto aa = sqr( a );
292 const auto y = ( aa - sqr( b ) + sqr( c ) ) / ( 2 * c );
293 const auto yy = sqr( y );
294 if ( yy > aa )
295 return {};
296 const auto x = std::sqrt( aa - yy );
297 return Vector2<T>{ x, y };
298}
299
303template <typename T>
304[[nodiscard]] std::optional<T> quadrangleOtherDiagonal( T a, T b, T c, T a1, T b1 )
305{
306 const auto p = posFromTriEdgeLengths( a, b, c );
307 if ( !p )
308 return {};
309 auto p1 = posFromTriEdgeLengths( a1, b1, c );
310 if ( !p1 )
311 return {};
312 p1->x = -p1->x;
313 //where the other diagonal crosses axis Oy
314 auto y = ( p->x * p1->y - p1->x * p->y ) / ( p->x - p1->x );
315 if ( y < 0 || y > c )
316 return {};
317 return ( *p - *p1 ).length();
318}
319
323template <typename T>
324[[nodiscard]] inline T tanSqOfHalfAngle( T a, T b, T c )
325{
326 const T den = ( a + b + c ) * ( b + c - a );
327 if ( den <= 0 )
328 return std::numeric_limits<T>::infinity();
329 const T num = ( a + c - b ) * ( a + b - c );
330 if ( num <= 0 )
331 return 0;
332 return num / den;
333}
334
337template <typename T>
338[[nodiscard]] inline T cotan( const Triangle3<T> & t, T absMaxVal = std::numeric_limits<T>::max() )
339{
340 auto a = t[0] - t[2];
341 auto b = t[1] - t[2];
342 auto nom = dot( a, b );
343 auto den = cross( a, b ).length();
344 if ( fabs( nom ) >= absMaxVal * den )
345 return absMaxVal * sgn( nom );
346 return nom / den;
347}
348
352template <typename T>
353[[nodiscard]] inline T cotan( T a, T b, T c )
354{
355 const T den = ( a + b + c ) * ( b + c - a );
356 if ( den <= 0 )
357 return -std::numeric_limits<T>::infinity();
358 const T num = ( a + c - b ) * ( a + b - c );
359 if ( num <= 0 )
360 return std::numeric_limits<T>::infinity();
361 const auto tanSq = num / den;
362 return ( 1 - tanSq ) / ( 2 * std::sqrt( tanSq ) );
363}
364
367template <typename T>
368[[nodiscard]] std::optional<Vector3<T>> gradientInTri( const Vector3<T> & b, const Vector3<T> & c, T vb, T vc )
369{
370 const auto bb = dot( b, b );
371 const auto bc = dot( b, c );
372 const auto cc = dot( c, c );
373 const auto det = bb * cc - bc * bc;
374 if ( det <= 0 )
375 return {};
376 const auto kb = ( 1 / det ) * ( cc * vb - bc * vc );
377 const auto kc = ( 1 / det ) * (-bc * vb + bb * vc );
378 return kb * b + kc * c;
379}
380
383template <typename T>
384[[nodiscard]] std::optional<Vector3<T>> gradientInTri( const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c, T va, T vb, T vc )
385{
386 return gradientInTri( b - a, c - a, vb - va, vc - va );
387}
388
389// consider triangle 0BC, where gradient of linear scalar field is given;
390// computes the intersection of the ray (org=0, dir=-grad) with the open line segment BC
391template <typename T>
392[[nodiscard]] std::optional<T> findTriExitPos( const Vector3<T> & b, const Vector3<T> & c, const Vector3<T> & grad )
393{
394 const auto gradSq = grad.lengthSq();
395 if ( gradSq <= 0 )
396 return {};
397 const auto d = c - b;
398 // gort is a vector in the triangle plane orthogonal to grad
399 const auto gort = d - ( dot( d, grad ) / gradSq ) * grad;
400 const auto god = dot( gort, d );
401 if ( god <= 0 )
402 return {};
403 const auto gob = -dot( gort, b );
404 if ( gob <= 0 || gob >= god )
405 return {};
406 const auto a = gob / god;
407 assert( a < std::numeric_limits<T>::max() );
408 const auto ip = a * c + ( 1 - a ) * b;
409 if ( dot( grad, ip ) >= 0 )
410 return {}; // (b,c) is intersected in the direction +grad
411 return a;
412}
413
420template <typename T>
421[[nodiscard]] std::optional<Vector3<T>> tangentPlaneNormalToSpheres( const Vector3<T> & b, const Vector3<T> & c, T rb, T rc )
422{
423 auto grad = gradientInTri( b, c, rb, rc );
424 if ( !grad.has_value() )
425 return {}; // degenerate triangle
426 auto gradSq = grad->lengthSq();
427 if ( gradSq >= 1 )
428 return {}; // the larger of the spheres contains the origin inside - no touch plane exists
429 return sqrt( 1 - gradSq ) * normal( b, c ) - *grad; // unit normal
430}
431
438template <typename T>
439[[nodiscard]] std::optional<Plane3<T>> tangentPlaneToSpheres( const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c, T ra, T rb, T rc )
440{
441 if ( auto n = tangentPlaneNormalToSpheres( b - a, c - a, rb - ra, rc - ra ) )
442 return Plane3<T>( *n, dot( *n, a ) + ra );
443 return {}; // the larger of the spheres contains another sphere inside - no touch plane exists
444}
445
446} // namespace MR
length
Definition MRObjectDimensionsEnum.h:14
T minTriangleAngleSin(const Vector3< T > &a, const Vector3< T > &b, const Vector3< T > &c)
Definition MRTriMath.h:92
Definition MRCameraOrientationPlugin.h:8
float area(const MeshTopology &topology, const VertCoords &points, FaceId f)
returns the area of given face
Definition MRMeshMath.h:165
constexpr T sqr(T x) noexcept
squared value
Definition MRMesh/MRMeshFwd.h:753
MRMESH_API float circumcircleDiameter(const MeshTopology &topology, const VertCoords &points, FaceId f)
returns circumcircle diameter of given mesh triangle
MRMESH_API float dihedralAngle(const MeshTopology &topology, const VertCoords &points, UndirectedEdgeId e)
Vector3< T > circumcircleCenter(const Vector3< T > &a, const Vector3< T > &b)
Computes the center of the the triangle's 0AB circumcircle.
Definition MRTriMath.h:45
MRMESH_API float dihedralAngleCos(const MeshTopology &topology, const VertCoords &points, UndirectedEdgeId e)
std::optional< T > findTriExitPos(const Vector3< T > &b, const Vector3< T > &c, const Vector3< T > &grad)
Definition MRTriMath.h:392
std::optional< T > quadrangleOtherDiagonal(T a, T b, T c, T a1, T b1)
Definition MRTriMath.h:304
Triangle3< T > makeDegenerate(const Triangle3< T > &t)
make degenerate triangle (all 3 points on a line) that maximally resembles the input one and has the ...
Definition MRTriMath.h:204
T tanSqOfHalfAngle(T a, T b, T c)
Definition MRTriMath.h:324
Triangle3< T > triangleWithNormal(const Triangle3< T > &t, const Vector3< T > &n)
Definition MRTriMath.h:229
Vector3f dirDblArea(const MeshTopology &topology, const VertCoords &points, FaceId f)
computes directed double area for a triangular face from its vertices
Definition MRMeshMath.h:150
MRMESH_API float triangleAspectRatio(const MeshTopology &topology, const VertCoords &points, FaceId f)
returns aspect ratio of given mesh triangle equal to the ratio of the circum-radius to twice its in-r...
constexpr int sgn(T x) noexcept
sign of given value in { -1, 0, 1 }
Definition MRMesh/MRMeshFwd.h:757
std::optional< Vector3< T > > gradientInTri(const Vector3< T > &b, const Vector3< T > &c, T vb, T vc)
Definition MRTriMath.h:368
std::optional< Plane3< T > > tangentPlaneToSpheres(const Vector3< T > &a, const Vector3< T > &b, const Vector3< T > &c, T ra, T rb, T rc)
Definition MRTriMath.h:439
MRMESH_API float circumcircleDiameterSq(const MeshTopology &topology, const VertCoords &points, FaceId f)
returns squared circumcircle diameter of given mesh triangle
bool circumballCenters(const Vector3< T > &a, const Vector3< T > &b, const Vector3< T > &c, T radius, Vector3< T > &centerPos, Vector3< T > &centerNeg)
Definition MRTriMath.h:71
std::optional< Vector3< T > > tangentPlaneNormalToSpheres(const Vector3< T > &b, const Vector3< T > &c, T rb, T rc)
Definition MRTriMath.h:421
float dblArea(const MeshTopology &topology, const VertCoords &points, FaceId f)
returns twice the area of given face
Definition MRMeshMath.h:159
T minTriangleAngle(const Vector3< T > &a, const Vector3< T > &b, const Vector3< T > &c)
Definition MRTriMath.h:104
std::optional< Vector2< T > > posFromTriEdgeLengths(T a, T b, T c)
Definition MRTriMath.h:282
float cotan(const MeshTopology &topology, const VertCoords &points, UndirectedEdgeId ue)
Definition MRMeshMath.h:309
T dblAreaSq(const Vector3< T > &p, const Vector3< T > &q, const Vector3< T > &r)
computes the square of double area of given triangle
Definition MRTriMath.h:162
MRMESH_API float dihedralAngleSin(const MeshTopology &topology, const VertCoords &points, UndirectedEdgeId e)
Definition MRPlane3.h:17
Definition MRVector2.h:27
Definition MRMesh/MRVector3.h:28
Vector3 normalized() const MR_REQUIRES_IF_SUPPORTED(std
Definition MRMesh/MRVector3.h:68
T lengthSq() const
Definition MRMesh/MRVector3.h:59