0030435: Improving performance of Approx_ComputeCLine
[occt.git] / src / IntTools / IntTools_EdgeEdge.cxx
1 // Created by: Eugeny MALTCHIKOV
2 // Copyright (c) 2013-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15
16 #include <Bnd_Box.hxx>
17 #include <BndLib_Add3dCurve.hxx>
18 #include <TColStd_MapOfInteger.hxx>
19 #include <BRep_Tool.hxx>
20 #include <BRepAdaptor_Curve.hxx>
21 #include <ElCLib.hxx>
22 #include <Geom_BezierCurve.hxx>
23 #include <Geom_BSplineCurve.hxx>
24 #include <Geom_Line.hxx>
25 #include <Geom_Circle.hxx>
26 #include <Geom_Curve.hxx>
27 #include <Geom_Ellipse.hxx>
28 #include <Geom_OffsetCurve.hxx>
29 #include <GeomAPI_ProjectPointOnCurve.hxx>
30 #include <gp_Dir.hxx>
31 #include <gp_Lin.hxx>
32 #include <IntTools_CommonPrt.hxx>
33 #include <IntTools_EdgeEdge.hxx>
34 #include <IntTools_Range.hxx>
35 #include <IntTools_Tools.hxx>
36 #include <TopoDS_Edge.hxx>
37 #include <TopoDS_Iterator.hxx>
38 #include <BRepExtrema_DistShapeShape.hxx>
39
40 static 
41   void BndBuildBox(const BRepAdaptor_Curve& theBAC,
42                    const Standard_Real aT1,
43                    const Standard_Real aT2,
44                    const Standard_Real theTol,
45                    Bnd_Box& theBox);
46 static
47   Standard_Real PointBoxDistance(const Bnd_Box& aB,
48                                  const gp_Pnt& aP);
49 static 
50   Standard_Integer SplitRangeOnSegments(const Standard_Real aT1, 
51                                         const Standard_Real aT2,
52                                         const Standard_Real theResolution,
53                                         const Standard_Integer theNbSeg,
54                                         IntTools_SequenceOfRanges& theSegments);
55 static
56  Standard_Integer DistPC(const Standard_Real aT1, 
57                          const Handle(Geom_Curve)& theC1,
58                          const Standard_Real theCriteria, 
59                          GeomAPI_ProjectPointOnCurve& theProjector,
60                          Standard_Real& aD, 
61                          Standard_Real& aT2,
62                          const Standard_Integer iC = 1);
63 static
64  Standard_Integer DistPC(const Standard_Real aT1, 
65                          const Handle(Geom_Curve)& theC1,
66                          const Standard_Real theCriteria,
67                          GeomAPI_ProjectPointOnCurve& theProjector, 
68                          Standard_Real& aD, 
69                          Standard_Real& aT2,
70                          Standard_Real& aDmax,
71                          Standard_Real& aT1max,
72                          Standard_Real& aT2max,
73                          const Standard_Integer iC = 1);
74 static
75   Standard_Integer FindDistPC(const Standard_Real aT1A, 
76                               const Standard_Real aT1B,
77                               const Handle(Geom_Curve)& theC1,
78                               const Standard_Real theCriteria,
79                               const Standard_Real theEps,
80                               GeomAPI_ProjectPointOnCurve& theProjector,
81                               Standard_Real& aDmax, 
82                               Standard_Real& aT1max,
83                               Standard_Real& aT2max,
84                               const Standard_Boolean bMaxDist = Standard_True);
85 static
86   Standard_Real ResolutionCoeff(const BRepAdaptor_Curve& theBAC,
87                                 const IntTools_Range& theRange);
88 static
89   Standard_Real Resolution(const Handle(Geom_Curve)& theCurve,
90                            const GeomAbs_CurveType theCurveType,
91                            const Standard_Real theResCoeff,
92                            const Standard_Real theR3D);
93 static
94   Standard_Real CurveDeflection(const BRepAdaptor_Curve& theBAC,
95                                 const IntTools_Range& theRange);
96 static 
97   Standard_Boolean IsClosed(const Handle(Geom_Curve)& theCurve,
98                             const Standard_Real aT1,
99                             const Standard_Real aT2,
100                             const Standard_Real theTol,
101                             const Standard_Real theRes);
102 static 
103   Standard_Integer TypeToInteger(const GeomAbs_CurveType theCType);
104
105 //=======================================================================
106 //function : Prepare
107 //purpose  : 
108 //=======================================================================
109 void IntTools_EdgeEdge::Prepare()
110 {
111   GeomAbs_CurveType aCT1, aCT2;
112   Standard_Integer iCT1, iCT2;
113   //
114   myCurve1.Initialize(myEdge1);
115   myCurve2.Initialize(myEdge2);
116   //
117   if (myRange1.First() == 0. && myRange1.Last() == 0.) {
118     myRange1.SetFirst(myCurve1.FirstParameter());
119     myRange1.SetLast (myCurve1.LastParameter());
120   }
121   //
122   if (myRange2.First() == 0. && myRange2.Last() == 0.) {
123     myRange2.SetFirst(myCurve2.FirstParameter());
124     myRange2.SetLast (myCurve2.LastParameter());
125   }
126   //
127   aCT1 = myCurve1.GetType();
128   aCT2 = myCurve2.GetType();
129   //
130   iCT1 = TypeToInteger(aCT1);
131   iCT2 = TypeToInteger(aCT2);
132   //
133   if (iCT1 == iCT2) {
134     if (iCT1 != 0) {
135       //compute deflection
136       Standard_Real aC1, aC2;
137       //
138       aC2 = CurveDeflection(myCurve2, myRange2);
139       aC1 = (aC2 > Precision::Confusion()) ? 
140         CurveDeflection(myCurve1, myRange1) : 1.;
141       //
142       if (aC1 < aC2) {
143         --iCT1;
144       }
145     }
146   }
147   //
148   if (iCT1 < iCT2) {
149     TopoDS_Edge tmpE = myEdge1;
150     myEdge1 = myEdge2;
151     myEdge2 = tmpE;
152     //
153     BRepAdaptor_Curve tmpC = myCurve1;
154     myCurve1 = myCurve2;
155     myCurve2 = tmpC;
156     //
157     IntTools_Range tmpR = myRange1;
158     myRange1 = myRange2;
159     myRange2 = tmpR;
160     //
161     mySwap = Standard_True;
162   }
163   //
164   Standard_Real aTolAdd = myFuzzyValue / 2.;
165   myTol1 = myCurve1.Tolerance() + aTolAdd;
166   myTol2 = myCurve2.Tolerance() + aTolAdd;
167   myTol = myTol1 + myTol2;
168   //
169   if (iCT1 != 0 || iCT2 != 0) {
170     Standard_Real f, l, aTM;
171     //
172     myGeom1 = BRep_Tool::Curve(myEdge1, f, l);
173     myGeom2 = BRep_Tool::Curve(myEdge2, f, l);
174     //
175     myResCoeff1 = ResolutionCoeff(myCurve1, myRange1);
176     myResCoeff2 = ResolutionCoeff(myCurve2, myRange2);
177     //
178     myRes1 = Resolution(myCurve1.Curve().Curve(), myCurve1.GetType(), myResCoeff1, myTol1);
179     myRes2 = Resolution(myCurve2.Curve().Curve(), myCurve2.GetType(), myResCoeff2, myTol2);
180     //
181     myPTol1 = 5.e-13;
182     aTM = Max(fabs(myRange1.First()), fabs(myRange1.Last()));
183     if (aTM > 999.) {
184       myPTol1 = 5.e-16 * aTM;
185     }
186     //
187     myPTol2 = 5.e-13;
188     aTM = Max(fabs(myRange2.First()), fabs(myRange2.Last()));
189     if (aTM > 999.) {
190       myPTol2 = 5.e-16 * aTM;
191     }
192   }
193 }
194
195 //=======================================================================
196 //function : Perform
197 //purpose  : 
198 //=======================================================================
199 void IntTools_EdgeEdge::Perform()
200 {
201   //1. Check data
202   CheckData();
203   if (myErrorStatus) {
204     return;
205   }
206   //
207   //2. Prepare Data
208   Prepare();
209   //
210   //3.1. Check Line/Line case
211   if (myCurve1.GetType() == GeomAbs_Line &&
212       myCurve2.GetType() == GeomAbs_Line) {
213     ComputeLineLine();
214     return;
215   }
216   //
217   if (myQuickCoincidenceCheck) {
218     if (IsCoincident()) {
219       Standard_Real aT11, aT12, aT21, aT22;
220       //
221       myRange1.Range(aT11, aT12);
222       myRange2.Range(aT21, aT22);
223       AddSolution(aT11, aT12, aT21, aT22, TopAbs_EDGE);
224       return;
225     }
226   }
227   //
228   if ((myCurve1.GetType() <= GeomAbs_Parabola && myCurve2.GetType() <= GeomAbs_Parabola) &&
229       (myCurve1.GetType() == GeomAbs_Line || myCurve2.GetType() == GeomAbs_Line))
230   {
231     //Improvement of performance for cases of searching common parts between line  
232     //and analytical curve. This code allows to define that edges have no
233     //common parts more fast, then regular algorithm (FindSolution(...))
234     //Check minimal distance between edges
235     BRepExtrema_DistShapeShape  aMinDist(myEdge1, myEdge2, Extrema_ExtFlag_MIN);
236     if (aMinDist.IsDone())
237     {
238       Standard_Real d = aMinDist.Value();
239       if (d > 1.1 * myTol)
240       {
241         return;
242       }
243     }
244   }
245
246   IntTools_SequenceOfRanges aRanges1, aRanges2;
247   //
248   //3.2. Find ranges containig solutions
249   Standard_Boolean bSplit2;
250   FindSolutions(aRanges1, aRanges2, bSplit2);
251   //
252   //4. Merge solutions and save common parts
253   MergeSolutions(aRanges1, aRanges2, bSplit2);
254 }
255
256 //=======================================================================
257 //function :  IsCoincident
258 //purpose  : 
259 //=======================================================================
260 Standard_Boolean IntTools_EdgeEdge::IsCoincident() 
261 {
262   Standard_Integer i, iCnt, aNbSeg, aNbP2;
263   Standard_Real dT, aT1, aCoeff, aTresh, aD;
264   Standard_Real aT11, aT12, aT21, aT22;
265   GeomAPI_ProjectPointOnCurve aProjPC;
266   gp_Pnt aP1;
267   //
268   aTresh=0.5;
269   aNbSeg=23;
270   myRange1.Range(aT11, aT12);
271   myRange2.Range(aT21, aT22);
272   //
273   aProjPC.Init(myGeom2, aT21, aT22);
274   //
275   dT=(aT12-aT11)/aNbSeg;
276   //
277   iCnt=0;
278   for(i=0; i <= aNbSeg; ++i) {
279     aT1 = aT11+i*dT;
280     myGeom1->D0(aT1, aP1);
281     //
282     aProjPC.Perform(aP1);
283     aNbP2=aProjPC.NbPoints();
284     if (!aNbP2) {
285       continue;
286     }
287     //
288     aD=aProjPC.LowerDistance();
289     if(aD < myTol) {
290       ++iCnt; 
291     }
292   }
293   //
294   aCoeff=(Standard_Real)iCnt/((Standard_Real)aNbSeg+1);
295   return aCoeff > aTresh;
296 }
297 //=======================================================================
298 //function : FindSolutions
299 //purpose  : 
300 //=======================================================================
301 void IntTools_EdgeEdge::FindSolutions(IntTools_SequenceOfRanges& theRanges1,
302                                       IntTools_SequenceOfRanges& theRanges2,
303                                       Standard_Boolean& bSplit2)
304 {
305   Standard_Boolean bIsClosed2;
306   Standard_Real aT11, aT12, aT21, aT22;
307   Bnd_Box aB2;
308   //
309   bSplit2 = Standard_False;
310   myRange1.Range(aT11, aT12);
311   myRange2.Range(aT21, aT22);
312   //
313   bIsClosed2 = IsClosed(myGeom2, aT21, aT22, myTol2, myRes2);
314   //
315   if (bIsClosed2) {
316     Bnd_Box aB1;
317     BndBuildBox(myCurve1, aT11, aT12, myTol1, aB1);
318     //
319     gp_Pnt aP = myGeom2->Value(aT21);
320     bIsClosed2 = !aB1.IsOut(aP);
321   }
322   //
323   if (!bIsClosed2) {
324     BndBuildBox(myCurve2, aT21, aT22, myTol2, aB2);
325     FindSolutions(myRange1, myRange2, aB2, theRanges1, theRanges2);
326     return;
327   }
328   //
329   if (!CheckCoincidence(aT11, aT12, aT21, aT22, myTol, myRes1)) {
330     theRanges1.Append(myRange1);
331     theRanges2.Append(myRange2);
332     return;
333   }
334   //
335   Standard_Integer i, j, aNb1, aNb2;
336   IntTools_SequenceOfRanges aSegments1, aSegments2;
337   //
338   aNb1 = IsClosed(myGeom1, aT11, aT12, myTol1, myRes1) ? 2 : 1;
339   aNb2 = 2;
340   //
341   aNb1 = SplitRangeOnSegments(aT11, aT12, myRes1, aNb1, aSegments1);
342   aNb2 = SplitRangeOnSegments(aT21, aT22, myRes2, aNb2, aSegments2);
343   //
344   for (i = 1; i <= aNb1; ++i) {
345     const IntTools_Range& aR1 = aSegments1(i);
346     for (j = 1; j <= aNb2; ++j) {
347       const IntTools_Range& aR2 = aSegments2(j);
348       BndBuildBox(myCurve2, aR2.First(), aR2.Last(), myTol2, aB2);
349       FindSolutions(aR1, aR2, aB2, theRanges1, theRanges2);
350     }
351   }
352   //
353   bSplit2 = aNb2 > 1;
354 }
355
356 //=======================================================================
357 //function : FindSolutions
358 //purpose  : 
359 //=======================================================================
360 void IntTools_EdgeEdge::FindSolutions(const IntTools_Range& theR1,
361                                       const IntTools_Range& theR2,
362                                       const Bnd_Box& theBox2,
363                                       IntTools_SequenceOfRanges& theRanges1,
364                                       IntTools_SequenceOfRanges& theRanges2)
365 {
366   Standard_Boolean bOut, bStop, bThin;
367   Standard_Real aT11, aT12, aT21, aT22;
368   Standard_Real aTB11, aTB12, aTB21, aTB22;
369   Standard_Real aSmallStep1, aSmallStep2;
370   Standard_Integer iCom;
371   Bnd_Box aB1, aB2;
372   //
373   theR1.Range(aT11, aT12);
374   theR2.Range(aT21, aT22);
375   //
376   aB2 = theBox2;
377   //
378   bThin = Standard_False;
379   bStop = Standard_False;
380   iCom  = 1;
381   //
382   do {
383     aTB11 = aT11;
384     aTB12 = aT12;
385     aTB21 = aT21;
386     aTB22 = aT22;
387     //
388     //1. Build box for first edge and find parameters 
389     //   of the second one in that box
390     BndBuildBox(myCurve1, aT11, aT12, myTol1, aB1);
391     bOut = aB1.IsOut(aB2);
392     if (bOut) {
393       break;
394     }
395     //
396     bThin = ((aT12 - aT11) < myRes1) ||
397       (aB1.IsXThin(myTol) && aB1.IsYThin(myTol) && aB1.IsZThin(myTol));
398     //
399     bOut = !FindParameters(myCurve2, aTB21, aTB22, myTol2, myRes2, myPTol2, 
400                            myResCoeff2, aB1, aT21, aT22);
401     if (bOut || bThin) {
402       break;
403     }
404     //
405     //2. Build box for second edge and find parameters 
406     //   of the first one in that box
407     BndBuildBox(myCurve2, aT21, aT22, myTol2, aB2);
408     bOut = aB1.IsOut(aB2);
409     if (bOut) {
410       break;
411     }
412     //
413     bThin = ((aT22 - aT21) < myRes2) ||
414       (aB2.IsXThin(myTol) && aB2.IsYThin(myTol) && aB2.IsZThin(myTol));
415     //
416     bOut = !FindParameters(myCurve1, aTB11, aTB12, myTol1, myRes1, myPTol1,
417                            myResCoeff1, aB2, aT11, aT12);
418     //
419     if (bOut || bThin) {
420       break;
421     }
422     //
423     //3. Check if it makes sense to continue
424     aSmallStep1 = (aTB12 - aTB11) / 250.;
425     aSmallStep2 = (aTB22 - aTB21) / 250.;
426     //
427     if (aSmallStep1 < myRes1) {
428       aSmallStep1 = myRes1;
429     }
430     if (aSmallStep2 < myRes2) {
431       aSmallStep2 = myRes2;
432     }
433     //
434     if (((aT11 - aTB11) < aSmallStep1) && ((aTB12 - aT12) < aSmallStep1) &&
435         ((aT21 - aTB21) < aSmallStep2) && ((aTB22 - aT22) < aSmallStep2)) {
436       bStop = Standard_True;
437     }
438     //
439   } while (!bStop);
440   //
441   if (bOut) {
442     //no intersection;
443     return;
444   }
445   //
446   if (!bThin) {
447     //check curves for coincidence on the ranges
448     iCom = CheckCoincidence(aT11, aT12, aT21, aT22, myTol, myRes1);
449     if (!iCom) {
450       bThin = Standard_True;
451     }
452   }
453   //
454   if (bThin) {
455     if (iCom != 0) {
456       //check intermediate points
457       Standard_Boolean bSol;
458       Standard_Real aT1;
459       gp_Pnt aP1;
460       GeomAPI_ProjectPointOnCurve aProjPC;
461       //
462       aT1 = (aT11 + aT12) * .5;
463       myGeom1->D0(aT1, aP1);
464       //
465       aProjPC.Init(myGeom2, aT21, aT22);
466       aProjPC.Perform(aP1);
467       //
468       if (aProjPC.NbPoints()) {
469         bSol = aProjPC.LowerDistance() <= myTol;
470       }
471       else {
472         Standard_Real aT2;
473         gp_Pnt aP2;
474         //
475         aT2 = (aT21 + aT22) * .5;
476         myGeom2->D0(aT2, aP2);
477         //
478         bSol = aP1.IsEqual(aP2, myTol);
479       }
480       //
481       if (!bSol) {
482         return;
483       }
484     }
485     //add common part
486     IntTools_Range aR1(aT11, aT12), aR2(aT21, aT22);
487     //
488     theRanges1.Append(aR1);
489     theRanges2.Append(aR2);
490     return;
491   }
492   //
493   if (!IsIntersection(aT11, aT12, aT21, aT22)) {
494     return;
495   }
496   //
497   //split ranges on segments and repeat
498   Standard_Integer i, aNb1;
499   IntTools_SequenceOfRanges aSegments1;
500   //
501   IntTools_Range aR2(aT21, aT22);
502   BndBuildBox(myCurve2, aT21, aT22, myTol2, aB2);
503   //
504   aNb1 = SplitRangeOnSegments(aT11, aT12, myRes1, 3, aSegments1);
505   for (i = 1; i <= aNb1; ++i) {
506     const IntTools_Range& aR1 = aSegments1(i);
507     FindSolutions(aR1, aR2, aB2, theRanges1, theRanges2);
508   }
509 }
510
511 //=======================================================================
512 //function : FindParameters
513 //purpose  : 
514 //=======================================================================
515 Standard_Boolean IntTools_EdgeEdge::FindParameters(const BRepAdaptor_Curve& theBAC,
516                                                    const Standard_Real aT1, 
517                                                    const Standard_Real aT2,
518                                                    const Standard_Real theTol,
519                                                    const Standard_Real theRes,
520                                                    const Standard_Real thePTol,
521                                                    const Standard_Real theResCoeff,
522                                                    const Bnd_Box& theCBox,
523                                                    Standard_Real& aTB1, 
524                                                    Standard_Real& aTB2)
525 {
526   Standard_Boolean bRet;
527   Standard_Integer aC, i;
528   Standard_Real aCf, aDiff, aDt, aT, aTB, aTOut, aTIn;
529   Standard_Real aDist, aDistP;
530   gp_Pnt aP;
531   Bnd_Box aCBx;
532   //
533   bRet = Standard_False;
534   aCf = 0.6180339887498948482045868343656;// =0.5*(1.+sqrt(5.))/2.;
535   aCBx = theCBox;
536   aCBx.SetGap(aCBx.GetGap() + theTol);
537   //
538   const Handle(Geom_Curve)& aCurve = theBAC.Curve().Curve();
539   const GeomAbs_CurveType aCurveType = theBAC.GetType();
540   Standard_Real aMaxDt = (aT2 - aT1) * 0.01;
541   //
542   for (i = 0; i < 2; ++i) {
543     aTB = !i ? aT1 : aT2;
544     aT = !i ? aT2 : aTB1;
545     aC = !i ? 1 : -1;
546     aDt = theRes;
547     aDistP = 0.;
548     bRet = Standard_False;
549     Standard_Real k = 1;
550     //looking for the point on the edge which is in the box;
551     while (aC*(aT-aTB) >= 0) {
552       theBAC.D0(aTB, aP);
553       aDist = PointBoxDistance(theCBox, aP);
554       if (aDist > theTol) {
555         if (aDistP > 0.) {
556           Standard_Boolean toGrow = Standard_False;
557           if (Abs(aDistP - aDist) / aDistP < 0.1) {
558             aDt = Resolution(aCurve, aCurveType, theResCoeff, k*aDist);
559             if (aDt < aMaxDt)
560             {
561               toGrow = Standard_True;
562               k *= 2;
563             }
564           }
565           if (!toGrow) {
566             k = 1;
567             aDt = Resolution(aCurve, aCurveType, theResCoeff, aDist);
568           }
569         }
570         aTB += aC*aDt;
571       } else {
572         bRet = Standard_True;
573         break;
574       }
575       aDistP = aDist;
576     }
577     //
578     if (!bRet) {
579       if (!i) {
580         //edge is out of the box;
581         return bRet;
582       } else {
583         bRet = !bRet;
584         aTB = aTB1;
585         aDt = aT2 - aTB1;
586       }
587     }
588     //
589     aT = !i ? aT1 : aT2;
590     if (aTB != aT) {
591       //one point IN, one point OUT; looking for the bounding point;
592       aTIn = aTB;
593       aTOut = aTB - aC*aDt;
594       aDiff = aTIn - aTOut;
595       while (fabs(aDiff) > thePTol) {
596         aTB = aTOut + aDiff*aCf;
597         theBAC.D0(aTB, aP);
598         if (aCBx.IsOut(aP)) {
599           aTOut = aTB;
600         } else {
601           aTIn = aTB;
602         }
603         aDiff = aTIn - aTOut;
604       }
605     }
606     if (!i) {
607       aTB1 = aTB;
608     } else {
609       aTB2 = aTB;
610     }
611   }
612   return bRet;
613 }
614
615 //=======================================================================
616 //function : MergeSolutions
617 //purpose  : 
618 //=======================================================================
619 void IntTools_EdgeEdge::MergeSolutions(const IntTools_SequenceOfRanges& theRanges1, 
620                                        const IntTools_SequenceOfRanges& theRanges2,
621                                        const Standard_Boolean bSplit2)
622 {
623   Standard_Integer aNbCP = theRanges1.Length();
624   if (aNbCP == 0) {
625     return;
626   }
627   //
628   IntTools_Range aRi1, aRi2, aRj1, aRj2;
629   Standard_Boolean bCond;
630   Standard_Integer i, j;
631   TopAbs_ShapeEnum aType;
632   Standard_Real aT11, aT12, aT21, aT22;
633   Standard_Real aTi11, aTi12, aTi21, aTi22;
634   Standard_Real aTj11, aTj12, aTj21, aTj22;
635   Standard_Real aRes1, aRes2, dTR1, dTR2;
636   TColStd_MapOfInteger aMI;
637   //
638   aRes1 = Resolution(myCurve1.Curve().Curve(), 
639                      myCurve1.GetType(), myResCoeff1, myTol);
640   aRes2 = Resolution(myCurve2.Curve().Curve(), 
641                      myCurve2.GetType(), myResCoeff2, myTol);
642   //
643   myRange1.Range(aT11, aT12);
644   myRange2.Range(aT21, aT22);
645   dTR1 = 20*aRes1;
646   dTR2 = 20*aRes2;
647   aType = TopAbs_VERTEX;
648   //
649   for (i = 1; i <= aNbCP;) {
650     if (aMI.Contains(i)) {
651       ++i;
652       continue;
653     }
654     //
655     aRi1 = theRanges1(i);
656     aRi2 = theRanges2(i);
657     //
658     aRi1.Range(aTi11, aTi12);
659     aRi2.Range(aTi21, aTi22);
660     //
661     aMI.Add(i);
662     //
663     for (j = i+1; j <= aNbCP; ++j) {
664       if (aMI.Contains(j)) {
665         continue;
666       }
667       //
668       aRj1 = theRanges1(j);
669       aRj2 = theRanges2(j);
670       //
671       aRj1.Range(aTj11, aTj12);
672       aRj2.Range(aTj21, aTj22);
673       //
674       bCond = (fabs(aTi12 - aTj11) < dTR1) ||
675         (bSplit2 && (fabs(aTj12 - aTi11) < dTR1));
676       if (bCond && bSplit2) {
677         bCond = (fabs((Max(aTi22, aTj22) - Min(aTi21, aTj21)) - 
678                       ((aTi22 - aTi21) + (aTj22 - aTj21))) < dTR2);
679       }
680       //
681       if (bCond) {
682         aTi11 = Min(aTi11, aTj11);
683         aTi12 = Max(aTi12, aTj12);
684         aTi21 = Min(aTi21, aTj21);
685         aTi22 = Max(aTi22, aTj22);
686         aMI.Add(j);
687       }
688       else if (!bSplit2) {
689         i = j;
690         break;
691       }
692     }
693     //
694     if (((fabs(aT11 - aTi11) < myRes1) && (fabs(aT12 - aTi12) < myRes1)) ||
695         ((fabs(aT21 - aTi21) < myRes2) && (fabs(aT22 - aTi22) < myRes2))) {
696       aType = TopAbs_EDGE;
697       myCommonParts.Clear();
698     }
699     //
700     AddSolution(aTi11, aTi12, aTi21, aTi22, aType);
701     if (aType == TopAbs_EDGE) {
702       break;
703     }
704     //
705     if (bSplit2) {
706       ++i;
707     }
708   }
709 }
710
711 //=======================================================================
712 //function : AddSolution
713 //purpose  : 
714 //=======================================================================
715 void IntTools_EdgeEdge::AddSolution(const Standard_Real aT11,
716                                     const Standard_Real aT12,
717                                     const Standard_Real aT21,
718                                     const Standard_Real aT22,
719                                     const TopAbs_ShapeEnum theType)
720 {
721   IntTools_CommonPrt aCPart;
722   //
723   aCPart.SetType(theType);
724   if (!mySwap) {
725     aCPart.SetEdge1(myEdge1);
726     aCPart.SetEdge2(myEdge2);
727     aCPart.SetRange1(aT11, aT12);
728     aCPart.AppendRange2(aT21, aT22);
729   } else {
730     aCPart.SetEdge1(myEdge2);
731     aCPart.SetEdge2(myEdge1);
732     aCPart.SetRange1(aT21, aT22);
733     aCPart.AppendRange2(aT11, aT12);
734   }
735   //
736   if (theType == TopAbs_VERTEX) {
737     Standard_Real aT1, aT2;
738     //
739     FindBestSolution(aT11, aT12, aT21, aT22, aT1, aT2);
740     //
741     if (!mySwap) {
742       aCPart.SetVertexParameter1(aT1);
743       aCPart.SetVertexParameter2(aT2);
744     } else {
745       aCPart.SetVertexParameter1(aT2);
746       aCPart.SetVertexParameter2(aT1);
747     }
748   }
749   myCommonParts.Append(aCPart);
750 }
751
752 //=======================================================================
753 //function : FindBestSolution
754 //purpose  : 
755 //=======================================================================
756 void IntTools_EdgeEdge::FindBestSolution(const Standard_Real aT11,
757                                          const Standard_Real aT12,
758                                          const Standard_Real aT21,
759                                          const Standard_Real aT22,
760                                          Standard_Real& aT1,
761                                          Standard_Real& aT2)
762 {
763   Standard_Integer i, aNbS, iErr;
764   Standard_Real aDMin, aD, aRes1, aSolCriteria, aTouchCriteria;
765   Standard_Real aT1A, aT1B, aT1Min, aT2Min;
766   GeomAPI_ProjectPointOnCurve aProjPC;
767   IntTools_SequenceOfRanges aRanges;
768   //
769   aDMin = Precision::Infinite();
770   aSolCriteria   = 5.e-16;
771   aTouchCriteria = 5.e-13;
772   Standard_Boolean bTouch = Standard_False;
773   Standard_Boolean bTouchConfirm = Standard_False;
774   //
775   aRes1 = Resolution(myCurve1.Curve().Curve(), 
776                      myCurve1.GetType(), myResCoeff1, myTol);
777   aNbS = 10;
778   aNbS = SplitRangeOnSegments(aT11, aT12, 3*aRes1, aNbS, aRanges);
779   //
780   aProjPC.Init(myGeom2, aT21, aT22);
781   //
782   Standard_Real aT11Touch = aT11, aT12Touch = aT12;
783   Standard_Real aT21Touch = aT21, aT22Touch = aT22;
784   Standard_Boolean isSolFound = Standard_False;
785   for (i = 1; i <= aNbS; ++i) {
786     const IntTools_Range& aR1 = aRanges(i);
787     aR1.Range(aT1A, aT1B);
788     //
789     aD = myTol;
790     iErr = FindDistPC(aT1A, aT1B, myGeom1, aSolCriteria, myPTol1,
791                       aProjPC, aD, aT1Min, aT2Min, Standard_False);
792     if (iErr != 1) {
793       if (aD < aDMin) {
794         aT1 = aT1Min;
795         aT2 = aT2Min;
796         aDMin = aD;
797         isSolFound = Standard_True;
798       }
799       //
800       if (aD < aTouchCriteria) {
801         if (bTouch) {
802           aT12Touch = aT1Min;
803           aT22Touch = aT2Min;
804           bTouchConfirm = Standard_True;
805         }
806         else {
807           aT11Touch = aT1Min;
808           aT21Touch = aT2Min;
809           bTouch = Standard_True;
810         }
811       }
812     }
813   }
814   if (!isSolFound || bTouchConfirm)
815   {
816     aT1 = (aT11Touch + aT12Touch) * 0.5;
817     iErr = DistPC(aT1, myGeom1, aSolCriteria, aProjPC, aD, aT2, -1);
818     if (iErr == 1) {
819       aT2 = (aT21Touch + aT22Touch) * 0.5;
820     }
821   }
822 }
823
824 //=======================================================================
825 //function : ComputeLineLine
826 //purpose  : 
827 //=======================================================================
828 void IntTools_EdgeEdge::ComputeLineLine()
829 {
830   Standard_Boolean IsParallel, IsCoincide;
831   Standard_Real aSin, aCos, aAng, aTol;
832   Standard_Real aT1, aT2, aT11, aT12, aT21, aT22;
833   gp_Pnt aP11, aP12;
834   gp_Lin aL1, aL2;
835   gp_Dir aD1, aD2;
836   IntTools_CommonPrt aCommonPrt;
837   //
838   IsParallel = Standard_False;
839   IsCoincide = Standard_False;
840   aTol = myTol*myTol;
841   aL1 = myCurve1.Line();
842   aL2 = myCurve2.Line();
843   aD1 = aL1.Position().Direction();
844   aD2 = aL2.Position().Direction();
845   myRange1.Range(aT11, aT12);
846   myRange2.Range(aT21, aT22);
847   //
848   aCommonPrt.SetEdge1(myEdge1);
849   aCommonPrt.SetEdge2(myEdge2);
850   //
851   aCos = aD1.Dot(aD2);
852   aAng = (aCos >= 0.) ? 2.*(1. - aCos) : 2.*(1. + aCos);
853   //
854   if(aAng <= Precision::Angular()) {
855     IsParallel = Standard_True;
856     if(aL1.SquareDistance(aL2.Location()) <= aTol) {
857       IsCoincide = Standard_True;
858       aP11 = ElCLib::Value(aT11, aL1);
859       aP12 = ElCLib::Value(aT12, aL1);
860     }
861   }
862   else {
863     aP11 = ElCLib::Value(aT11, aL1);
864     aP12 = ElCLib::Value(aT12, aL1);
865     if(aL2.SquareDistance(aP11) <= aTol && aL2.SquareDistance(aP12) <= aTol) {
866       IsCoincide = Standard_True;
867     }
868   }
869   //
870   if (IsCoincide) {
871     Standard_Real t21, t22;
872     //
873     t21 = ElCLib::Parameter(aL2, aP11);
874     t22 = ElCLib::Parameter(aL2, aP12);
875     if((t21 > aT22 && t22 > aT22) || (t21 < aT21 && t22 < aT21)) {
876       return;
877     }
878     //
879     Standard_Real temp;
880     if(t21 > t22) {
881       temp = t21;
882       t21 = t22;
883       t22 = temp;
884     }
885     //
886     if(t21 >= aT21) {
887       if(t22 <= aT22) {
888         aCommonPrt.SetRange1(aT11, aT12);
889         aCommonPrt.SetAllNullFlag(Standard_True);
890         aCommonPrt.AppendRange2(t21, t22);
891       }
892       else {
893         aCommonPrt.SetRange1(aT11, aT12 - (t22 - aT22));
894         aCommonPrt.AppendRange2(t21, aT22);
895       }
896     }
897     else {
898       aCommonPrt.SetRange1(aT11 + (aT21 - t21), aT12);
899       aCommonPrt.AppendRange2(aT21, t22);
900     }
901     aCommonPrt.SetType(TopAbs_EDGE);  
902     myCommonParts.Append(aCommonPrt);
903     return;
904   }
905   //
906   if (IsParallel) {
907     return;
908   }
909   //
910   {
911     TopoDS_Iterator aIt1, aIt2;
912     aIt1.Initialize(myEdge1);
913     for (; aIt1.More(); aIt1.Next()) {
914       const TopoDS_Shape& aV1 = aIt1.Value();
915       aIt2.Initialize(myEdge2);
916       for (; aIt2.More(); aIt2.Next()) {
917         const TopoDS_Shape& aV2 = aIt2.Value();
918         if (aV2.IsSame(aV1)) {
919           return;
920         }
921       }
922     }
923   }
924   //
925   aSin = 1. - aCos*aCos;
926   gp_Pnt O1 = aL1.Location();
927   gp_Pnt O2 = aL2.Location();
928   gp_Vec O1O2 (O1, O2);
929   //
930   aT2 = (aD1.XYZ()*(O1O2.Dot(aD1))-(O1O2.XYZ())).Dot(aD2.XYZ());
931   aT2 /= aSin;
932   //
933   if(aT2 < aT21 || aT2 > aT22) {
934     return;
935   }
936   //
937   gp_Pnt aP2(ElCLib::Value(aT2, aL2));
938   aT1 = (gp_Vec(O1, aP2)).Dot(aD1);
939   //
940   if(aT1 < aT11 || aT1 > aT12) {
941     return;
942   }
943   //
944   gp_Pnt aP1(ElCLib::Value(aT1, aL1));
945   Standard_Real aDist = aP1.SquareDistance(aP2);
946   //
947   if (aDist > aTol) {
948     return;
949   }
950   //
951   // compute correct range on the edges
952   Standard_Real anAngle, aDt1, aDt2;
953   //
954   anAngle = aD1.Angle(aD2);
955   //
956   aDt1 = IntTools_Tools::ComputeIntRange(myTol1, myTol2, anAngle);
957   aDt2 = IntTools_Tools::ComputeIntRange(myTol2, myTol1, anAngle);
958   //
959   aCommonPrt.SetRange1(aT1 - aDt1, aT1 + aDt1);
960   aCommonPrt.AppendRange2(aT2 - aDt2, aT2 + aDt2);
961   aCommonPrt.SetType(TopAbs_VERTEX);
962   aCommonPrt.SetVertexParameter1(aT1);
963   aCommonPrt.SetVertexParameter2(aT2);
964   myCommonParts.Append(aCommonPrt);
965 }
966
967 //=======================================================================
968 //function : IsIntersection
969 //purpose  : 
970 //=======================================================================
971 Standard_Boolean IntTools_EdgeEdge::IsIntersection(const Standard_Real aT11,
972                                                    const Standard_Real aT12,
973                                                    const Standard_Real aT21,
974                                                    const Standard_Real aT22)
975 {
976   Standard_Boolean bRet;
977   gp_Pnt aP11, aP12, aP21, aP22;
978   gp_Vec aV11, aV12, aV21, aV22;
979   Standard_Real aD11_21, aD11_22, aD12_21, aD12_22, aCriteria, aCoef;
980   Standard_Boolean bSmall_11_21, bSmall_11_22, bSmall_12_21, bSmall_12_22;
981   //
982   bRet = Standard_True;
983   aCoef = 1.e+5;
984   if (((aT12 - aT11) > aCoef*myRes1) && ((aT22 - aT21) > aCoef*myRes2)) {
985     aCoef = 5000;
986   } else {
987     Standard_Real aTRMin = Min((aT12 - aT11)/myRes1, (aT22 - aT21)/myRes2);
988     aCoef = aTRMin / 100.;
989     if (aCoef < 1.) {
990       aCoef = 1.;
991     }
992   }
993   aCriteria = aCoef * myTol;
994   aCriteria *= aCriteria;
995   //
996   myGeom1->D1(aT11, aP11, aV11);
997   myGeom1->D1(aT12, aP12, aV12);
998   myGeom2->D1(aT21, aP21, aV21);
999   myGeom2->D1(aT22, aP22, aV22);
1000   //
1001   aD11_21 = aP11.SquareDistance(aP21);
1002   aD11_22 = aP11.SquareDistance(aP22);
1003   aD12_21 = aP12.SquareDistance(aP21);
1004   aD12_22 = aP12.SquareDistance(aP22);
1005   //
1006   bSmall_11_21 = aD11_21 < aCriteria;
1007   bSmall_11_22 = aD11_22 < aCriteria;
1008   bSmall_12_21 = aD12_21 < aCriteria;
1009   bSmall_12_22 = aD12_22 < aCriteria;
1010   //
1011   if ((bSmall_11_21 && bSmall_12_22) ||
1012       (bSmall_11_22 && bSmall_12_21)) {
1013     if (aCoef == 1.) {
1014       return bRet;
1015     }
1016     //
1017     Standard_Real anAngleCriteria;
1018     Standard_Real anAngle1 = 0.0,
1019                   anAngle2 = 0.0;
1020     //
1021     anAngleCriteria = 5.e-3;
1022     if (aV11.SquareMagnitude() > Precision::SquareConfusion() &&
1023         aV12.SquareMagnitude() > Precision::SquareConfusion() &&
1024         aV21.SquareMagnitude() > Precision::SquareConfusion() &&
1025         aV22.SquareMagnitude() > Precision::SquareConfusion() )
1026     {
1027       if (bSmall_11_21 && bSmall_12_22) {
1028         anAngle1 = aV11.Angle(aV21);
1029         anAngle2 = aV12.Angle(aV22);
1030       } else {
1031         anAngle1 = aV11.Angle(aV22);
1032         anAngle2 = aV12.Angle(aV21);
1033       }
1034     }
1035     //
1036     if (((anAngle1 < anAngleCriteria) || ((M_PI - anAngle1) < anAngleCriteria)) ||
1037         ((anAngle2 < anAngleCriteria) || ((M_PI - anAngle2) < anAngleCriteria))) {
1038       GeomAPI_ProjectPointOnCurve aProjPC;
1039       Standard_Integer iErr;
1040       Standard_Real aD, aT1Min, aT2Min;
1041       //
1042       aD = Precision::Infinite();
1043       aProjPC.Init(myGeom2, aT21, aT22);
1044       iErr = FindDistPC(aT11, aT12, myGeom1, myTol, myRes1, 
1045                         aProjPC, aD, aT1Min, aT2Min, Standard_False);
1046       bRet = (iErr == 2);
1047     }
1048   }
1049   return bRet;
1050 }
1051
1052 //=======================================================================
1053 //function : CheckCoincidence
1054 //purpose  : 
1055 //=======================================================================
1056 Standard_Integer IntTools_EdgeEdge::CheckCoincidence(const Standard_Real aT11,
1057                                                      const Standard_Real aT12,
1058                                                      const Standard_Real aT21,
1059                                                      const Standard_Real aT22,
1060                                                      const Standard_Real theCriteria,
1061                                                      const Standard_Real theCurveRes1)
1062 {
1063   Standard_Integer iErr, aNb, aNb1, i;
1064   Standard_Real aT1A, aT1B, aT1max, aT2max, aDmax;
1065   GeomAPI_ProjectPointOnCurve aProjPC;
1066   IntTools_SequenceOfRanges aRanges;
1067   //
1068   iErr  = 0;
1069   aDmax = -1.;
1070   aProjPC.Init(myGeom2, aT21, aT22);
1071   //
1072   // 1. Express evaluation
1073   aNb = 10; // Number of intervals on the curve #1
1074   aNb1 = SplitRangeOnSegments(aT11, aT12, theCurveRes1, aNb, aRanges);
1075   for (i = 1; i < aNb1; ++i) {
1076     const IntTools_Range& aR1 = aRanges(i);
1077     aR1.Range(aT1A, aT1B);
1078     //
1079     iErr = DistPC(aT1B, myGeom1, theCriteria, aProjPC, aDmax, aT2max);
1080     if (iErr) {
1081       return iErr;
1082     }
1083   }
1084   //
1085   // if the ranges in aRanges are less than theCurveRes1,
1086   // there is no need to do step 2 (deep evaluation)
1087   if (aNb1 < aNb) {
1088     return iErr;
1089   }
1090   //
1091   // 2. Deep evaluation
1092   for (i = 2; i < aNb1; ++i) {
1093     const IntTools_Range& aR1 = aRanges(i);
1094     aR1.Range(aT1A, aT1B);
1095     //
1096     iErr = FindDistPC(aT1A, aT1B, myGeom1, theCriteria, theCurveRes1, 
1097                       aProjPC, aDmax, aT1max, aT2max);
1098     if (iErr) {
1099       return iErr;
1100     }
1101   }
1102   // Possible values:
1103   // iErr == 0 - the patches are coincided
1104   // iErr == 1 - a point from aC1 can not be projected on aC2
1105   // iErr == 2 - the distance is too big
1106   return iErr;
1107 }
1108
1109 //=======================================================================
1110 //function : FindDistPC
1111 //purpose  : 
1112 //=======================================================================
1113 Standard_Integer FindDistPC(const Standard_Real aT1A, 
1114                             const Standard_Real aT1B,
1115                             const Handle(Geom_Curve)& theC1,
1116                             const Standard_Real theCriteria,
1117                             const Standard_Real theEps,
1118                             GeomAPI_ProjectPointOnCurve& theProjPC,
1119                             Standard_Real& aDmax, 
1120                             Standard_Real& aT1max,
1121                             Standard_Real& aT2max,
1122                             const Standard_Boolean bMaxDist) 
1123 {
1124   Standard_Integer iErr, iC;
1125   Standard_Real aGS, aXP, aA, aB, aXL, aYP, aYL, aT2P, aT2L;
1126   //
1127   iC = bMaxDist ? 1 : -1;
1128   iErr = 0;
1129   aT1max = aT2max = 0.; // silence GCC warning
1130   //
1131   aGS = 0.6180339887498948482045868343656;// =0.5*(1.+sqrt(5.))-1.;
1132   aA = aT1A;
1133   aB = aT1B;
1134   //
1135   // check bounds
1136   iErr = DistPC(aA, theC1, theCriteria, theProjPC, 
1137                 aYP, aT2P, aDmax, aT1max, aT2max, iC);
1138   if (iErr == 2) {
1139     return iErr;
1140   }
1141   //
1142   iErr = DistPC(aB, theC1, theCriteria, theProjPC, 
1143                 aYL, aT2L, aDmax, aT1max, aT2max, iC);
1144   if (iErr == 2) {
1145     return iErr;
1146   }
1147   //
1148   aXP = aA + (aB - aA)*aGS;
1149   aXL = aB - (aB - aA)*aGS;
1150   //
1151   iErr = DistPC(aXP, theC1, theCriteria, theProjPC, 
1152                 aYP, aT2P, aDmax, aT1max, aT2max, iC);
1153   if (iErr) {
1154     return iErr;
1155   }
1156   //
1157   iErr = DistPC(aXL, theC1, theCriteria, theProjPC, 
1158                 aYL, aT2L, aDmax, aT1max, aT2max, iC);
1159   if (iErr) {
1160     return iErr;
1161   }
1162   //
1163   Standard_Real anEps = Max(theEps, Epsilon(Max(Abs(aA), Abs(aB))) * 10.);
1164   for (;;) {
1165     if (iC*(aYP - aYL) > 0) {
1166       aA = aXL;
1167       aXL = aXP;
1168       aYL = aYP;
1169       aXP = aA + (aB - aA)*aGS;
1170       iErr = DistPC(aXP, theC1, theCriteria, theProjPC, 
1171                     aYP, aT2P, aDmax, aT1max, aT2max, iC);
1172     }
1173     else {
1174       aB = aXP;
1175       aXP = aXL;
1176       aYP = aYL;
1177       aXL = aB - (aB - aA)*aGS;
1178       iErr = DistPC(aXL, theC1, theCriteria, theProjPC, 
1179                     aYL, aT2L, aDmax, aT1max, aT2max, iC);
1180     }
1181     //
1182     if (iErr) {
1183       if ((iErr == 2) && !bMaxDist) {
1184         aXP = (aA + aB) * 0.5;
1185         DistPC(aXP, theC1, theCriteria, theProjPC, 
1186                aYP, aT2P, aDmax, aT1max, aT2max, iC);
1187       }
1188       return iErr;
1189     }
1190     //
1191     if ((aB - aA) < anEps) {
1192       break;
1193     }
1194   }// for (;;) {
1195   //
1196   return iErr;
1197 }
1198 //=======================================================================
1199 //function : DistPC
1200 //purpose  : 
1201 //=======================================================================
1202 Standard_Integer DistPC(const Standard_Real aT1, 
1203                         const Handle(Geom_Curve)& theC1,
1204                         const Standard_Real theCriteria,
1205                         GeomAPI_ProjectPointOnCurve& theProjPC, 
1206                         Standard_Real& aD, 
1207                         Standard_Real& aT2,
1208                         Standard_Real& aDmax,
1209                         Standard_Real& aT1max,
1210                         Standard_Real& aT2max,
1211                         const Standard_Integer iC)
1212 {
1213   Standard_Integer iErr;
1214   //
1215   iErr = DistPC(aT1, theC1, theCriteria, theProjPC, aD, aT2, iC);
1216   if (iErr == 1) {
1217     return iErr;
1218   }
1219   //
1220   if (iC*(aD - aDmax) > 0) {
1221     aDmax = aD;
1222     aT1max = aT1;
1223     aT2max = aT2;
1224   }
1225   //
1226   return iErr;
1227 }
1228 //=======================================================================
1229 //function : DistPC
1230 //purpose  : 
1231 //=======================================================================
1232 Standard_Integer DistPC(const Standard_Real aT1, 
1233                         const Handle(Geom_Curve)& theC1,
1234                         const Standard_Real theCriteria, 
1235                         GeomAPI_ProjectPointOnCurve& theProjPC,
1236                         Standard_Real& aD, 
1237                         Standard_Real& aT2,
1238                         const Standard_Integer iC) 
1239 {
1240   Standard_Integer iErr, aNbP2;
1241   gp_Pnt aP1;
1242   //
1243   iErr = 0;
1244   theC1->D0(aT1, aP1);
1245   //
1246   theProjPC.Perform(aP1);
1247   aNbP2 = theProjPC.NbPoints();
1248   if (!aNbP2) {
1249     iErr = 1;// the point from aC1 can not be projected on aC2
1250     return iErr;
1251   }
1252   //
1253   aD  = theProjPC.LowerDistance();
1254   aT2 = theProjPC.LowerDistanceParameter();
1255   if (iC*(aD - theCriteria) > 0) {
1256     iErr = 2;// the distance is too big or small
1257   }
1258   //
1259   return iErr;
1260 }
1261
1262 //=======================================================================
1263 //function : SplitRangeOnSegments
1264 //purpose  : 
1265 //=======================================================================
1266 Standard_Integer SplitRangeOnSegments(const Standard_Real aT1, 
1267                                       const Standard_Real aT2,
1268                                       const Standard_Real theResolution,
1269                                       const Standard_Integer theNbSeg,
1270                                       IntTools_SequenceOfRanges& theSegments)
1271 {
1272   Standard_Real aDiff = aT2 - aT1;
1273   if (aDiff < theResolution || theNbSeg == 1) {
1274     theSegments.Append(IntTools_Range(aT1, aT2));
1275     return 1;
1276   }
1277   //
1278   Standard_Real aDt, aT1x, aT2x, aSeg;
1279   Standard_Integer aNbSegments, i;
1280   //
1281   aNbSegments = theNbSeg;
1282   aDt = aDiff / aNbSegments;
1283   if (aDt < theResolution) {
1284     aSeg = aDiff / theResolution;
1285     aNbSegments = Standard_Integer(aSeg) + 1;
1286     aDt = aDiff / aNbSegments;
1287   }
1288   //
1289   aT1x = aT1;
1290   for (i = 1; i < aNbSegments; ++i) {
1291     aT2x = aT1x + aDt;
1292     //
1293     IntTools_Range aR(aT1x, aT2x);
1294     theSegments.Append(aR);
1295     //
1296     aT1x = aT2x;
1297   }
1298   //
1299   IntTools_Range aR(aT1x, aT2);
1300   theSegments.Append(aR);
1301   //
1302   return aNbSegments;
1303 }
1304
1305 //=======================================================================
1306 //function : BndBuildBox
1307 //purpose  : 
1308 //=======================================================================
1309 void BndBuildBox(const BRepAdaptor_Curve& theBAC,
1310                  const Standard_Real aT1,
1311                  const Standard_Real aT2,
1312                  const Standard_Real theTol,
1313                  Bnd_Box& theBox)
1314 {
1315   Bnd_Box aB;
1316   BndLib_Add3dCurve::Add(theBAC, aT1, aT2, theTol, aB);
1317   theBox = aB;
1318 }
1319
1320 //=======================================================================
1321 //function : PointBoxDistance
1322 //purpose  : 
1323 //=======================================================================
1324 Standard_Real PointBoxDistance(const Bnd_Box& aB,
1325                                const gp_Pnt& aP)
1326 {
1327   Standard_Real aPCoord[3];
1328   Standard_Real aBMinCoord[3], aBMaxCoord[3];
1329   Standard_Real aDist, aR1, aR2;
1330   Standard_Integer i;
1331   //
1332   aP.Coord(aPCoord[0], aPCoord[1], aPCoord[2]);
1333   aB.Get(aBMinCoord[0], aBMinCoord[1], aBMinCoord[2], 
1334          aBMaxCoord[0], aBMaxCoord[1], aBMaxCoord[2]);
1335   //
1336   aDist = 0.;
1337   for (i = 0; i < 3; ++i) {
1338     aR1 = aBMinCoord[i] - aPCoord[i];
1339     if (aR1 > 0.) {
1340       aDist += aR1*aR1;
1341       continue;
1342     }
1343     //
1344     aR2 = aPCoord[i] - aBMaxCoord[i];
1345     if (aR2 > 0.) {
1346       aDist += aR2*aR2;
1347     }
1348   }
1349   //
1350   aDist = Sqrt(aDist);
1351   return aDist;
1352 }
1353
1354 //=======================================================================
1355 //function : TypeToInteger
1356 //purpose  : 
1357 //=======================================================================
1358 Standard_Integer TypeToInteger(const GeomAbs_CurveType theCType)
1359 {
1360   Standard_Integer iRet;
1361   //
1362   switch(theCType) {
1363   case GeomAbs_Line:
1364     iRet=0;
1365     break;
1366   case GeomAbs_Hyperbola:
1367   case GeomAbs_Parabola:
1368     iRet=1;
1369     break;
1370   case GeomAbs_Circle:
1371   case GeomAbs_Ellipse:
1372     iRet=2;
1373     break;
1374   case GeomAbs_BezierCurve:
1375   case GeomAbs_BSplineCurve:
1376     iRet=3;
1377     break;
1378   default:
1379     iRet=4;
1380     break;
1381   }
1382   return iRet;
1383 }
1384
1385 //=======================================================================
1386 //function : ResolutionCoeff
1387 //purpose  : 
1388 //=======================================================================
1389 Standard_Real ResolutionCoeff(const BRepAdaptor_Curve& theBAC,
1390                               const IntTools_Range& theRange)
1391 {
1392   Standard_Real aResCoeff = 0.;
1393   //
1394   const Handle(Geom_Curve)& aCurve = theBAC.Curve().Curve();
1395   const GeomAbs_CurveType aCurveType = theBAC.GetType();
1396   //
1397   switch (aCurveType) {
1398   case GeomAbs_Circle :
1399     aResCoeff = 1. / (2 * Handle(Geom_Circle)::DownCast (aCurve)->Circ().Radius());
1400     break;
1401   case GeomAbs_Ellipse :
1402     aResCoeff =  1. / Handle(Geom_Ellipse)::DownCast (aCurve)->MajorRadius();
1403     break;
1404   case GeomAbs_OffsetCurve : {
1405     const Handle(Geom_OffsetCurve)& anOffsetCurve = Handle(Geom_OffsetCurve)::DownCast(aCurve);
1406     const Handle(Geom_Curve)& aBasisCurve = anOffsetCurve->BasisCurve();
1407     GeomAdaptor_Curve aGBasisCurve(aBasisCurve);
1408     const GeomAbs_CurveType aBCType = aGBasisCurve.GetType();
1409     if (aBCType == GeomAbs_Line) {
1410       break;
1411     }
1412     else if (aBCType == GeomAbs_Circle) {
1413       aResCoeff = 1. / (2 * (anOffsetCurve->Offset() + aGBasisCurve.Circle().Radius()));
1414       break;
1415     }
1416     else if (aBCType == GeomAbs_Ellipse) {
1417       aResCoeff = 1. / (anOffsetCurve->Offset() + aGBasisCurve.Ellipse().MajorRadius());
1418       break;
1419     }
1420   }
1421   Standard_FALLTHROUGH
1422   case GeomAbs_Hyperbola :
1423   case GeomAbs_Parabola : 
1424   case GeomAbs_OtherCurve :{
1425     Standard_Real k, kMin, aDist, aDt, aT1, aT2, aT;
1426     Standard_Integer aNbP, i;
1427     gp_Pnt aP1, aP2;
1428     //
1429     aNbP = 30;
1430     theRange.Range(aT1, aT2);
1431     aDt = (aT2 - aT1) / aNbP;
1432     aT = aT1;
1433     kMin = 10.;
1434     //
1435     theBAC.D0(aT1, aP1);
1436     for (i = 1; i <= aNbP; ++i) {
1437       aT += aDt;
1438       theBAC.D0(aT, aP2);
1439       aDist = aP1.Distance(aP2);
1440       k = aDt / aDist;
1441       if (k < kMin) {
1442         kMin = k;
1443       }
1444       aP1 = aP2;
1445     }
1446     //
1447     aResCoeff = kMin;
1448     break;
1449   }
1450   default:
1451     break;
1452   }
1453   //
1454   return aResCoeff;
1455 }
1456
1457 //=======================================================================
1458 //function : Resolution
1459 //purpose  : 
1460 //=======================================================================
1461 Standard_Real Resolution(const Handle(Geom_Curve)& theCurve,
1462                          const GeomAbs_CurveType theCurveType,
1463                          const Standard_Real theResCoeff,
1464                          const Standard_Real theR3D)
1465 {
1466   Standard_Real aRes;
1467   //
1468   switch (theCurveType) {
1469   case GeomAbs_Line :
1470     aRes = theR3D;
1471     break;
1472   case GeomAbs_Circle: {
1473     Standard_Real aDt = theResCoeff * theR3D;
1474     aRes = (aDt <= 1.) ? 2*ASin(aDt) : 2*M_PI;
1475     break;
1476   }
1477   case GeomAbs_BezierCurve:
1478     Handle(Geom_BezierCurve)::DownCast (theCurve)->Resolution(theR3D, aRes);
1479     break;
1480   case GeomAbs_BSplineCurve:
1481     Handle(Geom_BSplineCurve)::DownCast (theCurve)->Resolution(theR3D, aRes);
1482     break;
1483   case GeomAbs_OffsetCurve: {
1484     const Handle(Geom_Curve)& aBasisCurve = 
1485       Handle(Geom_OffsetCurve)::DownCast(theCurve)->BasisCurve();
1486     const GeomAbs_CurveType aBCType = GeomAdaptor_Curve(aBasisCurve).GetType();
1487     if (aBCType == GeomAbs_Line) {
1488       aRes = theR3D;
1489       break;
1490     }
1491     else if (aBCType == GeomAbs_Circle) {
1492       Standard_Real aDt = theResCoeff * theR3D;
1493       aRes = (aDt <= 1.) ? 2*ASin(aDt) : 2*M_PI;
1494       break;
1495     }
1496   }
1497   Standard_FALLTHROUGH
1498   default:
1499     aRes = theResCoeff * theR3D;
1500     break;
1501   }
1502   //
1503   return aRes;
1504 }
1505
1506 //=======================================================================
1507 //function : CurveDeflection
1508 //purpose  : 
1509 //=======================================================================
1510 Standard_Real CurveDeflection(const BRepAdaptor_Curve& theBAC,
1511                               const IntTools_Range& theRange)
1512 {
1513   Standard_Real aDt, aT, aT1, aT2, aDefl;
1514   Standard_Integer i, aNbP;
1515   gp_Vec aV1, aV2;
1516   gp_Pnt aP;
1517   //
1518   aDefl = 0;
1519   aNbP = 10;
1520   theRange.Range(aT1, aT2);
1521   aDt = (aT2 - aT1) / aNbP;
1522   aT = aT1;
1523   //
1524   theBAC.D1(aT1, aP, aV1);
1525   for (i = 1; i <= aNbP; ++i) {
1526     aT += aDt;
1527     theBAC.D1(aT, aP, aV2);
1528     if (aV1.Magnitude() > gp::Resolution() &&
1529         aV2.Magnitude() > gp::Resolution()) {
1530       gp_Dir aD1(aV1), aD2(aV2);
1531       aDefl += aD1.Angle(aD2);
1532     }
1533     aV1 = aV2;
1534   }
1535   //
1536   return aDefl;
1537 }
1538
1539 //=======================================================================
1540 //function : IsClosed
1541 //purpose  : 
1542 //=======================================================================
1543 Standard_Boolean IsClosed(const Handle(Geom_Curve)& theCurve,
1544                           const Standard_Real aT1,
1545                           const Standard_Real aT2,
1546                           const Standard_Real theTol,
1547                           const Standard_Real theRes)
1548 {
1549   if (Abs(aT1 - aT2) < theRes)
1550   {
1551     return Standard_False;
1552   }
1553
1554   gp_Pnt aP1, aP2;
1555   theCurve->D0(aT1, aP1);
1556   theCurve->D0(aT2, aP2);
1557   //
1558   Standard_Real aD = aP1.Distance(aP2);
1559   return aD < theTol;
1560 }