From 78df4424b70769d4cbb48a4c02463260cf4e03de Mon Sep 17 00:00:00 2001 From: wmayer Date: Tue, 11 Dec 2018 18:42:47 +0100 Subject: [PATCH] fix MeshGeomFacet::IntersectWithPlane --- src/Mod/Mesh/App/Core/Elements.cpp | 174 ++++++++++++++++++----------- 1 file changed, 111 insertions(+), 63 deletions(-) diff --git a/src/Mod/Mesh/App/Core/Elements.cpp b/src/Mod/Mesh/App/Core/Elements.cpp index 016f283be5..6faf90af38 100644 --- a/src/Mod/Mesh/App/Core/Elements.cpp +++ b/src/Mod/Mesh/App/Core/Elements.cpp @@ -515,75 +515,123 @@ bool MeshGeomFacet::IntersectBoundingBox ( const Base::BoundBox3f &rclBB ) const bool MeshGeomFacet::IntersectWithPlane (const Base::Vector3f &rclBase, const Base::Vector3f &rclNormal, Base::Vector3f &rclP1, Base::Vector3f &rclP2) const { - // the triangle's corner points - const Base::Vector3f& v0 = _aclPoints[0]; - const Base::Vector3f& v1 = _aclPoints[1]; - const Base::Vector3f& v2 = _aclPoints[2]; + const float eps = 1e-06f; - // edge lengths - float len0 = (v0-v1).Length(); - float len1 = (v1-v2).Length(); - float len2 = (v2-v0).Length(); + // the triangle's corner points + const Base::Vector3f& v0 = _aclPoints[0]; + const Base::Vector3f& v1 = _aclPoints[1]; + const Base::Vector3f& v2 = _aclPoints[2]; - // Build up the line segments - Vector3 p0(0.5f*(v0.x+v1.x), 0.5f*(v0.y+v1.y), 0.5f*(v0.z+v1.z)); - Vector3 p1(0.5f*(v1.x+v2.x), 0.5f*(v1.y+v2.y), 0.5f*(v1.z+v2.z)); - Vector3 p2(0.5f*(v2.x+v0.x), 0.5f*(v2.y+v0.y), 0.5f*(v2.z+v0.z)); - - Vector3 d0(v1.x - v0.x, v1.y - v0.y, v1.z - v0.z); - d0.Normalize(); - Vector3 d1(v2.x - v1.x, v2.y - v1.y, v2.z - v1.z); - d1.Normalize(); - Vector3 d2(v0.x - v2.x, v0.y - v2.y, v0.z - v2.z); - d2.Normalize(); - - Segment3 akSeg0(p0, d0, len0/2.0f ); - Segment3 akSeg1(p1, d1, len1/2.0f); - Segment3 akSeg2(p2, d2, len2/2.0f); - - // Build up the plane - Vector3 p(rclBase.x, rclBase.y, rclBase.z); - Vector3 n(rclNormal.x, rclNormal.y, rclNormal.z); - Plane3 akPln(n, p); - - // Check for intersection with plane for each line segment - IntrSegment3Plane3 test0(akSeg0, akPln); - IntrSegment3Plane3 test1(akSeg1, akPln); - IntrSegment3Plane3 test2(akSeg2, akPln); - - Vector3 intr; - if ( test0.Find() ) - { - intr = p0 + test0.GetSegmentT() * d0; - rclP1.Set( intr[0], intr[1], intr[2]); - - if ( test1.Find() ) - { - intr = p1 + test1.GetSegmentT() * d1; - rclP2.Set( intr[0], intr[1], intr[2]); - return true; + // first check if a triangle's edge lies on the plane + float dist0 = fabs(v0.DistanceToPlane(rclBase, rclNormal)); + float dist1 = fabs(v1.DistanceToPlane(rclBase, rclNormal)); + float dist2 = fabs(v2.DistanceToPlane(rclBase, rclNormal)); + if (dist0 < eps && dist1 < eps) { + rclP1 = v0; + rclP2 = v1; + return true; } - else if ( test2.Find() ) - { - intr = p2 + test2.GetSegmentT() * d2; - rclP2.Set( intr[0], intr[1], intr[2]); - return true; + if (dist1 < eps && dist2 < eps) { + rclP1 = v1; + rclP2 = v2; + return true; } - } - else if ( test1.Find() ) - { - intr = p1 + test1.GetSegmentT() * d1; - rclP1.Set( intr[0], intr[1], intr[2]); - - if ( test2.Find() ) - { - intr = p2 + test2.GetSegmentT() * d2; - rclP2.Set( intr[0], intr[1], intr[2]); - return true; + if (dist2 < eps && dist0 < eps) { + rclP1 = v2; + rclP2 = v0; + return true; } - } - return false; + // edge lengths + float len0 = (v0-v1).Length(); + float len1 = (v1-v2).Length(); + float len2 = (v2-v0).Length(); + + // Build up the line segments + Vector3 p0(0.5f*(v0.x+v1.x), 0.5f*(v0.y+v1.y), 0.5f*(v0.z+v1.z)); + Vector3 p1(0.5f*(v1.x+v2.x), 0.5f*(v1.y+v2.y), 0.5f*(v1.z+v2.z)); + Vector3 p2(0.5f*(v2.x+v0.x), 0.5f*(v2.y+v0.y), 0.5f*(v2.z+v0.z)); + + Vector3 d0(v1.x - v0.x, v1.y - v0.y, v1.z - v0.z); + d0.Normalize(); + Vector3 d1(v2.x - v1.x, v2.y - v1.y, v2.z - v1.z); + d1.Normalize(); + Vector3 d2(v0.x - v2.x, v0.y - v2.y, v0.z - v2.z); + d2.Normalize(); + + Segment3 akSeg0(p0, d0, len0/2.0f ); + Segment3 akSeg1(p1, d1, len1/2.0f); + Segment3 akSeg2(p2, d2, len2/2.0f); + + // Build up the plane + Vector3 p(rclBase.x, rclBase.y, rclBase.z); + Vector3 n(rclNormal.x, rclNormal.y, rclNormal.z); + Plane3 akPln(n, p); + + // Check for intersection with plane for each line segment + IntrSegment3Plane3 test0(akSeg0, akPln); + IntrSegment3Plane3 test1(akSeg1, akPln); + IntrSegment3Plane3 test2(akSeg2, akPln); + + Vector3 intr; + + // now check if a triangle's corner lies on the plane + if (dist0 < eps) { + rclP1 = v0; + rclP2 = v0; + if (test1.Find()) { + intr = p1 + test1.GetSegmentT() * d1; + rclP2.Set(intr[0], intr[1], intr[2]); + } + return true; + } + else if (dist1 < eps) { + rclP1 = v1; + rclP2 = v1; + if (test2.Find()) { + intr = p2 + test2.GetSegmentT() * d2; + rclP2.Set(intr[0], intr[1], intr[2]); + } + return true; + } + else if (dist2 < eps) { + rclP1 = v2; + rclP2 = v2; + if (test0.Find()) { + intr = p0 + test0.GetSegmentT() * d0; + rclP2.Set(intr[0], intr[1], intr[2]); + } + return true; + } + + // check for arbitrary intersections + if (test0.Find()) { + intr = p0 + test0.GetSegmentT() * d0; + rclP1.Set( intr[0], intr[1], intr[2]); + + if (test1.Find()) { + intr = p1 + test1.GetSegmentT() * d1; + rclP2.Set( intr[0], intr[1], intr[2]); + return true; + } + else if (test2.Find()) { + intr = p2 + test2.GetSegmentT() * d2; + rclP2.Set( intr[0], intr[1], intr[2]); + return true; + } + } + else if (test1.Find()) { + intr = p1 + test1.GetSegmentT() * d1; + rclP1.Set( intr[0], intr[1], intr[2]); + + if (test2.Find()) { + intr = p2 + test2.GetSegmentT() * d2; + rclP2.Set( intr[0], intr[1], intr[2]); + return true; + } + } + + return false; } bool MeshGeomFacet::Foraminate (const Base::Vector3f &P, const Base::Vector3f &dir, Base::Vector3f &I, float fMaxAngle) const