Skip to content

Restrict steepest descent path to a region #1616

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Fedr merged 2 commits into from
Sep 15, 2023
Merged

Restrict steepest descent path to a region #1616

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
beda4f6
Restrict steepest descent path to a region
Fedr Sep 15, 2023
35cf6c2
more debug checks
Fedr Sep 15, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions source/MRMesh/MRMeshTopology.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -430,6 +430,14 @@ void MeshTopology::setLeft( EdgeId a, FaceId f )
}
}

bool MeshTopology::isInnerOrBdVertex( VertId v, const FaceBitSet * region ) const
{
for ( auto e : orgRing( *this, v ) )
if ( contains( region, left( e ) ) )
return true;
return false;
}

EdgeId MeshTopology::nextLeftBd( EdgeId e, const FaceBitSet * region ) const
{
assert( isLeftBdEdge( e, region ) );
Expand Down
2 changes: 2 additions & 0 deletions source/MRMesh/MRMeshTopology.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -212,6 +212,8 @@ class MeshTopology
[[nodiscard]] EdgeId bdEdgeWithOrigin( VertId v, const FaceBitSet * region = nullptr ) const { return bdEdgeSameOrigin( edgeWithOrg( v ), region ); }
/// returns true if given vertex is on (region) boundary
[[nodiscard]] bool isBdVertex( VertId v, const FaceBitSet * region = nullptr ) const { return isBdVertexInOrg( edgeWithOrg( v ), region ); }
/// returns true if one of incident faces of given vertex pertain to given region (or any such face exists if region is nullptr)
[[nodiscard]] MRMESH_API bool isInnerOrBdVertex( VertId v, const FaceBitSet * region = nullptr ) const;
/// returns true if left face of given edge belongs to given region (if provided) and right face either does not exist or does not belong
[[nodiscard]] bool isLeftBdEdge( EdgeId e, const FaceBitSet * region = nullptr ) const { return region ? ( isLeftInRegion( e, region ) && !isLeftInRegion( e.sym(), region ) ) : !right( e ); }
/// return true if given edge is inner or boundary for given region (or for whole mesh if region is nullptr), returns false for lone edges
Expand Down
92 changes: 49 additions & 43 deletions source/MRMesh/MRSurfacePath.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,16 +103,20 @@ static std::optional<float> computeExitPos( const Vector3f & b, const Vector3f &
return a;
}

MeshEdgePoint findSteepestDescentPoint( const Mesh & mesh, const VertScalars & field, VertId v )
MeshEdgePoint findSteepestDescentPoint( const MeshPart & mp, const VertScalars & field, VertId v )
{
assert( mp.mesh.topology.isInnerOrBdVertex( v, mp.region ) );

MeshEdgePoint res;
float maxGradSq = 0;
const auto vv = field[v];
const auto pv = mesh.points[v];
for ( EdgeId e : orgRing( mesh.topology, v ) )
const auto pv = mp.mesh.points[v];
for ( EdgeId e : orgRing( mp.mesh.topology, v ) )
{
const auto d = mesh.topology.dest( e );
const auto pd = mesh.points[d] - pv;
if ( mp.region && !mp.mesh.topology.isInnerOrBdEdge( e, mp.region ) )
continue;
const auto d = mp.mesh.topology.dest( e );
const auto pd = mp.mesh.points[d] - pv;
const auto pdSq = pd.lengthSq();
if ( field[d] == FLT_MAX )
continue;
Expand All @@ -131,11 +135,11 @@ MeshEdgePoint findSteepestDescentPoint( const Mesh & mesh, const VertScalars & f
}
}
}
if ( mesh.topology.left( e ) )
if ( auto f = mp.mesh.topology.left( e ); contains( mp.region, f ) )
{
const auto eBd = mesh.topology.prev( e.sym() );
const auto x = mesh.topology.dest( eBd );
const auto px = mesh.points[x] - pv;
const auto eBd = mp.mesh.topology.prev( e.sym() );
const auto x = mp.mesh.topology.dest( eBd );
const auto px = mp.mesh.points[x] - pv;
if ( auto fx = field[x]; fx < FLT_MAX )
{
const auto vx = fx - vv;
Expand All @@ -155,21 +159,22 @@ MeshEdgePoint findSteepestDescentPoint( const Mesh & mesh, const VertScalars & f
return res;
}

MeshEdgePoint findSteepestDescentPoint( const Mesh & mesh, const VertScalars & field, const MeshEdgePoint & ep )
MeshEdgePoint findSteepestDescentPoint( const MeshPart & mp, const VertScalars & field, const MeshEdgePoint & ep )
{
if ( auto v = ep.inVertex( mesh.topology ) )
return findSteepestDescentPoint( mesh, field, v );
if ( auto v = ep.inVertex( mp.mesh.topology ) )
return findSteepestDescentPoint( mp, field, v );
assert( mp.mesh.topology.isInnerOrBdEdge( ep.e, mp.region ) );

// point is not in vertex
const auto p = mesh.edgePoint( ep );
const auto p = mp.mesh.edgePoint( ep );

const auto o = mesh.topology.org( ep.e );
const auto d = mesh.topology.dest( ep.e );
const auto o = mp.mesh.topology.org( ep.e );
const auto d = mp.mesh.topology.dest( ep.e );
const auto fo = field[o];
const auto fd = field[d];
const auto v = ( 1 - ep.a ) * fo + ep.a * fd;
const auto po = mesh.points[o];
const auto pd = mesh.points[d];
const auto po = mp.mesh.points[o];
const auto pd = mp.mesh.points[d];

MeshEdgePoint result;
float maxGradSq = -FLT_MAX;
Expand All @@ -181,12 +186,12 @@ MeshEdgePoint findSteepestDescentPoint( const Mesh & mesh, const VertScalars & f
result = fo < fd ? MeshEdgePoint{ ep.e, 0 } : MeshEdgePoint{ ep.e.sym(), 0 };
}

if ( mesh.topology.left( ep.e ) )
if ( auto f = mp.mesh.topology.left( ep.e ); contains( mp.region, f ) )
{
const auto el = mesh.topology.next( ep.e );
const auto l = mesh.topology.dest( el );
const auto el = mp.mesh.topology.next( ep.e );
const auto l = mp.mesh.topology.dest( el );
const auto fl = field[l];
const auto pl = mesh.points[l];
const auto pl = mp.mesh.points[l];
const auto triGrad = computeGradient( pd - po, pl - po, fd - fo, fl - fo );
const auto triGradSq = triGrad.lengthSq();
bool moveL = true;
Expand All @@ -202,7 +207,7 @@ MeshEdgePoint findSteepestDescentPoint( const Mesh & mesh, const VertScalars & f
if ( a <= 1 )
{
moveL = false;
result = MeshEdgePoint{ mesh.topology.prev( ep.e.sym() ), a };
result = MeshEdgePoint{ mp.mesh.topology.prev( ep.e.sym() ), a };
maxGradSq = triGradSq;
}
else
Expand Down Expand Up @@ -233,12 +238,12 @@ MeshEdgePoint findSteepestDescentPoint( const Mesh & mesh, const VertScalars & f
}
}

if ( mesh.topology.right( ep.e ) )
if ( auto f = mp.mesh.topology.right( ep.e ); contains( mp.region, f ) )
{
const auto er = mesh.topology.prev( ep.e );
const auto r = mesh.topology.dest( er );
const auto er = mp.mesh.topology.prev( ep.e );
const auto r = mp.mesh.topology.dest( er );
const auto fr = field[r];
const auto pr = mesh.points[r];
const auto pr = mp.mesh.points[r];
const auto triGrad = computeGradient( pr - po, pd - po, fr - fo, fd - fo );
const auto triGradSq = triGrad.lengthSq();
bool moveR = true;
Expand All @@ -254,7 +259,7 @@ MeshEdgePoint findSteepestDescentPoint( const Mesh & mesh, const VertScalars & f
if ( a >= 0 )
{
moveR = false;
result = MeshEdgePoint{ mesh.topology.next( ep.e.sym() ).sym(), a };
result = MeshEdgePoint{ mp.mesh.topology.next( ep.e.sym() ).sym(), a };
maxGradSq = triGradSq;
}
else
Expand Down Expand Up @@ -296,28 +301,29 @@ MeshEdgePoint findSteepestDescentPoint( const Mesh & mesh, const VertScalars & f
return result;
}

MeshEdgePoint findSteepestDescentPoint( const Mesh & mesh, const VertScalars & field, const MeshTriPoint & tp )
MeshEdgePoint findSteepestDescentPoint( const MeshPart & mp, const VertScalars & field, const MeshTriPoint & tp )
{
if ( auto ep = tp.onEdge( mesh.topology ) )
return findSteepestDescentPoint( mesh, field, ep );
if ( auto ep = tp.onEdge( mp.mesh.topology ) )
return findSteepestDescentPoint( mp, field, ep );
assert( contains( mp.region, mp.mesh.topology.left( tp.e ) ) );

// point is not on edge
MeshEdgePoint res;
const auto p = mesh.triPoint( tp );
const auto p = mp.mesh.triPoint( tp );

VertId v[3];
mesh.topology.getLeftTriVerts( tp.e, v );
mp.mesh.topology.getLeftTriVerts( tp.e, v );

Vector3f pv[3];
float vv[3];
EdgeId e[3];
auto ei = tp.e;
for ( int i = 0; i < 3; ++i )
{
pv[i] = mesh.points[v[i]];
pv[i] = mp.mesh.points[v[i]];
vv[i] = field[v[i]];
e[i] = ei;
ei = mesh.topology.prev( ei.sym() );
ei = mp.mesh.topology.prev( ei.sym() );
}
if ( vv[0] == vv[1] && vv[1] == vv[2] )
return res; // the triangle is completely "flat"
Expand Down Expand Up @@ -429,51 +435,51 @@ Expected<SurfacePath, PathError> computeGeodesicPathApprox( const Mesh & mesh,
return res;
}

SurfacePath computeSteepestDescentPath( const Mesh & mesh, const VertScalars & field,
SurfacePath computeSteepestDescentPath( const MeshPart & mp, const VertScalars & field,
const MeshTriPoint & start, const ComputeSteepestDescentPathSettings & settings )
{
SurfacePath res;
computeSteepestDescentPath( mesh, field, start, &res, settings );
computeSteepestDescentPath( mp, field, start, &res, settings );
return res;
}

void computeSteepestDescentPath( const Mesh & mesh, const VertScalars & field,
void computeSteepestDescentPath( const MeshPart & mp, const VertScalars & field,
const MeshTriPoint & start, SurfacePath * outPath, const ComputeSteepestDescentPathSettings & settings )
{
assert( start );
assert( settings.outVertexReached || settings.outBdReached || outPath );
size_t iniPathSize = outPath ? outPath->size() : 0;
size_t edgesPassed = 0;
auto curr = findSteepestDescentPoint( mesh, field, start );
auto curr = findSteepestDescentPoint( mp, field, start );
while ( curr )
{
if ( settings.outVertexReached )
{
if ( auto v = curr.inVertex( mesh.topology ) )
if ( auto v = curr.inVertex( mp.mesh.topology ) )
{
*settings.outVertexReached = v;
return;
}
}
if ( settings.outBdReached && curr.isBd( mesh.topology ) )
if ( settings.outBdReached && curr.isBd( mp.mesh.topology ) )
{
*settings.outBdReached = curr;
return;
}
++edgesPassed;
if ( outPath )
outPath->push_back( curr );
if ( settings.end && fromSameTriangle( mesh.topology, MeshTriPoint( settings.end ), MeshTriPoint( curr ) ) )
if ( settings.end && fromSameTriangle( mp.mesh.topology, MeshTriPoint( settings.end ), MeshTriPoint( curr ) ) )
break; // reached triangle with end point
if ( edgesPassed > mesh.topology.numValidFaces() )
if ( edgesPassed > mp.mesh.topology.numValidFaces() )
{
// normal path cannot visit any triangle more than once
assert( false );
if ( outPath )
outPath->resize( iniPathSize );
return;
}
curr = findSteepestDescentPoint( mesh, field, curr );
curr = findSteepestDescentPoint( mp, field, curr );
}
}

Expand Down
10 changes: 5 additions & 5 deletions source/MRMesh/MRSurfacePath.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,24 +77,24 @@ struct ComputeSteepestDescentPathSettings
/// computes the path (edge points crossed by the path) staring in given point
/// and moving in each triangle in minus gradient direction of given field;
/// the path stops when it reaches a local minimum in the field or one of stop conditions in settings
[[nodiscard]] MRMESH_API SurfacePath computeSteepestDescentPath( const Mesh & mesh, const VertScalars & field,
[[nodiscard]] MRMESH_API SurfacePath computeSteepestDescentPath( const MeshPart & mp, const VertScalars & field,
const MeshTriPoint & start, const ComputeSteepestDescentPathSettings & settings = {} );

/// computes the path (edge points crossed by the path) staring in given point
/// and moving in each triangle in minus gradient direction of given field,
/// and outputs the path in \param outPath if requested;
/// the path stops when it reaches a local minimum in the field or one of stop conditions in settings
MRMESH_API void computeSteepestDescentPath( const Mesh & mesh, const VertScalars & field,
MRMESH_API void computeSteepestDescentPath( const MeshPart & mp, const VertScalars & field,
const MeshTriPoint & start, SurfacePath * outPath, const ComputeSteepestDescentPathSettings & settings = {} );

/// finds the point along minus maximal gradient on the boundary of first ring boundary around given vertex
[[nodiscard]] MRMESH_API MeshEdgePoint findSteepestDescentPoint( const Mesh & mesh, const VertScalars & field, VertId v );
[[nodiscard]] MRMESH_API MeshEdgePoint findSteepestDescentPoint( const MeshPart & mp, const VertScalars & field, VertId v );

/// finds the point along minus maximal gradient on the boundary of triangles around given point (the boundary of left and right edge triangles' union in case (ep) is inner edge point)
[[nodiscard]] MRMESH_API MeshEdgePoint findSteepestDescentPoint( const Mesh & mesh, const VertScalars & field, const MeshEdgePoint & ep );
[[nodiscard]] MRMESH_API MeshEdgePoint findSteepestDescentPoint( const MeshPart & mp, const VertScalars & field, const MeshEdgePoint & ep );

/// finds the point along minus maximal gradient on the boundary of triangles around given point (the boundary of the triangle itself in case (tp) is inner triangle point)
[[nodiscard]] MRMESH_API MeshEdgePoint findSteepestDescentPoint( const Mesh & mesh, const VertScalars & field, const MeshTriPoint & tp );
[[nodiscard]] MRMESH_API MeshEdgePoint findSteepestDescentPoint( const MeshPart & mp, const VertScalars & field, const MeshTriPoint & tp );

enum class ExtremeEdgeType
{
Expand Down