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