From fe070903d3ce587b0e4f855b321ad33ad44746a1 Mon Sep 17 00:00:00 2001 From: wmayer Date: Tue, 19 Oct 2021 19:22:21 +0200 Subject: [PATCH] Mesh: add method to get intersection of edges --- src/Mod/Mesh/App/Core/Elements.cpp | 99 +++++++++++++++++++++++++++++- src/Mod/Mesh/App/Core/Elements.h | 16 ++++- 2 files changed, 109 insertions(+), 6 deletions(-) diff --git a/src/Mod/Mesh/App/Core/Elements.cpp b/src/Mod/Mesh/App/Core/Elements.cpp index 488effb791..b3d41952bc 100644 --- a/src/Mod/Mesh/App/Core/Elements.cpp +++ b/src/Mod/Mesh/App/Core/Elements.cpp @@ -280,6 +280,61 @@ bool MeshGeomEdge::IntersectWithLine (const Base::Vector3f &rclPt, return dist2 + dist3 <= dist1 + eps; } +bool MeshGeomEdge::IntersectWithEdge (const MeshGeomEdge &edge, Base::Vector3f &res) const +{ + const float eps = 1e-06f; + Base::Vector3f p(_aclPoints[0]); + Base::Vector3f r(_aclPoints[1] - _aclPoints[0]); + Base::Vector3f q(edge._aclPoints[0]); + Base::Vector3f s(edge._aclPoints[1] - edge._aclPoints[0]); + Base::Vector3f n = r.Cross(s); + Base::Vector3f d = q - p; + + // lines are collinear or parallel + if (n.IsNull()) { + if (d.Cross(r).IsNull()) { + // Collinear + if (IsProjectionPointOf(edge._aclPoints[0])) { + res = edge._aclPoints[0]; + return true; + } + if (IsProjectionPointOf(edge._aclPoints[1])) { + res = edge._aclPoints[1]; + return true; + } + + return false; + } + else { + // Parallel + return false; + } + } + else { + // Get the distance of q to the plane defined by p and n + float distance = q.DistanceToPlane(p, n); + + // lines are warped + if (fabs(distance) > eps) + return false; + + float t = d.Cross(s).Dot(n) / n.Sqr(); + float u = d.Cross(r).Dot(n) / n.Sqr(); + + auto is_in_range = [](float v) { + return v >= 0.0f && v <= 1.0f; + }; + + if (is_in_range(t) && is_in_range(u)) { + res = p + t * r; + res = q + u * s; + return true; + } + + return false; + } +} + bool MeshGeomEdge::IntersectWithPlane (const Base::Vector3f &rclPt, const Base::Vector3f &rclDir, Base::Vector3f &rclRes) const @@ -370,6 +425,14 @@ bool MeshGeomEdge::IsPointOf (const Base::Vector3f &rclPoint, float fDistance) c return Base::Distance(ptEdge, rclPoint) <= fDistance; } +bool MeshGeomEdge::IsProjectionPointOf(const Base::Vector3f& point) const +{ + Base::Vector3f fromStartToPoint = point - _aclPoints[0]; + Base::Vector3f fromPointToEnd = _aclPoints[1] - point; + float dot = fromStartToPoint * fromPointToEnd; + return dot >= 0.0f; +} + // ----------------------------------------------------------------- MeshGeomFacet::MeshGeomFacet () @@ -1010,9 +1073,31 @@ int MeshGeomFacet::IntersectWithFacet (const MeshGeomFacet& rclFacet, // Note: tri_tri_intersect_with_isection() does not return line of // intersection when triangles are coplanar. See tritritest.h:18 and 658. - // So rclPt* may be garbage values and we cannot continue. - if (coplanar) - return 2; // equivalent to rclPt0 != rclPt1 + if (coplanar) { + // Since tri_tri_intersect_with_isection may return garbage values try to get + // sensible values with edge/edge intersections + std::vector intersections; + for (short i=0; i<3; i++) { + MeshGeomEdge edge1 = GetEdge(i); + for (short j=0; j<3; j++) { + MeshGeomEdge edge2 = rclFacet.GetEdge(j); + Base::Vector3f point; + if (edge1.IntersectWithEdge(edge2, point)) { + intersections.push_back(point); + } + } + } + + if (intersections.size() == 2) { + rclPt0 = intersections[0]; + rclPt1 = intersections[1]; + } + else if (intersections.size() == 1) { + rclPt0 = intersections[0]; + rclPt1 = intersections[0]; + } + return 2; + } // With extremely acute-angled triangles it may happen that the algorithm // claims an intersection but the intersection points are far outside the @@ -1281,6 +1366,14 @@ void MeshGeomFacet::NearestEdgeToPoint(const Base::Vector3f& rclPt, float& fDist } } +MeshGeomEdge MeshGeomFacet::GetEdge(short side) const +{ + MeshGeomEdge edge; + edge._aclPoints[0] = this->_aclPoints[side %3]; + edge._aclPoints[1] = this->_aclPoints[(side+1)%3]; + return edge; +} + float MeshGeomFacet::VolumeOfPrism (const MeshGeomFacet& rclF1) const { Base::Vector3f P1 = this->_aclPoints[0]; diff --git a/src/Mod/Mesh/App/Core/Elements.h b/src/Mod/Mesh/App/Core/Elements.h index b2c1f62372..db9fa791ca 100644 --- a/src/Mod/Mesh/App/Core/Elements.h +++ b/src/Mod/Mesh/App/Core/Elements.h @@ -168,6 +168,10 @@ public: * with the edge. The intersection must be inside the edge. If there is no intersection false is returned. */ bool IntersectWithLine (const Base::Vector3f &rclPt, const Base::Vector3f &rclDir, Base::Vector3f &rclRes) const; + /** Calculates the intersection point of an edge with this edge. + * The intersection must be inside both edges. If there is no intersection false is returned. + */ + bool IntersectWithEdge (const MeshGeomEdge &edge, Base::Vector3f &res) const; /** Calculates the intersection point of the plane defined by the base \a rclPt and the direction \a rclDir * with the edge. The intersection must be inside the edge. If there is no intersection false is returned. */ @@ -185,11 +189,15 @@ public: void ClosestPointsToLine(const Base::Vector3f &linePt, const Base::Vector3f &lineDir, Base::Vector3f& rclPnt1, Base::Vector3f& rclPnt2) const; /** - * Checks if the point is part of the facet. A point is regarded as part - * of a facet if the distance is lower \a fDistance and the projected point - * in the facet normal direction is inside the triangle. + * Checks if the point is part of the edge. A point is regarded as part + * of an edge if the distance is lower than \a fDistance to the projected point + * of \a rclPoint on the edge. */ bool IsPointOf (const Base::Vector3f &rclPoint, float fDistance) const; + /** + * Checks if the projection point of \a point lies on the edge. + */ + bool IsProjectionPointOf(const Base::Vector3f& point) const; public: Base::Vector3f _aclPoints[2]; /**< Corner points */ @@ -506,6 +514,8 @@ public: unsigned short NearestEdgeToPoint(const Base::Vector3f& rclPt) const; /** Returns the edge number \a side of the facet and the distance to the edge that is nearest to the point \a rclPt. */ void NearestEdgeToPoint(const Base::Vector3f& rclPt, float& fDistance, unsigned short& side) const; + /** Returns the edge for \a side. */ + MeshGeomEdge GetEdge(short side) const; /** The center and radius of the circum circle define a sphere in 3D. If the point \a rP is part of this sphere true is * returned, otherwise false. */