From 4cdee629a64fba845983fed120bccaef57db47e9 Mon Sep 17 00:00:00 2001 From: Pasukhin Dmitry Date: Tue, 16 Dec 2025 15:04:00 +0000 Subject: [PATCH] Mesh - Fix point-in-polygon check for CCW polygons in BRepMesh_Delaun (#920) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The isVertexInsidePolygon method uses the winding number algorithm, which computes cumulative angles from the test point to polygon vertices. For a point inside a polygon, this sum should be ±2π. The gp_Vec2d::Angle method uses an inverted sign convention where CCW rotation gives negative angles. For CCW polygons, the total angle is -2π, but the original check only accepted +2π. Changed the condition from: std::abs(Angle2PI - aTotalAng) > Precision::Angular() to: std::abs(std::abs(aTotalAng) - Angle2PI) > Precision::Angular() This correctly handles both CCW (-2π) and CW (+2π) polygon orientations. --- .../TKMesh/BRepMesh/BRepMesh_Delaun.cxx | 2 +- .../TKMesh/GTests/BRepMesh_Delaun_Test.cxx | 237 ++++++++++++++++++ .../TKMesh/GTests/FILES.cmake | 1 + 3 files changed, 239 insertions(+), 1 deletion(-) create mode 100644 src/ModelingAlgorithms/TKMesh/GTests/BRepMesh_Delaun_Test.cxx diff --git a/src/ModelingAlgorithms/TKMesh/BRepMesh/BRepMesh_Delaun.cxx b/src/ModelingAlgorithms/TKMesh/BRepMesh/BRepMesh_Delaun.cxx index 703d9cc729..227a340975 100644 --- a/src/ModelingAlgorithms/TKMesh/BRepMesh/BRepMesh_Delaun.cxx +++ b/src/ModelingAlgorithms/TKMesh/BRepMesh/BRepMesh_Delaun.cxx @@ -1550,7 +1550,7 @@ Standard_Boolean BRepMesh_Delaun::isVertexInsidePolygon( aPrevVertexDir = aCurVertexDir; } - if (std::abs(Angle2PI - aTotalAng) > Precision::Angular()) + if (std::abs(std::abs(aTotalAng) - Angle2PI) > Precision::Angular()) return Standard_False; return Standard_True; diff --git a/src/ModelingAlgorithms/TKMesh/GTests/BRepMesh_Delaun_Test.cxx b/src/ModelingAlgorithms/TKMesh/GTests/BRepMesh_Delaun_Test.cxx new file mode 100644 index 0000000000..f8e31cabd2 --- /dev/null +++ b/src/ModelingAlgorithms/TKMesh/GTests/BRepMesh_Delaun_Test.cxx @@ -0,0 +1,237 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +// Test for bug 0032395: BRepMesh_Delaun::isVertexInsidePolygon +// The method uses winding number algorithm to check if a point is inside a polygon. +// For CCW polygons, the cumulative angle is -2*PI, but the original code only +// checked for +2*PI. This test verifies that both orientations are handled correctly. + +TEST(BRepMesh_DelaunTest, Vec2dAngleSignConvention) +{ + // Verify gp_Vec2d::Angle sign convention used in isVertexInsidePolygon + // The angle from aPrev to aCur is computed as aCur.Angle(aPrev) + + // Test CCW rotation: from (1,0) to (0,1) is 90 degrees CCW + gp_Vec2d aVecRight(1.0, 0.0); + gp_Vec2d aVecUp(0.0, 1.0); + + // aCur.Angle(aPrev) where aCur=up, aPrev=right + // Going from right to up is CCW, angle should be negative in OCCT convention + const double anAngleCCW = aVecUp.Angle(aVecRight); + EXPECT_LT(anAngleCCW, 0.0) << "CCW rotation should give negative angle in OCCT convention"; + EXPECT_NEAR(anAngleCCW, -M_PI_2, Precision::Angular()); + + // Test CW rotation: from (0,1) to (1,0) is 90 degrees CW + const double anAngleCW = aVecRight.Angle(aVecUp); + EXPECT_GT(anAngleCW, 0.0) << "CW rotation should give positive angle in OCCT convention"; + EXPECT_NEAR(anAngleCW, M_PI_2, Precision::Angular()); +} + +TEST(BRepMesh_DelaunTest, WindingAngleCCWPolygon) +{ + // Simulate the winding angle calculation for a CCW square polygon + // with a point at the center + // Vertices: (1,0), (0,1), (-1,0), (0,-1) - CCW order + // Center point: (0,0) + + gp_XY aCenter(0.0, 0.0); + gp_XY aVertices[5] = { + gp_XY(1.0, 0.0), + gp_XY(0.0, 1.0), + gp_XY(-1.0, 0.0), + gp_XY(0.0, -1.0), + gp_XY(1.0, 0.0) // Close the polygon + }; + + gp_Vec2d aPrevDir(aVertices[0] - aCenter); + double aTotalAngle = 0.0; + + for (int i = 1; i < 5; ++i) + { + gp_Vec2d aCurDir(aVertices[i] - aCenter); + aTotalAngle += aCurDir.Angle(aPrevDir); + aPrevDir = aCurDir; + } + + // For CCW polygon with point inside, total angle should be -2*PI + EXPECT_NEAR(std::abs(aTotalAngle), 2.0 * M_PI, Precision::Angular()) + << "Winding angle magnitude should be 2*PI for point inside polygon"; + + // The fix checks std::abs(std::abs(aTotalAng) - Angle2PI), so both +2PI and -2PI work + const double Angle2PI = 2.0 * M_PI; + const bool isInsideOldFix = std::abs(Angle2PI - aTotalAngle) <= Precision::Angular(); + const bool isInsideNewFix = std::abs(std::abs(aTotalAngle) - Angle2PI) <= Precision::Angular(); + + EXPECT_FALSE(isInsideOldFix) << "Old check fails for CCW polygon (angle is negative)"; + EXPECT_TRUE(isInsideNewFix) << "New check handles both CCW and CW polygons"; +} + +TEST(BRepMesh_DelaunTest, WindingAngleCWPolygon) +{ + // Simulate the winding angle calculation for a CW square polygon + // Vertices in CW order: (1,0), (0,-1), (-1,0), (0,1) + // Center point: (0,0) + + gp_XY aCenter(0.0, 0.0); + gp_XY aVertices[5] = { + gp_XY(1.0, 0.0), + gp_XY(0.0, -1.0), + gp_XY(-1.0, 0.0), + gp_XY(0.0, 1.0), + gp_XY(1.0, 0.0) // Close the polygon + }; + + gp_Vec2d aPrevDir(aVertices[0] - aCenter); + double aTotalAngle = 0.0; + + for (int i = 1; i < 5; ++i) + { + gp_Vec2d aCurDir(aVertices[i] - aCenter); + aTotalAngle += aCurDir.Angle(aPrevDir); + aPrevDir = aCurDir; + } + + // For CW polygon with point inside, total angle should be +2*PI + EXPECT_NEAR(aTotalAngle, 2.0 * M_PI, Precision::Angular()) + << "Winding angle should be +2*PI for CW polygon"; + + const double Angle2PI = 2.0 * M_PI; + const bool isInsideOldFix = std::abs(Angle2PI - aTotalAngle) <= Precision::Angular(); + const bool isInsideNewFix = std::abs(std::abs(aTotalAngle) - Angle2PI) <= Precision::Angular(); + + EXPECT_TRUE(isInsideOldFix) << "Old check works for CW polygon"; + EXPECT_TRUE(isInsideNewFix) << "New check also works for CW polygon"; +} + +TEST(BRepMesh_DelaunTest, MeshPlanarFaceWithHole) +{ + // Create a planar face with a hole - this exercises the Delaunay triangulation + // code path that uses isVertexInsidePolygon during mesh refinement + + // Outer wire: square 10x10 centered at origin + gp_Pnt aP1(-5.0, -5.0, 0.0); + gp_Pnt aP2(5.0, -5.0, 0.0); + gp_Pnt aP3(5.0, 5.0, 0.0); + gp_Pnt aP4(-5.0, 5.0, 0.0); + + BRepBuilderAPI_MakeWire anOuterWireMaker; + anOuterWireMaker.Add(BRepBuilderAPI_MakeEdge(aP1, aP2)); + anOuterWireMaker.Add(BRepBuilderAPI_MakeEdge(aP2, aP3)); + anOuterWireMaker.Add(BRepBuilderAPI_MakeEdge(aP3, aP4)); + anOuterWireMaker.Add(BRepBuilderAPI_MakeEdge(aP4, aP1)); + ASSERT_TRUE(anOuterWireMaker.IsDone()); + + // Inner wire (hole): smaller square 2x2 centered at origin + gp_Pnt aH1(-1.0, -1.0, 0.0); + gp_Pnt aH2(1.0, -1.0, 0.0); + gp_Pnt aH3(1.0, 1.0, 0.0); + gp_Pnt aH4(-1.0, 1.0, 0.0); + + BRepBuilderAPI_MakeWire anInnerWireMaker; + anInnerWireMaker.Add(BRepBuilderAPI_MakeEdge(aH1, aH2)); + anInnerWireMaker.Add(BRepBuilderAPI_MakeEdge(aH2, aH3)); + anInnerWireMaker.Add(BRepBuilderAPI_MakeEdge(aH3, aH4)); + anInnerWireMaker.Add(BRepBuilderAPI_MakeEdge(aH4, aH1)); + ASSERT_TRUE(anInnerWireMaker.IsDone()); + + // Create face with hole + Handle(Geom_Plane) aPlane = new Geom_Plane(gp_Pln(gp::Origin(), gp::DZ())); + BRepBuilderAPI_MakeFace aFaceMaker(aPlane, anOuterWireMaker.Wire()); + aFaceMaker.Add(anInnerWireMaker.Wire()); + ASSERT_TRUE(aFaceMaker.IsDone()); + + TopoDS_Face aFace = aFaceMaker.Face(); + + // Mesh the face + BRepMesh_IncrementalMesh aMesher(aFace, 0.1); + EXPECT_TRUE(aMesher.IsDone()) << "Meshing should succeed"; + + // Verify triangulation exists + TopLoc_Location aLoc; + const Handle(Poly_Triangulation) aTri = BRep_Tool::Triangulation(aFace, aLoc); + ASSERT_FALSE(aTri.IsNull()) << "Triangulation should be created"; + EXPECT_GT(aTri->NbTriangles(), 0) << "Should have triangles"; + EXPECT_GT(aTri->NbNodes(), 0) << "Should have nodes"; +} + +TEST(BRepMesh_DelaunTest, MeshBoxAllFaces) +{ + // Create and mesh a box - tests triangulation on faces with different orientations + TopoDS_Shape aBox = BRepPrimAPI_MakeBox(10.0, 10.0, 10.0).Shape(); + + BRepMesh_IncrementalMesh aMesher(aBox, 0.5); + EXPECT_TRUE(aMesher.IsDone()) << "Meshing should succeed"; + + // Verify all faces are triangulated + int aFaceCount = 0; + for (TopExp_Explorer anExp(aBox, TopAbs_FACE); anExp.More(); anExp.Next()) + { + const TopoDS_Face& aFace = TopoDS::Face(anExp.Current()); + TopLoc_Location aLoc; + const Handle(Poly_Triangulation) aTri = BRep_Tool::Triangulation(aFace, aLoc); + + EXPECT_FALSE(aTri.IsNull()) << "Face " << aFaceCount << " should have triangulation"; + if (!aTri.IsNull()) + { + EXPECT_GT(aTri->NbTriangles(), 0) << "Face " << aFaceCount << " should have triangles"; + } + ++aFaceCount; + } + + EXPECT_EQ(aFaceCount, 6) << "Box should have 6 faces"; +} + +TEST(BRepMesh_DelaunTest, MeshCylinderCurvedFaces) +{ + // Create and mesh a cylinder - curved faces may have different polygon orientations + TopoDS_Shape aCylinder = BRepPrimAPI_MakeCylinder(5.0, 10.0).Shape(); + + BRepMesh_IncrementalMesh aMesher(aCylinder, 0.5); + EXPECT_TRUE(aMesher.IsDone()) << "Meshing should succeed"; + + // Verify faces are triangulated + int aTotalTriangles = 0; + for (TopExp_Explorer anExp(aCylinder, TopAbs_FACE); anExp.More(); anExp.Next()) + { + const TopoDS_Face& aFace = TopoDS::Face(anExp.Current()); + TopLoc_Location aLoc; + const Handle(Poly_Triangulation) aTri = BRep_Tool::Triangulation(aFace, aLoc); + + EXPECT_FALSE(aTri.IsNull()) << "Face should have triangulation"; + if (!aTri.IsNull()) + { + aTotalTriangles += aTri->NbTriangles(); + } + } + + EXPECT_GT(aTotalTriangles, 0) << "Cylinder should have triangles"; +} diff --git a/src/ModelingAlgorithms/TKMesh/GTests/FILES.cmake b/src/ModelingAlgorithms/TKMesh/GTests/FILES.cmake index df0a8f7a27..c6dd438425 100644 --- a/src/ModelingAlgorithms/TKMesh/GTests/FILES.cmake +++ b/src/ModelingAlgorithms/TKMesh/GTests/FILES.cmake @@ -2,5 +2,6 @@ set(OCCT_TKMesh_GTests_FILES_LOCATION "${CMAKE_CURRENT_LIST_DIR}") set(OCCT_TKMesh_GTests_FILES + BRepMesh_Delaun_Test.cxx BRepMesh_GeomTool_Test.cxx )