From 7d511b2bd90d64a48cd726d3791e9b1a11bd7ed8 Mon Sep 17 00:00:00 2001 From: wandererfan Date: Fri, 21 Mar 2025 08:45:48 -0400 Subject: [PATCH] [TD]use Ptolemy's Theorem for bspline to circle conversion --- src/Mod/TechDraw/App/Geometry.cpp | 209 ++++++++++++++++++------------ src/Mod/TechDraw/App/Geometry.h | 34 ++--- 2 files changed, 144 insertions(+), 99 deletions(-) diff --git a/src/Mod/TechDraw/App/Geometry.cpp b/src/Mod/TechDraw/App/Geometry.cpp index 8519cf623f..0cb09d530b 100644 --- a/src/Mod/TechDraw/App/Geometry.cpp +++ b/src/Mod/TechDraw/App/Geometry.cpp @@ -1493,87 +1493,124 @@ TopoDS_Edge GeometryUtils::edgeFromCircleArc(TechDraw::AOCPtr c) } //used by DVDim for approximate dims -bool GeometryUtils::isCircle(TopoDS_Edge occEdge) +bool GeometryUtils::isCircle(const TopoDS_Edge& occEdge) { - double radius; + double radius{0}; Base::Vector3d center; bool isArc = false; return GeometryUtils::getCircleParms(occEdge, radius, center, isArc); } //! tries to interpret a B-spline edge as a circle. Used by DVDim for approximate dimensions. -//! calculates the curvature of the spline at a number of places and measures the deviation from the average -//! a true circle has constant curvature and would have no deviation from the average. -bool GeometryUtils::getCircleParms(TopoDS_Edge occEdge, double& radius, Base::Vector3d& center, bool& isArc) +//! calculates the radius and center of circles using groups of 4 points on the b-spline. if the +//! groups of 4 points all lie on a circle, we use that circle to get the radius and center. +bool GeometryUtils::getCircleParms(const TopoDS_Edge& occEdge, double& radius, Base::Vector3d& center, bool& isArc) { - int testCount = 5; - double curveLimit = EWTOLERANCE; - BRepAdaptor_Curve c(occEdge); - Handle(Geom_BSplineCurve) spline = c.BSpline(); - double f, l; - f = c.FirstParameter(); - l = c.LastParameter(); - double parmRange = fabs(l - f); - double parmStep = parmRange/testCount; - std::vector curvatures; - std::vector centers; - gp_Pnt curveCenter; - double sumCurvature = 0; - Base::Vector3d sumCenter, valueAt; - try { - GeomLProp_CLProps prop(spline, f, 3, Precision::Confusion()); + constexpr int PointCount{8}; // number of points on the edge to examine (>= 8) + constexpr int TestCount{3}; // number of candidate circles to test - // check only the interior points of the edge - for (int i = 1; i < (testCount - 1); i++) { - prop.SetParameter(parmStep * i); - curvatures.push_back(prop.Curvature()); - sumCurvature += prop.Curvature(); - prop.CentreOfCurvature(curveCenter); - centers.push_back(curveCenter); - sumCenter += Base::convertTo(curveCenter); - } + BRepAdaptor_Curve curveAdapt(occEdge); + double firstParam = curveAdapt.FirstParameter(); + auto firstPoint = Base::convertTo(curveAdapt.Value(firstParam)); + double lastParam = curveAdapt.LastParameter(); + auto lastPoint = Base::convertTo(curveAdapt.Value(lastParam)); + double parmRange = fabs(lastParam - firstParam); + double parmStep = parmRange / PointCount; + + std::vector pointsOnCurve; + for (size_t iPoint = 0; iPoint < PointCount; iPoint++) { + auto iPointMath = static_cast(iPoint); + auto newpoint = curveAdapt.Value(firstParam + iPointMath * parmStep); + pointsOnCurve.push_back(Base::convertTo(newpoint)); } - catch (Standard_Failure&) { + + auto edgeLong = edgeLength(occEdge); + constexpr double LimitFactor{0.001}; // 0.1% not sure about this value + double tolerance = edgeLong * LimitFactor; + + isArc = true; + if (firstPoint.IsEqual(lastPoint, tolerance)) { + isArc = false; + } + + int passCount{0}; + int firstIndex{0}; + for (int iTest = 0; iTest < TestCount; iTest++) { + firstIndex++; + auto A = pointsOnCurve.at(firstIndex); + auto B = pointsOnCurve.at(firstIndex + 1); + auto C = pointsOnCurve.at(firstIndex + 2); + auto D = pointsOnCurve.at(firstIndex + 3); + if (pointsAreOnCircle(A, B, C, D, tolerance)) { + passCount++; + } + } + + if (passCount != TestCount) { + // at least 1 test failed. return false; } - Base::Vector3d avgCenter = sumCenter/ centers.size(); - double avgCurve = sumCurvature/ centers.size(); - double errorCurve = 0; - // sum the errors in curvature - for (auto& cv: curvatures) { - errorCurve += avgCurve - cv; + // each group of 4 points lies on a circle. Since the groups of 4 overlap, all the points lie + // on the same circle. https://en.wikipedia.org/wiki/Ptolemy%27s_theorem and + // https://math.stackexchange.com/questions/3130053/how-to-check-if-a-set-of-points-in-cartesian-space-could-lie-on-the-circumferenc + // so we can use any three points to make our circle. + + auto gPoint0 = Base::convertTo(pointsOnCurve.at(1)); + auto gPoint1 = Base::convertTo(pointsOnCurve.at(3)); + auto gPoint2 = Base::convertTo(pointsOnCurve.at(5)); //NOLINT readability-magic-numbers + try { + GC_MakeCircle mkCircle(gPoint0, gPoint1, gPoint2); + if (!mkCircle.IsDone()) { + return false; + } + + const Handle(Geom_Circle) circleFromParms = mkCircle.Value(); + radius = circleFromParms->Circ().Radius(); + center = Base::convertTo(circleFromParms->Circ().Location()); + return true; + } + catch (Standard_Failure& err) { + Base::Console().Message("Geo::getCircleParms - failed to make a circle\n"); } - double errorCenter{0}; - for (auto& observe : centers) { - auto error = (Base::convertTo(observe)- avgCenter).Length(); - errorCenter += error; + return false; +} + + +//! returns true if the A, B, C and D all lie on the same circle according to Ptolemy's theorem +//! we can skip the test for same plane, since the points are all on the XY plane(?not true +//! for 3d dims?). +bool GeometryUtils::pointsAreOnCircle(Base::Vector3d A, + Base::Vector3d B, + Base::Vector3d C, + Base::Vector3d D, + double tolerance) +{ + auto AB = (B-A).Length(); + auto AC = (C-A).Length(); + auto AD = (D-A).Length(); + auto BC = (C-B).Length(); + auto BD = (D-B).Length(); + auto CD = (D-C).Length(); + + auto pieceLength = AB + BC + CD; + auto wholeLength = AD; + if (DU::fpCompare(pieceLength, wholeLength, tolerance)) { + // these points are colinear + return false; } - // calculate average error in curvature. we are only interested in the magnitude of the error - errorCurve = fabs(errorCurve / curvatures.size()); - // calculate the average error in center of curvature - errorCenter = errorCenter / curvatures.size(); - auto edgeLong = edgeLength(occEdge); - double centerLimit = edgeLong * 0.01; - - isArc = !c.IsClosed(); - bool isCircle(false); - if ( errorCurve <= curveLimit && - errorCenter <= centerLimit) { - isCircle = true; - radius = 1.0/avgCurve; - center = avgCenter; - } - - return isCircle; + bool eq1 = DU::fpCompare(AB*CD + AC*BD, AD*BC, tolerance); + bool eq2 = DU::fpCompare(AB*CD + AD*BC, AC*BD, tolerance); + bool eq3 = DU::fpCompare(AC*BD + AD*BC, AB*CD, tolerance); + return eq1 || eq2 || eq3; } //! make a circle or arc of circle Edge from BSpline Edge // Note that the input edge has been inverted by GeometryObject, so +Y points down. -TopoDS_Edge GeometryUtils::asCircle(TopoDS_Edge splineEdge, bool& arc) +TopoDS_Edge GeometryUtils::asCircle(const TopoDS_Edge& splineEdge, bool& arc) { double radius{0}; Base::Vector3d center; @@ -1582,11 +1619,16 @@ TopoDS_Edge GeometryUtils::asCircle(TopoDS_Edge splineEdge, bool& arc) if (!canMakeCircle) { throw Base::RuntimeError("GU::asCircle received non-circular edge!"); } + arc = isArc; gp_Pnt gCenter = Base::convertTo(center); gp_Dir gNormal{0, 0, 1}; Handle(Geom_Circle) circleFromParms = GC_MakeCircle(gCenter, gNormal, radius); + if (!isArc) { + return BRepBuilderAPI_MakeEdge(circleFromParms); + } + // find the ends of the edge from the underlying curve BRepAdaptor_Curve curveAdapt(splineEdge); double firstParam = curveAdapt.FirstParameter(); @@ -1594,40 +1636,42 @@ TopoDS_Edge GeometryUtils::asCircle(TopoDS_Edge splineEdge, bool& arc) gp_Pnt startPoint = curveAdapt.Value(firstParam); gp_Pnt endPoint = curveAdapt.Value(lastParam); - if (startPoint.IsEqual(endPoint, EWTOLERANCE)) { //more reliable than IsClosed flag - arc = false; - return BRepBuilderAPI_MakeEdge(circleFromParms); - } - - arc = true; double midRange = (lastParam + firstParam) / 2; gp_Pnt midPoint = curveAdapt.Value(midRange); + // this should be using circleFromParms as a base instead of points on the original spline. + // could be problems with very small errors?? other versions of GC_MakeArcOfCircle have + // poorly explained parameters. GC_MakeArcOfCircle mkArc(startPoint, midPoint, endPoint); auto circleArc = mkArc.Value(); - + if (!mkArc.IsDone()) { + throw Base::RuntimeError("GU::asCircle failed to create arc"); + } return BRepBuilderAPI_MakeEdge(circleArc); } - -bool GeometryUtils::isLine(TopoDS_Edge occEdge) +bool GeometryUtils::isLine(const TopoDS_Edge& occEdge) { - BRepAdaptor_Curve c(occEdge); + BRepAdaptor_Curve adapt(occEdge); - Handle(Geom_BSplineCurve) spline = c.BSpline(); - double f = c.FirstParameter(); - double l = c.LastParameter(); - gp_Pnt s = c.Value(f); - gp_Pnt e = c.Value(l); + Handle(Geom_BSplineCurve) spline = adapt.BSpline(); + double firstParm = adapt.FirstParameter(); + double lastParm = adapt.LastParameter(); + auto startPoint = Base::convertTo(adapt.Value(firstParm)); + auto endPoint = Base::convertTo(adapt.Value(lastParm)); + auto edgeLong = edgeLength(occEdge); - bool samePnt = s.IsEqual(e, FLT_EPSILON); - if (samePnt) { + constexpr double LimitFactor{0.001}; // 0.1% not sure about this value + double tolerance = edgeLong * LimitFactor; + if (startPoint.IsEqual(endPoint, tolerance)) { + // either not a line or a zero length line? return false; } - Base::Vector3d vs = Base::convertTo(s); - Base::Vector3d ve = Base::convertTo(e); - double endLength = (vs - ve).Length(); + // in a line the sum of the lengths of the segments should equal the distance + // from start to end + double endPointLength = (endPoint - startPoint).Length(); + int low = 0; int high = spline->NbPoles() - 1; TColgp_Array1OfPnt poles(low, high); @@ -1641,15 +1685,12 @@ bool GeometryUtils::isLine(TopoDS_Edge occEdge) lenTotal += (v2-v1).Length(); } - if (DrawUtil::fpCompare(lenTotal, endLength)) { - return true; - } - return false; + return DrawUtil::fpCompare(lenTotal, endPointLength, tolerance); } //! make a line Edge from B-spline Edge -TopoDS_Edge GeometryUtils::asLine(TopoDS_Edge occEdge) +TopoDS_Edge GeometryUtils::asLine(const TopoDS_Edge& occEdge) { BRepAdaptor_Curve c(occEdge); diff --git a/src/Mod/TechDraw/App/Geometry.h b/src/Mod/TechDraw/App/Geometry.h index 8373881d12..81a316a1e4 100644 --- a/src/Mod/TechDraw/App/Geometry.h +++ b/src/Mod/TechDraw/App/Geometry.h @@ -315,7 +315,6 @@ class TechDrawExport BSpline: public BaseGeom bool isLine(); bool isCircle(); TopoDS_Edge asCircle(bool& isArc); -// void getCircleParms(bool& isCircle, double& radius, Base::Vector3d& center, bool& isArc); bool intersectsArc(Base::Vector3d p1, Base::Vector3d p2); std::vector segments; }; @@ -382,23 +381,23 @@ class TechDrawExport Vertex : public TechDraw::Tag Base::Vector3d point() const { return Base::Vector3d(pnt.x, pnt.y, 0.0); } void point(Base::Vector3d v){ pnt = Base::Vector3d(v.x, v.y); } - double x() {return pnt.x;} - double y() {return pnt.y;} + double x() const {return pnt.x;} + double y() const {return pnt.y;} // attribute setters and getters - bool getHlrVisible() { return hlrVisible; } + bool getHlrVisible() const { return hlrVisible; } void setHlrVisible(bool state) { hlrVisible = state; } - int getRef3d() { return ref3D; } + int getRef3d() const { return ref3D; } void setRef3d(int ref) { ref3D = ref; } TopoDS_Vertex getOCCVertex() { return occVertex; } - void setOCCVertex(TopoDS_Vertex newVertex) { occVertex = newVertex; } - bool getCosmetic() { return cosmetic; } + void setOCCVertex(const TopoDS_Vertex& newVertex) { occVertex = newVertex; } + bool getCosmetic() const { return cosmetic; } void setCosmetic (bool state) { cosmetic = state; } std::string getCosmeticTag() { return cosmeticTag; } - void setCosmeticTag(std::string t) { cosmeticTag = t; } - bool isCenter() {return m_center;} + void setCosmeticTag(const std::string& t) { cosmeticTag = t; } + bool isCenter() const {return m_center;} void isCenter(bool state) { m_center = state; } - bool isReference() { return m_reference; } + bool isReference() const { return m_reference; } void isReference(bool state) { m_reference = state; } Part::TopoShape asTopoShape(double scale = 1.0); @@ -447,17 +446,22 @@ class TechDrawExport GeometryUtils static TopoDS_Edge edgeFromCircle(TechDraw::CirclePtr c); static TopoDS_Edge edgeFromCircleArc(TechDraw::AOCPtr c); - static bool isCircle(TopoDS_Edge occEdge); - static bool getCircleParms(TopoDS_Edge occEdge, double& radius, Base::Vector3d& center, bool& isArc); - static TopoDS_Edge asCircle(TopoDS_Edge splineEdge, bool& arc); - static bool isLine(TopoDS_Edge occEdge); - static TopoDS_Edge asLine(TopoDS_Edge occEdge); + static bool isCircle(const TopoDS_Edge& occEdge); + static bool getCircleParms(const TopoDS_Edge& occEdge, double& radius, Base::Vector3d& center, bool& isArc); + static TopoDS_Edge asCircle(const TopoDS_Edge& splineEdge, bool& arc); + static bool isLine(const TopoDS_Edge& occEdge); + static TopoDS_Edge asLine(const TopoDS_Edge& occEdge); static double edgeLength(TopoDS_Edge occEdge); static TopoDS_Face makePerforatedFace(FacePtr bigCheese, const std::vector& holesAll); static std::vector findHolesInFace(const DrawViewPart* dvp, const std::string& bigCheeseSubRef); + static bool pointsAreOnCircle(Base::Vector3d A, + Base::Vector3d B, + Base::Vector3d C, + Base::Vector3d D, + double tolerance); };