Introduce plane tracked geodesic path

source/MRMesh/MRContoursCut.cpp

@@ -896,9 +896,8 @@ OneMeshIntersection intersectionFromMeshTriPoint( const Mesh& mesh, const MeshTr
     return res;
 }
 
-
-Expected<OneMeshContour, PathError> convertMeshTriPointsToMeshContour( const Mesh& mesh, const std::vector<MeshTriPoint>& meshTriPointsOrg,
-    SearchPathSettings searchSettings, std::vector<int>* pivotIndices )
+Expected<OneMeshContour> convertMeshTriPointsToMeshContour( const Mesh& mesh, const std::vector<MeshTriPoint>& meshTriPointsOrg, 
+    MeshTriPointsConnector connectorFn /*= {}*/, std::vector<int>* pivotIndices /*= nullptr */ )
 {
     MR_TIMER;
     if ( meshTriPointsOrg.size() < 2 )
@@ -950,6 +949,17 @@ Expected<OneMeshContour, PathError> convertMeshTriPointsToMeshContour( const Mes
     if ( meshTriPoints.size() < 2 )
         return {};
 
+    if ( !connectorFn )
+    {
+        connectorFn = [&] ( const MeshTriPoint& start, const MeshTriPoint& end, int, int )->Expected<SurfacePath>
+        {
+            auto res = computeGeodesicPath( mesh, start, end );
+            if ( !res.has_value() )
+                return unexpected( toString( res.error() ) );
+            return *res;
+        };
+    }
+
     int sameMtpsNavigator = 0;
     int pivotNavigator = 0;
     sizeMTP = closed ? meshTriPoints.size() : meshTriPoints.size() - 1;
@@ -958,8 +968,27 @@ Expected<OneMeshContour, PathError> convertMeshTriPointsToMeshContour( const Mes
     std::vector<OneMeshContour> surfacePaths( sizeMTP );
     for ( int i = 0; i < sizeMTP; ++i )
     {
+        while ( sameMtpsNavigator < sameEdgeMTPs.size() && pivotNavigator == sameEdgeMTPs[sameMtpsNavigator] )
+        {
+            ++pivotNavigator;
+            ++sameMtpsNavigator;
+        }
+        int firstIndex = pivotNavigator;
+        ++pivotNavigator;
+        int secondIndex = pivotNavigator;
+        int secSameNav = sameMtpsNavigator;
+        while ( secSameNav < sameEdgeMTPs.size() && secondIndex == sameEdgeMTPs[secSameNav] )
+        {
+            ++secondIndex;
+            ++secSameNav;
+            if ( secondIndex >= meshTriPointsOrg.size() )
+            {
+                secondIndex = 0;
+                secSameNav = 0;
+            }
+        }
         // using DijkstraAStar here might be faster, in most case points are close to each other
-        auto sp = computeGeodesicPath( mesh, meshTriPoints[i], meshTriPoints[( i + 1 ) % meshTriPoints.size()], searchSettings.geodesicPathApprox, searchSettings.maxReduceIters );
+        auto sp = connectorFn( meshTriPoints[i], meshTriPoints[( i + 1 ) % meshTriPoints.size()], firstIndex, secondIndex );
         if ( !sp.has_value() )
             return unexpected( sp.error() );
         auto partContours = convertSurfacePathsToMeshContours( mesh, { std::move( sp.value() ) } );
@@ -999,6 +1028,8 @@ Expected<OneMeshContour, PathError> convertMeshTriPointsToMeshContour( const Mes
         }
     }
 
+    sameMtpsNavigator = 0;
+    pivotNavigator = 0;
     const float closeEdgeEps = std::numeric_limits<float>::epsilon() * box.diagonal();
     // add interjacent
     for ( int i = 0; i < meshTriPoints.size(); ++i )
@@ -1094,11 +1125,24 @@ Expected<OneMeshContour, PathError> convertMeshTriPointsToMeshContour( const Mes
     return res;
 }
 
+Expected<OneMeshContour> convertMeshTriPointsToMeshContour( const Mesh& mesh, const std::vector<MeshTriPoint>& meshTriPointsOrg,
+    SearchPathSettings searchSettings, std::vector<int>* pivotIndices )
+{
+    MeshTriPointsConnector conFn = [&] ( const MeshTriPoint& start, const MeshTriPoint& end, int, int )->Expected<SurfacePath>
+    {
+        auto res = computeGeodesicPath( mesh, start, end, searchSettings.geodesicPathApprox, searchSettings.maxReduceIters );
+        if ( !res.has_value() )
+            return unexpected( toString( res.error() ) );
+        return *res;
+    };
+    return convertMeshTriPointsToMeshContour( mesh, meshTriPointsOrg, conFn, pivotIndices );
+}
+
 Expected<OneMeshContours> convertMeshTriPointsIsoLineToMeshContour( const Mesh& mesh, const std::vector<MeshTriPoint>& meshTriPoints, float isoValue, SearchPathSettings searchSettings /*= {} */ )
 {
     auto initCutContours = convertMeshTriPointsToMeshContour( mesh, meshTriPoints, searchSettings );
     if ( !initCutContours.has_value() )
-        return unexpected( toString( initCutContours.error() ) );
+        return unexpected( initCutContours.error() );
 
     if ( isoValue == 0.0f )
         return OneMeshContours{ std::move( *initCutContours ) };
@@ -1157,7 +1201,7 @@ Expected<OneMeshContours> convertMeshTriPointsIsoLineToMeshContour( const Mesh&
     return res;
 }
 
-Expected<OneMeshContour, PathError> convertMeshTriPointsToClosedContour( const Mesh& mesh, const std::vector<MeshTriPoint>& meshTriPointsOrg,
+Expected<OneMeshContour> convertMeshTriPointsToClosedContour( const Mesh& mesh, const std::vector<MeshTriPoint>& meshTriPointsOrg,
     SearchPathSettings searchSettings, std::vector<int>* pivotIndices )
 {
     auto conts = meshTriPointsOrg;

source/MRMesh/MRContoursCut.h

@@ -58,6 +58,19 @@ MRMESH_API void subdivideLoneContours( Mesh& mesh, const OneMeshContours& contou
 MRMESH_API OneMeshContours getOneMeshIntersectionContours( const Mesh& meshA, const Mesh& meshB, const ContinuousContours& contours, bool getMeshAIntersections,
     const CoordinateConverters& converters, const AffineXf3f* rigidB2A = nullptr );
 
+
+using MeshTriPointsConnector = std::function<Expected<SurfacePath>( const MeshTriPoint& start, const MeshTriPoint& end, int startIndex, int endIndex )>;
+/** \ingroup BooleanGroup
+  * \brief Makes continuous contour by mesh tri points, if first and last meshTriPoint is the same, makes closed contour
+  *
+  * Finds paths between neighbor \p meshTriPoints with MeshTriPointsConnector function and build contour MR::cutMesh input
+  * \param connectorFn function to build path between neighbor meshTriPoints, if not present simple geodesic path function is used
+  * \param pivotIndices optional output indices of given meshTriPoints in result OneMeshContour
+  */
+[[nodiscard]]
+MRMESH_API Expected<OneMeshContour> convertMeshTriPointsToMeshContour( const Mesh& mesh, const std::vector<MeshTriPoint>& meshTriPoints,
+    MeshTriPointsConnector connectorFn, std::vector<int>* pivotIndices = nullptr );
+
 /// Geo path search settings
 struct SearchPathSettings
 {
@@ -73,7 +86,7 @@ struct SearchPathSettings
   * \param pivotIndices optional output indices of given meshTriPoints in result OneMeshContour
   */
 [[nodiscard]]
-MRMESH_API Expected<OneMeshContour, PathError> convertMeshTriPointsToMeshContour( const Mesh& mesh, const std::vector<MeshTriPoint>& meshTriPoints,
+MRMESH_API Expected<OneMeshContour> convertMeshTriPointsToMeshContour( const Mesh& mesh, const std::vector<MeshTriPoint>& meshTriPoints,
     SearchPathSettings searchSettings = {}, std::vector<int>* pivotIndices = nullptr );
 
 /** \ingroup BooleanGroup
@@ -95,7 +108,7 @@ MRMESH_API Expected<OneMeshContours> convertMeshTriPointsIsoLineToMeshContour( c
   * \note better use convertMeshTriPointsToMeshContour(...) instead, note that it requires same front and back MeshTriPoints for closed contour
   */
 [[nodiscard]]
-MRMESH_API Expected<OneMeshContour, PathError> convertMeshTriPointsToClosedContour( const Mesh& mesh, const std::vector<MeshTriPoint>& meshTriPoints,
+MRMESH_API Expected<OneMeshContour> convertMeshTriPointsToClosedContour( const Mesh& mesh, const std::vector<MeshTriPoint>& meshTriPoints,
     SearchPathSettings searchSettings = {}, std::vector<int>* pivotIndices = nullptr );
 
 /** \ingroup BooleanGroup

source/MRMesh/MREmbedTerrainStructure.cpp

@@ -17,6 +17,7 @@
 #include "MRMeshIntersect.h"
 #include "MRPolyline.h"
 #include "MRMapEdge.h"
+#include "MRSurfacePath.h"
 
 namespace MR
 {
@@ -297,9 +298,44 @@ Expected<TerrainEmbedder::MappedMeshContours> TerrainEmbedder::prepareTerrainCut
         for ( int i = 0; i < res.contours.size(); ++i )
         {
             bool lone = true;
-            auto contourRes = convertMeshTriPointsToClosedContour( result_, noBowtiesMtps[i], {}, &res.map[i] );
+            auto cont = noBowtiesMtps[i];
+            cont.push_back( cont.front() );
+            auto contourRes = convertMeshTriPointsToMeshContour( result_, cont,
+                [&] ( const MeshTriPoint& start, const MeshTriPoint& end, int startInd, int endInd )->Expected<SurfacePath>
+            {
+                auto initSMtpIndex = res.filtBowTiesMap[i][startInd];
+                auto initEMtpIndex = res.filtBowTiesMap[i][endInd];
+                auto baseSIndex = findOffsetContourIndex_( initSMtpIndex, offCont.idsShifts ) % bounds_[0].size();
+                auto baseEIndex = findOffsetContourIndex_( initEMtpIndex, offCont.idsShifts ) % bounds_[0].size();
+                auto planePoint = ( cutStructure_.orgPnt( bounds_[0][baseSIndex] ) + cutStructure_.orgPnt( bounds_[0][baseEIndex] ) ) * 0.5f;
+                auto sPoint = result_.triPoint( start );
+                auto ePoint = result_.triPoint( end );
+                auto midPoint = ( sPoint + ePoint ) * 0.5f;
+                auto ccwPath = trackSection( result_, start, end, planePoint, true );
+                auto cwPath = trackSection( result_, start, end, planePoint, false );
+                if ( ccwPath.has_value() && cwPath.has_value() )
+                {
+                    auto ccwL = surfacePathLength( result_, *ccwPath );
+                    auto cwL = surfacePathLength( result_, *cwPath );
+                    if ( ccwL < cwL )
+                        return ccwPath;
+                    else
+                        return cwPath;
+                }
+                else if ( ccwPath.has_value() )
+                    return ccwPath;
+                else if ( cwPath.has_value() )
+                    return cwPath;
+                else
+                {
+                    auto locRes = computeGeodesicPath( result_, start, end );
+                    if ( !locRes.has_value() )
+                        return unexpected( toString( locRes.error() ) );
+                    return *locRes;
+                }
+            }, &res.map[i] );
             if ( !contourRes.has_value() )
-                return unexpected( toString( contourRes.error() ) );
+                return unexpected( contourRes.error() );
             res.contours[i] = std::move( *contourRes );
             for ( int j = 0; j < res.contours[i].intersections.size(); ++j )
             {

source/MRMesh/MRExtractIsolines.cpp

@@ -12,6 +12,7 @@
 #include "MRBitSetParallelFor.h"
 #include "MRParallelFor.h"
 #include "MRTimer.h"
+#include "MRMeshTriPoint.h"
 #include <atomic>
 
 namespace MR
@@ -254,8 +255,10 @@ IsoLine Isoliner::extractOneLine_( EdgeId first, ContinueTrack continueTrack )
 
     if ( !addCrossedEdge( first ) )
         return res;
-    assert( activeEdges_.test( first.undirected() ) );
-        activeEdges_.reset( first.undirected() );
+
+    // looks like this assert is excessive for trackSection functions
+    //assert( activeEdges_.test( first.undirected() ) );
+    activeEdges_.reset( first.undirected() );
 
     bool closed = false;
     while ( auto next = findNextEdge_( res.back().e ) )
@@ -437,6 +440,39 @@ PlaneSection trackSection( const MeshPart& mp,
     return res;
 }
 
+Expected<PlaneSection> trackSection( const MeshPart& mp, const MeshTriPoint& start, const MeshTriPoint& end, const Vector3f& planePoint, bool ccw )
+{
+    MR_TIMER;
+
+    if ( fromSameTriangle( mp.mesh.topology, MeshTriPoint( start ), MeshTriPoint( end ) ) )
+        return PlaneSection();
+
+    auto startPoint = mp.mesh.triPoint( start );
+    auto endPoint = mp.mesh.triPoint( end );
+    auto crossDir = cross( startPoint - planePoint, endPoint - planePoint ).normalized();
+    auto plane = Plane3f::fromDirAndPt( ccw ? crossDir : -crossDir, planePoint );
+    VertMetric valueInVertex = [&] ( VertId v )
+    {
+        return plane.distance( mp.mesh.points[v] );
+    };
+    ContinueTrack continueTrack = [&] ( const MeshEdgePoint& next )->bool
+    {
+        return !fromSameTriangle( mp.mesh.topology, MeshTriPoint( next ), MeshTriPoint( end ) );
+    };
+    Isoliner s( mp.mesh.topology, valueInVertex, mp.region );
+    auto res = s.track( start, continueTrack );
+    if ( res.empty() )
+    {
+        assert( false );
+        return unexpected( "Empty section" );
+    }
+    if ( res.size() > 1 && res.front() == res.back() )
+        return unexpected( "Looped section" );
+    if ( !fromSameTriangle( mp.mesh.topology, MeshTriPoint( res.back() ), MeshTriPoint( end ) ) )
+        return unexpected( "Interrupted section" );
+    return res;
+}
+
 Contour2f planeSectionToContour2f( const Mesh& mesh, const PlaneSection& section, const AffineXf3f& meshToPlane )
 {
     MR_TIMER;

source/MRMesh/MRExtractIsolines.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include "MRMeshFwd.h"
+#include "MRExpected.h"
 
 namespace MR
 {
@@ -39,6 +40,13 @@ namespace MR
 [[nodiscard]] MRMESH_API PlaneSection trackSection( const MeshPart& mp,
     const MeshTriPoint& start, MeshTriPoint& end, const Vector3f& direction, float distance );
 
+/// track section of plane set by start point, end point and planePoint
+/// from start to end
+/// \param ccw - if true use start->end->planePoint plane, otherwise use start->planePoint->end (changes direction of plane tracking)
+/// returns track on surface without end point (return error if path was looped or reached boundary)
+[[nodiscard]] MRMESH_API Expected<PlaneSection> trackSection( const MeshPart& mp,
+    const MeshTriPoint& start, const MeshTriPoint& end, const Vector3f& planePoint, bool ccw );
+
 /// converts PlaneSections in 2D contours by computing coordinate of each point, applying given xf to it, and retaining only x and y
 [[nodiscard]] MRMESH_API Contour2f planeSectionToContour2f( const Mesh & mesh, const PlaneSection & section, const AffineXf3f & meshToPlane );
 [[nodiscard]] MRMESH_API Contours2f planeSectionsToContours2f( const Mesh & mesh, const PlaneSections & sections, const AffineXf3f & meshToPlane );