0028605: Improve the algorithm of calculation of valid intersection range of an edge
[occt.git] / src / BRepLib / BRepLib_1.cxx
1 // Created on: 2017-03-24
2 // Created by: Mikhail Sazonov
3 // Copyright (c) 2017 OPEN CASCADE SAS
4 //
5 // This file is part of Open CASCADE Technology software library.
6 //
7 // This library is free software; you can redistribute it and/or modify it under
8 // the terms of the GNU Lesser General Public License version 2.1 as published
9 // by the Free Software Foundation, with special exception defined in the file
10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11 // distribution for complete text of the license and disclaimer of any warranty.
12 //
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
15
16 #include <BRepLib.hxx>
17 #include <BRep_Tool.hxx>
18 #include <BRepAdaptor_Curve.hxx>
19 #include <Geom_OffsetCurve.hxx>
20 #include <Precision.hxx>
21 #include <TopExp.hxx>
22 #include <TopoDS_Vertex.hxx>
23
24 //=======================================================================
25 // function: findNearestValidPoint
26 // purpose : Starting from the appointed end of the curve, find the nearest
27 //           point on the curve that is an intersection with the sphere with
28 //           center theVertPnt and radius theTol.
29 //=======================================================================
30 static Standard_Boolean findNearestValidPoint(
31   const Adaptor3d_Curve& theCurve,
32   const Standard_Real theFirst, const Standard_Real theLast,
33   const Standard_Boolean isFirst,
34   const gp_Pnt& theVertPnt,
35   const Standard_Real theTol,
36   const Standard_Real theEps,
37   Standard_Real& thePar)
38 {
39   // 1. Check that the needed end is inside the sphere
40
41   Standard_Real aStartU = theFirst;
42   Standard_Real anEndU = theLast;
43   if (!isFirst)
44     std::swap(aStartU, anEndU);
45   gp_Pnt aP = theCurve.Value(aStartU);
46   const Standard_Real aSqTol = theTol * theTol;
47   if (aP.SquareDistance(theVertPnt) > aSqTol)
48     // the vertex does not cover the corresponding to this vertex end of the curve
49     return Standard_False;
50
51   // 2. Find a nearest point that is outside
52
53   // stepping along the curve by theTol till go out
54   //
55   // the general step is computed using general curve resolution
56   Standard_Real aStep = theCurve.Resolution(theTol) * 1.01;
57   // aD1Mag is a threshold to consider local derivative magnitude too small
58   // and to accelerate going out of sphere
59   // (inverse of resolution is the maximal derivative);
60   // this is actual for bezier and b-spline types only
61   Standard_Real aD1Mag = 0.;
62   GeomAbs_CurveType aType = theCurve.GetType();
63   if (aType == GeomAbs_OffsetCurve)
64   {
65     Handle(Geom_OffsetCurve) anOffsetCurve = theCurve.OffsetCurve();
66     Handle(Geom_Curve) aBaseCurve = anOffsetCurve->BasisCurve();
67     aType = GeomAdaptor_Curve(aBaseCurve).GetType();
68   }
69   if (aType == GeomAbs_BezierCurve || aType == GeomAbs_BSplineCurve)
70   {
71     aD1Mag = 1. / theCurve.Resolution(1.) * 0.01;
72     aD1Mag *= aD1Mag;
73   }
74   if (!isFirst)
75     aStep = -aStep;
76   Standard_Boolean isOut = Standard_False;
77   Standard_Real anUIn = aStartU;
78   Standard_Real anUOut = anUIn;
79   while (!isOut)
80   {
81     anUIn = anUOut;
82     anUOut += aStep;
83     if ((isFirst && anUOut > anEndU) || (!isFirst && anUOut < anEndU))
84     {
85       // step is too big and we go out of bounds,
86       // check if the opposite bound is outside
87       aP = theCurve.Value(anEndU);
88       isOut = (aP.SquareDistance(theVertPnt) > aSqTol);
89       if (!isOut)
90         // all range is inside sphere
91         return Standard_False;
92       anUOut = anEndU;
93       break;
94     }
95     if (aD1Mag > 0.)
96     {
97       Standard_Real aStepLocal = aStep;
98       for (;;)
99       {
100         // cycle to go out of local singularity
101         gp_Vec aD1;
102         theCurve.D1(anUOut, aP, aD1);
103         if (aD1.SquareMagnitude() < aD1Mag)
104         {
105           aStepLocal *= 2.;
106           anUOut += aStepLocal;
107           if ((isFirst && anUOut < anEndU) || (!isFirst && anUOut > anEndU))
108             // still in range
109             continue;
110           // went out of range, so check if the end point has out state
111           anUOut = anEndU;
112           aP = theCurve.Value(anUOut);
113           isOut = (aP.SquareDistance(theVertPnt) > aSqTol);
114           if (!isOut)
115             // all range is inside sphere
116             return Standard_False;
117         }
118         break;
119       }
120     }
121     else
122     {
123       aP = theCurve.Value(anUOut);
124     }
125     if (!isOut)
126       isOut = (aP.SquareDistance(theVertPnt) > aSqTol);
127   }
128
129   // 3. Precise solution with binary search
130
131   Standard_Real aDelta = Abs(anUOut - anUIn);
132   while (aDelta > theEps)
133   {
134     Standard_Real aMidU = (anUIn + anUOut) * 0.5;
135     aP = theCurve.Value(aMidU);
136     isOut = (aP.SquareDistance(theVertPnt) > aSqTol);
137     if (isOut)
138       anUOut = aMidU;
139     else
140       anUIn = aMidU;
141     aDelta = Abs(anUOut - anUIn);
142   }
143   thePar = (anUIn + anUOut) * 0.5;
144   return Standard_True;
145 }
146
147 //=======================================================================
148 // function: FindValidRange
149 // purpose : 
150 //=======================================================================
151 Standard_Boolean BRepLib::FindValidRange
152   (const Adaptor3d_Curve& theCurve, const Standard_Real theTolE,
153    const Standard_Real theParV1, const gp_Pnt& thePntV1, const Standard_Real theTolV1,
154    const Standard_Real theParV2, const gp_Pnt& thePntV2, const Standard_Real theTolV2,
155    Standard_Real& theFirst, Standard_Real& theLast)
156 {
157   if (theParV2 - theParV1 < Precision::PConfusion())
158     return Standard_False;
159   
160   Standard_Real anEps = Max(theCurve.Resolution(theTolE) * 0.1, Precision::PConfusion());
161
162   if (Precision::IsInfinite(theParV1))
163     theFirst = theParV1;
164   else
165   {
166     if (!findNearestValidPoint(theCurve, theParV1, theParV2, Standard_True,
167                                thePntV1, theTolV1, anEps, theFirst))
168       return Standard_False;
169     if (theParV2 - theFirst < anEps)
170       return Standard_False;
171   }
172
173   if (Precision::IsInfinite(theParV2))
174     theLast = theParV2;
175   else
176   {
177     if (!findNearestValidPoint(theCurve, theParV1, theParV2, Standard_False,
178                                thePntV2, theTolV2, anEps, theLast))
179       return Standard_False;
180     if (theLast - theParV1 < anEps)
181       return Standard_False;
182   }
183
184   // check found parameters
185   if (theFirst > theLast)
186   {
187     // overlapping, not valid range
188     return Standard_False;
189   }
190
191   return Standard_True;
192 }
193
194 //=======================================================================
195 // function: FindValidRange
196 // purpose : 
197 //=======================================================================
198 Standard_Boolean BRepLib::FindValidRange
199   (const TopoDS_Edge& theEdge, Standard_Real& theFirst, Standard_Real& theLast)
200 {
201   TopLoc_Location aLoc;
202   Standard_Real f, l;
203   if (BRep_Tool::Curve(theEdge, aLoc, f, l).IsNull())
204     return Standard_False;
205   BRepAdaptor_Curve anAC(theEdge);
206   Standard_Real aParV[2] = { anAC.FirstParameter(), anAC.LastParameter() };
207   if (aParV[1] - aParV[0] < Precision::PConfusion())
208     return Standard_False;
209
210   // get vertices
211   TopoDS_Vertex aV[2];
212   TopExp::Vertices(theEdge, aV[0], aV[1]);
213
214   Standard_Real aTolE = BRep_Tool::Tolerance(theEdge);
215   // to have correspondence with intersection precision
216   // the tolerances of vertices are increased on Precision::Confusion()
217   Standard_Real aTolV[2] = { Precision::Confusion(), Precision::Confusion() };
218   gp_Pnt aPntV[2];
219   for (Standard_Integer i = 0; i < 2; i++)
220   {
221     if (!aV[i].IsNull())
222     {
223       aTolV[i] += BRep_Tool::Tolerance(aV[i]);
224       aPntV[i] = BRep_Tool::Pnt(aV[i]);
225     }
226     else if (!Precision::IsInfinite(aParV[i]))
227     {
228       aTolV[i] += aTolE;
229       aPntV[i] = anAC.Value(aParV[i]);
230     }
231   }
232   return FindValidRange(anAC, aTolE, 
233                         aParV[0], aPntV[0], aTolV[0],
234                         aParV[1], aPntV[1], aTolV[1],
235                         theFirst, theLast);
236 }