diff --git a/source/MRMesh/MR2DContoursTriangulation.cpp b/source/MRMesh/MR2DContoursTriangulation.cpp index 43c5cf35eff1..33a9221c5b52 100644 --- a/source/MRMesh/MR2DContoursTriangulation.cpp +++ b/source/MRMesh/MR2DContoursTriangulation.cpp @@ -18,6 +18,7 @@ #include "MRBitSetParallelFor.h" #include "MRPrecisePredicates2.h" #include "MRPrecisePredicates3.h" +#include "MRCutTrianglePredicates.h" #include "MRPch/MRTBB.h" #include #include @@ -51,6 +52,9 @@ struct SweepLinePredicates std::function addIntersectionPoint; // position of vertex v; tp is passed at call time because the queue's topology is moved into the output mesh std::function point; + // whether the predicates met points they cannot tell apart, after which their answers are arbitrary (empty if they cannot fail); + // findIntersections and makeMonotone check it after every event, the former also before the first one, and then stop with a failure + std::function failed; }; // the sweep-line predicates that depend only on the projected integer coordinates `pts2`; @@ -256,8 +260,10 @@ int findClosestToFront( const MeshTopology& tp, const SweepLinePredicates& predi break; i = nextI; } - assert( false ); - return 0; + // only failed predicates get here (CutTrianglePredicates on invalid points), whose run stops at the end of the event; + // 1 is a valid index for every caller, while 0 is the queried edge itself (a position of -1 for the start vertex) + assert( predicates.failed && predicates.failed() ); + return 1; } // same as above, reading the reference ray target from edges[0]'s dest (so edges[0] must be a valid @@ -304,6 +310,9 @@ struct SweepLineParams /// one can disable merge for identical vertices, merge is useful on symbol contours bool allowMerge{ true }; + /// if set, the loops of 2 edges (a slit: one edge taken both ways) are swept too, otherwise skipped; shorter loops are always skipped + bool allowTwoEdgeLoops{ false }; + /// optional out per-face winding numbers Vector* outFaceWinding{ nullptr }; @@ -334,9 +343,12 @@ class SweepLineQueue bool findIntersections(); void injectIntersections( IntersectionsMap* interMap ); - void makeMonotone(); + bool makeMonotone(); void triangulate(); private: + // whether the predicates failed, see SweepLinePredicates::failed + bool predicatesFailed_() const { return predicates_.failed && predicates_.failed(); } + // clears the cached buffers of the previous run, keeping their capacity void resetCache_(); @@ -458,6 +470,8 @@ class SweepLineQueue bool doSegmentSegmentIntersect_( VertId aOrg, VertId aDest, VertId bOrg, VertId bDest ) const; void checkIntersection_( int index, bool lower ); void checkIntersection_( int indexLower ); + // an intersection was found while params_.abortWhenIntersect, so the run fails + bool intersectionFound_{ false }; public: // all the buffers the sweep-line triangulation reuses between runs sharing the cache; @@ -580,14 +594,15 @@ MeshTopology* SweepLineQueue::runTopology( IntersectionsMap* interMap ) if ( !findIntersections() ) return nullptr; injectIntersections( interMap ); - makeMonotone(); + if ( !makeMonotone() ) + return nullptr; triangulate(); return &tp_; } bool SweepLineQueue::findIntersections() { - if ( multipleEdge_ ) + if ( multipleEdge_ || predicatesFailed_() ) return false; stage_ = Stage::Intersections; events_.clear(); @@ -599,11 +614,9 @@ bool SweepLineQueue::findIntersections() else if ( event.type == EventType::Destination ) processDestenationEvent_( event.index ); else if ( event.type == EventType::Intersection ) - { - if ( params_.abortWhenIntersect ) - return false; processIntersectionEvent_( event.index ); - } + if ( intersectionFound_ || predicatesFailed_() ) + return false; events_.push_back( event ); } return true; @@ -717,7 +730,7 @@ void SweepLineQueue::injectIntersections( IntersectionsMap* interMap ) e = newE; } -void SweepLineQueue::makeMonotone() +bool SweepLineQueue::makeMonotone() { stage_ = Stage::Monotonation; startVertIndex_ = 0; @@ -729,7 +742,10 @@ void SweepLineQueue::makeMonotone() else processDestenationEvent_( event.index ); calculateWinding_(); + if ( predicatesFailed_() ) + return false; } + return true; } void SweepLineQueue::triangulate() @@ -1145,6 +1161,11 @@ void SweepLineQueue::checkIntersection_( int i ) if ( !doSegmentSegmentIntersect_( org1, dest1, org2, dest2 ) ) return; + if ( params_.abortWhenIntersect ) + { + intersectionFound_ = true; // the run fails, so the point of the intersection is not needed + return; + } auto minEdgeId = std::min( activeSweepEdges_[i].edgeId, activeSweepEdges_[i + 1].edgeId ); auto maxEdgeId = std::max( activeSweepEdges_[i].edgeId, activeSweepEdges_[i + 1].edgeId ); @@ -1349,7 +1370,7 @@ void SweepLineQueue::initMeshByLoops_( const MeshTopology& inTp, const EdgeLoops for ( int loopId = 0; loopId < int( loops.size() ); ++loopId ) { const auto& loop = loops[loopId]; - if ( loop.size() < 3 ) + if ( loop.size() < ( params_.allowTwoEdgeLoops ? 2 : 3 ) ) continue; for ( int lId = 0; lId < loop.size() && !multipleEdge_; ++lId ) addNewEdge( loop[lId], loopId, lId ); @@ -1809,25 +1830,25 @@ MeshTopology* triangulateDisjointContoursTopology( const Mesh& mesh, const EdgeL return triangulator.runTopology(); } -std::optional getMonotonePlan( const Mesh& mesh, const EdgeLoops& loops, const Vector3f& normal ) +// the monotone plan of the region bounded by loops, with the predicates that makePredicates builds for the given patch-to-input edge map +template +static std::optional getMonotonePlanImpl( const Mesh& mesh, const EdgeLoops& loops, MakePredicates&& makePredicates ) { HoleFillPlan res; if ( loops.empty() ) return res; - for ( const auto& loop : loops ) - if ( loop.size() < 3 ) - return {}; // the sweep skips such loops entirely, so their holes would stay unaddressed SweepLineQueue::Cache cache; WholeEdgeMap patchToInEdges; // initMeshByLoops_ reserves it from the loop sizes std::vector chords; - SweepLineQueue triangulator( cache, mesh.topology, meshSpacePredicates( mesh, loops, normal, patchToInEdges, &cache.pts2Buffer ), loops, - { .abortWhenIntersect = true, .outPatchMap = &patchToInEdges, .outChords = &chords } ); + SweepLineQueue triangulator( cache, mesh.topology, makePredicates( patchToInEdges, cache ), loops, + { .abortWhenIntersect = true, .allowTwoEdgeLoops = true, .outPatchMap = &patchToInEdges, .outChords = &chords } ); if ( !triangulator.findIntersections() ) return {}; // passing findIntersections here means there is nothing to inject, so monotonation can run at // once; and it is where to stop - monotone parts are what the caller fills - triangulator.makeMonotone(); + if ( !triangulator.makeMonotone() ) + return {}; // an anchor is a patch edge: either a loop edge copied from the mesh (patchToInEdges maps it back), // or, past those, a chord referenced by the item that will create it @@ -1861,6 +1882,76 @@ std::optional getMonotonePlan( const Mesh& mesh, const EdgeLoops& return res; } +std::optional getMonotonePlan( const Mesh& mesh, const EdgeLoops& loops, const Vector3f& normal ) +{ + return getMonotonePlanImpl( mesh, loops, [&] ( const WholeEdgeMap& patchToInEdges, SweepLineQueue::Cache& cache ) + { + return meshSpacePredicates( mesh, loops, normal, patchToInEdges, &cache.pts2Buffer ); + } ); +} + +std::optional getMonotonePlan( const Mesh& mesh, const EdgeLoops& loops, + const std::array& tri, const std::function& pointOf ) +{ + // the distinct vertices of the loops (a vertex can be on several loops, or twice on a slit), sorted; vertex verts[i] is point i + // of the predicates, all added before the sweep, since no predicate may run after a failed addPoint + std::vector verts; + for ( const auto& loop : loops ) + for ( auto e : loop ) + verts.push_back( mesh.topology.org( e ) ); + std::sort( verts.begin(), verts.end() ); + verts.erase( std::unique( verts.begin(), verts.end() ), verts.end() ); + CutTrianglePredicates preds; + preds.reset( tri ); + // two vertices defined alike make the predicates fail only if the sweep compares them, unless the shifts of their edges tell them apart + // (as the crossings of the edges of a multiple edge of the other mesh with their wedges), so every two of them are compared before + // the sweep, and the loops are rejected if that fails: a vertex point by its id, a crossing by the ids of its segment and triangle + // in any order, with the index of the point + std::vector, int>> keys( verts.size() ); + for ( int i = 0; i < verts.size(); ++i ) + { + const auto d = pointOf( verts[i] ); + preds.addPoint( VertId( i ), d ); + auto& [k, pos] = keys[i]; + pos = i; + if ( d.isVertex() ) + { + k = { d.seg[0].id, -1, -1, -1, -1 }; + continue; + } + k = { d.seg[0].id, d.seg[1].id, d.tri[0].id, d.tri[1].id, d.tri[2].id }; + std::sort( k.begin(), k.begin() + 2 ); + std::sort( k.begin() + 2, k.end() ); + } + if ( preds.failed() ) + return {}; + std::sort( keys.begin(), keys.end() ); + for ( size_t i = 0; i < keys.size(); ++i ) + for ( size_t j = i + 1; j < keys.size() && keys[j].first == keys[i].first; ++j ) + (void)preds.less( VertId( keys[i].second ), VertId( keys[j].second ) ); + if ( preds.failed() ) + return {}; + + Vector pointIds; // the point of every vertex of the sweep + return getMonotonePlanImpl( mesh, loops, [&] ( const WholeEdgeMap&, SweepLineQueue::Cache& ) + { + // samePos, addIntersectionPoint and point stay empty: the sweep on loops merges no points, stops at the first crossing + // and triangulates nothing here, so it never calls them + SweepLinePredicates p; + p.less = [&] ( VertId l, VertId r ) { return preds.less( pointIds[l], pointIds[r] ); }; + p.ccw = [&] ( VertId a, VertId b, VertId c ) { return preds.ccw( pointIds[a], pointIds[b], pointIds[c] ); }; + p.ccwFromBehind = [&] ( VertId b, VertId pivot ) { return preds.ccwFromBehind( pointIds[b], pointIds[pivot] ); }; + p.addInputPoint = [&] ( VertId v, int contourId, int pointId ) + { + const auto it = std::lower_bound( verts.begin(), verts.end(), mesh.topology.org( loops[contourId][pointId] ) ); + pointIds.autoResizeSet( v, VertId( int( it - verts.begin() ) ) ); + }; + // the predicates can still fail on invalid points that the checks above miss, as three crossings of one cut on one line + p.failed = [&] { return preds.failed(); }; + return p; + } ); +} + } } diff --git a/source/MRMesh/MR2DContoursTriangulation.h b/source/MRMesh/MR2DContoursTriangulation.h index 7eb5837ef7de..4d22e402f996 100644 --- a/source/MRMesh/MR2DContoursTriangulation.h +++ b/source/MRMesh/MR2DContoursTriangulation.h @@ -8,6 +8,9 @@ namespace MR { +struct PreciseVertCoords; +struct PreciseCutPoint; + namespace PlanarTriangulation { @@ -139,8 +142,11 @@ MR_BIND_IGNORE MRMESH_API MeshTopology* triangulateDisjointContoursTopology( con /** * @brief splits the near-planar region bounded by \p loops in monotone parts, without triangulating them * this is the sweep-line triangulation stopped right after monotonation - * @param loops one closed EdgeLoop of at least 3 edges per boundary of one region, bounding holes of - * \p mesh as trackRightBoundaryLoop returns them; which way round a loop runs does not matter + * @param loops one closed EdgeLoop per boundary of one region, bounding holes of \p mesh as trackRightBoundaryLoop + * returns them (a slit of one edge is a loop of 2 edges, and a loop of one edge is skipped); which way round a loop runs + * does not matter; the holes can be several parts of one region split by edges between them, as all the holes of a face + * cut along contours (the parts between the contours that cross its sides, and both sides of every island), since an edge + * taken both ways adds zero winding * @param normal the orientation the region is monotone around, see \ref triangulateDisjointContours * @return std::nullopt if the loops intersect or a chord would land in a face instead of a hole, so that * the caller can fall back on filling every loop separately; otherwise a plan with numTris == 0 @@ -149,5 +155,15 @@ MR_BIND_IGNORE MRMESH_API MeshTopology* triangulateDisjointContoursTopology( con */ MRMESH_API std::optional getMonotonePlan( const Mesh& mesh, const EdgeLoops& loops, const Vector3f& normal ); +/// same as above for the region of one triangle of a mesh cut by the boolean (the loops of all the holes of the triangle), +/// on the exact predicates of \ref CutTrianglePredicates, which never contradict the boolean's ones, so the loops are disjoint +/// whenever the boolean's contours are +/// \param tri the triangle, counterclockwise around its normal +/// \param pointOf the definition of every vertex on the loops, valid as \ref CutTrianglePredicates::addPoint requires +/// \return std::nullopt also if two vertices on the loops are defined alike and their wedges do not tell them apart, or the predicates fail +/// on the points +MR_BIND_IGNORE MRMESH_API std::optional getMonotonePlan( const Mesh& mesh, const EdgeLoops& loops, + const std::array& tri, const std::function& pointOf ); + } } \ No newline at end of file diff --git a/source/MRMesh/MRContoursCut.cpp b/source/MRMesh/MRContoursCut.cpp index dd7fdb16ce76..9a49a51d4d9b 100644 --- a/source/MRMesh/MRContoursCut.cpp +++ b/source/MRMesh/MRContoursCut.cpp @@ -2,8 +2,11 @@ #include "MRMesh.h" #include "MRPolyline.h" #include "MRTriangleIntersection.h" +#include "MRTriMath.h" #include "MRMeshTopology.h" #include "MRMeshDelone.h" +#include "MR2DContoursTriangulation.h" +#include "MRCutTrianglePredicates.h" #include "MRRingIterator.h" #include "MRMeshFillHole.h" #include "MRTimer.h" @@ -53,6 +56,7 @@ struct PathsEdgeIndex { bool hasLeft{true}; bool hasRight{true}; + FaceId right; // the face on the right of a path edge of the old mesh, while removedFaces has the one on its left }; struct RemovedFaceInfo @@ -71,16 +75,20 @@ struct PreCutResult std::vector> oldEdgesInfo; }; +/// prepares precise-coordinates for a vertex of mesh A or B as the boolean's predicates take them: +/// the ids of B shifted by meshAVertsNum, its points moved by rigidB2A in float before the conversion +PreciseVertCoords preciseVert( VertId v, const Mesh& mesh, bool isA, const SortIntersectionsData& sortData ) +{ + if ( isA ) + return { v, sortData.converter( mesh.points[v] ) }; + const auto& p = mesh.points[v]; + return { v + int( sortData.meshAVertsNum ), sortData.converter( sortData.rigidB2A ? ( *sortData.rigidB2A )( p ) : p ) }; +} + /// prepares precise-coordinates for a vertex from the other mesh PreciseVertCoords preciseOtherVert( VertId v, const SortIntersectionsData& sortData ) { - if ( sortData.isOtherA ) - return { v, sortData.converter( sortData.otherMesh.points[v] ) }; - - if ( !sortData.rigidB2A ) - return { v + int( sortData.meshAVertsNum ), sortData.converter( sortData.otherMesh.points[v] )}; - - return { v + int( sortData.meshAVertsNum ),sortData.converter( ( *sortData.rigidB2A )( sortData.otherMesh.points[v] ) )}; + return preciseVert( v, sortData.otherMesh, sortData.isOtherA, sortData ); } // sets left face of given edge invalid and saves info about its old left ring @@ -236,6 +244,7 @@ PreCutResult doPreCutMesh( Mesh& mesh, const OneMeshContours& contours ) newEdgeId = e; oldEdgesInfo[intersectionId].hasLeft = mesh.topology.left( newEdgeId ).valid(); oldEdgesInfo[intersectionId].hasRight = mesh.topology.right( newEdgeId ).valid(); + oldEdgesInfo[intersectionId].right = mesh.topology.right( newEdgeId ); break; } } @@ -458,6 +467,179 @@ void fixOrphans( Mesh& mesh, const std::vector& paths, const FullRemov } } +// a lone contour is inside one face crossing none of its edges (e.g. where a tip of the other mesh pokes through it), +// so it is an island in the hole of its face, or a slit if it is open +struct LoneFace +{ + FaceId f; + std::array edges; // counterclockwise, before the cut + ThreeVertIds verts; // verts[i] = org( edges[i] ) +}; + +// finds the faces of lone contours, before the cut changes their topology +std::vector findLoneFaces( const MeshTopology& topology, const OneMeshContours& contours ) +{ + std::vector faces; + for ( const auto& c : contours ) + { + if ( c.intersections.size() > 1 && std::all_of( c.intersections.begin(), c.intersections.end(), + [] ( const OneMeshIntersection& i ) { return std::holds_alternative( i.primitiveId ); } ) ) + faces.push_back( std::get( c.intersections.front().primitiveId ) ); + } + std::sort( faces.begin(), faces.end() ); + faces.erase( std::unique( faces.begin(), faces.end() ), faces.end() ); + std::vector res( faces.size() ); + for ( int i = 0; i < faces.size(); ++i ) + { + auto& lf = res[i]; + lf.f = faces[i]; + topology.getTriEdges( lf.f, lf.edges[0], lf.edges[1], lf.edges[2] ); + for ( int k = 0; k < 3; ++k ) + lf.verts[k] = topology.org( lf.edges[k] ); + } + return res; +} + +// connects the islands of lone contours with the rest of the holes of their faces by the chords that split the holes in monotone parts, +// found on the exact predicates consistent with the boolean if sortData is given; the faces where that fails are added to badFaces +// and stay in the mesh unless another contour removed them +void connectLoneContours( Mesh& mesh, const OneMeshContours& contours, const std::vector& loneFaces, const PreCutResult& preRes, + UndirectedEdgeId firstNewEdge, const SortIntersectionsData* sortData, NewEdgesMap* new2OldEdgeMap, FaceBitSet& badFaces ) +{ + if ( loneFaces.empty() ) + return; + MR_TIMER; + + // the path edges inside every lone face, with the face on the left, and the contour points at their ends for the exact definitions + std::vector faceEdges( loneFaces.size() ); + HashMap> vertInters; + for ( int i = 0; i < preRes.paths.size(); ++i ) + { + for ( int j = 0; j < preRes.paths[i].size(); ++j ) + { + const auto f = preRes.removedFaces[i][j].f; + const auto it = std::lower_bound( loneFaces.begin(), loneFaces.end(), f, [] ( const LoneFace& lf, FaceId f ) { return lf.f < f; } ); + if ( it == loneFaces.end() || it->f != f ) + continue; + auto& edges = faceEdges[it - loneFaces.begin()]; + const auto e = preRes.paths[i][j]; + edges.push_back( e ); + if ( e.undirected() >= firstNewEdge ) // a new edge is inside the face, while an old one has it on the left only + edges.push_back( e.sym() ); + if ( !sortData ) + continue; + vertInters[mesh.topology.org( e )] = { i, j }; + vertInters[mesh.topology.dest( e )] = { i, j + 1 }; + } + } + + auto ownVert = [&] ( VertId v ) { return preciseVert( v, mesh, !sortData->isOtherA, *sortData ); }; + // the exact definition of a vertex on the loops of a lone face + auto pointOf = [&] ( const LoneFace& lf, VertId v ) + { + PreciseCutPoint res; + const auto it = vertInters.find( v ); + if ( it == vertInters.end() ) + { + res.seg[0] = ownVert( v ); // a corner of the face + return res; + } + const auto [c, i] = it->second; + const auto& prim = contours[c].intersections[i].primitiveId; + const auto& et = sortData->contours[c][i]; + const auto& otherTopology = sortData->otherMesh.topology; + if ( const auto pe = std::get_if( &prim ) ) + { + // a side of the face crosses a triangle of the other mesh + for ( int k = 0; k < 3; ++k ) + if ( lf.edges[k].undirected() == pe->undirected() ) + res.seg = { ownVert( lf.verts[k] ), ownVert( lf.verts[( k + 1 ) % 3] ) }; + const auto t = otherTopology.getTriVerts( et.tri() ); + res.tri = { preciseOtherVert( t[0], *sortData ), preciseOtherVert( t[1], *sortData ), preciseOtherVert( t[2], *sortData ) }; + } + else if ( std::holds_alternative( prim ) ) + { + // an edge of the other mesh crosses the face; the third vertices of its triangles tell apart the crossings of a multiple edge + const EdgeId e = et.edge; + res.seg = { preciseOtherVert( otherTopology.org( e ), *sortData ), preciseOtherVert( otherTopology.dest( e ), *sortData ) }; + res.tri = { ownVert( lf.verts[0] ), ownVert( lf.verts[1] ), ownVert( lf.verts[2] ) }; + if ( otherTopology.left( e ) && otherTopology.right( e ) ) + res.wedge = { preciseOtherVert( otherTopology.getLeftTriVerts( e )[2], *sortData ), + preciseOtherVert( otherTopology.getLeftTriVerts( e.sym() )[2], *sortData ) }; + } + else + res.seg[0] = ownVert( std::get( prim ) ); + return res; + }; + + // the boundaries of the hole of every lone face + std::vector faceLoops( loneFaces.size() ); + EdgeBitSet inLoops( mesh.topology.edgeSize() ); + BitSet removedHere( loneFaces.size() ); // the faces that no other contour removed, to put back if their plans fail + for ( int i = 0; i < loneFaces.size(); ++i ) + { + const auto& lf = loneFaces[i]; + auto& edges = faceEdges[i]; + if ( mesh.topology.hasFace( lf.f ) ) // no other contour removed it + { + mesh.topology.setLeft( lf.edges[0], {} ); + edges.push_back( lf.edges[0] ); + removedHere.set( i ); + } + for ( auto e : edges ) + { + if ( mesh.topology.left( e ) || inLoops.test( e ) ) + continue; + auto& loop = faceLoops[i].emplace_back(); + for ( auto le : leftRing( mesh.topology, e ) ) + { + inLoops.set( le ); + loop.push_back( le ); + } + } + } + + // the faces are independent, so plan them in parallel + std::vector> plans( loneFaces.size() ); + ParallelFor( plans, [&] ( size_t i ) + { + const auto& lf = loneFaces[i]; + const auto& loops = faceLoops[i]; + if ( sortData ) + { + plans[i] = PlanarTriangulation::getMonotonePlan( mesh, loops, { ownVert( lf.verts[0] ), ownVert( lf.verts[1] ), ownVert( lf.verts[2] ) }, + [&] ( VertId v ) { return pointOf( lf, v ); } ); + return; + } + // the loops go counterclockwise around the normal of the face, of which only the dominant component and its sign matter; + // the corners in doubles give a normal free of the rounding of the islands' points, which a sum over the loops would add; + // normalized before the conversion to float, which could flush it to zero or overflow at extreme scales + const auto& p = mesh.points; + const auto normal = dirDblArea( Triangle3d{ Vector3d( p[lf.verts[0]] ), Vector3d( p[lf.verts[1]] ), Vector3d( p[lf.verts[2]] ) } ); + // corners on one line give no normal, so no plane to connect the islands in: the plan stays empty, and the face is reported + if ( normal == Vector3d() ) + return; + plans[i] = PlanarTriangulation::getMonotonePlan( mesh, loops, Vector3f( normal.normalized() ) ); + } ); + + for ( int i = 0; i < loneFaces.size(); ++i ) + { + auto& plan = plans[i]; + // every island or slit needs a chord to the rest of its face, so a plan without any leaves them unconnected + // (as where the cut of a side has already filled the face again, ending an open contour on that side) + if ( plan && !plan->items.empty() && isFillingMultipleEdgeFree( mesh.topology, *plan ) ) + executeTriangulateContourPlan( mesh, faceLoops[i].front().front(), *plan, loneFaces[i].f, nullptr, new2OldEdgeMap ); + else + { + badFaces.autoResizeSet( loneFaces[i].f ); + // the holes to fill are found by the cut paths, and the old sides of the face have none, so the face is put back; + // no chord of it was added, so its old sides still make its ring + if ( removedHere.test( i ) ) + mesh.topology.setLeft( loneFaces[i].edges[0], loneFaces[i].f ); + } + } +} + [[maybe_unused]] void debugSortingInfo( EdgeId baseE, const EdgeData& edgeData, const std::vector& res, @@ -510,27 +692,10 @@ CompareResult compareEdgeTrisIntersections( EdgeId intersectedEdge, const Mesh& if ( vsl == vsr ) return CompareResult::Unknown; - VertId vo = mesh.topology.org( intersectedEdge ); - VertId vd = mesh.topology.dest( intersectedEdge ); - - auto po = mesh.points[vo]; - auto pd = mesh.points[vd]; - - if ( sortData.isOtherA ) - { - vo += (int)sortData.meshAVertsNum; - vd += (int)sortData.meshAVertsNum; - if ( sortData.rigidB2A ) - { - po = (*sortData.rigidB2A)( po ); - pd = (*sortData.rigidB2A)( pd ); - } - } - std::array preciseVerts { - PreciseVertCoords{ vo, sortData.converter( po ) }, - PreciseVertCoords{ vd, sortData.converter( pd ) }, + preciseVert( mesh.topology.org( intersectedEdge ), mesh, !sortData.isOtherA, sortData ), + preciseVert( mesh.topology.dest( intersectedEdge ), mesh, !sortData.isOtherA, sortData ), preciseOtherVert( vsl[0], sortData ), preciseOtherVert( vsl[1], sortData ), preciseOtherVert( vsl[2], sortData ), @@ -1117,6 +1282,8 @@ CutMeshResult cutMesh( Mesh& mesh, const OneMeshContours& contours, const CutMes if ( params.new2OldMap ) prepareFacesMap( mesh.topology, *params.new2OldMap ); + const auto loneFaces = findLoneFaces( mesh.topology, contours ); + const auto firstNewEdge = UndirectedEdgeId( mesh.topology.undirectedEdgeSize() ); auto preRes = doPreCutMesh( mesh, contours ); if ( params.new2oldEdgesMap ) @@ -1132,9 +1299,11 @@ CutMeshResult cutMesh( Mesh& mesh, const OneMeshContours& contours, const CutMes } cutEdgesIntoPieces( mesh, std::move( preRes.edgeData ), contours, params.sortData, params.new2OldMap, params.new2oldEdgesMap ); + FaceBitSet badLoneFaces; + connectLoneContours( mesh, contours, loneFaces, preRes, firstNewEdge, params.sortData, params.new2oldEdgesMap, badLoneFaces ); fixOrphans( mesh, preRes.paths, preRes.removedFaces, params.new2OldMap, params.new2oldEdgesMap ); - res.fbsWithContourIntersections = getBadFacesAfterCut( mesh.topology, preRes, preRes.removedFaces ); + res.fbsWithContourIntersections = getBadFacesAfterCut( mesh.topology, preRes, preRes.removedFaces ) | badLoneFaces; if ( params.forceFillMode == CutMeshParameters::ForceFill::None && res.fbsWithContourIntersections.any() ) return res; @@ -1182,16 +1351,18 @@ CutMeshResult cutMesh( Mesh& mesh, const OneMeshContours& contours, const CutMes continue; bool addLeftDesc = oldEdgesInfo[edgeId].hasLeft && !mesh.topology.left( path[edgeId] ).valid(); bool addRightDesc = oldEdgesInfo[edgeId].hasRight && !mesh.topology.right( path[edgeId] ).valid(); + // the hole on the right of an edge of the old mesh is in its right face, which another contour or a lone one removed later + const FaceId oldRightF = oldEdgesInfo[edgeId].right ? oldEdgesInfo[edgeId].right : oldf; // order does matter because bad holes have tunnels and passes through both left and right if ( !fillLeft ) { - if ( addRightDesc ) addHoleDesc( path[edgeId].sym(), oldf, !fillRight ); + if ( addRightDesc ) addHoleDesc( path[edgeId].sym(), oldRightF, !fillRight ); if ( addLeftDesc ) addHoleDesc( path[edgeId], oldf, true ); // pseudo hole should always be second } else { if ( addLeftDesc ) addHoleDesc( path[edgeId], oldf, false ); // real hole should alwauys be first - if ( addRightDesc ) addHoleDesc( path[edgeId].sym(), oldf, !fillRight ); + if ( addRightDesc ) addHoleDesc( path[edgeId].sym(), oldRightF, !fillRight ); } } } diff --git a/source/MRMesh/MRContoursCut.h b/source/MRMesh/MRContoursCut.h index 46c77f0ffe74..1627e47b625e 100644 --- a/source/MRMesh/MRContoursCut.h +++ b/source/MRMesh/MRContoursCut.h @@ -78,7 +78,8 @@ struct CutMeshResult /** \ingroup BooleanGroup * \brief Cuts mesh by given contours * - * This function cuts mesh making new edges paths on place of input contours + * This function cuts mesh making new edges paths on place of input contours; + * a contour inside one face is connected with the rest of the face by new edges * \param mesh Input mesh that will be cut * \param contours Input contours to cut mesh with, find more \ref MR::OneMeshContours * \param params Parameters describing some cut options, find more \ref MR::CutMeshParameters diff --git a/source/MRMesh/MRCutTrianglePredicates.cpp b/source/MRMesh/MRCutTrianglePredicates.cpp index ce02ac968ecc..42da6ede8767 100644 --- a/source/MRMesh/MRCutTrianglePredicates.cpp +++ b/source/MRMesh/MRCutTrianglePredicates.cpp @@ -96,6 +96,45 @@ bool sameTriangle( const std::array & a, const std::array< return hasVert( b, a[0].id ) && hasVert( b, a[1].id ) && hasVert( b, a[2].id ); } +/// whether p and q are defined alike: crossings of one segment, in either direction, with one triangle, with its vertices in any order +bool definedAlike( const PreciseCutPoint & p, const PreciseCutPoint & q ) +{ + if ( p.isVertex() || q.isVertex() || !sameTriangle( p.tri, q.tri ) ) + return false; + const auto & s = p.seg; + const auto & t = q.seg; + return ( s[0].id == t[0].id && s[1].id == t[1].id ) || ( s[0].id == t[1].id && s[1].id == t[0].id ); +} + +/// whether the segment of p has a wedge: two distinct third vertices of its triangles +bool hasWedge( const PreciseCutPoint & p ) +{ + return p.wedge[0].id.valid() && p.wedge[1].id.valid() && p.wedge[0].id != p.wedge[1].id; +} + +/// whether the half-plane from line ab through y is in the angle below 180 degrees between the half-planes through x and z, +/// on the sides of planes abx and abz where the boolean's orient3d puts it +bool insideAngle( const std::array & ab, const PreciseVertCoords & x, const PreciseVertCoords & z, const PreciseVertCoords & y ) +{ + const auto & [a, b] = ab; + return orient3d( { a, b, x, y } ) == orient3d( { a, b, x, z } ) && orient3d( { a, b, z, y } ) == orient3d( { a, b, z, x } ); +} + +/// whether the triangles of the edge of q, defined alike as p, are in the angle below 180 degrees between the triangles of the edge of p: +/// each of them shared with p or in that angle, not both shared +bool nestedIn( const PreciseCutPoint & p, const PreciseCutPoint & q ) +{ + int numShared = 0; + for ( const auto & y : q.wedge ) + { + if ( y.id == p.wedge[0].id || y.id == p.wedge[1].id ) + ++numShared; + else if ( !insideAngle( p.seg, p.wedge[0], p.wedge[1], y ) ) + return false; + } + return numShared < 2; +} + #ifndef NDEBUG /// whether every coordinate of v is below 2^30 by absolute value, for which the integer widths and the error bounds here are derived bool inRange( const PreciseVertCoords & v ) @@ -516,11 +555,15 @@ struct CutTrianglePredicates::Impl /// whether the maps above may have entries, since clearing the empty submaps of the parallel map takes a lock each std::atomic cached{ false }; + /// whether the maps may have entries found by the shifts of the points defined alike, which depend on the other points defined alike + std::atomic shifted{ false }; + std::atomic failed{ false }; /// forgets the orders and orientations found by the slower rules void clearCaches() { + shifted.store( false, std::memory_order_relaxed ); if ( !cached.exchange( false, std::memory_order_relaxed ) ) return; dirs.clear(); @@ -528,6 +571,14 @@ struct CutTrianglePredicates::Impl } }; +/// the shift of the edge of a point defined alike with others (see CutTrianglePredicates): x and y are the projected coordinates +/// of the direction in which the point moves, times scale > 0, which is the same for all points defined alike +struct CutTrianglePredicates::Shift +{ + FastInt256 x, y; + FastInt128 scale; +}; + CutTrianglePredicates::CutTrianglePredicates() : impl_( std::make_unique() ) { } @@ -605,6 +656,7 @@ void CutTrianglePredicates::reset( const std::array & tri void CutTrianglePredicates::addPoint( VertId v, const PreciseCutPoint & p ) { assert( inRange( p.seg[0] ) && ( p.isVertex() || ( inRange( p.seg[1] ) && std::all_of( p.tri.begin(), p.tri.end(), inRange ) ) ) ); + assert( !hasWedge( p ) || ( inRange( p.wedge[0] ) && inRange( p.wedge[1] ) ) ); auto & points = impl_->points; if ( points.size() <= v ) { @@ -612,8 +664,10 @@ void CutTrianglePredicates::addPoint( VertId v, const PreciseCutPoint & p ) approx_.resize( v + 1 ); } auto & pt = points[v]; - if ( pt.def.seg[0].id.valid() ) - impl_->clearCaches(); // the point is set again: forget the orders found for its old definition + // the point is set again: forget the orders found for its old definition; or a point with a wedge may change the shifts + // of the points defined alike as it, and so the orders found by them + if ( pt.def.seg[0].id.valid() || ( impl_->shifted.load( std::memory_order_relaxed ) && hasWedge( p ) ) ) + impl_->clearCaches(); pt.def = p; auto corner = [&] ( const PreciseVertCoords & x ) { @@ -724,25 +778,16 @@ bool CutTrianglePredicates::less( VertId p, VertId q ) const return c < 0; if ( const auto dir = coincidentDirection_( p, q ) ) return sweepForward( *dir ); - return p < q; // the failure: points defined alike, ordered by their ids + return p < q; // the failure: points defined alike, not told apart by their shifts, ordered by their ids } /// the signs of the projected coordinates of perturbed q - p for points p and q that coincide before the perturbation, -/// as the sweep sees them (see sweepSigns); nullopt and the failure if they coincide after the perturbation too +/// as the sweep sees them (see sweepSigns), or of the difference of their shifts for points defined alike; +/// nullopt and the failure if they coincide after the perturbation and the shifts too std::optional CutTrianglePredicates::coincidentDirection_( VertId p, VertId q ) const { const auto & a = impl_->points[p]; const auto & b = impl_->points[q]; - if ( !a.def.isVertex() && !b.def.isVertex() && sameTriangle( a.def.tri, b.def.tri ) ) - { - const auto & s = a.def.seg; - const auto & t = b.def.seg; - if ( ( s[0].id == t[0].id && s[1].id == t[1].id ) || ( s[0].id == t[1].id && s[1].id == t[0].id ) ) - { - fail_(); // the same point twice - return {}; - } - } // a vertex and another point of a side from it, a crossing strictly inside the side or the vertex at its other end: q - p goes // along side i in the direction that reset found if it goes from vertex i or to vertex i+1; the same vertex twice has two sides in common if ( const int common = a.sides & b.sides; std::has_single_bit( unsigned( common ) ) && ( a.def.isVertex() || b.def.isVertex() ) ) @@ -752,10 +797,10 @@ std::optional CutTrianglePredicates::coincidentDirection_( VertId p, V return forward ? sideDirs_[i] : -sideDirs_[i]; } // the other pairs by the direction of the perturbed points, which the slower rules find and the cache keeps: two crossings of one side - // as the boolean's cut orders them along it (see directionFromSide_) - const auto dir = perturbedDirection_( p, q ); + // as the boolean's cut orders them along it (see directionFromSide_), and the same point twice by the shifts of its edges + const auto dir = perturbedDirection_( p, q, definedAlike( a.def, b.def ) ); if ( !dir ) - fail_(); // the points coincide after the perturbation too + fail_(); // the points coincide after the perturbation and the shifts too return dir; } @@ -803,17 +848,23 @@ std::optional CutTrianglePredicates::directionFromSide_( VertId p, Ver return {}; } -std::optional CutTrianglePredicates::perturbedDirection_( VertId p, VertId q ) const +std::optional CutTrianglePredicates::perturbedDirection_( VertId p, VertId q, bool alike ) const { const bool swapped = q < p; const auto key = swapped ? ( std::uint64_t( q ) << 32 ) | std::uint64_t( p ) : ( std::uint64_t( p ) << 32 ) | std::uint64_t( q ); Vector2i dir; if ( !impl_->dirs.if_contains( key, [&] ( const auto & x ) { dir = x.second; } ) ) { - // a side gives the direction faster than the expansion of the perturbation, if it applies - auto res = swapped ? directionFromSide_( q, p ) : directionFromSide_( p, q ); - if ( !res ) - res = swapped ? directionByPerturbation( q, p ) : directionByPerturbation( p, q ); + std::optional res; + if ( alike ) + res = swapped ? shiftDirection_( q, p ) : shiftDirection_( p, q ); // the perturbation moves both points alike + else + { + // a side gives the direction faster than the expansion of the perturbation, if it applies + res = swapped ? directionFromSide_( q, p ) : directionFromSide_( p, q ); + if ( !res ) + res = swapped ? directionByPerturbation( q, p ) : directionByPerturbation( p, q ); + } dir = res.value_or( Vector2i() ); impl_->dirs.emplace( key, dir ); impl_->cached.store( true, std::memory_order_relaxed ); @@ -828,6 +879,61 @@ std::optional CutTrianglePredicates::directionByPerturbation( VertId p return sosDirection( impl_->points[p].def, impl_->points[q].def, kx_, ky_ ); } +/// the signs of the projected coordinates of q - p for points p and q defined alike: of the difference of their shifts, with the same scale; +/// nullopt if they have no shifts, or equal ones +std::optional CutTrianglePredicates::shiftDirection_( VertId p, VertId q ) const +{ + const auto sp = shift_( p ), sq = shift_( q ); + if ( !sp || !sq ) + return {}; + const Vector2i dir( ( sq->x - sp->x ).sign(), ( sq->y - sp->y ).sign() ); + if ( dir == Vector2i() ) + return {}; + return dir; +} + +/// the shift of crossing p of an edge of another mesh with the triangle, by which the predicates tell apart the points defined alike: +/// the edge ab moves by delta * u, u = w0 + w1 - a - b toward its wedge w0, w1, by delta infinitely smaller than every term of the perturbation +/// of the vertices, and so its crossing with the plane of normal n by delta * ( u - ( n.u ) / ( n.d ) * d ), d = b - a, along the cut of the plane +/// by the plane of ab and u, into the angle between the cuts by the triangles of the edge (the same for either direction of the segment and +/// either orientation of the triangle); or by -delta * u if the triangles of another edge defined alike are in that angle, so that the sheets +/// of the edges that touch along the segment give contours that do not cross; nullopt without a wedge, or if the segment is parallel to the plane +std::optional CutTrianglePredicates::shift_( VertId p ) const +{ + const auto & points = impl_->points; + const auto & pt = points[p]; + const auto & def = pt.def; + // the segment of a crossing of a side is that side + if ( def.isVertex() || pt.sides != 0 || !hasWedge( def ) ) + return {}; + const Vector3i64 a( def.seg[0].pt ), b( def.seg[1].pt ); + const auto d = b - a; + const auto u = Vector3i64( def.wedge[0].pt ) + Vector3i64( def.wedge[1].pt ) - a - b; + // with the components of the normal below 2^63 (see normalOf), those of d below 2^31 and those of u below 2^32, + // n.d and n.u are below 3 * 2^95, and the projected coordinates below 6 * 2^126 + auto dot = [&] ( const Vector3i64 & v ) + { + return Int64Mul128( normal_.x ) * Int64Mul128( v.x ) + Int64Mul128( normal_.y ) * Int64Mul128( v.y ) + + Int64Mul128( normal_.z ) * Int64Mul128( v.z ); + }; + const FastInt128 nd = dot( d ), nu = dot( u ); + if ( nd == 0 ) + return {}; // the segment is in the plane before the perturbation, or the triangle has a zero normal + auto coord = [&] ( int k ) { return Int128Mul256( nd ) * Int128Mul256( u[k] ) - Int128Mul256( nu ) * Int128Mul256( d[k] ); }; + Shift res{ coord( kx_ ), coord( ky_ ), nd < 0 ? -nd : nd }; + // away from the wedge if the triangles of another point defined alike are in its angle; found only here, since points defined alike are rare + bool away = false; + for ( auto x = points.beginId(); x < points.endId() && !away; ++x ) + away = x != p && points[x].sides == 0 && hasWedge( points[x].def ) && definedAlike( def, points[x].def ) && nestedIn( def, points[x].def ); + if ( ( nd < 0 ) != away ) + { + res.x = -res.x; + res.y = -res.y; + } + impl_->shifted.store( true, std::memory_order_relaxed ); + return res; +} + std::optional CutTrianglePredicates::commonLine_( VertId p, VertId q, bool withVertex ) const { const auto & a = impl_->points[p]; @@ -889,43 +995,32 @@ bool CutTrianglePredicates::along_( const Vector2i & dir, VertId p, VertId q ) c return less( p, q ) == sweepForward( dir ); } -bool CutTrianglePredicates::sideOfPlane_( const std::array & plane, VertId x ) const +/// 1 or -1 as point x is on the side of the plane where orient3d is true or false, or 0 if it is on the plane, so on the line +int CutTrianglePredicates::sideOfPlane_( const std::array & plane, VertId x ) const { const auto & p = impl_->points[x].def; auto orient = [&] ( const PreciseVertCoords & v ) { return orient3d( { plane[0], plane[1], plane[2], v } ); }; + auto side = [] ( bool positive ) { return positive ? 1 : -1; }; if ( p.isVertex() ) - { - if ( hasVert( plane, p.seg[0].id ) ) - { - fail_(); // a vertex of the plane, so on the line - return false; - } - return orient( p.seg[0] ); - } + return hasVert( plane, p.seg[0].id ) ? 0 : side( orient( p.seg[0] ) ); // a vertex of the plane is on it // an end of the segment in the plane: the segment meets the plane only there, so the point is on the side of the other end const bool in0 = hasVert( plane, p.seg[0].id ); const bool in1 = hasVert( plane, p.seg[1].id ); if ( in0 && in1 ) - { - fail_(); // the segment is in the plane, so the point is on the line - return false; - } + return 0; // the segment is in the plane if ( in0 ) - return orient( p.seg[1] ); + return side( orient( p.seg[1] ) ); if ( in1 ) - return orient( p.seg[0] ); + return side( orient( p.seg[0] ) ); // the point is on the segment, so on the side of both ends if they are on one side const bool o0 = orient( p.seg[0] ); if ( o0 == orient( p.seg[1] ) ) - return o0; + return side( o0 ); if ( sameTriangle( p.tri, plane ) ) - { - fail_(); // the point is on the plane, so on the line - return false; - } + return 0; // the point is on the plane // the segment crosses the plane, so the point is on the side of the end it meets first if it meets its triangle before the plane; // the boolean's predicate takes any plane through three vertices with distinct ids, as the made-up planes vcd and dce of commonLine_ are, // that shares at most two of them with the triangle and none with the segment, as the checks above and in addPoint ensure @@ -935,7 +1030,7 @@ bool CutTrianglePredicates::sideOfPlane_( const std::array if ( doTriangleSegmentIntersect( { plane[0], plane[1], plane[2], p.seg[0], p.seg[1] } ) ) assert( res == ( segmentIntersectionOrder( { p.seg[0], p.seg[1], p.tri[0], p.tri[1], p.tri[2], plane[0], plane[1], plane[2] } ) == o0 ) ); #endif - return res; + return side( res ); } int CutTrianglePredicates::exactOrient_( VertId a, VertId b, VertId c ) const @@ -992,7 +1087,10 @@ int CutTrianglePredicates::orientSlow_( VertId a, VertId b, VertId c ) const return *res ? 1 : -1; if ( const int s = exactOrient_( a, b, c ) ) return s; - // the points are on one line before the perturbation, and no two of them on a known one: the expansion of the perturbation decides + // the points are on one line before the perturbation, and no two of them on a known one: the expansion of the perturbation decides, + // unless two of them are defined alike, which it moves alike, and which only their shifts tell apart + if ( const auto res = ccwByShifts_( a, b, c ) ) + return *res ? 1 : -1; const auto & points = impl_->points; if ( const int s = sosOrient( points[a].def, points[b].def, points[c].def, kx_, ky_ ) ) return s; @@ -1053,22 +1151,94 @@ std::optional CutTrianglePredicates::ccwByPlanes_( VertId a, VertId b, Ver { assert( l->side < 0 ); // no two of the points on one side, see ccwBySides_ // the left of the line in its direction is the side of the plane opposite to the plane's normal, which is where orient3d is true - return along_( lineDirection_( *l ), p, q ) == sideOfPlane_( l->plane, r ); + if ( const int side = sideOfPlane_( l->plane, r ) ) + return along_( lineDirection_( *l ), p, q ) == ( side > 0 ); + // r is on the line too, which in valid input only a point defined alike as p or q is: their shifts decide + if ( const auto res = ccwByShifts_( a, b, c ) ) + return res; + fail_(); // three points on one line also after the perturbation + return false; } } return {}; // no two of the points are on a known line } +/// the orientation of triangle abc with two points defined alike, by the shifts of their edges (see shift_), or nullopt if no two of them are; +/// the failure if the shifts do not tell it either +std::optional CutTrianglePredicates::ccwByShifts_( VertId a, VertId b, VertId c ) const +{ + const auto & points = impl_->points; + // rotating the vertices keeps the orientation of the triangle, so the two defined alike can come first + const std::array v{ a, b, c }; + int i = 0; + while ( i < 3 && !definedAlike( points[v[i]].def, points[v[( i + 1 ) % 3]].def ) ) + ++i; + if ( i == 3 ) + return {}; + const VertId p = v[i], q = v[( i + 1 ) % 3], r = v[( i + 2 ) % 3]; + const auto sp = shift_( p ), sq = shift_( q ); + if ( !sp || !sq ) + { + fail_(); // nothing tells apart p and q + return false; + } + // p and q at the common point P are shifted by delta * wp and delta * wq, and r by delta * wr, so the doubled area of the triangle is + // delta * cross( wq - wp, r - P ) + delta^2 * cross( wq - wp, wr - wp ), where r - P is perturbed, and all its terms come before delta + const FastInt256 dx = sq->x - sp->x, dy = sq->y - sp->y; + if ( !definedAlike( points[p].def, points[r].def ) ) + { + // r - P = ( f.nx * e.d - e.nx * f.d, f.ny * e.d - e.ny * f.d ) / ( e.d * f.d ) with positive denominators, and the numerators + // below 2^223 (see exactCoords), times the coordinates of the difference of the shifts below 2^130 + const auto & e = points[p].exact; + const auto & f = points[r].exact; + using M = Int128Mul256; + const FastInt256 ex = M( f.nx ) * M( e.d ) - M( e.nx ) * M( f.d ), ey = M( f.ny ) * M( e.d ) - M( e.ny ) * M( f.d ); + if ( const int s = ( dx * ey - dy * ex ).sign() ) + return s > 0; + // r - P is parallel to the difference of the shifts before the perturbation, so the expansion of the perturbation decides + auto constant = [] ( const FastInt256 & x ) { return x.sign() ? SosPoly256( std::vector{ { 0, x } } ) : SosPoly256(); }; + const SosPoly256 wx = constant( dx ), wy = constant( dy ); + const int sos = sosFirstStage<2>( { &points[p].def, &points[r].def }, kx_, ky_, [&] ( const auto & pts, std::int64_t from ) -> std::optional + { + const SosFraction & pp = *pts[0]; + const SosFraction & rr = *pts[1]; + const SosPoly512 nx = mulAs( rr.nx, pp.d ) - mulAs( pp.nx, rr.d ); + const SosPoly512 ny = mulAs( rr.ny, pp.d ) - mulAs( pp.ny, rr.d ); + const SignedProduct products[2] = { { &wx, &ny, false }, { &wy, &nx, true } }; + if ( const int sign = lowestTermOfSumOfProducts( products, 2, from ).second ) + return sign; + return {}; + } ).value_or( 0 ); + if ( sos ) + return sos > 0; + } + // the second-order term, with the shift of r of its own scale (zero for a point without a shift): + // wr - wp = ( Wr * scale - Wp * rScale ) / ( scale * rScale ), the products below 2^129 * 2^96 + FastInt256 rx, ry; + FastInt128 rScale = 1; + if ( const auto sr = shift_( r ) ) + { + rx = sr->x; + ry = sr->y; + rScale = sr->scale; + } + const FastInt256 fx( rx * sp->scale - sp->x * rScale ), fy( ry * sp->scale - sp->y * rScale ); + if ( const int s = ( dx * fy - dy * fx ).sign() ) + return s > 0; + fail_(); // the shifts are on one line too + return false; +} + bool CutTrianglePredicates::ccwFromBehind( VertId b, VertId pivot ) const { // the triangle ( pivot - R * ( 1, eta ), b, pivot ) for the sweep direction ( 1, eta ) has the doubled area // R * ( eta * d.x - d.y ) for d = b - pivot, so d.y decides, and only an exactly zero d.y leaves d.x to decide, since eps << eta if ( const int c = compareCoord( approx_[b].y, approx_[pivot].y, [&] { return compareExact_( b, pivot, true ); } ) ) return std::abs( c ) == 1 ? c < 0 : c > 0; // below pivot, or level with it and ahead - // the points coincide before the perturbation: the direction of the perturbed points, by which less orders them, gives d + // the points coincide before the perturbation: the direction of the perturbed (or shifted) points, by which less orders them, gives d if ( const auto d = coincidentDirection_( pivot, b ) ) return d->y != 0 ? d->y < 0 : d->x > 0; - return false; // the failure: points defined alike + return false; // the failure: points defined alike, not told apart by their shifts } #if __GNUC__ >= 12 diff --git a/source/MRMesh/MRCutTrianglePredicates.h b/source/MRMesh/MRCutTrianglePredicates.h index 3f5105596b4f..7bf9b67c3067 100644 --- a/source/MRMesh/MRCutTrianglePredicates.h +++ b/source/MRMesh/MRCutTrianglePredicates.h @@ -25,6 +25,11 @@ struct MR_BIND_IGNORE PreciseCutPoint /// the triangle whose plane the segment crosses; not used for a vertex std::array tri; + /// optional for a segment that is an edge of another mesh crossing the triangle: the third vertices of the two triangles of that edge, + /// in any order, with invalid ids if not given (as for a boundary edge); they tell apart the crossings of the edges of a multiple edge, + /// which are defined alike otherwise (see \ref CutTrianglePredicates) + std::array wedge; + [[nodiscard]] bool isVertex() const { return !seg[1].id.valid(); } }; @@ -32,6 +37,11 @@ struct MR_BIND_IGNORE PreciseCutPoint /// line over the triangle; they never contradict the boolean's exact predicates, because every point gets the same /// simulation-of-simplicity perturbation by vertex ids as in \ref orient3d, and a point is on the side of a cut where /// the boolean's predicates put it; +/// the crossings of a segment with the triangle defined alike (the edges of a multiple edge of another mesh), which that perturbation +/// cannot tell apart and the boolean's predicates never compare, are told apart by the triangles of their edges: each edge is shifted by +/// an amount infinitely smaller than every term of the perturbation, toward the third vertices of its triangles (PreciseCutPoint::wedge), +/// or away from them if the triangles of another edge defined alike are in the angle between its own ones, so that the contours of two edges +/// whose sheets touch along the segment without crossing do not cross either, while those of crossing sheets do; /// the sweep goes over the projection of the perturbed triangle on the plane of the two axes least aligned with its normal /// (dropping the axis of the largest component of the normal, the first one of equal components, or the axis of the leading /// term of the perturbed normal if the exact one is zero), meeting the points in the order of one linear functional, @@ -56,7 +66,8 @@ class MR_BIND_IGNORE CutTrianglePredicates /// to the largest v + 1, so the points take small local ids, not those of the mesh vertices; a crossing is inside the triangle it crosses, /// as \ref doTriangleSegmentIntersect finds them, which only Debug checks, since it takes several orient3d per point; /// a segment in the plane of its triangle crosses it after the perturbation, at a point that tends to an end of the segment - /// or to a point inside it; setting a point again forgets the orders found for its old definition, but not a failure; + /// or to a point inside it; setting a point again forgets the orders found for its old definition, but not a failure, and so does + /// setting one with a wedge after points defined alike were compared, since it may change their shifts; /// a vertex point that is not a corner of the triangle, repeated ids, a segment with an end on a vertex of its own triangle, /// a segment off the sides with a triangle other than this one (as one with a single end at a corner), or a segment that does not /// cross the plane of its triangle are failures (see failed) @@ -67,7 +78,8 @@ class MR_BIND_IGNORE CutTrianglePredicates /// the signs of the projected coordinates of q - p found by expanding the perturbation alone, as the sweep sees them: /// a coordinate whose leading term has a larger degree than that of the other one is zero; this is slow, - /// while less and ccwFromBehind decide most pairs by faster rules; nullopt if the points coincide after the perturbation too + /// while less and ccwFromBehind decide most pairs by faster rules; nullopt if the points coincide after the perturbation too, + /// as points defined alike do (only the shifts of their edges tell them apart) [[nodiscard]] MRMESH_API std::optional directionByPerturbation( VertId p, VertId q ) const; /// whether triangle abc of distinct points is counterclockwise in the projection, for any three points: the first rule that applies @@ -83,23 +95,27 @@ class MR_BIND_IGNORE CutTrianglePredicates /// with a common vertex (for their crossings with this triangle), which ccw needs only for points on one line before the perturbation; /// the points on one side stay on one line after the perturbation, and they are taken /// as bulging inward a little, so three of them make a counterclockwise triangle if they come in the order - /// from vertex i+1 to vertex i of side i, up to a cyclic shift; + /// from vertex i+1 to vertex i of side i, up to a cyclic shift; where the line has the third point too, two of the points + /// are defined alike in valid input, and the shifts of their edges decide; /// nullopt if no two of the points are on a known line [[nodiscard]] MRMESH_API std::optional ccwByLines( VertId a, VertId b, VertId c ) const; /// the same as ccw found by expanding the perturbation alone, which is slow, while ccw decides most triangles by faster rules; - /// nullopt if the points are on one line after the perturbation too, as the points of one side are + /// nullopt if the points are on one line after the perturbation too, as the points of one side are, and two points defined alike [[nodiscard]] MRMESH_API std::optional ccwByPerturbation( VertId a, VertId b, VertId c ) const; /// whether the triangle of a point far behind pivot against the sweep, b and pivot is counterclockwise, for distinct b and pivot [[nodiscard]] MRMESH_API bool ccwFromBehind( VertId b, VertId pivot ) const; - /// whether the input broke a precondition that reset or addPoint checks, or a predicate met points it cannot tell apart, which valid - /// input never gives either: two points defined alike, or three points off the sides on one line also after the perturbation (as three - /// crossings of one cut, while valid input has two at most); it stays set till reset; the caller checks it after the addPoint calls + /// whether the input broke a precondition that reset or addPoint checks, or a predicate met points it cannot tell apart: two points + /// defined alike that the shifts of their edges do not tell apart (a vertex twice, crossings of a side, crossings without wedges or + /// with the same one, or with the segment parallel to the triangle before the perturbation, as with a zero normal), three points defined alike with + /// their shifts on one line, or three points off the sides on one line also after the perturbation (as three crossings of one cut, + /// while valid input has two at most, apart from points defined alike); it stays set till reset; the caller checks it after the addPoint calls /// and calls no predicate if it is set, since they take the points as valid (less may be cyclic then, and Debug asserts may fire); /// a failure that a predicate finds leaves the predicates callable, with arbitrary answers that the caller discards, and less a strict - /// total order (points defined alike compare alike with every other point, and by their ids with one another) + /// total order (points defined alike compare alike with every other point, and by their ids with one another if their shifts do not + /// tell them apart) [[nodiscard]] MRMESH_API bool failed() const; private: @@ -114,19 +130,25 @@ class MR_BIND_IGNORE CutTrianglePredicates /// defined in the .cpp file struct Impl; + /// the shift of a point defined alike with others; defined in the .cpp file + struct Shift; + int compareExact_( VertId p, VertId q, bool yFirst ) const; int exactOrient_( VertId a, VertId b, VertId c ) const; bool ccwExact_( VertId a, VertId b, VertId c ) const; int orientSlow_( VertId a, VertId b, VertId c ) const; std::optional ccwBySides_( VertId a, VertId b, VertId c ) const; std::optional ccwByPlanes_( VertId a, VertId b, VertId c ) const; + std::optional ccwByShifts_( VertId a, VertId b, VertId c ) const; std::optional coincidentDirection_( VertId p, VertId q ) const; std::optional directionFromSide_( VertId p, VertId q ) const; - std::optional perturbedDirection_( VertId p, VertId q ) const; + std::optional perturbedDirection_( VertId p, VertId q, bool alike ) const; + std::optional shiftDirection_( VertId p, VertId q ) const; + std::optional shift_( VertId p ) const; std::optional commonLine_( VertId p, VertId q, bool withVertex ) const; Vector2i lineDirection_( const Line & l ) const; bool along_( const Vector2i & dir, VertId p, VertId q ) const; - bool sideOfPlane_( const std::array & plane, VertId x ) const; + int sideOfPlane_( const std::array & plane, VertId x ) const; void fail_() const; std::array tri_; diff --git a/source/MRMesh/MREmbedTerrainStructure.cpp b/source/MRMesh/MREmbedTerrainStructure.cpp index 3912e73f9aa4..d275a67de78a 100644 --- a/source/MRMesh/MREmbedTerrainStructure.cpp +++ b/source/MRMesh/MREmbedTerrainStructure.cpp @@ -176,7 +176,9 @@ Expected TerrainEmbedder::createCutStructure_() VertBitSet cutVerts; if ( !structPrecutRes.contours.empty() ) { - auto structCutRes = cutMesh( structPrecutRes.mesh, structPrecutRes.contours ); + // the boolean's exact sorting data, which also connects the contours inside one face of the structure exactly + const auto sortData = structPrecutRes.sortData( result_ ); + auto structCutRes = cutMesh( structPrecutRes.mesh, structPrecutRes.contours, { .sortData = &sortData } ); if ( structCutRes.fbsWithContourIntersections.any() ) return unexpected( "Intersection contour of structure and terrain has self-intersections" ); cutVerts = getIncidentVerts( structPrecutRes.mesh.topology, @@ -228,135 +230,116 @@ Expected TerrainEmbedder::prepareTerrainCut auto offCont = offsetContour_( mc, cutOffset, fillOffset ); - for ( int loneResolveTries = 0; loneResolveTries < 5; ++loneResolveTries ) + std::vector mtps( offCont.contour.size() - 1 ); + tbb::task_group_context ctx; + bool canceled = false; + ParallelFor( mtps, [&] ( size_t i ) { - std::vector mtps( offCont.contour.size() - 1 ); - tbb::task_group_context ctx; - bool canceled = false; - ParallelFor( mtps, [&] ( size_t i ) + auto index = findOffsetContourIndex_( int( i ), offCont.idsShifts ); + const auto& startPoint = mc.contour[index]; + if ( mc.intBitSet.test( VertId( index ) ) ) { - auto index = findOffsetContourIndex_( int( i ), offCont.idsShifts ); - const auto& startPoint = mc.contour[index]; - if ( mc.intBitSet.test( VertId( index ) ) ) - { - mtps[i] = findProjection( startPoint, result_ ).mtp; - return; - } - auto line = - Line3f( startPoint, - to3dim( offCont.contour[i] ) + Vector3f( 0, 0, startPoint.z ) + - ( mc.cutBitSet.test( VertId( index ) ) ? Vector3f::plusZ() : Vector3f::minusZ() ) - startPoint ); - auto interRes = rayMeshIntersect( result_, line, -FLT_MIN ); // - FLT_MIN here to handle vertex lying in the same plane - if ( !interRes ) - { - if ( ctx.cancel_group_execution() ) - canceled = true; - return; - } - mtps[i] = interRes.mtp; - } ); + mtps[i] = findProjection( startPoint, result_ ).mtp; + return; + } + auto line = + Line3f( startPoint, + to3dim( offCont.contour[i] ) + Vector3f( 0, 0, startPoint.z ) + + ( mc.cutBitSet.test( VertId( index ) ) ? Vector3f::plusZ() : Vector3f::minusZ() ) - startPoint ); + auto interRes = rayMeshIntersect( result_, line, -FLT_MIN ); // - FLT_MIN here to handle vertex lying in the same plane + if ( !interRes ) + { + if ( ctx.cancel_group_execution() ) + canceled = true; + return; + } + mtps[i] = interRes.mtp; + } ); - if ( canceled ) - return unexpected( "Cannot embed structure beyond terrain" ); + if ( canceled ) + return unexpected( "Cannot embed structure beyond terrain" ); - filterDuplicates( mtps, offCont.idsShifts ); + filterDuplicates( mtps, offCont.idsShifts ); - Contour2f planarCont( mtps.size() + 1 ); - ParallelFor( mtps, [&] ( size_t i ) - { - planarCont[i] = to2dim( result_.triPoint( mtps[i] ) ); - } ); - planarCont.back() = planarCont.front(); + Contour2f planarCont( mtps.size() + 1 ); + ParallelFor( mtps, [&] ( size_t i ) + { + planarCont[i] = to2dim( result_.triPoint( mtps[i] ) ); + } ); + planarCont.back() = planarCont.front(); - auto filterBt = filterBowties( planarCont ); - std::vector> noBowtiesMtps( filterBt.initIndices.size() ); - for ( int i = 0; i < noBowtiesMtps.size(); ++i ) + auto filterBt = filterBowties( planarCont ); + std::vector> noBowtiesMtps( filterBt.initIndices.size() ); + for ( int i = 0; i < noBowtiesMtps.size(); ++i ) + { + const auto& initInds = filterBt.initIndices[i]; + const auto& noBTCont = filterBt.contours[i]; + auto& noBowtiesMtp = noBowtiesMtps[i]; + noBowtiesMtp.resize( initInds.size() ); + for ( int j = 0; j < initInds.size(); ++j ) { - const auto& initInds = filterBt.initIndices[i]; - const auto& noBTCont = filterBt.contours[i]; - auto& noBowtiesMtp = noBowtiesMtps[i]; - noBowtiesMtp.resize( initInds.size() ); - for ( int j = 0; j < initInds.size(); ++j ) + if ( initInds[j] != -1 ) + noBowtiesMtp[j] = mtps[initInds[j]]; + else { - if ( initInds[j] != -1 ) - noBowtiesMtp[j] = mtps[initInds[j]]; - else - { - auto line = Line3f( to3dim( noBTCont[j] ), Vector3f::plusZ() ); - auto interRes = rayMeshIntersect( result_, line, -FLT_MAX, FLT_MAX ); - if ( !interRes ) - return unexpected( "Cannot resolve bow ties on embedded structure wall" ); - noBowtiesMtp[j] = interRes.mtp; - } + auto line = Line3f( to3dim( noBTCont[j] ), Vector3f::plusZ() ); + auto interRes = rayMeshIntersect( result_, line, -FLT_MAX, FLT_MAX ); + if ( !interRes ) + return unexpected( "Cannot resolve bow ties on embedded structure wall" ); + noBowtiesMtp[j] = interRes.mtp; } } + } - MappedMeshContours res; - res.filtBowTiesMap = std::move( filterBt.initIndices ); - res.contours.resize( noBowtiesMtps.size() ); - res.map.resize( noBowtiesMtps.size() ); - OneMeshContours loneContours; - for ( int i = 0; i < res.contours.size(); ++i ) + MappedMeshContours res; + res.filtBowTiesMap = std::move( filterBt.initIndices ); + res.contours.resize( noBowtiesMtps.size() ); + res.map.resize( noBowtiesMtps.size() ); + for ( int i = 0; i < res.contours.size(); ++i ) + { + auto cont = noBowtiesMtps[i]; + cont.push_back( cont.front() ); + auto contourRes = surfaceLineToMeshContour( result_, cont, + [&] ( const MeshTriPoint& start, const MeshTriPoint& end, int startInd, int endInd )->Expected { - bool lone = true; - auto cont = noBowtiesMtps[i]; - cont.push_back( cont.front() ); - auto contourRes = surfaceLineToMeshContour( result_, cont, - [&] ( const MeshTriPoint& start, const MeshTriPoint& end, int startInd, int endInd )->Expected + auto initSMtpIndex = res.filtBowTiesMap[i][startInd]; + // the end of the last section is the closing point of cont, which repeats its first point + auto initEMtpIndex = res.filtBowTiesMap[i][endInd % res.filtBowTiesMap[i].size()]; + 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 ccwPath = trackSection( result_, start, end, planePoint, true ); + auto cwPath = trackSection( result_, start, end, planePoint, false ); + if ( ccwPath.has_value() && cwPath.has_value() ) { - 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 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() ) + auto ccwL = surfacePathLength( result_, *ccwPath ); + auto cwL = surfacePathLength( result_, *cwPath ); + if ( ccwL < cwL ) 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( contourRes.error() ); - res.contours[i] = std::move( *contourRes ); - for ( int j = 0; j < res.contours[i].intersections.size(); ++j ) + return cwPath; + } + else if ( ccwPath.has_value() ) + return ccwPath; + else if ( cwPath.has_value() ) + return cwPath; + else { - if ( std::holds_alternative( res.contours[i].intersections[j].primitiveId ) ) - { - lone = false; - break; - } + auto locRes = computeGeodesicPath( result_, start, end ); + if ( !locRes.has_value() ) + return unexpected( toString( locRes.error() ) ); + return *locRes; } - if ( !lone ) - continue; - loneContours.emplace_back( std::move( res.contours[i] ) ); - } - if ( loneContours.empty() ) - { - res.offsetMap = std::move( offCont.idsShifts ); - res.intBitSet = std::move( mc.intBitSet ); - res.cutBitSet = std::move( mc.cutBitSet ); - return res; - } - subdivideLoneContours( result_, loneContours ); + }, &res.map[i] ); + if ( !contourRes.has_value() ) + return unexpected( contourRes.error() ); + res.contours[i] = std::move( *contourRes ); } - return unexpected( "Cannot resolve lone cut on terrain" ); + // a contour inside one triangle of the terrain, crossing none of its edges, is connected with the rest of the triangle by cutMesh in cutTerrain + res.offsetMap = std::move( offCont.idsShifts ); + res.intBitSet = std::move( mc.intBitSet ); + res.cutBitSet = std::move( mc.cutBitSet ); + return res; } Expected> TerrainEmbedder::cutTerrain( const MappedMeshContours& mmc ) diff --git a/source/MRMesh/MRIntersectionContour.cpp b/source/MRMesh/MRIntersectionContour.cpp index c16dbd135053..9bd92fa727fb 100644 --- a/source/MRMesh/MRIntersectionContour.cpp +++ b/source/MRMesh/MRIntersectionContour.cpp @@ -7,6 +7,7 @@ #include "MRAffineXf3.h" #include "MRParallelFor.h" #include "MRphmap.h" +#include "MRPch/MRSuppressWarning.h" #include namespace MR @@ -353,7 +354,10 @@ void removeLoneDegeneratedContours( const MeshTopology& edgesTopology, OneMeshCo void removeLoneContours( ContinuousContours& contours, bool ignoreOpen ) { MR_TIMER; + MR_SUPPRESS_WARNING_PUSH + MR_SUPPRESS_WARNING( "-Wdeprecated-declarations", 4996 ) // deprecated together with this function auto loneContours = detectLoneContours( contours, ignoreOpen ); + MR_SUPPRESS_WARNING_POP for ( int i = int( loneContours.size() ) - 1; i >= 0; --i ) { contours.erase( contours.begin() + loneContours[i] ); diff --git a/source/MRMesh/MRIntersectionContour.h b/source/MRMesh/MRIntersectionContour.h index ca2fe2145d60..217a0f9a5b40 100644 --- a/source/MRMesh/MRIntersectionContour.h +++ b/source/MRMesh/MRIntersectionContour.h @@ -55,17 +55,20 @@ MRMESH_API bool isClosed( const ContinuousContour& contour ); /// Detects contours that fully lay inside one triangle /// if `ignoreOpen` then do not mark non-closed contours as lone, even if they really are /// returns they indices in contours +[[deprecated( "cutMesh connects lone contours itself, and the boolean no longer uses this" )]] MRMESH_API std::vector detectLoneContours( const ContinuousContours& contours, bool ignoreOpen = false ); /// Removes contours with zero area (do not remove if contour is handle on topology) /// edgesTopology - topology on which contours are represented with edges /// faceContours - lone contours represented by faces (all intersections are in same mesh A face) /// edgeContours - lone contours represented by edges (all intersections are in mesh B edges, edgesTopology: meshB.topology) +[[deprecated( "cutMesh connects lone contours itself, and the boolean no longer uses this" )]] MRMESH_API void removeLoneDegeneratedContours( const MeshTopology& edgesTopology, OneMeshContours& faceContours, OneMeshContours& edgeContours ); /// Removes contours that fully lay inside one triangle from the contours /// if `ignoreOpen` then do not consider non-closed contours as lone, even if they really are +[[deprecated( "cutMesh connects lone contours itself, and the boolean no longer uses this" )]] MRMESH_API void removeLoneContours( ContinuousContours& contours, bool ignoreOpen = false ); } diff --git a/source/MRMesh/MRMesh.h b/source/MRMesh/MRMesh.h index cb5d6720ab56..d90160143ad5 100644 --- a/source/MRMesh/MRMesh.h +++ b/source/MRMesh/MRMesh.h @@ -494,8 +494,8 @@ struct [[nodiscard]] Mesh /// updates existing caches after some faces were split into parts (e.g. by splitFace() or splitEdge()) with no other changes in the mesh; /// it shall be considered as a faster alternative to invalidateCaches() and following rebuild of AABB tree /// \param new2Old maps every face appeared after the splits to the face it was split from (as filled by splitFace() and splitEdge()); - /// it must contain only the splits made since the AABB tree was built or last updated: a map accumulated over several calls - /// (like new2oldMap of subdivideLoneContours in the iterations of boolean) would add the earlier faces to the tree once more + /// it must contain only the splits made since the AABB tree was built or last updated: a map accumulated over several updates + /// would add the earlier faces to the tree once more MRMESH_API void updateCachesAfterSplits( const FaceHashMap & new2Old ); // returns the amount of memory this object occupies on heap diff --git a/source/MRMesh/MRMeshBoolean.cpp b/source/MRMesh/MRMeshBoolean.cpp index a1b87351a5b7..00b8465daf75 100644 --- a/source/MRMesh/MRMeshBoolean.cpp +++ b/source/MRMesh/MRMeshBoolean.cpp @@ -76,6 +76,12 @@ void gatherEdgeInfo( const MeshTopology& topology, EdgeId e, FaceBitSet& faces, dests.set( topology.dest( e ) ); } +// removes the closed contours of two points, which have zero area, as their two segments coincide +void removeTwoPointLoops( ContinuousContours& contours ) +{ + std::erase_if( contours, [] ( const ContinuousContour& c ) { return c.size() == 3 && isClosed( c ); } ); +} + } //anonymous namespace void convertIntFloatAllVerts( Mesh & mesh, const CoordinateConverters& conv ) @@ -168,75 +174,6 @@ std::vector findSelfContoursMapping( const ContinuousContours& } return holePairs; } - -enum LoneProccessingState -{ - None, - DegenerateOnly, - Processed -}; - -LoneProccessingState subdivideSelfLone( Mesh& mesh, const CoordinateConverters& converters, const ContinuousContours& contours ) -{ - auto loneContoursIds = detectLoneContours( contours, true ); - if ( loneContoursIds.empty() ) - return LoneProccessingState::None; // none left, go on - - auto loneHolePairs = detail::findSelfContoursMapping( contours ); - assert( loneContoursIds.size() % 2 == 0 ); - ContinuousContours fLone( loneContoursIds.size() / 2 ), joinedELone( loneContoursIds.size() / 2 ); - BitSet visitedLone( loneContoursIds.size(), true ); - for ( int i = 0; i < fLone.size(); ++i ) - { - int loneIndex = -1; - int otherLoneIndex = -1; - for ( auto l : visitedLone ) - { - if ( loneIndex == -1 ) - { - loneIndex = loneContoursIds[l]; - visitedLone.reset( l ); - } - else - { - for ( auto [f, s, _] : loneHolePairs ) - { - if ( f == loneIndex ) - { - otherLoneIndex = s; - break; - } - else if ( s == loneIndex ) - { - otherLoneIndex = f; - break; - } - } - } - if ( otherLoneIndex != -1 ) - { - auto loIt = std::find( loneContoursIds.begin(), loneContoursIds.end(), otherLoneIndex ); - assert( loIt != loneContoursIds.end() ); - visitedLone.reset( size_t( std::distance( loneContoursIds.begin(), loIt ) ) ); - break; - } - } - assert( loneIndex != -1 ); - assert( otherLoneIndex != -1 ); - if ( contours[loneIndex].front().isEdgeATriB() ) - std::swap( loneIndex, otherLoneIndex ); - fLone[i] = contours[loneIndex]; - joinedELone[i] = contours[otherLoneIndex]; - } - auto loneFInts = getOneMeshSelfIntersectionContours( mesh, fLone, converters ); - auto joinedLoneEInts = getOneMeshSelfIntersectionContours( mesh, joinedELone, converters ); - - removeLoneDegeneratedContours( mesh.topology, loneFInts, joinedLoneEInts ); - if ( loneFInts.empty() ) - return LoneProccessingState::DegenerateOnly; - subdivideLoneContours( mesh, loneFInts ); - return LoneProccessingState::Processed; -} } BooleanResult forceBoolean( const Mesh& meshA, const Mesh& meshB, BooleanOperation operation, const BooleanParameters& inParams /*= {} */ ) @@ -255,35 +192,13 @@ Expected selfBoolean( const Mesh& inMesh ) { auto box = Box3d( inMesh.computeBoundingBox() ); CoordinateConverters converters{ .toInt = getToIntConverter( box ),.toFloat = getToFloatConverter( box ) }; - std::vector intersections; - ContinuousContours contours; - Mesh mesh = inMesh; - int iters = 0; - const int cMaxFixLoneIterations = 3; // lets not do it many times - for ( ;; ++iters ) - { - // find intersections - intersections = findSelfCollidingEdgeTrisPrecise( mesh, converters.toInt ); - // order intersections - contours = orderSelfIntersectionContours( mesh.topology, intersections ); - - // find lone - auto loneProcessing = detail::subdivideSelfLone( mesh, converters, contours ); - if ( loneProcessing == detail::LoneProccessingState::None ) - break; - else if ( loneProcessing == detail::LoneProccessingState::DegenerateOnly || iters == cMaxFixLoneIterations ) - { - // in some rare cases there are lone contours with zero area that cannot be resolved - // they lead to infinite loop, so just try to remove them - removeLoneContours( contours, true ); - break; - } - } + // find and order intersections + auto contours = orderSelfIntersectionContours( mesh.topology, findSelfCollidingEdgeTrisPrecise( mesh, converters.toInt ) ); + removeTwoPointLoops( contours ); // both contours of a pair together, as they have the same points auto holePairs = detail::findSelfContoursMapping( contours ); - auto meshCpy = mesh; // for now lets do copy each time, TODO: do copy only if needed auto intContours = getOneMeshSelfIntersectionContours( mesh, contours, converters ); // update non-closed @@ -340,7 +255,7 @@ Expected selfBoolean( const Mesh& inMesh ) } ); // meshAVertsNum = 0 because both meshes are actually the same and no need to change the IDs of any mesh vertices - auto sortData = std::make_unique( SortIntersectionsData{ meshCpy, contours, converters.toInt, nullptr, 0, false } ); + auto sortData = std::make_unique( SortIntersectionsData{ inMesh, contours, converters.toInt, nullptr, 0, false } ); auto cutRes = cutMesh( mesh, intContours, { .sortData = sortData.get() } ); if ( cutRes.fbsWithContourIntersections.any() ) return unexpected( fmt::format( "Self-Boolean has nested intersections on {} faces", cutRes.fbsWithContourIntersections.count() ) ); @@ -400,97 +315,28 @@ BooleanResult booleanImpl( Mesh&& meshA, Mesh&& meshB, BooleanOperation operatio MR_TIMER; BooleanResult result; CoordinateConverters converters; - PreciseCollisionResult intersections; ContinuousContours contours; bool needCutMeshA = operation != BooleanOperation::InsideB && operation != BooleanOperation::OutsideB; bool needCutMeshB = operation != BooleanOperation::InsideA && operation != BooleanOperation::OutsideA; converters = getVectorConverters( meshA, meshB, params.rigidB2A ); + const size_t meshAVertsNum = meshA.topology.vertSize(); // before any cut of A - auto [loneCb, mainCb] = splitProgress( params.cb, 0.8f ); - - FaceHashMap new2orgSubdivideMapA; - FaceHashMap new2orgSubdivideMapB; - std::vector prevLoneContoursIds; - int iters = 0; - const int cMaxFixLoneIterations = 100; - bool aSubdivided = false; - bool bSubdivided = false; - // find intersections - intersections = findCollidingEdgeTrisPrecise( meshA, meshB, converters.toInt, params.rigidB2A ); - for ( ;; iters++ ) - { - // order intersections - contours = orderIntersectionContours( meshA.topology, meshB.topology, intersections ); - // find lone - auto loneContoursIds = detectLoneContours( contours ); - - if ( loneCb && !loneCb( ( std::log10( float( iters + 1 ) * 0.1f ) + 2.0f ) / 3.0f ) ) - return { .errorString = stringOperationCanceled() }; + // find and order intersections + contours = orderIntersectionContours( meshA.topology, meshB.topology, + findCollidingEdgeTrisPrecise( meshA, meshB, converters.toInt, params.rigidB2A ) ); + removeTwoPointLoops( contours ); - if ( !loneContoursIds.empty() && ( loneContoursIds == prevLoneContoursIds || iters == cMaxFixLoneIterations ) ) - { - // in some rare cases there are lone contours with zero area that cannot be resolved - // they lead to infinite loop, so just try to remove them - removeLoneContours( contours ); - break; - } - - prevLoneContoursIds = loneContoursIds; - - // separate A lone from B lone - ContinuousContours loneA; - ContinuousContours loneB; - for ( int i = 0; i < loneContoursIds.size(); ++i ) - { - const auto& contour = contours[loneContoursIds[i]]; - if ( contour[0].isEdgeATriB() ) - loneB.push_back( contour ); - else - loneA.push_back( contour ); - } - if ( loneContoursIds.empty() || - ( loneA.empty() && !needCutMeshB ) || - ( loneB.empty() && !needCutMeshA ) ) - break; - // subdivide owners of lone - const VertId aFirstNewVert( meshA.topology.vertSize() ); - const VertId bFirstNewVert( meshB.topology.vertSize() ); - if ( !loneA.empty() && needCutMeshA ) - { - aSubdivided = true; - OneMeshContours loneIntsA, loneIntsAonB; - getOneMeshIntersectionContours( meshA, meshB, loneA, &loneIntsA, &loneIntsAonB, converters, params.rigidB2A ); - removeLoneDegeneratedContours( meshB.topology, loneIntsA, loneIntsAonB ); - subdivideLoneContours( meshA, loneIntsA, &new2orgSubdivideMapA ); - } - if ( !loneB.empty() && needCutMeshB ) - { - bSubdivided = true; - OneMeshContours loneIntsB, loneIntsBonA; - getOneMeshIntersectionContours( meshA, meshB, loneB, &loneIntsBonA, &loneIntsB, converters, params.rigidB2A ); - removeLoneDegeneratedContours( meshA.topology, loneIntsB, loneIntsBonA ); - subdivideLoneContours( meshB, loneIntsB, &new2orgSubdivideMapB ); - } - // update intersections: only the ones near new vertices could change - updateCollidingEdgeTrisPrecise( intersections, meshA, aFirstNewVert, meshB, bFirstNewVert, converters.toInt, params.rigidB2A ); - } - if ( iters == cMaxFixLoneIterations ) - { - spdlog::warn( "Boolean: fix lone contours iteration limit reached." ); - result.errorString = "Fix lone contours iteration limit reached."; - return result; - } - // clear intersections - intersections = {}; - - if ( mainCb && !mainCb( 0.0f ) ) + if ( params.cb && !params.cb( 0.0f ) ) return { .errorString = stringOperationCanceled() }; std::vector cutA, cutB; OneMeshContours meshAContours; OneMeshContours meshBContours; + // the meshes cut here, which need the sorting data: a precut output takes its mesh uncut + const bool doCutA = needCutMeshA && !params.outPreCutA; + const bool doCutB = needCutMeshB && !params.outPreCutB; // original mesh is needed to properly sort edges intersections while cutting meshes // as far as meshes is going to be change in the same time we need to take copies before start cutting Mesh meshACopyBuffer; // (only used if meshes were not copied before, otherwise old copy will be used) @@ -499,47 +345,41 @@ BooleanResult booleanImpl( Mesh&& meshA, Mesh&& meshB, BooleanOperation operatio std::unique_ptr dataForB; tbb::task_group taskGroup; - if ( needCutMeshA ) + if ( doCutA ) { taskGroup.run( [&] () { - if ( needCutMeshB ) + if ( doCutB ) { - if ( bSubdivided || !intParams.originalMeshB ) + if ( !intParams.originalMeshB ) { Timer t( "meshBCopyBuffer" ); meshBCopyBuffer = meshB; - if ( !intParams.originalMeshB ) - intParams.originalMeshB = &meshBCopyBuffer; + intParams.originalMeshB = &meshBCopyBuffer; } - assert( intParams.originalMeshB ); - const Mesh* sortMeshPtr = bSubdivided ? &meshBCopyBuffer : intParams.originalMeshB; - dataForA = std::make_unique( SortIntersectionsData{ *sortMeshPtr, contours, converters.toInt, params.rigidB2A, meshA.topology.vertSize(), false } ); + dataForA = std::make_unique( SortIntersectionsData{ *intParams.originalMeshB, contours, converters.toInt, params.rigidB2A, meshAVertsNum, false } ); } else // B is stable so no need copy - dataForA = std::make_unique( SortIntersectionsData{ meshB, contours, converters.toInt, params.rigidB2A, meshA.topology.vertSize(), false } ); + dataForA = std::make_unique( SortIntersectionsData{ meshB, contours, converters.toInt, params.rigidB2A, meshAVertsNum, false } ); } ); } - if ( needCutMeshB ) + if ( doCutB ) { - if ( needCutMeshA ) + if ( doCutA ) { - if ( aSubdivided || !intParams.originalMeshA ) + if ( !intParams.originalMeshA ) { Timer t( "meshACopyBuffer" ); meshACopyBuffer = meshA; - if ( !intParams.originalMeshA ) - intParams.originalMeshA = &meshACopyBuffer; + intParams.originalMeshA = &meshACopyBuffer; } - assert( intParams.originalMeshA ); - const Mesh* sortMeshPtr = aSubdivided ? &meshACopyBuffer : intParams.originalMeshA; - dataForB = std::make_unique( SortIntersectionsData{ *sortMeshPtr, contours, converters.toInt, params.rigidB2A, meshA.topology.vertSize(), true } ); + dataForB = std::make_unique( SortIntersectionsData{ *intParams.originalMeshA, contours, converters.toInt, params.rigidB2A, meshAVertsNum, true } ); } else // A is stable so no need to copy - dataForB = std::make_unique( SortIntersectionsData{ meshA, contours, converters.toInt, params.rigidB2A, meshA.topology.vertSize(), true } ); + dataForB = std::make_unique( SortIntersectionsData{ meshA, contours, converters.toInt, params.rigidB2A, meshAVertsNum, true } ); } taskGroup.wait(); @@ -549,18 +389,12 @@ BooleanResult booleanImpl( Mesh&& meshA, Mesh&& meshB, BooleanOperation operatio needCutMeshB ? &meshBContours : nullptr, converters, params.rigidB2A ); - if ( mainCb && !mainCb( 0.33f ) ) + if ( params.cb && !params.cb( 0.33f ) ) return { .errorString = stringOperationCanceled() }; - if ( params.outPreCutA ) - { - params.outPreCutA->contours = std::move( meshAContours ); - params.outPreCutA->mesh = std::move( meshA ); - } - bool needInsideA = operation == BooleanOperation::InsideA || operation == BooleanOperation::Intersection || operation == BooleanOperation::DifferenceBA; bool needInsideB = operation == BooleanOperation::InsideB || operation == BooleanOperation::Intersection || operation == BooleanOperation::DifferenceAB; - if ( needCutMeshA && !params.outPreCutA ) + if ( doCutA ) { taskGroup.run( [&] () { @@ -577,29 +411,12 @@ BooleanResult booleanImpl( Mesh&& meshA, Mesh&& meshB, BooleanOperation operatio auto res = cutMesh( meshA, meshAContours, cmParams ); meshAContours.clear(); meshAContours.shrink_to_fit(); // free memory - if ( cut2oldAPtr && !new2orgSubdivideMapA.empty() ) - { - ParallelFor( *cut2oldAPtr, [&] ( FaceId i ) - { - if ( !( *cut2oldAPtr )[i] ) - return; - FaceId refFace = ( *cut2oldAPtr )[i]; - auto it = new2orgSubdivideMapA.find( refFace ); - if ( it != new2orgSubdivideMapA.end() && it->second ) - ( *cut2oldAPtr )[i] = it->second; - } ); - } result.meshABadContourFaces = std::move( res.fbsWithContourIntersections ); cutA = std::move( res.resultCut ); } ); } - if ( params.outPreCutB ) - { - params.outPreCutB->contours = std::move( meshBContours ); - params.outPreCutB->mesh = std::move( meshB ); - } - if ( needCutMeshB && !params.outPreCutB ) + if ( doCutB ) { Timer t( "CutMeshB" ); FaceMap* cut2oldBPtr = params.mapper ? ¶ms.mapper->maps[int( BooleanResultMapper::MapObject::B )].cut2origin : nullptr; @@ -614,23 +431,28 @@ BooleanResult booleanImpl( Mesh&& meshA, Mesh&& meshB, BooleanOperation operatio auto res = cutMesh( meshB, meshBContours, cmParams ); meshBContours.clear(); meshBContours.shrink_to_fit(); // free memory - if ( cut2oldBPtr && !new2orgSubdivideMapB.empty() ) - { - ParallelFor( *cut2oldBPtr, [&] ( FaceId i ) - { - if ( !( *cut2oldBPtr )[i] ) - return; - FaceId refFace = ( *cut2oldBPtr )[i]; - auto it = new2orgSubdivideMapB.find( refFace ); - if ( it != new2orgSubdivideMapB.end() && it->second ) - ( *cut2oldBPtr )[i] = it->second; - } ); - } result.meshBBadContourFaces = std::move( res.fbsWithContourIntersections ); cutB = std::move( res.resultCut ); } taskGroup.wait(); + // the precut outputs take the meshes only after the cuts: where one mesh is not cut (by the operation or for its precut output), + // the cut of the other one sorts its intersections by that mesh itself, not by a copy; + // nothing reads the intersections after that, so the last output takes them + auto fillPreCut = [&] ( BooleanPreCutResult& out, Mesh& mesh, OneMeshContours& meshContours, bool isMeshA, bool last ) + { + out.mesh = std::move( mesh ); + out.contours = std::move( meshContours ); + out.intersections = last ? std::move( contours ) : contours; + out.toInt = converters.toInt; + out.rigidB2A = params.rigidB2A ? std::optional( *params.rigidB2A ) : std::nullopt; + out.meshAVertsNum = meshAVertsNum; + out.isMeshA = isMeshA; + }; + if ( params.outPreCutA ) + fillPreCut( *params.outPreCutA, meshA, meshAContours, true, !params.outPreCutB ); + if ( params.outPreCutB ) + fillPreCut( *params.outPreCutB, meshB, meshBContours, false, true ); if ( !params.forceCut && result.meshABadContourFaces.any() ) { @@ -646,7 +468,7 @@ BooleanResult booleanImpl( Mesh&& meshA, Mesh&& meshB, BooleanOperation operatio } - if ( mainCb && !mainCb( 0.66f ) ) + if ( params.cb && !params.cb( 0.66f ) ) return { .errorString = stringOperationCanceled() }; if ( params.outPreCutA || params.outPreCutB ) @@ -659,7 +481,7 @@ BooleanResult booleanImpl( Mesh&& meshA, Mesh&& meshB, BooleanOperation operatio auto res = doBooleanOperation( std::move( meshA ), std::move( meshB ), std::move( cutA ), std::move( cutB ), operation, params.rigidB2A, params.mapper, params.mergeAllNonIntersectingComponents, intParams ); - if ( mainCb && !mainCb( 1.0f ) ) + if ( params.cb && !params.cb( 1.0f ) ) return { .errorString = stringOperationCanceled() }; if ( res.has_value() ) diff --git a/source/MRMesh/MRMeshBoolean.h b/source/MRMesh/MRMeshBoolean.h index 3385d0ede13d..85b915f58b83 100644 --- a/source/MRMesh/MRMeshBoolean.h +++ b/source/MRMesh/MRMeshBoolean.h @@ -5,6 +5,7 @@ #include "MRMesh.h" #include "MRBitSet.h" #include "MRExpected.h" +#include #include namespace MR @@ -64,10 +65,33 @@ MR_BIND_OVERLOAD_SUFFIX(withXf) MRMESH_API BooleanResult boolean( const Mesh& me // no version of the function with rvalue-referenced meshes, because meshes' copies are currently used for sorting intersections +/// a mesh of the boolean and its contours of intersection with the other mesh before the cut, with the boolean's data for an exact cut; +/// the contours can include lone ones (inside one face, crossing none of its edges), which cutMesh connects exactly only given +/// the boolean's sorting data, see sortData struct BooleanPreCutResult { Mesh mesh; OneMeshContours contours; + + /// the boolean's ordered intersections of both meshes, without the closed contours of two points + ContinuousContours intersections; + + /// the boolean's conversion of the coordinates of both meshes to integers + ConvertToIntVector toInt; + + /// a copy of the boolean's transformation from mesh B to mesh A + std::optional rigidB2A; + + /// the number of vertices of mesh A before any cut + size_t meshAVertsNum = 0; + + /// whether mesh is mesh A of the boolean + bool isMeshA = false; + + /// the data to pass cutMesh( mesh, contours ) in CutMeshParameters::sortData for the boolean's exact cut, where otherMesh is the other mesh + /// of the boolean as it was passed to it; the result points into this struct and otherMesh, so both must outlive its use + [[nodiscard]] SortIntersectionsData sortData( const Mesh& otherMesh ) const + { return { otherMesh, intersections, toInt, rigidB2A ? &*rigidB2A : nullptr, meshAVertsNum, !isMeshA }; } }; /** \struct MR::BooleanResult diff --git a/source/MRMesh/MRMeshCollidePrecise.h b/source/MRMesh/MRMeshCollidePrecise.h index faf6b1d86360..3890c34bbb6b 100644 --- a/source/MRMesh/MRMeshCollidePrecise.h +++ b/source/MRMesh/MRMeshCollidePrecise.h @@ -81,6 +81,7 @@ MR_BIND_OVERLOAD_SUFFIX() MRMESH_API PreciseCollisionResult findCollidingEdgeTri * \param res the result of findCollidingEdgeTrisPrecise with the same conv and rigidB2A before the splits * \param aFirstNewVert (bFirstNewVert) all vertices of mesh A (B) with this or larger ids appeared in the splits */ +[[deprecated( "cutMesh connects lone contours itself, and the boolean no longer uses this" )]] MRMESH_API void updateCollidingEdgeTrisPrecise( PreciseCollisionResult & res, const Mesh & a, VertId aFirstNewVert, const Mesh & b, VertId bFirstNewVert, ConvertToIntVector conv, const AffineXf3f* rigidB2A = nullptr ); diff --git a/source/MRMesh/MROneMeshContours.h b/source/MRMesh/MROneMeshContours.h index 706d33b15f7e..a069a9f8495c 100644 --- a/source/MRMesh/MROneMeshContours.h +++ b/source/MRMesh/MROneMeshContours.h @@ -46,6 +46,7 @@ using OneMeshContours = std::vector; // Divides faces that fully own contours into 3 parts with center in center mass of one of the face contours // if there is more than one contour on face it guarantee to subdivide at least one lone contour on this face // the AABB tree of the mesh (if any) is updated instead of being invalidated +[[deprecated( "cutMesh connects lone contours itself, and the boolean no longer uses this" )]] MRMESH_API void subdivideLoneContours( Mesh& mesh, const OneMeshContours& contours, FaceHashMap* new2oldMap = nullptr ); /// Converts contours given in topological terms as the intersections of one mesh's edge and another mesh's triangle (ContinuousContours), diff --git a/source/MRTest/MR2DContoursTriangulationTests.cpp b/source/MRTest/MR2DContoursTriangulationTests.cpp index 8d89076f2641..8015b8f3a5f1 100644 --- a/source/MRTest/MR2DContoursTriangulationTests.cpp +++ b/source/MRTest/MR2DContoursTriangulationTests.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include @@ -8,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -412,6 +414,99 @@ TEST( MRMesh, PlanarTriangulationMonotonePlan ) EXPECT_EQ( parts.size(), size_t( 2 ) ); } +// on the exact predicates of a cut triangle, two vertices of the loops defined alike (as where both edges of a multiple edge +// of the other mesh cross the triangle) are told apart only by different wedges (the third vertices of the triangles of their edges), +// and an invalid definition makes addPoint fail, so without such wedges, or with an invalid point, the plan is rejected without running the sweep +TEST( MRMesh, PlanarTriangulationMonotonePlanSameDefinitions ) +{ + // the hole of triangle 0_f with island 1_f inside it + VertCoords pts{ Vector3f( 0, 0, 0 ), Vector3f( 1000, 0, 0 ), Vector3f( 0, 1000, 0 ), + Vector3f( 200, 200, 0 ), Vector3f( 400, 200, 0 ), Vector3f( 200, 400, 0 ) }; + Mesh mesh = Mesh::fromTriangles( std::move( pts ), { { 0_v, 1_v, 2_v }, { 3_v, 4_v, 5_v } } ); + const auto outer = mesh.topology.edgeWithLeft( 0_f ); + mesh.topology.setLeft( outer, {} ); + auto ringOf = [&] ( EdgeId e ) + { + EdgeLoop res; + for ( auto le : leftRing( mesh.topology, e ) ) + res.push_back( le ); + return res; + }; + const EdgeLoops loops{ ringOf( outer ), ringOf( mesh.topology.edgeWithLeft( 1_f ).sym() ) }; + + // the corners of the triangle are vertex points, the island vertices are the crossings of vertical segments with the triangle, + // the segment of vertex v having the ends of ids 2v and 2v + 1 (larger than those of the triangle) + const std::array tri{ PreciseVertCoords{ 0_v, { 0, 0, 0 } }, { 1_v, { 1000, 0, 0 } }, { 2_v, { 0, 1000, 0 } } }; + auto pointOf = [&] ( VertId v, VertId sameAs ) + { + if ( v < 3 ) + return PreciseCutPoint{ .seg = { tri[v] } }; + const auto u = v == 4_v && sameAs ? sameAs : v; + const auto p = Vector3i( mesh.points[u] ); + return PreciseCutPoint{ .seg = { PreciseVertCoords{ VertId( 2 * u ), p - Vector3i( 0, 0, 100 ) }, { VertId( 2 * u + 1 ), p + Vector3i( 0, 0, 100 ) } }, + .tri = tri }; + }; + + auto plan = PlanarTriangulation::getMonotonePlan( mesh, loops, tri, [&] ( VertId v ) { return pointOf( v, {} ); } ); + ASSERT_TRUE( plan.has_value() ); + EXPECT_EQ( plan->items.size(), size_t( 2 ) ); // where the sweep reaches the island and where it leaves it + + // vertex 4 defined as vertex 3, also with the ends of the segment and the vertices of the triangle in another order + EXPECT_FALSE( PlanarTriangulation::getMonotonePlan( mesh, loops, tri, [&] ( VertId v ) { return pointOf( v, 3_v ); } ).has_value() ); + auto reordered = [&] ( VertId v ) + { + auto res = pointOf( v, 3_v ); + if ( v == 4_v ) + { + std::swap( res.seg[0], res.seg[1] ); + std::rotate( res.tri.begin(), res.tri.begin() + 1, res.tri.end() ); + } + return res; + }; + EXPECT_FALSE( PlanarTriangulation::getMonotonePlan( mesh, loops, tri, reordered ).has_value() ); + + // with wedges, of the triangles of two edges between the ends of the segment of vertex 3 at ( 200, 200 ), one to the left and one + // to the right of the segment: the shifts toward them keep vertex 4 to the right of vertex 3, so the island keeps its orientation, + // and the plan is the same as with vertex 4 at its own place; with the same wedge, in either order, or without the wedge of one of them, + // the points stay defined alike + const Vector3i p3( mesh.points[3_v] ); + auto wedgeVert = [&] ( int id, int dx, int dy, int dz ) { return PreciseVertCoords{ VertId( id ), p3 + Vector3i( dx, dy, dz ) }; }; + const std::array left{ wedgeVert( 20, -10, 5, 50 ), wedgeVert( 21, -10, -5, -50 ) }; + const std::array right{ wedgeVert( 22, 10, 5, 50 ), wedgeVert( 23, 10, -5, -50 ) }; + auto withWedges = [&] ( const std::array& w3, const std::array& w4 ) + { + return [&reordered, w3, w4] ( VertId v ) + { + auto res = reordered( v ); + if ( v == 3_v ) + res.wedge = w3; + else if ( v == 4_v ) + res.wedge = w4; + return res; + }; + }; + plan = PlanarTriangulation::getMonotonePlan( mesh, loops, tri, withWedges( left, right ) ); + ASSERT_TRUE( plan.has_value() ); + EXPECT_EQ( plan->items.size(), size_t( 2 ) ); + EXPECT_FALSE( PlanarTriangulation::getMonotonePlan( mesh, loops, tri, withWedges( left, left ) ).has_value() ); + EXPECT_FALSE( PlanarTriangulation::getMonotonePlan( mesh, loops, tri, withWedges( left, { left[1], left[0] } ) ).has_value() ); + EXPECT_FALSE( PlanarTriangulation::getMonotonePlan( mesh, loops, tri, withWedges( left, {} ) ).has_value() ); + + // vertex 4 as a vertex point off the corners of the triangle, or as a segment not reaching the plane of the triangle + for ( bool asVertex : { true, false } ) + { + EXPECT_FALSE( PlanarTriangulation::getMonotonePlan( mesh, loops, tri, [&] ( VertId v ) + { + auto res = pointOf( v, {} ); + if ( v == 4_v && asVertex ) + res = PreciseCutPoint{ .seg = { res.seg[0] } }; + else if ( v == 4_v ) + res.seg[0].pt.z += 150; + return res; + } ).has_value() ); + } +} + namespace { diff --git a/source/MRTest/MRContoursCutTests.cpp b/source/MRTest/MRContoursCutTests.cpp index 44a02b7c337f..dff4a78d664a 100644 --- a/source/MRTest/MRContoursCutTests.cpp +++ b/source/MRTest/MRContoursCutTests.cpp @@ -3,10 +3,16 @@ #include #include #include +#include +#include #include #include #include +#include +#include +#include #include +#include namespace MR { @@ -432,6 +438,292 @@ TEST( MRMesh, CutMeshContourTurnsInsideFace ) } } +// a contour inside one triangle (a lone contour) is connected with the rest of the triangle by the cut, near the origin and far from it, +// also where the contour fills most of the triangle; a closed contour is an island in the hole of the triangle, an open one a slit +TEST( MRMesh, CutMeshLoneContour ) +{ + // uniform values in [-1, 1] straight from the generator, whose output the standard specifies (unlike that of the distributions), + // so they are the same on every platform; the points made of them with std::cos and std::sin can differ in their last bits there, + // which the tolerances absorb + std::mt19937 gen( 1 ); + auto random = [&] { return 2.0 * gen() / std::mt19937::max() - 1; }; + auto randomVec = [&] + { + Vector3d res; + for ( int i = 0; i < 3; ++i ) + res[i] = random(); + return res; + }; + for ( double dist : { 0.0, 1e5 } ) + { + for ( float fill : { 0.5f, 0.9f } ) + { + for ( int trial = 0; trial < 50; ++trial ) + { + // an equilateral triangle in a random plane, so that every axis is dropped by the projection in some trials + const auto base = dist * randomVec(); + const auto u = randomVec().normalized(); + const auto v = cross( u, randomVec() ).normalized(); + VertCoords points; + for ( int i = 0; i < 3; ++i ) + points.push_back( Vector3f( base + 10.0 * ( std::cos( 2 * PI * i / 3 ) * u + std::sin( 2 * PI * i / 3 ) * v ) ) ); + const auto mesh = Mesh::fromTriangles( std::move( points ), { { 0_v, 1_v, 2_v } } ); + + // the triangle shrunk around its centroid to the given part of its area, with four points on every side + const auto e = mesh.topology.edgeWithLeft( 0_f ); + const float k = std::sqrt( fill ), m = ( 1 - k ) / 3; + const std::array corners{ Vector2f( m, m ), Vector2f( m + k, m ), Vector2f( m, m + k ) }; + std::vector mtps; + for ( int i = 0; i < 3; ++i ) + for ( int j = 0; j < 4; ++j ) + { + const auto b = corners[i] + ( corners[( i + 1 ) % 3] - corners[i] ) * ( j / 4.0f ); + mtps.push_back( MeshTriPoint( e, TriPointf( b.x, b.y ) ) ); + } + auto closed = mtps; + closed.push_back( mtps.front() ); + // the closed contour, the same without its last segment, and a slit of one segment (a hole loop of 2 edges) + for ( const auto& line : { closed, mtps, std::vector{ mtps[0], mtps[6] } } ) + { + auto contour = surfaceLineToMeshContour( mesh, line ); + ASSERT_TRUE( contour.has_value() ); + + Mesh cut = mesh; + NewEdgesMap new2OldEdges; + auto res = cutMesh( cut, { *contour }, { .new2oldEdgesMap = &new2OldEdges } ); + EXPECT_TRUE( cut.topology.checkValidity() ); + EXPECT_TRUE( res.fbsWithContourIntersections.none() ) + << "distance " << dist << ", fill " << fill << ", trial " << trial << ", points " << line.size(); + EXPECT_EQ( cut.topology.findNumHoles(), 1 ); // the island or the slit is connected and filled on both sides + EXPECT_NEAR( cut.area(), mesh.area(), 1e-3 * mesh.area() ); + // every new edge, as a chord connecting the contour with the sides, maps to the triangle (surface offset relies on it) + for ( auto ue = UndirectedEdgeId( mesh.topology.undirectedEdgeSize() ); ue < cut.topology.undirectedEdgeSize(); ++ue ) + { + const auto it = new2OldEdges.map.find( ue ); + ASSERT_NE( it, new2OldEdges.map.end() ); + EXPECT_FALSE( new2OldEdges.splitEdges.test( ue ) ); + EXPECT_EQ( it->second, 0 ); // 0_f + } + // the cut goes along the contour + ASSERT_EQ( res.resultCut.size(), 1 ); + const auto& path = res.resultCut[0]; + ASSERT_EQ( path.size() + 1, contour->intersections.size() ); + EXPECT_TRUE( isEdgePath( cut.topology, path ) ); + for ( int i = 0; i < path.size(); ++i ) + EXPECT_EQ( cut.orgPnt( path[i] ), contour->intersections[i].coordinate ); + } + } + } + } +} + +// a lone contour crossing itself (a bowtie inside one triangle of a closed mesh) cannot be connected with the triangle; +// the cut reports the triangle as bad and keeps it, so that its sides make no hole whichever faces are filled +TEST( MRMesh, CutMeshSelfCrossingLoneContour ) +{ + const auto mesh = makeCube(); + const auto f = 0_f; + const auto e = mesh.topology.edgeWithLeft( f ); + // side 0-1 crosses side 2-3 + std::vector mtps; + for ( const auto& b : { Vector2f( 0.2f, 0.2f ), Vector2f( 0.5f, 0.4f ), Vector2f( 0.5f, 0.2f ), Vector2f( 0.2f, 0.4f ), Vector2f( 0.2f, 0.2f ) } ) + mtps.push_back( MeshTriPoint( e, TriPointf( b.x, b.y ) ) ); + const auto contour = surfaceLineToMeshContour( mesh, mtps ); + ASSERT_TRUE( contour.has_value() ); + + for ( auto mode : { CutMeshParameters::ForceFill::None, CutMeshParameters::ForceFill::Good, CutMeshParameters::ForceFill::All } ) + { + Mesh cut = mesh; + const auto res = cutMesh( cut, { *contour }, { .forceFillMode = mode } ); + EXPECT_TRUE( cut.topology.checkValidity() ); + EXPECT_TRUE( res.fbsWithContourIntersections.test( f ) ); + EXPECT_TRUE( cut.topology.hasFace( f ) ); + EXPECT_EQ( cut.topology.left( e ), f ); + // the two sides of the island are left without faces (Good fills no hole of a bad face), or filled + EXPECT_EQ( cut.topology.findNumHoles(), mode == CutMeshParameters::ForceFill::All ? 0 : 2 ); + } +} + +// a lone contour inside a triangle of zero area (its corners on one line) has no plane to connect it with the triangle in; +// without the exact sorting data, the cut reports the triangle as bad and keeps it +TEST( MRMesh, CutMeshLoneContourInZeroAreaFace ) +{ + // corner 2_v of 1_f is on its side from 1_v to 0_v + const auto mesh = Mesh::fromTriangles( { Vector3f( 0, 0, 0 ), Vector3f( 10, 10, 0 ), Vector3f( 4, 4, 0 ), Vector3f( 0, 10, 0 ) }, + { { 0_v, 1_v, 3_v }, { 1_v, 0_v, 2_v } } ); + const auto e = mesh.topology.edgeWithLeft( 1_f ); + OneMeshContour island; + for ( float t : { 1.f, 3.f, 7.f, 1.f } ) + island.intersections.push_back( { .primitiveId = 1_f, .coordinate = Vector3f( t, t, 0 ) } ); + island.closed = true; + + for ( auto mode : { CutMeshParameters::ForceFill::None, CutMeshParameters::ForceFill::All } ) + { + Mesh cut = mesh; + const auto res = cutMesh( cut, { island }, { .forceFillMode = mode } ); + EXPECT_TRUE( cut.topology.checkValidity() ); + EXPECT_TRUE( res.fbsWithContourIntersections.test( 1_f ) ); + EXPECT_TRUE( cut.topology.hasFace( 1_f ) ); + EXPECT_EQ( cut.topology.left( e ), 1_f ); + // the border of the mesh, and the two sides of the island unless they are filled + EXPECT_EQ( cut.topology.findNumHoles(), mode == CutMeshParameters::ForceFill::All ? 1 : 3 ); + } +} + +// a small square pierces a big triangle of a cube without crossing its edges, so the cube is cut along an open lone contour (a slit) +// of three points on one line, where the sides and the diagonal of the square cross the triangle; the steps of the boolean before its cut +// (as in the terrain embedding) give the exact sorting data, with the square as mesh A, and as mesh B, whose vertex ids the predicates shift +TEST( MRMesh, CutMeshLoneSlitExact ) +{ + const auto cube = makeCube( Vector3f::diagonal( 4.f ), Vector3f::diagonal( -2.f ) ); + // the square of side 0.4 in the plane y = -0.4, centered on the top z = 2 of the cube, in its triangle below the diagonal y = x + auto square = makePlane(); + square.transform( AffineXf3f( Matrix3f( { 0.4f, 0, 0 }, { 0, 0, -0.4f }, { 0, 0.4f, 0 } ), Vector3f( 1.f, -0.4f, 2.f ) ) ); + + for ( bool squareIsA : { true, false } ) + { + const auto& a = squareIsA ? square : cube; + const auto& b = squareIsA ? cube : square; + const auto conv = getVectorConverters( a, b ); + const auto contours = orderIntersectionContours( a.topology, b.topology, findCollidingEdgeTrisPrecise( a, b, conv.toInt ) ); + OneMeshContours cubeContours; + getOneMeshIntersectionContours( a, b, contours, squareIsA ? nullptr : &cubeContours, squareIsA ? &cubeContours : nullptr, conv ); + ASSERT_EQ( cubeContours.size(), 1 ); + ASSERT_EQ( cubeContours[0].intersections.size(), 3 ); + EXPECT_FALSE( cubeContours[0].closed ); + const SortIntersectionsData sortData{ square, contours, conv.toInt, nullptr, a.topology.vertSize(), squareIsA }; + + Mesh cut = cube; + const auto res = cutMesh( cut, cubeContours, { .sortData = &sortData } ); + EXPECT_TRUE( res.fbsWithContourIntersections.none() ) << "square is " << ( squareIsA ? "A" : "B" ); + EXPECT_TRUE( cut.topology.checkValidity() ); + EXPECT_EQ( cut.topology.findNumHoles(), 0 ); // the slit is connected and filled on both sides + EXPECT_NEAR( cut.area(), cube.area(), 1e-6 * cube.area() ); + ASSERT_EQ( res.resultCut.size(), 1 ); + EXPECT_EQ( res.resultCut[0].size(), 2 ); + EXPECT_TRUE( isEdgePath( cut.topology, res.resultCut[0] ) ); + } +} + +namespace +{ + +// the square [0, 10]x[0, 10] in the plane z = 0 of two triangles: 0_f below its diagonal from 0_v to 2_v, 1_f above it +Mesh makeSquareOfTwoTriangles() +{ + return Mesh::fromTriangles( { Vector3f( 0, 0, 0 ), Vector3f( 10, 0, 0 ), Vector3f( 10, 10, 0 ), Vector3f( 0, 10, 0 ) }, + { { 0_v, 1_v, 2_v }, { 0_v, 2_v, 3_v } } ); +} + +// an island inside triangle 0_f of the square: the closed contour around the triangle (6, 2), (8, 2), (8, 4) +OneMeshContour makeIslandInSquare( const Mesh& square ) +{ + std::vector mtps; + for ( const auto& p : { Vector3f( 6, 2, 0 ), Vector3f( 8, 2, 0 ), Vector3f( 8, 4, 0 ), Vector3f( 6, 2, 0 ) } ) + mtps.push_back( square.toTriPoint( 0_f, p ) ); + auto res = surfaceLineToMeshContour( square, mtps ); + EXPECT_TRUE( res.has_value() && res->closed && res->intersections.size() == 4 ); + return res.value_or( OneMeshContour{} ); +} + +// the point of a contour on side ab of the square, which the contour crosses from the right of ab to its left +OneMeshIntersection onSquareSide( const Mesh& square, VertId a, VertId b, const Vector3f& p ) +{ + return { .primitiveId = square.topology.findEdge( a, b ), .coordinate = p }; +} + +// cuts the square by the contours, which cross none of one another, and checks that the cut is valid without bad faces and keeps +// the area, every island is connected, and the maps take every new face, and every new edge inside a triangle, to the triangle it is in +void checkSquareCut( const Mesh& square, const OneMeshContours& contours ) +{ + Mesh cut = square; + FaceMap new2Old; + NewEdgesMap new2OldEdges; + const auto res = cutMesh( cut, contours, { .new2OldMap = &new2Old, .new2oldEdgesMap = &new2OldEdges } ); + const auto& tp = cut.topology; + EXPECT_TRUE( tp.checkValidity() ); + EXPECT_TRUE( res.fbsWithContourIntersections.none() ); + EXPECT_EQ( tp.findNumHoles(), 1 ); // the border of the square + EXPECT_EQ( MeshComponents::getNumComponents( cut ), 1 ); // no island is left unconnected + EXPECT_NEAR( cut.area(), square.area(), 1e-5 * square.area() ); + auto triangleOf = [] ( const Vector3f& p ) { return p.y < p.x ? 0_f : 1_f; }; + for ( auto f : tp.getValidFaces() ) + { + EXPECT_EQ( f < new2Old.size() ? new2Old[f] : FaceId(), triangleOf( cut.triCenter( f ) ) ) << "face " << f; + } + for ( const auto& [ue, old] : new2OldEdges.map ) + { + const EdgeId e = ue; + const auto mid = 0.5f * ( cut.orgPnt( e ) + cut.destPnt( e ) ); + if ( new2OldEdges.splitEdges.test( ue ) || std::abs( mid.x - mid.y ) < 1e-5f ) + continue; // a part of an old edge, or on the diagonal + EXPECT_EQ( FaceId( old ), triangleOf( mid ) ) << "edge " << e; + } +} + +} //anonymous namespace + +// an island inside triangle 0_f of the square, and an open contour from the side of the square in triangle 1_f to the middle of +// the diagonal: the cut of the diagonal fills 0_f again with the end of the open contour as a new corner before the island is +// connected, so the island has nothing to connect with there, and 0_f is reported as bad instead of the island being left unconnected +TEST( MRMesh, CutMeshLoneContourInRefilledFace ) +{ + const auto square = makeSquareOfTwoTriangles(); + const auto island = makeIslandInSquare( square ); + const OneMeshContour open{ .intersections = { + onSquareSide( square, 3_v, 0_v, Vector3f( 0, 5, 0 ) ), onSquareSide( square, 2_v, 0_v, Vector3f( 5, 5, 0 ) ) } }; + for ( bool islandFirst : { false, true } ) + { + Mesh cut = square; + const auto res = cutMesh( cut, islandFirst ? OneMeshContours{ island, open } : OneMeshContours{ open, island } ); + EXPECT_TRUE( cut.topology.checkValidity() ); + EXPECT_TRUE( res.fbsWithContourIntersections.test( 0_f ) ) << "island first " << islandFirst; + } +} + +// a contour along the diagonal from 0_v to 2_v (0_f on its right) or back (0_f on its left), with an island inside 0_f, which the cut +// connects with the rest of 0_f after removing it, or with a contour cutting off the corner of 0_f at 1_v, in either order: the faces +// filled in 0_f map to it, also on the right of the diagonal, where they took the face on its left (1_f) before +TEST( MRMesh, CutMeshContourAlongSideMaps ) +{ + const auto square = makeSquareOfTwoTriangles(); + auto along = [&] ( VertId a, VertId b ) + { + return OneMeshContour{ .intersections = { + { .primitiveId = a, .coordinate = square.points[a] }, { .primitiveId = b, .coordinate = square.points[b] } } }; + }; + const OneMeshContour corner{ .intersections = { + onSquareSide( square, 0_v, 1_v, Vector3f( 8, 0, 0 ) ), onSquareSide( square, 2_v, 1_v, Vector3f( 10, 2, 0 ) ) } }; + const std::pair alongs[] = { { "along 0-2", along( 0_v, 2_v ) }, { "along 2-0", along( 2_v, 0_v ) } }; + const std::pair others[] = { { "island", makeIslandInSquare( square ) }, { "corner", corner } }; + for ( const auto& [alongName, alongContour] : alongs ) + { + for ( const auto& [otherName, other] : others ) + { + for ( bool alongFirst : { true, false } ) + { + SCOPED_TRACE( std::string( alongName ) + ( alongFirst ? " before " : " after " ) + otherName ); + checkSquareCut( square, alongFirst ? OneMeshContours{ alongContour, other } : OneMeshContours{ other, alongContour } ); + } + } + } +} + +// an island inside triangle 0_f of the square, and a contour from the side of the square in triangle 1_f through 0_f to its side +// on the border of the square, in either order: the cut connects the island with the part of 0_f around it on the float path +TEST( MRMesh, CutMeshLoneContourWithCrossingContour ) +{ + const auto square = makeSquareOfTwoTriangles(); + const auto island = makeIslandInSquare( square ); + const OneMeshContour crossing{ .intersections = { onSquareSide( square, 3_v, 0_v, Vector3f( 0, 2, 0 ) ), + onSquareSide( square, 2_v, 0_v, Vector3f( 1.5f, 1.5f, 0 ) ), onSquareSide( square, 1_v, 0_v, Vector3f( 6, 0, 0 ) ) } }; + for ( bool islandFirst : { true, false } ) + { + SCOPED_TRACE( islandFirst ? "island first" : "island second" ); + checkSquareCut( square, islandFirst ? OneMeshContours{ island, crossing } : OneMeshContours{ crossing, island } ); + } +} + TEST( MRMesh, MeshCollidePrecise ) { const auto meshA = makeTorus( 1.1f, 0.5f, 8, 8 ); diff --git a/source/MRTest/MRCutTrianglePredicatesTests.cpp b/source/MRTest/MRCutTrianglePredicatesTests.cpp index 7ee8e4980429..0b31d3240dd8 100644 --- a/source/MRTest/MRCutTrianglePredicatesTests.cpp +++ b/source/MRTest/MRCutTrianglePredicatesTests.cpp @@ -962,7 +962,7 @@ TEST( MRMesh, CutTrianglePredicatesFailed ) EXPECT_TRUE( preds.failed() ); } - // two points defined alike, the second one with the segment backward, which less cannot order: it orders them by their ids, + // two points defined alike without wedges, the second one with the segment backward, which less cannot order: it orders them by their ids, // and alike with the other points, also with a third point at their position, so it stays a strict total order auto alike = [&] { @@ -1009,4 +1009,299 @@ TEST( MRMesh, CutTrianglePredicatesFailed ) EXPECT_TRUE( preds.failed() ); } +namespace +{ + +// the direction in which the crossing of segment ab of p with the plane of f moves when the segment moves toward the wedge w0, w1 of p, +// by u = w0 + w1 - a - b: u - ( n.u ) / ( n.d ) * d for d = b - a and the normal n of f, times |n.d|, projected on axes kx, ky +Vector2i64 shiftOf( const Tri& f, const PreciseCutPoint& p, int kx, int ky ) +{ + const auto n = exactNormal( f ); + const Vector3i64 a( p.seg[0].pt ), b( p.seg[1].pt ); + const auto d = b - a; + const auto u = Vector3i64( p.wedge[0].pt ) + Vector3i64( p.wedge[1].pt ) - a - b; + const auto nd = dot( n, d ); + const auto res = ( nd < 0 ? -nd : nd ) * u - ( nd < 0 ? -dot( n, u ) : dot( n, u ) ) * d; + return { res[kx], res[ky] }; +} + +std::int64_t cross2( const Vector2i64& a, const Vector2i64& b ) +{ + return a.x * b.y - a.y * b.x; +} + +// whether direction y is strictly inside the angle below 180 degrees between directions x and z +bool insideAngle2( const Vector2i64& x, const Vector2i64& z, const Vector2i64& y ) +{ + return ( cross2( x, y ) > 0 ) == ( cross2( x, z ) > 0 ) && ( cross2( z, y ) > 0 ) == ( cross2( z, x ) > 0 ); +} + +} // anonymous namespace + +// Crossings of one segment ab with f defined alike, as of the edges of a multiple edge of another mesh, with different wedges (the third vertices +// of the triangles of the edges, which may share one of them): the predicates must agree with the exact model, where every such point is moved +// off their common position P by a tiny delta along the projection of its shift u (see shiftOf), reversed if the triangles of another of them +// are in the angle between its own ones, which this test finds from the directions of the cuts of f by the triangles; less must be a strict +// total order, and ccw must make a CC system, with the triangles where the first-order terms of the model vanish left to the expansion +// of the perturbation; also, two crossings of a contour from each point along the cuts by its triangles, as the boolean's contour +// goes there, must cross the ones from the other point exactly if the triangles of the two points interleave around ab (the sheets of +// the edges cross); with equal wedges, or without one, the predicates must fail +TEST( MRMesh, CutTrianglePredicatesDefinedAlike ) +{ + RandomVerts rv( 11, 4 ); + CutTrianglePredicates preds; + OrderStats orderStats; + CcwStats ccwStats; + int numSets = 0, numSideBySide = 0, numNested = 0, numInterleaved = 0, numSharedRay = 0, numThree = 0, numFirstOrderZero = 0; + const FastInt256 scale( FastInt128( 1 ) << 100 ); // 1 / delta + for ( int trial = 0; numSets < 400 / cTrialDivisor && trial < 20000; ++trial ) + { + rv.newTrial(); + const Tri f = rv.tri(); + const auto n = exactNormal( f ); + const auto a = rv.vert(), b = rv.vert(); + if ( n == Vector3i64() || dot( n, Vector3i64( b.pt - a.pt ) ) == 0 || !doTriangleSegmentIntersect( { f[0], f[1], f[2], a, b } ) ) + continue; + const auto [kx, ky] = projectionAxes( f ); + // two points defined alike, or three in every fourth set, the second one with a wedge vertex of the first one in every third set + const int numAlike = numSets % 4 == 3 ? 3 : 2; + std::vector opp; + for ( int i = 0; i < 2 * numAlike; ++i ) + opp.push_back( rv.vert() ); + const bool sharedRay = numSets % 3 == 1; + if ( sharedRay ) + opp[2] = opp[0]; + // the direction of the cut of f by triangle abo from P; all of them must differ and not be opposite, so that no exact ties remain + std::vector rays; + for ( const auto& o : opp ) + rays.push_back( shiftOf( f, { .seg = { a, b }, .wedge = { o, o } }, kx, ky ) ); + bool generic = true; + for ( int i = 0; i < int( opp.size() ); ++i ) + for ( int j = i + 1; j < int( opp.size() ); ++j ) + generic = generic && ( opp[i].id == opp[j].id || cross2( rays[i], rays[j] ) != 0 ); + for ( const auto& r : rays ) + generic = generic && r != Vector2i64(); + if ( !generic ) + continue; + + // the points defined alike, with the segment backward, the triangle rotated or reversed, or the wedge swapped in some of them + std::vector alike; + for ( int k = 0; k < numAlike; ++k ) + { + PreciseCutPoint p{ .seg = { a, b }, .tri = f, .wedge = { opp[2 * k], opp[2 * k + 1] } }; + if ( k % 2 == 1 ) + std::swap( p.seg[0], p.seg[1] ); + if ( k == 1 ) + std::rotate( p.tri.begin(), p.tri.begin() + 1, p.tri.end() ); + if ( k == 2 ) + { + std::swap( p.tri[1], p.tri[2] ); + std::swap( p.wedge[0], p.wedge[1] ); + } + alike.push_back( p ); + } + // the shift of every point, reversed if the triangles of another one are in the angle between its own ones (shared or strictly inside) + auto nestedIn = [&] ( int k, int j ) + { + int numShared = 0; + for ( int i : { 2 * j, 2 * j + 1 } ) + { + if ( opp[i].id == opp[2 * k].id || opp[i].id == opp[2 * k + 1].id ) + ++numShared; + else if ( !insideAngle2( rays[2 * k], rays[2 * k + 1], rays[i] ) ) + return false; + } + return numShared < 2; + }; + std::vector shifts; + for ( int k = 0; k < numAlike; ++k ) + { + bool away = false; + for ( int j = 0; j < numAlike; ++j ) + away = away || ( j != k && nestedIn( k, j ) ); + const auto s = shiftOf( f, alike[k], kx, ky ); + shifts.push_back( away ? -s : s ); + } + bool distinct = true; + for ( int k = 0; k < numAlike; ++k ) + for ( int j = k + 1; j < numAlike; ++j ) + distinct = distinct && shifts[k] != shifts[j]; + if ( !distinct || ( numAlike == 3 && cross2( shifts[1] - shifts[0], shifts[2] - shifts[0] ) == 0 ) ) + continue; + if ( numAlike == 2 ) + { + const bool in2 = insideAngle2( rays[0], rays[1], rays[2] ), in3 = insideAngle2( rays[0], rays[1], rays[3] ); + const bool interleaved = !sharedRay && in2 != in3; + numInterleaved += interleaved; + numNested += nestedIn( 0, 1 ) || nestedIn( 1, 0 ); + numSideBySide += !interleaved && !nestedIn( 0, 1 ) && !nestedIn( 1, 0 ); + } + numSharedRay += sharedRay; + numThree += numAlike == 3; + ++numSets; + + // the cuts of f by the triangles abo and two other triangles; the crossing P of ab is replaced by the points defined alike, + // marked on no cut, since the shifts move them off the cuts by the triangles of the others, which ccwByLines takes + std::vector tris; + std::vector rayTri; // the index in tris of the triangle of every wedge vertex + for ( const auto& o : opp ) + { + const auto it = std::find_if( tris.begin(), tris.end(), [&] ( const Tri& t ) { return t[2].id == o.id; } ); + rayTri.push_back( int( it - tris.begin() ) ); + if ( it == tris.end() ) + tris.push_back( { a, b, o } ); + } + for ( int i = 0; i < 2; ++i ) + tris.push_back( rv.tri() ); + auto base = cutPoints( f, tris ); + const auto itP = std::find_if( base.defs.begin(), base.defs.end(), [&] ( const PreciseCutPoint& p ) + { + return !p.isVertex() && ( ( p.seg[0].id == a.id && p.seg[1].id == b.id ) || ( p.seg[0].id == b.id && p.seg[1].id == a.id ) ); + } ); + ASSERT_TRUE( itP != base.defs.end() ); + const int iP = int( itP - base.defs.begin() ); + + // the far end of the cut by the triangle of every wedge vertex + std::vector farEnd; + for ( int i = 0; i < int( opp.size() ); ++i ) + { + int end = -1; + for ( int j = 0; j < int( base.defs.size() ); ++j ) + if ( j != iP && ( base.lines[j] & ( 8 << rayTri[i] ) ) ) + end = j; + EXPECT_GE( end, 0 ); + farEnd.push_back( end ); + } + base.lines[iP] = 0; + + // the points with those defined alike at P and after the others, the first two of them in either order + auto layout = [&] ( int order, std::vector& ids ) + { + auto res = base; + ids = { iP }; + for ( int k = 1; k < numAlike; ++k ) + { + ids.push_back( int( res.defs.size() ) ); + res.add( alike[k], 0 ); + } + if ( order == 1 ) + std::swap( ids[0], ids[1] ); + for ( int k = 0; k < numAlike; ++k ) + res.defs[ids[k]] = alike[k]; + return res; + }; + for ( int order = 0; order < 2; ++order ) + { + std::vector ids; + const auto pts = layout( order, ids ); + const int num = int( pts.defs.size() ); + preds.reset( f ); + for ( int i = 0; i < num; ++i ) + preds.addPoint( VertId( i ), pts.defs[i] ); + ASSERT_FALSE( preds.failed() ); + + // the exact positions, and the moved ones of the points defined alike: P + shift / ( |n.d| * scale ) + std::vector>> pos; + for ( const auto& d : pts.defs ) + pos.push_back( fractionOf( d ) ); + const auto& common = *pos[iP]; + const auto nd = dot( n, Vector3i64( b.pt - a.pt ) ); + const FastInt256 den( common.d * FastInt256( nd < 0 ? -nd : nd ) * scale ); + std::vector alikeOf( num, -1 ); + std::vector> moved( numAlike ); + for ( int k = 0; k < numAlike; ++k ) + { + alikeOf[ids[k]] = k; + auto& m = moved[k]; + m.d = den; + for ( int c = 0; c < 3; ++c ) + m.n[c] = FastInt256( common.n[c] * FastInt256( nd < 0 ? -nd : nd ) * scale ); + m.n[kx] = m.n[kx] + FastInt256( FastInt128( shifts[k].x ) * common.d ); + m.n[ky] = m.n[ky] + FastInt256( FastInt128( shifts[k].y ) * common.d ); + } + auto compare = [&] ( int p, int q, int c ) -> std::optional + { + if ( alikeOf[p] >= 0 && alikeOf[q] >= 0 ) + return compareFractions( moved[alikeOf[p]], moved[alikeOf[q]], c ); + if ( !pos[p] || !pos[q] ) + return {}; + return compareFractions( *pos[p], *pos[q], c ); + }; + auto orient = [&] ( int p, int q, int r ) + { + const int numMoved = ( alikeOf[p] >= 0 ) + ( alikeOf[q] >= 0 ) + ( alikeOf[r] >= 0 ); + if ( numMoved < 2 ) + return pos[p] && pos[q] && pos[r] ? orientFractions( *pos[p], *pos[q], *pos[r], kx, ky ) : 0; + if ( numMoved == 2 ) + { + // the first-order term: cross( wy - wx, z - P ) for the points x, y defined alike and the other one z; + // zero if z - P is parallel to the difference of the shifts, where the expansion of the perturbation decides + const int z = alikeOf[p] < 0 ? p : ( alikeOf[q] < 0 ? q : r ); + const int x = z == p ? q : p, y = z == r ? q : r; + if ( !pos[z] ) + return 0; + auto toward = common; + toward.n[kx] = toward.n[kx] + FastInt256( FastInt128( shifts[alikeOf[y]].x - shifts[alikeOf[x]].x ) * common.d ); + toward.n[ky] = toward.n[ky] + FastInt256( FastInt128( shifts[alikeOf[y]].y - shifts[alikeOf[x]].y ) * common.d ); + if ( orientFractions( common, toward, *pos[z], kx, ky ) == 0 ) + { + ++numFirstOrderZero; + return 0; + } + } + auto at = [&] ( int i ) { return alikeOf[i] >= 0 ? moved[alikeOf[i]] : *pos[i]; }; + return orientFractions( at( p ), at( q ), at( r ), kx, ky ); + }; + checkOrder( preds, f, pts, kx, ky, false, compare, orderStats ); + checkCcw( preds, pts, orient, 0, ccwStats ); + + // the contour from the first two points along the cuts by their triangles, as the boolean's contour goes: + // its pieces from one point cross those from the other one exactly if the triangles interleave + if ( numAlike == 2 ) + { + auto crosses = [&] ( int p, int pEnd, int q, int qEnd ) + { + auto ccw = [&] ( int x, int y, int z ) { return preds.ccw( VertId( x ), VertId( y ), VertId( z ) ); }; + return ccw( p, pEnd, q ) != ccw( p, pEnd, qEnd ) && ccw( q, qEnd, p ) != ccw( q, qEnd, pEnd ); + }; + bool anyCross = false; + for ( int i : { 0, 1 } ) + for ( int j : { 2, 3 } ) + if ( farEnd[i] >= 0 && farEnd[j] >= 0 && farEnd[i] != farEnd[j] ) + anyCross = anyCross || crosses( ids[0], farEnd[i], ids[1], farEnd[j] ); + const bool in2 = insideAngle2( rays[0], rays[1], rays[2] ), in3 = insideAngle2( rays[0], rays[1], rays[3] ); + EXPECT_EQ( anyCross, !sharedRay && in2 != in3 ) << "trial " << trial; + EXPECT_FALSE( preds.failed() ); + } + } + + // with the same wedge, in either order, or without the wedge of one of the points, the predicates cannot tell them apart + if ( numAlike == 2 ) + { + std::vector ids; + const auto pts = layout( 0, ids ); + for ( int k = 0; k < 3; ++k ) + { + auto p = alike[1]; + p.wedge = k == 0 ? alike[0].wedge : ( k == 1 ? std::array{ alike[0].wedge[1], alike[0].wedge[0] } : std::array{} ); + preds.reset( f ); + for ( int i = 0; i < int( pts.defs.size() ); ++i ) + preds.addPoint( VertId( i ), i == ids[1] ? p : pts.defs[i] ); + EXPECT_FALSE( preds.failed() ); + (void)preds.less( VertId( ids[0] ), VertId( ids[1] ) ); + EXPECT_TRUE( preds.failed() ); + } + } + } + EXPECT_EQ( numSets, 400 / cTrialDivisor ); + EXPECT_EQ( numThree, numSets / 4 ); + EXPECT_GE( numSideBySide, atLeast( 110, 17 ) ); + EXPECT_GE( numNested, atLeast( 85, 18 ) ); + EXPECT_GE( numInterleaved, atLeast( 55, 12 ) ); + EXPECT_GE( numSharedRay, atLeast( 110, 22 ) ); + EXPECT_GE( numFirstOrderZero, atLeast( 200, 45 ) ); + EXPECT_GE( orderStats.perturbed, atLeast( 320, 60 ) ); + EXPECT_GE( ccwStats.exact, atLeast( 110000, 21000 ) ); +} + } // namespace MR diff --git a/source/MRTest/MRMeshBooleanTests.cpp b/source/MRTest/MRMeshBooleanTests.cpp index 3c855b6b5019..18692847ebc9 100644 --- a/source/MRTest/MRMeshBooleanTests.cpp +++ b/source/MRTest/MRMeshBooleanTests.cpp @@ -15,7 +15,13 @@ #include #include #include +#include #include +#include +#include +#include +#include +#include #include #include #include @@ -142,6 +148,222 @@ TEST( MRMesh, BooleanDisjointMeshes ) makeBox( { 5.f, 5.f, 5.f }, { 6.f, 6.f, 6.f } ) ); } +// a small sphere pierces a big triangle of the cube without crossing its edges (a lone contour), alone, so that nothing else +// cuts the triangle, or with another sphere crossing the edges of the same triangle, or in the hole of a torus lying in the plane +// of the triangle, which cuts it along two nested circles around the sphere's one; the cut connects the lone contours with the rest +// of the triangle, and the mapper maps the faces made of the triangle back to it +TEST( MRMesh, BooleanLoneContours ) +{ + const Mesh cube = makeCube( Vector3f::diagonal( 4.f ), Vector3f::diagonal( -2.f ) ); + Mesh lone = makeSphere( { .radius = 0.2f, .numMeshVertices = 200 } ); + lone.transform( AffineXf3f::translation( Vector3f( 1.f, -0.4f, 2.f ) ) ); + Mesh crossing = makeSphere( { .radius = 0.3f, .numMeshVertices = 300 } ); + crossing.transform( AffineXf3f::translation( Vector3f( 0.f, 0.f, 2.f ) ) ); + Mesh both = lone; + both.addMesh( crossing ); + // the torus cuts the triangle along circles of radii about 0.4 and 0.8, and the small sphere in its hole along one of radius 0.15, + // around a point 1.13 or more from the sides; the torus is a little above the plane z = 2, which its inner circle of vertices would be in otherwise + Mesh islands = makeTorus( 0.6f, 0.2f, 16, 8 ); + islands.transform( AffineXf3f::translation( Vector3f( 0.8f, -0.8f, 2.01f ) ) ); + Mesh inner = makeSphere( { .radius = 0.15f, .numMeshVertices = 100 } ); + inner.transform( AffineXf3f::translation( Vector3f( 0.8f, -0.8f, 2.f ) ) ); + islands.addMesh( inner ); + const auto xf = AffineXf3f::translation( Vector3f( 0.3f, -0.2f, 0.25f ) ) * + AffineXf3f::linear( Matrix3f::rotation( Vector3f( 1.f, 2.f, 3.f ).normalized(), 0.7f ) ); + // the triangle of the lone contours on the top z = 2 of the cube, below its diagonal y = x + const FaceId loneFace = findProjection( Vector3f( 1.f, -0.4f, 2.f ), cube ).proj.face; + ASSERT_TRUE( loneFace.valid() ); + + // every face of a result comes from a valid face of one input; the faces in the plane of the triangle strictly inside it come from it, + // and those from it are inside it or on its sides + auto checkMapper = [&] ( const BooleanResult& r, const BooleanResultMapper& mapper, const Mesh& a, const Mesh& b, bool cubeIsA ) + { + const auto mapA = mapper.getNew2OldFaceMap( BooleanResultMapper::MapObject::A ); + const auto mapB = mapper.getNew2OldFaceMap( BooleanResultMapper::MapObject::B ); + const float eps = 1e-4f; + int numInside = 0; + for ( auto f : r.mesh.topology.getValidFaces() ) + { + const auto fa = f < mapA.size() ? mapA[f] : FaceId(); + const auto fb = f < mapB.size() ? mapB[f] : FaceId(); + EXPECT_NE( fa.valid(), fb.valid() ); + EXPECT_TRUE( fa ? a.topology.hasFace( fa ) : b.topology.hasFace( fb ) ); + const bool fromLone = ( cubeIsA ? fa : fb ) == loneFace; + const auto c = r.mesh.triCenter( f ); + // the slivers along the contours have no reliable normal, while the other faces of the spheres there are nearly vertical + if ( std::abs( c.z - 2.f ) < eps && r.mesh.area( f ) > 1e-6 && std::abs( r.mesh.normal( f ).z ) > 0.99f + && c.y < c.x - eps && c.x < 2.f - eps && c.y > -2.f + eps ) + { + EXPECT_TRUE( fromLone ); + ++numInside; + } + if ( fromLone ) + { + EXPECT_TRUE( std::abs( c.z - 2.f ) < eps && c.y < c.x + eps && c.x < 2.f + eps && c.y > -2.f - eps ); + } + } + EXPECT_GT( numInside, 0 ); + }; + + for ( const Mesh* spheres : { &lone, &both, &islands } ) + { + const double tol = 1e-5 * ( cube.volume() + spheres->volume() ); + for ( bool swap : { false, true } ) + { + for ( const AffineXf3f* rigidB2A : { (const AffineXf3f*)nullptr, &xf } ) + { + // rigidB2A moves the second mesh back in place + const Mesh& a = swap ? *spheres : cube; + Mesh b = swap ? cube : *spheres; + if ( rigidB2A ) + b.transform( rigidB2A->inverse() ); + BooleanResultMapper interMapper, uniMapper, diffMapper; + const auto inter = boolean( a, b, BooleanOperation::Intersection, rigidB2A, &interMapper ); + const auto uni = boolean( a, b, BooleanOperation::Union, rigidB2A, &uniMapper ); + const auto diff = boolean( a, b, BooleanOperation::DifferenceAB, rigidB2A, &diffMapper ); + ASSERT_TRUE( inter.valid() ) << inter.errorString; + ASSERT_TRUE( uni.valid() ) << uni.errorString; + ASSERT_TRUE( diff.valid() ) << diff.errorString; + for ( const auto* r : { &inter, &uni, &diff } ) + EXPECT_EQ( r->mesh.topology.findNumHoles(), 0 ); + // an island left in its triangle would make a separate component; the components of the intersection are the lower halves + // of the spheres and of the torus, which do not touch one another + EXPECT_EQ( MeshComponents::getNumComponents( uni.mesh ), 1 ); + EXPECT_EQ( MeshComponents::getNumComponents( inter.mesh ), spheres == &lone ? 1 : 2 ); + EXPECT_NEAR( uni.mesh.volume() + inter.mesh.volume(), a.volume() + b.volume(), tol ); + EXPECT_NEAR( diff.mesh.volume() + inter.mesh.volume(), a.volume(), tol ); + checkMapper( inter, interMapper, a, b, !swap ); + checkMapper( uni, uniMapper, a, b, !swap ); + checkMapper( diff, diffMapper, a, b, !swap ); + } + } + } +} + +// the precut output of the mesh that the operation does not cut (A for InsideB, B for InsideA) takes it only after the cut +// of the other mesh, which sorts its intersections by that mesh itself, not by a copy; here the cut of the cube connects +// a lone contour, whose exact points are made of the vertices of the sphere +TEST( MRMesh, BooleanPreCutOfUncutMesh ) +{ + const Mesh cube = makeCube( Vector3f::diagonal( 4.f ), Vector3f::diagonal( -2.f ) ); + Mesh sphere = makeSphere( { .radius = 0.2f, .numMeshVertices = 200 } ); + sphere.transform( AffineXf3f::translation( Vector3f( 1.f, -0.4f, 2.f ) ) ); + + for ( bool swap : { false, true } ) + { + const Mesh& a = swap ? sphere : cube; + const Mesh& b = swap ? cube : sphere; + const auto op = swap ? BooleanOperation::InsideB : BooleanOperation::InsideA; + const auto cubeObj = swap ? BooleanResultMapper::MapObject::B : BooleanResultMapper::MapObject::A; + BooleanResultMapper fullMapper; + const auto full = boolean( a, b, op, { .mapper = &fullMapper } ); + ASSERT_TRUE( full.valid() ) << full.errorString; + EXPECT_TRUE( full.mesh.topology.checkValidity() ); + EXPECT_EQ( full.mesh.topology.findNumHoles(), 1 ); // the disc of the cube inside the sphere + EXPECT_GT( full.mesh.topology.numValidFaces(), 0 ); + + // with the precut output, through the const& overload and the Mesh&& one, the cube is cut as without it + for ( bool rvalue : { false, true } ) + { + BooleanPreCutResult preCut; + BooleanResultMapper mapper; + const BooleanParameters params{ .mapper = &mapper, .outPreCutA = swap ? &preCut : nullptr, .outPreCutB = swap ? nullptr : &preCut }; + const auto res = rvalue ? boolean( Mesh( a ), Mesh( b ), op, params ) : boolean( a, b, op, params ); + ASSERT_TRUE( res.valid() ) << res.errorString; + EXPECT_TRUE( res.meshABadContourFaces.none() && res.meshBBadContourFaces.none() ); + EXPECT_TRUE( preCut.mesh == sphere ); + EXPECT_TRUE( mapper.getMaps( cubeObj ).cut2origin == fullMapper.getMaps( cubeObj ).cut2origin ); + } + } +} + +// the precut output of either mesh of an operation cutting both, or of both meshes at once, has what cutMesh needs for the boolean's +// exact cut: the cut of the precut mesh by its contours with sortData( the other mesh ) is the boolean's own one, here of a cube +// with a lone contour of a small sphere, as mesh A and as mesh B, also with a rotated rigidB2A, through both overloads (where the Mesh&& +// one copies no mesh: the cut of the other mesh sorts its intersections by the uncut one itself, and with both outputs nothing is cut) +TEST( MRMesh, BooleanPreCutExactCut ) +{ + const Mesh cube = makeCube( Vector3f::diagonal( 4.f ), Vector3f::diagonal( -2.f ) ); + Mesh sphere = makeSphere( { .radius = 0.2f, .numMeshVertices = 200 } ); + sphere.transform( AffineXf3f::translation( Vector3f( 1.f, -0.4f, 2.f ) ) ); + const auto xf = AffineXf3f::translation( Vector3f( 0.3f, -0.2f, 0.25f ) ) * + AffineXf3f::linear( Matrix3f::rotation( Vector3f( 1.f, 2.f, 3.f ).normalized(), 0.7f ) ); + const auto op = BooleanOperation::Intersection; // A is filled on the left of its cut, B on the right + + for ( bool swap : { false, true } ) + { + for ( const AffineXf3f* rigidB2A : { (const AffineXf3f*)nullptr, &xf } ) + { + // rigidB2A moves the second mesh back in place + const Mesh& a = swap ? sphere : cube; + Mesh b = swap ? cube : sphere; + if ( rigidB2A ) + b.transform( rigidB2A->inverse() ); + BooleanResultMapper mapper; + const auto full = boolean( a, b, op, { .rigidB2A = rigidB2A, .mapper = &mapper } ); + ASSERT_TRUE( full.valid() ) << full.errorString; + + for ( bool wantA : { true, false } ) + { + for ( bool wantB : { true, false } ) + { + if ( !wantA && !wantB ) + continue; + for ( bool rvalue : { false, true } ) + { + BooleanPreCutResult preA, preB; + const BooleanParameters params{ .rigidB2A = rigidB2A, .outPreCutA = wantA ? &preA : nullptr, .outPreCutB = wantB ? &preB : nullptr }; + const auto res = rvalue ? boolean( Mesh( a ), Mesh( b ), op, params ) : boolean( a, b, op, params ); + ASSERT_TRUE( res.valid() ) << res.errorString; + for ( bool isA : { true, false } ) + { + if ( !( isA ? wantA : wantB ) ) + continue; + SCOPED_TRACE( std::string( swap ? "sphere is A" : "cube is A" ) + ( rigidB2A ? ", rigidB2A" : "" ) + + ( wantA && wantB ? ", both outputs" : "" ) + ( isA ? ", output A" : ", output B" ) + ( rvalue ? ", Mesh&&" : "" ) ); + const auto& pre = isA ? preA : preB; + EXPECT_TRUE( pre.mesh == ( isA ? a : b ) ); + EXPECT_EQ( pre.isMeshA, isA ); + EXPECT_EQ( pre.meshAVertsNum, a.topology.vertSize() ); + EXPECT_FALSE( pre.intersections.empty() ); + EXPECT_EQ( pre.rigidB2A.has_value(), rigidB2A != nullptr ); + if ( rigidB2A && pre.rigidB2A ) + { + EXPECT_TRUE( *pre.rigidB2A == *rigidB2A ); + } + + Mesh cut = pre.mesh; + FaceMap new2Old; + const auto sortData = pre.sortData( isA ? b : a ); + const auto cutRes = cutMesh( cut, pre.contours, { .sortData = &sortData, .new2OldMap = &new2Old, + .fillPart = isA ? CutMeshParameters::FillPart::Left : CutMeshParameters::FillPart::Right } ); + EXPECT_TRUE( cutRes.fbsWithContourIntersections.none() ); + EXPECT_TRUE( new2Old == mapper.getMaps( isA ? BooleanResultMapper::MapObject::A : BooleanResultMapper::MapObject::B ).cut2origin ); + } + } + } + } + } + } +} + +// selfBoolean of a mesh of two parts, a cube and a small sphere piercing a big triangle of it without crossing its edges +// (a lone contour); the parts are only reconnected along their intersection, so the result keeps the sum of their volumes +TEST( MRMesh, SelfBooleanLoneContour ) +{ + Mesh mesh = makeCube( Vector3f::diagonal( 4.f ), Vector3f::diagonal( -2.f ) ); + Mesh sphere = makeSphere( { .radius = 0.2f, .numMeshVertices = 200 } ); + sphere.transform( AffineXf3f::translation( Vector3f( 1.f, -0.4f, 2.f ) ) ); + const double volume = mesh.volume() + sphere.volume(); + mesh.addMesh( sphere ); + + const auto res = selfBoolean( mesh ); + ASSERT_TRUE( res.has_value() ) << res.error(); + EXPECT_TRUE( res->topology.checkValidity() ); + EXPECT_EQ( res->topology.findNumHoles(), 0 ); + EXPECT_NEAR( res->volume(), volume, 1e-5 * volume ); +} + // closed mesh of two parts (apexes 3_v and 6_v) stitched along their common boundary, which passes through all vertices of triangle 0_v, 2_v, 1_v // of each part, so these two triangles share edges 0_v-2_v and 0_v-1_v, while their edges 1_v-2_v make a multiple edge static Mesh makeMeshWithMultipleEdge() @@ -188,6 +410,167 @@ TEST( MRMesh, BooleanMultipleEdgePropogationSort ) } } +// both edges of the multiple edge (see makeMeshWithMultipleEdge) cross a big triangle of the cube at one point inside it, on a lone contour +// that goes through this point twice, from the upper half of the lens to the lower one; the exact predicates of the cut tell apart the two +// crossings defined alike by the third vertices of the triangles of each edge, so every boolean and selfBoolean give valid closed meshes +// with the volumes of the halves of the lens on either side of the cube's side x = 0.1 +TEST( MRMesh, BooleanMultipleEdgeThroughLoneFace ) +{ + const Mesh lens = makeMeshWithMultipleEdge(); + const Mesh cube = makeCube( Vector3f::diagonal( 20.f ), Vector3f( 0.1f, -9.f, -0.6f ) ); + const double volCube = cube.volume(), volLens = lens.volume(); + EXPECT_NEAR( volCube, 8000., 1e-6 * 8000. ); + EXPECT_NEAR( volLens, 0.25, 1e-6 * 0.25 ); + for ( bool swap : { false, true } ) + { + const Mesh& a = swap ? lens : cube; + const Mesh& b = swap ? cube : lens; + double vols[3] = {}; + // the volumes of the intersection and the union, and of the difference of the cube and the lens or the other way round + const double expected[3] = { 0.100166672, 8000.149845, swap ? 0.1498333365 : 7999.899845 }; + for ( int i = 0; i < 3; ++i ) + { + const auto op = std::array{ BooleanOperation::Intersection, BooleanOperation::Union, BooleanOperation::DifferenceAB }[i]; + const auto res = boolean( a, b, op ); + ASSERT_TRUE( res.valid() ) << res.errorString; + EXPECT_TRUE( res.mesh.topology.checkValidity() ); + EXPECT_EQ( res.mesh.topology.findNumHoles(), 0 ); + vols[i] = res.mesh.volume(); + EXPECT_NEAR( vols[i], expected[i], 1e-6 * expected[i] ); + } + const double volA = swap ? volLens : volCube, volB = swap ? volCube : volLens; + EXPECT_NEAR( vols[1] + vols[0], volA + volB, 1e-6 * ( volA + volB ) ); + EXPECT_NEAR( vols[2], volA - vols[0], 1e-6 * ( volA + volB ) ); + } + + Mesh both = cube; + both.addMesh( lens ); + const auto res = selfBoolean( both ); + ASSERT_TRUE( res.has_value() ) << res.error(); + EXPECT_TRUE( res->topology.checkValidity() ); + EXPECT_EQ( res->topology.findNumHoles(), 0 ); +} + +// the tip at 0_v of two triangles with the same vertices (see makeMeshWithMultipleEdge) pokes through a big triangle of the cube: +// edges 0_v-1_v and 0_v-2_v, which both triangles share, cross it at two points, on a closed lone contour of two points; it has zero area +// and is dropped, so the boolean and selfBoolean give valid closed meshes +TEST( MRMesh, BooleanTwoPointLoneContour ) +{ + const Mesh lens = makeMeshWithMultipleEdge(); + const Mesh cube = makeCube( Vector3f::diagonal( 20.f ), Vector3f( -3.f, 0.5f, -15.f ) ); // only its side y = 0.5 crosses the lens + for ( bool swap : { false, true } ) + { + const Mesh& a = swap ? lens : cube; + const Mesh& b = swap ? cube : lens; + for ( auto op : { BooleanOperation::Intersection, BooleanOperation::Union, BooleanOperation::DifferenceAB } ) + { + const auto res = boolean( a, b, op ); + ASSERT_TRUE( res.valid() ) << res.errorString; + EXPECT_TRUE( res.mesh.topology.checkValidity() ); + EXPECT_EQ( res.mesh.topology.findNumHoles(), 0 ); + } + } + + Mesh both = cube; + both.addMesh( lens ); + const auto res = selfBoolean( both ); + ASSERT_TRUE( res.has_value() ) << res.error(); + EXPECT_TRUE( res->topology.checkValidity() ); + EXPECT_EQ( res->topology.findNumHoles(), 0 ); +} + +// a corner of a small cube, or the lowest vertex of a small sphere, lies exactly on a big triangle of a cube, the rest of it outside; +// the edges from that vertex meet the triangle at one point, so simulation of simplicity makes a lone contour of points that are all +// the same in float (here with the big cube as mesh A), or no contact (as mesh B); either way the results are valid closed meshes +// with the volumes of meshes touching at a point +TEST( MRMesh, BooleanVertexOnFace ) +{ + const Mesh cube = makeCube( Vector3f::diagonal( 4.f ), Vector3f::diagonal( -2.f ) ); + const Vector3f p( 1.f, -0.4f, 2.f ); // on the top z = 2 of the cube, in its triangle below the diagonal y = x + Mesh corner = makeCube( Vector3f::diagonal( 0.4f ), Vector3f::diagonal( -0.2f ) ); + corner.transform( AffineXf3f::linear( Matrix3f::rotation( Vector3f::diagonal( -1.f ), Vector3f::minusZ() ) ) ); // a corner down + Mesh sphere = makeSphere( { .radius = 0.2f, .numMeshVertices = 200 } ); + for ( Mesh* m : { &corner, &sphere } ) + { + VertId low; + for ( auto v : m->topology.getValidVerts() ) + if ( !low || m->points[v].z < m->points[low].z ) + low = v; + m->transform( AffineXf3f::translation( p - m->points[low] ) ); + m->points[low] = p; // exactly + m->invalidateCaches(); + } + + for ( const Mesh* tip : { &corner, &sphere } ) + { + for ( bool swap : { false, true } ) + { + const Mesh& a = swap ? *tip : cube; + const Mesh& b = swap ? cube : *tip; + for ( auto op : { BooleanOperation::Intersection, BooleanOperation::Union, BooleanOperation::DifferenceAB } ) + { + const auto res = boolean( a, b, op ); + ASSERT_TRUE( res.valid() ) << res.errorString; + EXPECT_TRUE( res.mesh.topology.checkValidity() ); + EXPECT_EQ( res.mesh.topology.findNumHoles(), 0 ); + const double volume = op == BooleanOperation::Union ? a.volume() + b.volume() : op == BooleanOperation::Intersection ? 0 : a.volume(); + EXPECT_NEAR( res.mesh.volume(), volume, 1e-6 * cube.volume() ); + } + } + } +} + +// a peak of the terrain pokes through a big triangle of the structure without crossing its edges, so the structure is cut along +// a lone contour; far from the origin with the peak just above the structure, all points of the contour are the same in float, +// so only the exact predicates tell how to connect it with the corners of the triangle +TEST( MRMesh, EmbedStructureLoneContour ) +{ + for ( float shift : { 0.f, 1000.f } ) + { + const float peak = shift == 0 ? 2.f : std::nextafter( 1.f, 2.f ); + const auto terrain = makeRegularGridMesh( 21, 21, [] ( size_t, size_t ) { return true; }, [&] ( size_t x, size_t y ) + { + return Vector3f( shift + float( x ) - 10.f, shift + float( y ) - 10.f, x == 10 && y == 10 ? peak : 0.f ); + } ); + ASSERT_TRUE( terrain.has_value() ); + Mesh structure = makePlane(); // the peak is inside one of its two triangles + structure.transform( AffineXf3f( Matrix3f::scale( 8.f ), Vector3f( shift + 1.f, shift - 1.f, 1.f ) ) ); + const auto res = embedStructureToTerrain( *terrain, structure, { .fillAngle = PI_F / 4, .cutAngle = PI_F / 4 } ); + ASSERT_TRUE( res.has_value() ) << res.error(); + EXPECT_TRUE( res->topology.checkValidity() ); + EXPECT_EQ( res->topology.findNumHoles(), 1 ); // the border of the terrain + } +} + +// the wall of a small structure above a terrain of two big triangles reaches the terrain inside one of them, so the terrain is cut +// along a lone contour, which the cut connects with the rest of the triangle; every face of the terrain in the result maps to +// the triangle of the input terrain it is in +TEST( MRMesh, EmbedStructureLoneTerrainContour ) +{ + Mesh terrain = makePlane(); + terrain.transform( AffineXf3f::linear( Matrix3f::scale( 40.f ) ) ); + Mesh structure = makePlane(); + structure.transform( AffineXf3f::translation( Vector3f( 8.f, -8.f, 0.5f ) ) ); // above the triangle below the diagonal y = x + const FaceId lower = findProjection( Vector3f( 8.f, -8.f, 0.f ), terrain ).proj.face; + ASSERT_TRUE( lower.valid() ); + FaceMap new2Old; + const auto res = embedStructureToTerrain( terrain, structure, { .fillAngle = PI_F / 4, .cutAngle = PI_F / 4, .new2oldFaces = &new2Old } ); + ASSERT_TRUE( res.has_value() ) << res.error(); + EXPECT_TRUE( res->topology.checkValidity() ); + EXPECT_EQ( res->topology.findNumHoles(), 1 ); // the border of the terrain + int numTerrainFaces = 0; + for ( auto f : res->topology.getValidFaces() ) + { + if ( f >= new2Old.size() ) + continue; // the faces of the structure and its walls come after those of the cut terrain + ++numTerrainFaces; + const auto c = res->triCenter( f ); + EXPECT_TRUE( terrain.topology.hasFace( new2Old[f] ) ) << "face " << f; + EXPECT_EQ( new2Old[f] == lower, c.y < c.x ) << "face " << f; + } + EXPECT_GT( numTerrainFaces, 0 ); +} + // two triangles sharing two vertices but not the edge between them (multiple edges connect the vertices) can meet only along // the segment between these vertices, so findSelfCollidingEdgeTrisPrecise skips them as the triangles with a shared edge, // but still finds all crossings of other edges and triangles @@ -284,6 +667,10 @@ TEST( MRMesh, BooleanResultMapper ) EXPECT_EQ( mapsB.cut2origin.size(), 320 ); } +// the tests of the functions deprecated since cutMesh connects lone contours itself +MR_SUPPRESS_WARNING_PUSH +MR_SUPPRESS_WARNING( "-Wdeprecated-declarations", 4996 ) + // the spheres of the Boolean benchmark (B rotated around Z by n times 0.1 degree) have lone contours (each inside one triangle) for some n; // subdivideLoneContours splits those triangles and updates the AABB tree instead of its rebuild TEST( MRMesh, SubdivideLoneContoursUpdatesAABBTree ) @@ -388,4 +775,6 @@ TEST( MRMesh, UpdateCollidingEdgeTrisPrecise ) } } +MR_SUPPRESS_WARNING_POP + } //namespace MR