Skip to content

Commit e085fff

Browse files
FedrFedr
authored
gradientInTri (#4283)
1 parent 28fc1bf commit e085fff

9 files changed

+186
-157
lines changed

‎source/MRMesh/MRMesh.vcxproj

-1
Original file line numberDiff line numberDiff line change
@@ -615,7 +615,6 @@
615615
<ClCompile Include="MRSurfaceDistanceBuilder.cpp" />
616616
<ClCompile Include="MRSurroundingContour.cpp" />
617617
<ClCompile Include="MRTerrainTriangulation.cpp" />
618-
<ClCompile Include="MRTriMath.cpp" />
619618
<ClCompile Include="MRTunnelDetector.cpp" />
620619
<ClCompile Include="MRTupleBindings.cpp" />
621620
<ClCompile Include="MRUniformSampling.cpp" />

‎source/MRMesh/MRMesh.vcxproj.filters

-3
Original file line numberDiff line numberDiff line change
@@ -1724,9 +1724,6 @@
17241724
<ClCompile Include="MROrder.cpp">
17251725
<Filter>Source Files\Mesh</Filter>
17261726
</ClCompile>
1727-
<ClCompile Include="MRTriMath.cpp">
1728-
<Filter>Source Files\Math</Filter>
1729-
</ClCompile>
17301727
<ClCompile Include="MRExtractIsolines.cpp">
17311728
<Filter>Source Files\SurfacePath</Filter>
17321729
</ClCompile>

‎source/MRMesh/MRSurfacePath.cpp

+3-120
Original file line numberDiff line numberDiff line change
@@ -2,42 +2,23 @@
22
#include "MRBitSetParallelFor.h"
33
#include "MREdgePaths.h"
44
#include "MRExtractIsolines.h"
5-
#include "MRGTest.h"
65
#include "MRLaplacian.h"
76
#include "MRMesh.h"
8-
#include "MRMeshBuilder.h"
97
#include "MRMeshComponents.h"
108
#include "MRMeshPart.h"
119
#include "MRGeodesicPath.h"
1210
#include "MRRegionBoundary.h"
1311
#include "MRRingIterator.h"
1412
#include "MRSurfaceDistance.h"
1513
#include "MRTimer.h"
14+
#include "MRTriMath.h"
1615

1716
namespace MR
1817
{
1918

20-
// consider triangle 0bc, where a linear scalar field is defined in all vertices: v(0) = 0, v(b) = vb, v(c) = vc;
21-
// computes field gradient in the triangle
22-
static Vector3d computeGradient( const Vector3d & b, const Vector3d & c, double vb, double vc )
23-
{
24-
const auto bb = dot( b, b );
25-
const auto bc = dot( b, c );
26-
const auto cc = dot( c, c );
27-
const auto det = bb * cc - bc * bc;
28-
if ( det <= 0 )
29-
{
30-
// degenerate triangle
31-
return {};
32-
}
33-
const auto kb = ( 1 / det ) * ( cc * vb - bc * vc );
34-
const auto kc = ( 1 / det ) * (-bc * vb + bb * vc );
35-
return kb * b + kc * c;
36-
}
37-
3819
static Vector3f computeGradient( const Vector3f & b, const Vector3f & c, float vb, float vc )
3920
{
40-
return Vector3f{ computeGradient( Vector3d( b ), Vector3d( c ), double( vb ), double( vc ) ) };
21+
return Vector3f{ gradientInTri( Vector3d( b ), Vector3d( c ), double( vb ), double( vc ) ) };
4122
}
4223

4324
/// given triangle with scalar field increasing in the direction \param dir;
@@ -78,30 +59,6 @@ static bool computeEnter01Cross( const Triangle3f & t, const Vector3f & unitDir,
7859
return computeLineLineCross( t[0] - p, t[1] - p, unitDir, a );
7960
}
8061

81-
// consider triangle 0bc, where gradient is given;
82-
// computes the intersection of the ray (org=0, dir=-grad) with the open segment (b,c)
83-
static std::optional<float> computeExitPos( const Vector3f & b, const Vector3f & c, const Vector3f & grad )
84-
{
85-
const auto gradSq = grad.lengthSq();
86-
if ( gradSq <= 0 )
87-
return {};
88-
const auto d = c - b;
89-
// gort is a vector in the triangle plane orthogonal to grad
90-
const auto gort = d - ( dot( d, grad ) / gradSq ) * grad;
91-
const auto god = dot( gort, d );
92-
if ( god <= 0 )
93-
return {};
94-
const auto gob = -dot( gort, b );
95-
if ( gob <= 0 || gob >= god )
96-
return {};
97-
const auto a = gob / god;
98-
assert( a < FLT_MAX );
99-
const auto ip = a * c + ( 1 - a ) * b;
100-
if ( dot( grad, ip ) >= 0 )
101-
return {}; // (b,c) is intersected in the direction +grad
102-
return a;
103-
}
104-
10562
MeshEdgePoint findSteepestDescentPoint( const MeshPart & mp, const VertScalars & field, VertId v )
10663
{
10764
assert( mp.mesh.topology.isInnerOrBdVertex( v, mp.region ) );
@@ -146,7 +103,7 @@ MeshEdgePoint findSteepestDescentPoint( const MeshPart & mp, const VertScalars &
146103
const auto triGradSq = triGrad.lengthSq();
147104
if ( triGradSq > maxGradSq )
148105
{
149-
if ( auto a = computeExitPos( pd, px, triGrad ) )
106+
if ( auto a = findTriExitPos( pd, px, triGrad ) )
150107
{
151108
maxGradSq = triGradSq;
152109
res = MeshEdgePoint{ eBd, *a };
@@ -694,78 +651,4 @@ Contours3f surfacePathsToContours3f( const Mesh & mesh, const SurfacePaths & lin
694651
return res;
695652
}
696653

697-
TEST(MRMesh, SurfacePath)
698-
{
699-
Vector3f g;
700-
701-
g = computeGradient( Vector3f{ 1, 0, 0 }, Vector3f{ 0.5f, 1, 0 }, 0, 1 );
702-
EXPECT_NEAR( ( g - Vector3f{ 0, 1, 0 } ).length(), 0, 1e-5f );
703-
704-
g = computeGradient( Vector3f{ 1, 0, 0 }, Vector3f{ 0.1f, 1, 0 }, 0, 1 );
705-
EXPECT_NEAR( ( g - Vector3f{ 0, 1, 0 } ).length(), 0, 1e-5f );
706-
707-
g = computeGradient( Vector3f{ 1, 0, 0 }, Vector3f{ 0.9f, 1, 0 }, 0, 1 );
708-
EXPECT_NEAR( ( g - Vector3f{ 0, 1, 0 } ).length(), 0, 1e-5f );
709-
710-
std::optional<float> e;
711-
g = computeGradient( Vector3f{ 1, 0, 0 }, Vector3f{ 0, 1, 0 }, 1, 1 );
712-
EXPECT_NEAR( ( g - Vector3f{ 1, 1, 0 } ).length(), 0, 1e-5f );
713-
e = computeExitPos ( Vector3f{ 1, 0, 0 }, Vector3f{ 0, 1, 0 }, g );
714-
EXPECT_FALSE( e.has_value() );
715-
e = computeExitPos ( Vector3f{ 1, 0, 0 }, Vector3f{ 0, 1, 0 }, -g );
716-
EXPECT_NEAR( *e, 0.5f, 1e-5f );
717-
718-
g = computeGradient( Vector3f{ 1, -1, 0 }, Vector3f{ 1, 1, 0 }, -1, -1 );
719-
EXPECT_NEAR( ( g - Vector3f{ -1, 0, 0 } ).length(), 0, 1e-5f );
720-
e = computeExitPos ( Vector3f{ 1, -1, 0 }, Vector3f{ 1, 1, 0 }, g );
721-
EXPECT_NEAR( *e, 0.5f, 1e-5f );
722-
e = computeExitPos ( Vector3f{ 1, -1, 0 }, Vector3f{ 1, 1, 0 }, -g );
723-
EXPECT_FALSE( e.has_value() );
724-
725-
g = computeGradient( Vector3f{ 1, -0.1f, 0 }, Vector3f{ 1, 0.9f, 0 }, -1, -1 );
726-
EXPECT_NEAR( ( g - Vector3f{ -1, 0, 0 } ).length(), 0, 1e-5f );
727-
e = computeExitPos ( Vector3f{ 1, -0.1f, 0 }, Vector3f{ 1, 0.9f, 0 }, g );
728-
EXPECT_NEAR( *e, 0.1f, 1e-5f );
729-
730-
g = computeGradient( Vector3f{ 1, -0.9f, 0 }, Vector3f{ 1, 0.1f, 0 }, -1, -1 );
731-
EXPECT_NEAR( ( g - Vector3f{ -1, 0, 0 } ).length(), 0, 1e-5f );
732-
e = computeExitPos ( Vector3f{ 1, -0.9f, 0 }, Vector3f{ 1, 0.1f, 0 }, g );
733-
EXPECT_NEAR( *e, 0.9f, 1e-5f );
734-
735-
g = computeGradient( Vector3f{ 1, 0.1f, 0 }, Vector3f{ 1, 0.9f, 0 }, -1, -1 );
736-
EXPECT_NEAR( ( g - Vector3f{ -1, 0, 0 } ).length(), 0, 1e-5f );
737-
e = computeExitPos ( Vector3f{ 1, 0.1f, 0 }, Vector3f{ 1, 0.9f, 0 }, g );
738-
EXPECT_FALSE( e.has_value() );
739-
e = computeExitPos ( Vector3f{ 1, 0.1f, 0 }, Vector3f{ 1, 0.9f, 0 }, -g );
740-
EXPECT_FALSE( e.has_value() );
741-
}
742-
743-
TEST( MRMesh, SurfacePathTargets )
744-
{
745-
Triangulation t{
746-
{ 0_v, 1_v, 2_v }
747-
};
748-
Mesh mesh;
749-
mesh.topology = MeshBuilder::fromTriangles( t );
750-
751-
mesh.points.emplace_back( 0.f, 0.f, 0.f ); // 0_v
752-
mesh.points.emplace_back( 1.f, 0.f, 0.f ); // 1_v
753-
mesh.points.emplace_back( 0.f, 1.f, 0.f ); // 2_v
754-
755-
VertBitSet starts(3);
756-
starts.set( 1_v );
757-
starts.set( 2_v );
758-
759-
VertBitSet ends(3);
760-
ends.set( 0_v );
761-
762-
const auto map = computeClosestSurfacePathTargets( mesh, starts, ends );
763-
EXPECT_EQ( map.size(), starts.count() );
764-
for ( const auto & [start, end] : map )
765-
{
766-
EXPECT_TRUE( starts.test( start ) );
767-
EXPECT_TRUE( ends.test( end ) );
768-
}
769-
}
770-
771654
} //namespace MR

‎source/MRMesh/MRTriMath.cpp

-33
This file was deleted.

‎source/MRMesh/MRTriMath.h

+53
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33

44
#include "MRVector3.h"
55
#include <algorithm>
6+
#include <cassert>
67
#include <cmath>
78
#include <limits>
89
#include <optional>
@@ -293,4 +294,56 @@ template <typename T>
293294
return ( *p - *p1 ).length();
294295
}
295296

297+
/// Consider triangle 0BC, where a linear scalar field is defined in all 3 vertices: v(0) = 0, v(b) = vb, v(c) = vc;
298+
/// computes and returns field gradient in the triangle
299+
template <typename T>
300+
[[nodiscard]] Vector3<T> gradientInTri( const Vector3<T> & b, const Vector3<T> & c, T vb, T vc )
301+
{
302+
const auto bb = dot( b, b );
303+
const auto bc = dot( b, c );
304+
const auto cc = dot( c, c );
305+
const auto det = bb * cc - bc * bc;
306+
if ( det <= 0 )
307+
{
308+
// degenerate triangle
309+
return {};
310+
}
311+
const auto kb = ( 1 / det ) * ( cc * vb - bc * vc );
312+
const auto kc = ( 1 / det ) * (-bc * vb + bb * vc );
313+
return kb * b + kc * c;
314+
}
315+
316+
/// Consider triangle ABC, where a linear scalar field is defined in all 3 vertices: v(a) = va, v(b) = vb, v(c) = vc;
317+
/// computes and returns field gradient in the triangle
318+
template <typename T>
319+
[[nodiscard]] Vector3<T> gradientInTri( const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c, T va, T vb, T vc )
320+
{
321+
return gradientInTri( b - a, c - a, vb - va, vc - va );
322+
}
323+
324+
// consider triangle 0BC, where gradient of linear scalar field is given;
325+
// computes the intersection of the ray (org=0, dir=-grad) with the open line segment BC
326+
template <typename T>
327+
[[nodiscard]] std::optional<T> findTriExitPos( const Vector3<T> & b, const Vector3<T> & c, const Vector3<T> & grad )
328+
{
329+
const auto gradSq = grad.lengthSq();
330+
if ( gradSq <= 0 )
331+
return {};
332+
const auto d = c - b;
333+
// gort is a vector in the triangle plane orthogonal to grad
334+
const auto gort = d - ( dot( d, grad ) / gradSq ) * grad;
335+
const auto god = dot( gort, d );
336+
if ( god <= 0 )
337+
return {};
338+
const auto gob = -dot( gort, b );
339+
if ( gob <= 0 || gob >= god )
340+
return {};
341+
const auto a = gob / god;
342+
assert( a < std::numeric_limits<T>::max() );
343+
const auto ip = a * c + ( 1 - a ) * b;
344+
if ( dot( grad, ip ) >= 0 )
345+
return {}; // (b,c) is intersected in the direction +grad
346+
return a;
347+
}
348+
296349
} // namespace MR

‎source/MRTest/MRSurfacePathTests.cpp

+36
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
#include <MRMesh/MRSurfacePath.h>
2+
#include <MRMesh/MRGTest.h>
3+
#include <MRMesh/MRMeshBuilder.h>
4+
5+
namespace MR
6+
{
7+
8+
TEST( MRMesh, SurfacePathTargets )
9+
{
10+
Triangulation t{
11+
{ 0_v, 1_v, 2_v }
12+
};
13+
Mesh mesh;
14+
mesh.topology = MeshBuilder::fromTriangles( t );
15+
16+
mesh.points.emplace_back( 0.f, 0.f, 0.f ); // 0_v
17+
mesh.points.emplace_back( 1.f, 0.f, 0.f ); // 1_v
18+
mesh.points.emplace_back( 0.f, 1.f, 0.f ); // 2_v
19+
20+
VertBitSet starts(3);
21+
starts.set( 1_v );
22+
starts.set( 2_v );
23+
24+
VertBitSet ends(3);
25+
ends.set( 0_v );
26+
27+
const auto map = computeClosestSurfacePathTargets( mesh, starts, ends );
28+
EXPECT_EQ( map.size(), starts.count() );
29+
for ( const auto & [start, end] : map )
30+
{
31+
EXPECT_TRUE( starts.test( start ) );
32+
EXPECT_TRUE( ends.test( end ) );
33+
}
34+
}
35+
36+
} //namespace MR

‎source/MRTest/MRTest.vcxproj

+2
Original file line numberDiff line numberDiff line change
@@ -24,9 +24,11 @@
2424
<ClCompile Include="MRPdfTests.cpp" />
2525
<ClCompile Include="MRPolylineTrimWithPlane.cpp" />
2626
<ClCompile Include="MRSpdlog.cpp" />
27+
<ClCompile Include="MRSurfacePathTests.cpp" />
2728
<ClCompile Include="MRTestApp.cpp" />
2829
<ClCompile Include="MRTriangleIntersection.cpp" />
2930
<ClCompile Include="MRMeshVoxelsConverter.cpp" />
31+
<ClCompile Include="MRTriMathTests.cpp" />
3032
<ClCompile Include="MRVolumeToMeshByPartsTests.cpp" />
3133
<ClCompile Include="MRZlib.cpp" />
3234
</ItemGroup>

‎source/MRTest/MRTest.vcxproj.filters

+6
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,12 @@
6161
<ClCompile Include="MRMeshIntersectTests.cpp">
6262
<Filter>Source Files</Filter>
6363
</ClCompile>
64+
<ClCompile Include="MRTriMathTests.cpp">
65+
<Filter>Source Files</Filter>
66+
</ClCompile>
67+
<ClCompile Include="MRSurfacePathTests.cpp">
68+
<Filter>Source Files</Filter>
69+
</ClCompile>
6470
</ItemGroup>
6571
<ItemGroup>
6672
<None Include="..\.editorconfig" />

0 commit comments

Comments
 (0)