// Created on: 1994-07-22 // Created by: Remi LEQUETTE // Copyright (c) 1994-1999 Matra Datavision // Copyright (c) 1999-2014 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 #include #include #include #include #include #include #include #include #include #include #include #include #include #include //================================================================================================= static Standard_Real Controle(const TColgp_SequenceOfPnt& thePoints, const Handle(Geom_Plane)& thePlane) { Standard_Real dfMaxDist = 0.; Standard_Real a, b, c, d, dist; Standard_Integer ii; thePlane->Coefficients(a, b, c, d); for (ii = 1; ii <= thePoints.Length(); ii++) { const gp_XYZ& xyz = thePoints(ii).XYZ(); dist = Abs(a * xyz.X() + b * xyz.Y() + c * xyz.Z() + d); if (dist > dfMaxDist) dfMaxDist = dist; } return dfMaxDist; } //======================================================================= // function : Is2DConnected // purpose : Return true if the last vertex of theEdge1 coincides with // the first vertex of theEdge2 in parametric space of theFace //======================================================================= inline static Standard_Boolean Is2DConnected(const TopoDS_Edge& theEdge1, const TopoDS_Edge& theEdge2, const Handle(Geom_Surface)& theSurface, const TopLoc_Location& theLocation) { Standard_Real f, l; // TopLoc_Location aLoc; Handle(Geom2d_Curve) aCurve; gp_Pnt2d p1, p2; // get 2D points aCurve = BRep_Tool::CurveOnSurface(theEdge1, theSurface, theLocation, f, l); p1 = aCurve->Value(theEdge1.Orientation() == TopAbs_FORWARD ? l : f); aCurve = BRep_Tool::CurveOnSurface(theEdge2, theSurface, theLocation, f, l); p2 = aCurve->Value(theEdge2.Orientation() == TopAbs_FORWARD ? f : l); // compare 2D points GeomAdaptor_Surface aSurface(theSurface); TopoDS_Vertex aV = TopExp::FirstVertex(theEdge2, Standard_True); Standard_Real tol3D = BRep_Tool::Tolerance(aV); Standard_Real tol2D = aSurface.UResolution(tol3D) + aSurface.VResolution(tol3D); Standard_Real dist2 = p1.SquareDistance(p2); return dist2 < tol2D * tol2D; } //======================================================================= // function : Is2DClosed // purpose : Return true if edges of theShape form a closed wire in // parametric space of theSurface //======================================================================= static Standard_Boolean Is2DClosed(const TopoDS_Shape& theShape, const Handle(Geom_Surface)& theSurface, const TopLoc_Location& theLocation) { try { OCC_CATCH_SIGNALS // get a wire theShape TopExp_Explorer aWireExp(theShape, TopAbs_WIRE); if (!aWireExp.More()) { return Standard_False; } TopoDS_Wire aWire = TopoDS::Wire(aWireExp.Current()); // a tmp face TopoDS_Face aTmpFace = BRepLib_MakeFace(theSurface, Precision::PConfusion()); // check topological closeness using wire explorer, if the wire is not closed // the 1st and the last vertices of wire are different BRepTools_WireExplorer aWireExplorer(aWire, aTmpFace); if (!aWireExplorer.More()) { return Standard_False; } // remember the 1st and the last edges of aWire TopoDS_Edge aFisrtEdge = aWireExplorer.Current(), aLastEdge = aFisrtEdge; // check if edges connected topologically (that is assured by BRepTools_WireExplorer) // are connected in 2D TopoDS_Edge aPrevEdge = aFisrtEdge; for (aWireExplorer.Next(); aWireExplorer.More(); aWireExplorer.Next()) { aLastEdge = aWireExplorer.Current(); if (!Is2DConnected(aPrevEdge, aLastEdge, theSurface, theLocation)) { return false; } aPrevEdge = aLastEdge; } // wire is closed if ( 1st vertex of aFisrtEdge ) == // ( last vertex of aLastEdge ) in 2D TopoDS_Vertex aV1 = TopExp::FirstVertex(aFisrtEdge, Standard_True); TopoDS_Vertex aV2 = TopExp::LastVertex(aLastEdge, Standard_True); return (aV1.IsSame(aV2) && Is2DConnected(aLastEdge, aFisrtEdge, theSurface, theLocation)); } catch (Standard_Failure const&) { return Standard_False; } } //================================================================================================= BRepLib_FindSurface::BRepLib_FindSurface() : myTolerance(0.0), myTolReached(0.0), isExisted(Standard_False) { } //================================================================================================= BRepLib_FindSurface::BRepLib_FindSurface(const TopoDS_Shape& S, const Standard_Real Tol, const Standard_Boolean OnlyPlane, const Standard_Boolean OnlyClosed) { Init(S, Tol, OnlyPlane, OnlyClosed); } namespace { static void fillParams(const TColStd_Array1OfReal& theKnots, Standard_Integer theDegree, Standard_Real theParMin, Standard_Real theParMax, NCollection_Vector& theParams) { Standard_Real aPrevPar = theParMin; theParams.Append(aPrevPar); Standard_Integer aNbP = Max(theDegree, 1); for (Standard_Integer i = 1; (i < theKnots.Length()) && (theKnots(i) < (theParMax - Precision::PConfusion())); ++i) { if (theKnots(i + 1) < theParMin + Precision::PConfusion()) continue; Standard_Real aStep = (theKnots(i + 1) - theKnots(i)) / aNbP; for (Standard_Integer k = 1; k <= aNbP; ++k) { Standard_Real aPar = theKnots(i) + k * aStep; if (aPar > theParMax - Precision::PConfusion()) break; if (aPar > aPrevPar + Precision::PConfusion()) { theParams.Append(aPar); aPrevPar = aPar; } } } theParams.Append(theParMax); } static void fillPoints(const BRepAdaptor_Curve& theCurve, const NCollection_Vector& theParams, TColgp_SequenceOfPnt& thePoints, TColStd_SequenceOfReal& theWeights) { Standard_Real aDistPrev = 0., aDistNext; gp_Pnt aPPrev(theCurve.Value(theParams(0))), aPNext; for (Standard_Integer iP = 1; iP <= theParams.Length(); ++iP) { if (iP < theParams.Length()) { Standard_Real aParam = theParams(iP); aPNext = theCurve.Value(aParam); aDistNext = aPPrev.Distance(aPNext); } else aDistNext = 0.0; thePoints.Append(aPPrev); theWeights.Append(aDistPrev + aDistNext); aDistPrev = aDistNext; aPPrev = aPNext; } } } // namespace //================================================================================================= void BRepLib_FindSurface::Init(const TopoDS_Shape& S, const Standard_Real Tol, const Standard_Boolean OnlyPlane, const Standard_Boolean OnlyClosed) { myTolerance = Tol; myTolReached = 0.; isExisted = Standard_False; myLocation.Identity(); mySurface.Nullify(); // compute the tolerance TopExp_Explorer ex; for (ex.Init(S, TopAbs_EDGE); ex.More(); ex.Next()) { Standard_Real t = BRep_Tool::Tolerance(TopoDS::Edge(ex.Current())); if (t > myTolerance) myTolerance = t; } // search an existing surface ex.Init(S, TopAbs_EDGE); if (!ex.More()) return; // no edges .... TopoDS_Edge E = TopoDS::Edge(ex.Current()); Standard_Real f, l, ff, ll; Handle(Geom2d_Curve) PC, aPPC; Handle(Geom_Surface) SS; TopLoc_Location L; Standard_Integer i = 0, j; // iterate on the surfaces of the first edge for (;;) { i++; BRep_Tool::CurveOnSurface(E, PC, mySurface, myLocation, f, l, i); if (mySurface.IsNull()) { break; } // check the other edges for (ex.Init(S, TopAbs_EDGE); ex.More(); ex.Next()) { if (!E.IsSame(ex.Current())) { j = 0; for (;;) { j++; BRep_Tool::CurveOnSurface(TopoDS::Edge(ex.Current()), aPPC, SS, L, ff, ll, j); if (SS.IsNull()) { break; } if ((SS == mySurface) && (L.IsEqual(myLocation))) { break; } SS.Nullify(); } if (SS.IsNull()) { mySurface.Nullify(); break; } } } // if OnlyPlane, eval if mySurface is a plane. if (OnlyPlane && !mySurface.IsNull()) { if (mySurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) mySurface = Handle(Geom_RectangularTrimmedSurface)::DownCast(mySurface)->BasisSurface(); mySurface = Handle(Geom_Plane)::DownCast(mySurface); } if (!mySurface.IsNull()) // if S is e.g. the bottom face of a cylinder, mySurface can be the // lateral (cylindrical) face of the cylinder; reject an improper mySurface if (!OnlyClosed || Is2DClosed(S, mySurface, myLocation)) break; } if (!mySurface.IsNull()) { isExisted = Standard_True; return; } // // no existing surface, search a plane // 07/02/02 akm vvv : (OCC157) changed algorithm // 1. Collect the points along all edges of the shape // For each point calculate the WEIGHT = sum of // distances from neighboring points (_only_ same edge) // 2. Minimizing the weighed sum of squared deviations // compute coefficients of the sought plane. TColgp_SequenceOfPnt aPoints; TColStd_SequenceOfReal aWeight; // ======================= Step #1 for (ex.Init(S, TopAbs_EDGE); ex.More(); ex.Next()) { BRepAdaptor_Curve c(TopoDS::Edge(ex.Current())); Standard_Real dfUf = c.FirstParameter(); Standard_Real dfUl = c.LastParameter(); if (IsEqual(dfUf, dfUl)) { // Degenerate continue; } Standard_Integer iNbPoints = 0; // Fill the parameters of the sampling points NCollection_Vector aParams; switch (c.GetType()) { case GeomAbs_BezierCurve: { Handle(Geom_BezierCurve) GC = c.Bezier(); TColStd_Array1OfReal aKnots(1, 2); aKnots.SetValue(1, GC->FirstParameter()); aKnots.SetValue(2, GC->LastParameter()); fillParams(aKnots, GC->Degree(), dfUf, dfUl, aParams); break; } case GeomAbs_BSplineCurve: { Handle(Geom_BSplineCurve) GC = c.BSpline(); fillParams(GC->Knots(), GC->Degree(), dfUf, dfUl, aParams); break; } case GeomAbs_Line: { // Two points on a straight segment aParams.Append(dfUf); aParams.Append(dfUl); break; } case GeomAbs_Circle: case GeomAbs_Ellipse: case GeomAbs_Hyperbola: case GeomAbs_Parabola: // Four points on other analytical curves iNbPoints = 4; Standard_FALLTHROUGH default: { // Put some points on other curves if (iNbPoints == 0) iNbPoints = 15 + c.NbIntervals(GeomAbs_C3); TColStd_Array1OfReal aBounds(1, 2); aBounds.SetValue(1, dfUf); aBounds.SetValue(2, dfUl); fillParams(aBounds, iNbPoints - 1, dfUf, dfUl, aParams); } } // Add the points with weights to the sequences fillPoints(c, aParams, aPoints, aWeight); } if (aPoints.Length() < 3) { return; } // ======================= Step #2 myLocation.Identity(); Standard_Integer iPoint; math_Matrix aMat(1, 3, 1, 3, 0.); math_Vector aVec(1, 3, 0.); // Find the barycenter and normalize weights Standard_Real dfMaxWeight = 0.; gp_XYZ aBaryCenter(0., 0., 0.); Standard_Real dfSumWeight = 0.; for (iPoint = 1; iPoint <= aPoints.Length(); iPoint++) { Standard_Real dfW = aWeight(iPoint); aBaryCenter += dfW * aPoints(iPoint).XYZ(); dfSumWeight += dfW; if (dfW > dfMaxWeight) { dfMaxWeight = dfW; } } aBaryCenter /= dfSumWeight; // Fill the matrix and the right vector for (iPoint = 1; iPoint <= aPoints.Length(); iPoint++) { gp_XYZ p = aPoints(iPoint).XYZ() - aBaryCenter; Standard_Real w = aWeight(iPoint) / dfMaxWeight; aMat(1, 1) += w * p.X() * p.X(); aMat(1, 2) += w * p.X() * p.Y(); aMat(1, 3) += w * p.X() * p.Z(); // aMat(2, 2) += w * p.Y() * p.Y(); aMat(2, 3) += w * p.Y() * p.Z(); // aMat(3, 3) += w * p.Z() * p.Z(); } aMat(2, 1) = aMat(1, 2); aMat(3, 1) = aMat(1, 3); aMat(3, 2) = aMat(2, 3); // math_Jacobi anEignval(aMat); math_Vector anEVals(1, 3); Standard_Boolean isSolved = anEignval.IsDone(); Standard_Integer isol = 0; if (isSolved) { anEVals = anEignval.Values(); // We need vector with eigenvalue ~ 0. Standard_Real anEMin = RealLast(); Standard_Real anEMax = -anEMin; for (i = 1; i <= 3; ++i) { Standard_Real anE = Abs(anEVals(i)); if (anEMin > anE) { anEMin = anE; isol = i; } if (anEMax < anE) { anEMax = anE; } } if (isol == 0) { isSolved = Standard_False; } else { Standard_Real eps = Epsilon(anEMax); if (anEMin <= eps) { anEignval.Vector(isol, aVec); } else { // try using vector product of other axes Standard_Integer ind[2] = {0, 0}; for (i = 1; i <= 3; ++i) { if (i == isol) { continue; } if (ind[0] == 0) { ind[0] = i; continue; } if (ind[1] == 0) { ind[1] = i; continue; } } math_Vector aVec1(1, 3, 0.), aVec2(1, 3, 0.); anEignval.Vector(ind[0], aVec1); anEignval.Vector(ind[1], aVec2); gp_Vec aV1(aVec1(1), aVec1(2), aVec1(3)); gp_Vec aV2(aVec2(1), aVec2(2), aVec2(3)); gp_Vec aN = aV1 ^ aV2; aVec(1) = aN.X(); aVec(2) = aN.Y(); aVec(3) = aN.Z(); } if (aVec.Norm2() < gp::Resolution()) { isSolved = Standard_False; } } } if (!isSolved) return; // Removing very small values Standard_Real aMaxV = Max(Abs(aVec(1)), Max(Abs(aVec(2)), Abs(aVec(3)))); Standard_Real eps = Epsilon(aMaxV); for (i = 1; i <= 3; ++i) { if (Abs(aVec(i)) <= eps) aVec(i) = 0.; } gp_Vec aN(aVec(1), aVec(2), aVec(3)); Handle(Geom_Plane) aPlane = new Geom_Plane(aBaryCenter, aN); myTolReached = Controle(aPoints, aPlane); const Standard_Real aWeakness = 5.0; if (myTolReached <= myTolerance || (Tol < 0 && myTolReached < myTolerance * aWeakness)) { mySurface = aPlane; // If S is wire, try to orient surface according to orientation of wire. if (S.ShapeType() == TopAbs_WIRE && S.Closed()) { TopoDS_Wire aW = TopoDS::Wire(S); TopoDS_Face aTmpFace = BRepLib_MakeFace(mySurface, Precision::Confusion()); BRep_Builder BB; BB.Add(aTmpFace, aW); BRepTopAdaptor_FClass2d FClass(aTmpFace, 0.); if (FClass.PerformInfinitePoint() == TopAbs_IN) { gp_Dir aNorm = aPlane->Position().Direction(); aNorm.Reverse(); mySurface = new Geom_Plane(aPlane->Position().Location(), aNorm); } } } } //================================================================================================= Standard_Boolean BRepLib_FindSurface::Found() const { return !mySurface.IsNull(); } //================================================================================================= Handle(Geom_Surface) BRepLib_FindSurface::Surface() const { return mySurface; } //================================================================================================= Standard_Real BRepLib_FindSurface::Tolerance() const { return myTolerance; } //================================================================================================= Standard_Real BRepLib_FindSurface::ToleranceReached() const { return myTolReached; } //================================================================================================= Standard_Boolean BRepLib_FindSurface::Existed() const { return isExisted; } //================================================================================================= TopLoc_Location BRepLib_FindSurface::Location() const { return myLocation; }