ed7dd9010473da679f639bf63be727f00a581679
[occt.git] / src / Extrema / Extrema_ExtCC.cxx
1 // Created on: 1994-07-06
2 // Created by: Laurent PAINNOT
3 // Copyright (c) 1994-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 // Modified by MPS (june 96) : correction du trap dans le cas droite/Bezier 
18 // Modified by MPS (mai 97) : PRO 7598 
19 //                            tri des solutions pour eviter de rendre plusieurs
20 //                            fois la meme solution 
21
22 #include <Adaptor3d_Curve.hxx>
23 #include <Bnd_Range.hxx>
24 #include <ElCLib.hxx>
25 #include <Extrema_CurveTool.hxx>
26 #include <Extrema_ECC.hxx>
27 #include <Extrema_ExtCC.hxx>
28 #include <Extrema_ExtElC.hxx>
29 #include <Extrema_ExtPElC.hxx>
30 #include <Extrema_POnCurv.hxx>
31 #include <Extrema_SequenceOfPOnCurv.hxx>
32 #include <Geom_Circle.hxx>
33 #include <Geom_Curve.hxx>
34 #include <Geom_Ellipse.hxx>
35 #include <Geom_Hyperbola.hxx>
36 #include <Geom_Line.hxx>
37 #include <Geom_Parabola.hxx>
38 #include <Geom_TrimmedCurve.hxx>
39 #include <GeomAbs_CurveType.hxx>
40 #include <gp_Pnt.hxx>
41 #include <Precision.hxx>
42 #include <Standard_Failure.hxx>
43 #include <Standard_NotImplemented.hxx>
44 #include <Standard_NullObject.hxx>
45 #include <Standard_OutOfRange.hxx>
46 #include <StdFail_NotDone.hxx>
47 #include <TColStd_Array1OfReal.hxx>
48 #include <TColStd_ListIteratorOfListOfTransient.hxx>
49 #include <TColStd_SequenceOfReal.hxx>
50
51 //=======================================================================
52 //function : Extrema_ExtCC
53 //purpose  : 
54 //=======================================================================
55 Extrema_ExtCC::Extrema_ExtCC (const Standard_Real TolC1,
56                               const Standard_Real TolC2)
57 : myIsFindSingleSolution(Standard_False),
58   myDone (Standard_False)
59 {
60   myC[0] = 0; myC[1] = 0;
61   myInf[0] = myInf[1] = -Precision::Infinite();
62   mySup[0] = mySup[1] = Precision::Infinite();
63   myTol[0] = TolC1; myTol[1] = TolC2;
64   mydist11 = mydist12 = mydist21 = mydist22 = RealFirst();
65 }
66
67 //=======================================================================
68 //function : Extrema_ExtCC
69 //purpose  : 
70 //=======================================================================
71
72 Extrema_ExtCC::Extrema_ExtCC(const Adaptor3d_Curve& C1,
73                              const Adaptor3d_Curve& C2,
74                              const Standard_Real      U1,
75                              const Standard_Real      U2,
76                              const Standard_Real      V1,
77                              const Standard_Real      V2,
78                              const Standard_Real      TolC1,
79                              const Standard_Real      TolC2)
80 : myIsFindSingleSolution(Standard_False),
81   myECC(C1, C2, U1, U2, V1, V2),
82   myDone (Standard_False)
83 {
84   SetCurve (1, C1, U1, U2);
85   SetCurve (2, C2, V1, V2);
86   SetTolerance (1, TolC1);
87   SetTolerance (2, TolC2);
88   mydist11 = mydist12 = mydist21 = mydist22 = RealFirst();
89   Perform();
90 }
91
92
93 //=======================================================================
94 //function : Extrema_ExtCC
95 //purpose  : 
96 //=======================================================================
97
98 Extrema_ExtCC::Extrema_ExtCC(const Adaptor3d_Curve& C1, 
99                              const Adaptor3d_Curve& C2,
100                              const Standard_Real      TolC1,
101                              const Standard_Real      TolC2)
102 : myIsFindSingleSolution(Standard_False),
103   myECC(C1, C2),
104   myDone (Standard_False)
105 {
106   SetCurve (1, C1, C1.FirstParameter(), C1.LastParameter());
107   SetCurve (2, C2, C2.FirstParameter(), C2.LastParameter());
108   SetTolerance (1, TolC1);
109   SetTolerance (2, TolC2);
110   mydist11 = mydist12 = mydist21 = mydist22 = RealFirst();
111   Perform();
112 }
113
114 //=======================================================================
115 //function : SetCurve
116 //purpose  : 
117 //=======================================================================
118
119 void Extrema_ExtCC::SetCurve (const Standard_Integer theRank, const Adaptor3d_Curve& C)
120 {
121   Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_ExtCC::SetCurve()")
122   Standard_Integer anInd = theRank - 1;
123   myC[anInd] = (Standard_Address)&C;
124 }
125
126 //=======================================================================
127 //function : SetCurve
128 //purpose  : 
129 //=======================================================================
130
131 void Extrema_ExtCC::SetCurve (const Standard_Integer theRank, const Adaptor3d_Curve& C,
132                                const Standard_Real Uinf, const Standard_Real Usup)
133 {
134   SetCurve (theRank, C);
135   SetRange (theRank, Uinf, Usup);
136 }
137
138 //=======================================================================
139 //function : SetRange
140 //purpose  : 
141 //=======================================================================
142
143 void Extrema_ExtCC::SetRange (const Standard_Integer theRank, 
144                                const Standard_Real Uinf, const Standard_Real Usup)
145 {
146   Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_ExtCC::SetRange()")
147   Standard_Integer anInd = theRank - 1;
148   myInf[anInd] = Uinf;
149   mySup[anInd] = Usup;
150 }
151
152 //=======================================================================
153 //function : SetTolerance
154 //purpose  : 
155 //=======================================================================
156
157 void Extrema_ExtCC::SetTolerance (const Standard_Integer theRank, const Standard_Real theTol)
158 {
159   Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_ExtCC::SetTolerance()")
160   Standard_Integer anInd = theRank - 1;
161   myTol[anInd] = theTol;
162 }
163
164
165 //=======================================================================
166 //function : Perform
167 //purpose  : 
168 //=======================================================================
169
170 void Extrema_ExtCC::Perform()
171 {  
172   Standard_NullObject_Raise_if (!myC[0] || !myC[1], "Extrema_ExtCC::Perform()")
173   myECC.SetParams(*((Adaptor3d_Curve*)myC[0]), 
174                   *((Adaptor3d_Curve*)myC[1]), myInf[0], mySup[0], myInf[1], mySup[1]);
175   myECC.SetTolerance(Min(myTol[0], myTol[1]));
176   myECC.SetSingleSolutionFlag(GetSingleSolutionFlag());
177   myDone = Standard_False;
178   mypoints.Clear();
179   mySqDist.Clear();
180   myIsPar = Standard_False;
181
182   GeomAbs_CurveType type1 = (*((Adaptor3d_Curve*)myC[0])).GetType();
183   GeomAbs_CurveType type2 = (*((Adaptor3d_Curve*)myC[1])).GetType();
184   Standard_Real U11, U12, U21, U22, Tol = Min(myTol[0], myTol[1]);
185
186   U11 = myInf[0];
187   U12 = mySup[0];
188   U21 = myInf[1];
189   U22 = mySup[1];
190
191   if (!Precision::IsInfinite(U11)) P1f = Extrema_CurveTool::Value(*((Adaptor3d_Curve*)myC[0]),  U11); 
192   if (!Precision::IsInfinite(U12)) P1l = Extrema_CurveTool::Value(*((Adaptor3d_Curve*)myC[0]),  U12);
193   if (!Precision::IsInfinite(U21)) P2f = Extrema_CurveTool::Value(*((Adaptor3d_Curve*)myC[1]), U21);
194   if (!Precision::IsInfinite(U22)) P2l = Extrema_CurveTool::Value(*((Adaptor3d_Curve*)myC[1]), U22);
195   
196
197   if (Precision::IsInfinite(U11) || Precision::IsInfinite(U21)) mydist11 = RealLast();
198   else mydist11 = P1f.SquareDistance(P2f);
199   if (Precision::IsInfinite(U11) || Precision::IsInfinite(U22)) mydist12 = RealLast();
200   else mydist12 = P1f.SquareDistance(P2l);
201   if (Precision::IsInfinite(U12) || Precision::IsInfinite(U21)) mydist21 = RealLast();
202   else mydist21 = P1l.SquareDistance(P2f);
203   if (Precision::IsInfinite(U12) || Precision::IsInfinite(U22)) mydist22 = RealLast();
204   else mydist22 = P1l.SquareDistance(P2l);
205
206   //Depending on the types of curves, the algorithm is chosen:
207   //- _ExtElC, when one of the curves is a line and the other is elementary,
208   //   or there are two circles;
209   //- _GenExtCC, in all other cases
210   if ( (type1 == GeomAbs_Line && type2 <= GeomAbs_Parabola) ||
211        (type2 == GeomAbs_Line && type1 <= GeomAbs_Parabola) ) {
212     //analytical case - one curve is always a line
213     Standard_Integer anInd1 = 0, anInd2 = 1;
214     GeomAbs_CurveType aType2 = type2;
215     Standard_Boolean isInverse = (type1 > type2);
216     if (isInverse)
217     {
218       //algorithm uses inverse order of arguments
219       anInd1 = 1;
220       anInd2 = 0;
221       aType2 = type1;
222     }
223     switch (aType2) {
224     case GeomAbs_Line: {
225       Extrema_ExtElC Xtrem((*((Adaptor3d_Curve*)myC[anInd1])).Line(), (*((Adaptor3d_Curve*)myC[anInd2])).Line(), Tol);
226       PrepareResults(Xtrem, isInverse, U11, U12, U21, U22);
227       break;
228     }
229     case GeomAbs_Circle: {
230       Extrema_ExtElC Xtrem((*((Adaptor3d_Curve*)myC[anInd1])).Line(), (*((Adaptor3d_Curve*)myC[anInd2])).Circle(), Tol);
231       PrepareResults(Xtrem, isInverse, U11, U12, U21, U22);
232       break;
233     }
234     case GeomAbs_Ellipse: {
235       Extrema_ExtElC Xtrem((*((Adaptor3d_Curve*)myC[anInd1])).Line(), (*((Adaptor3d_Curve*)myC[anInd2])).Ellipse());
236       PrepareResults(Xtrem, isInverse, U11, U12, U21, U22);
237       break;
238     }
239     case GeomAbs_Hyperbola: {
240       Extrema_ExtElC Xtrem((*((Adaptor3d_Curve*)myC[anInd1])).Line(), (*((Adaptor3d_Curve*)myC[anInd2])).Hyperbola());
241       PrepareResults(Xtrem, isInverse, U11, U12, U21, U22);
242       break;
243     }
244     case GeomAbs_Parabola: {
245       Extrema_ExtElC Xtrem((*((Adaptor3d_Curve*)myC[anInd1])).Line(), (*((Adaptor3d_Curve*)myC[anInd2])).Parabola());
246       PrepareResults(Xtrem, isInverse, U11, U12, U21, U22);
247       break;
248     }
249     default: break;
250     }
251   } else if (type1 == GeomAbs_Circle && type2 == GeomAbs_Circle) {
252     //analytical case - two circles
253     Standard_Boolean bIsDone;
254     Extrema_ExtElC CCXtrem ((*((Adaptor3d_Curve*)myC[0])).Circle(), (*((Adaptor3d_Curve*)myC[1])).Circle());
255     bIsDone = CCXtrem.IsDone();
256     if(bIsDone) {
257       PrepareResults(CCXtrem, Standard_False, U11, U12, U21, U22);
258     }
259     else {
260       myECC.Perform();
261       PrepareResults(myECC, U11, U12, U21, U22);
262     }
263   } else {
264     myECC.Perform();
265     PrepareResults(myECC, U11, U12, U21, U22);
266   }
267 }
268
269
270 //=======================================================================
271 //function : IsDone
272 //purpose  : 
273 //=======================================================================
274
275 Standard_Boolean Extrema_ExtCC::IsDone() const
276 {
277   return myDone;
278 }
279
280 //=======================================================================
281 //function : IsParallel
282 //purpose  : 
283 //=======================================================================
284
285 Standard_Boolean Extrema_ExtCC::IsParallel() const
286 {
287   if (!IsDone())
288   {
289     throw StdFail_NotDone();
290   }
291
292   return myIsPar;
293 }
294
295
296 //=======================================================================
297 //function : Value
298 //purpose  : 
299 //=======================================================================
300
301 Standard_Real Extrema_ExtCC::SquareDistance(const Standard_Integer N) const 
302 {
303   if ((N < 1) || (N > NbExt())) throw Standard_OutOfRange();
304   return mySqDist.Value(N);
305 }
306
307
308 //=======================================================================
309 //function : NbExt
310 //purpose  : 
311 //=======================================================================
312
313 Standard_Integer Extrema_ExtCC::NbExt() const
314 {
315   if(!myDone) throw StdFail_NotDone();
316   return mySqDist.Length();
317 }
318
319
320 //=======================================================================
321 //function : Points
322 //purpose  : 
323 //=======================================================================
324
325 void Extrema_ExtCC::Points(const Standard_Integer N, 
326                             Extrema_POnCurv& P1,
327                             Extrema_POnCurv& P2) const
328 {
329   if (N < 1 || N > NbExt())
330   {
331     throw Standard_OutOfRange();
332   }
333
334   P1 = mypoints.Value(2 * N - 1);
335   P2 = mypoints.Value(2 * N);
336 }
337
338
339
340 //=======================================================================
341 //function : TrimmedDistances
342 //purpose  : 
343 //=======================================================================
344
345 void Extrema_ExtCC::TrimmedSquareDistances(Standard_Real& dist11,
346                                       Standard_Real& dist12,
347                                       Standard_Real& dist21,
348                                       Standard_Real& dist22,
349                                       gp_Pnt&        P11   ,
350                                       gp_Pnt&        P12   ,
351                                       gp_Pnt&        P21   ,
352                                       gp_Pnt&        P22    ) const {
353                                         
354   dist11 = mydist11;
355   dist12 = mydist12;
356   dist21 = mydist21;
357   dist22 = mydist22;
358   P11 = P1f;
359   P12 = P1l;
360   P21 = P2f;
361   P22 = P2l;
362 }
363
364 //=======================================================================
365 //function : ParallelResult
366 //purpose  : 
367 //=======================================================================
368 void Extrema_ExtCC::PrepareParallelResult(const Standard_Real theUt11,
369                                           const Standard_Real theUt12,
370                                           const Standard_Real theUt21,
371                                           const Standard_Real theUt22,
372                                           const Standard_Real theSqDist)
373 {
374   if (!myIsPar)
375     return;
376
377   const GeomAbs_CurveType aType1 = Extrema_CurveTool::GetType(*((Adaptor3d_Curve*) myC[0]));
378   const GeomAbs_CurveType aType2 = Extrema_CurveTool::GetType(*((Adaptor3d_Curve*) myC[1]));
379   
380   if (((aType1 != GeomAbs_Line) && (aType1 != GeomAbs_Circle)) ||
381       ((aType2 != GeomAbs_Line) && (aType2 != GeomAbs_Circle)))
382   {
383     mySqDist.Append(theSqDist);
384     myDone = Standard_True;
385     myIsPar = Standard_True;
386     return;
387   }
388   
389   // Parallel case is only for line-line, circle-circle and circle-line!!!
390   // But really for trimmed curves extremas can not exist!
391   if (aType1 != aType2)
392   {
393     //The projection of the circle's location to the trimmed line must exist.
394     const Standard_Boolean isReversed = (aType1 != GeomAbs_Circle);
395     const gp_Pnt aPonC = !isReversed ?
396                       Extrema_CurveTool::Value(*((Adaptor3d_Curve*) myC[0]), theUt11) :
397                       Extrema_CurveTool::Value(*((Adaptor3d_Curve*) myC[1]), theUt21);
398
399     const gp_Lin aL = !isReversed ? ((Adaptor3d_Curve*) myC[1])->Line() :
400                                     ((Adaptor3d_Curve*) myC[0])->Line();
401     const Extrema_ExtPElC ExtPLin(aPonC, aL, Precision::Confusion(),
402                                   !isReversed ? theUt21 : theUt11,
403                                   !isReversed ? theUt22 : theUt12);
404
405     if (ExtPLin.IsDone())
406     {
407       mySqDist.Append(theSqDist);
408     }
409     else
410     {
411       myIsPar = Standard_False;
412     }
413
414     return;
415   }
416
417   if (aType1 == GeomAbs_Line)
418   {
419     // Line - Line
420
421     const Standard_Real isFirstInfinite = (Precision::IsInfinite(theUt11) &&
422                                            Precision::IsInfinite(theUt12));
423     const Standard_Real isLastInfinite = (Precision::IsInfinite(theUt21) &&
424                                           Precision::IsInfinite(theUt22));
425
426     if (isFirstInfinite || isLastInfinite)
427     {
428       // Infinite number of solution
429
430       mySqDist.Append(theSqDist);
431     }
432     else
433     {
434       // The range created by projection of both ends of the 1st line
435       // to the 2nd one must intersect the (native) trimmed range of
436       // the 2nd line.
437
438       myIsPar = Standard_False;
439
440       const gp_Lin aLin1 = ((Adaptor3d_Curve*) myC[0])->Line();
441       const gp_Lin aLin2 = ((Adaptor3d_Curve*) myC[1])->Line();
442       const Standard_Boolean isOpposite(aLin1.Direction().Dot(aLin2.Direction()) < 0.0);
443
444       Bnd_Range aRange2(theUt21, theUt22);
445       Bnd_Range aProjRng12;
446
447       if (Precision::IsInfinite(theUt11))
448       {
449         if (isOpposite)
450           aProjRng12.Add(Precision::Infinite());
451         else
452           aProjRng12.Add(-Precision::Infinite());
453       }
454       else
455       {
456         const gp_Pnt aPonC1 = ElCLib::Value(theUt11, aLin1);
457         const Standard_Real aPar = ElCLib::Parameter(aLin2, aPonC1);
458         aProjRng12.Add(aPar);
459       }
460
461       if (Precision::IsInfinite(theUt12))
462       {
463         if (isOpposite)
464           aProjRng12.Add(-Precision::Infinite());
465         else
466           aProjRng12.Add(Precision::Infinite());
467       }
468       else
469       {
470         const gp_Pnt aPonC1 = ElCLib::Value(theUt12, aLin1);
471         const Standard_Real aPar = ElCLib::Parameter(aLin2, aPonC1);
472         aProjRng12.Add(aPar);
473       }
474
475       aRange2.Common(aProjRng12);
476       if (aRange2.Delta() > Precision::Confusion())
477       {
478         ClearSolutions();
479         mySqDist.Append(theSqDist);
480         myIsPar = Standard_True;
481       }
482       else if (!aRange2.IsVoid())
483       {
484         //Case like this:
485
486         //  **************     aLin1
487         //               o
488         //               o
489         //               ***************  aLin2
490
491         ClearSolutions();
492         Standard_Real aPar1 = 0.0, aPar2 = 0.0;
493         aRange2.GetBounds(aPar1, aPar2);
494         aPar2 = 0.5*(aPar1 + aPar2);
495         gp_Pnt aP = ElCLib::Value(aPar2, aLin2);
496         const Extrema_POnCurv aP2(aPar2, aP);
497         aPar1 = ElCLib::Parameter(aLin1, aP);
498         aP = ElCLib::Value(aPar1, aLin1);
499         const Extrema_POnCurv aP1(aPar1, aP);
500         mypoints.Append(aP1);
501         mypoints.Append(aP2);
502         mySqDist.Append(theSqDist);
503       }
504     }
505   }
506   else
507   {
508     // Circle - Circle
509     myIsPar = Standard_False;
510
511     //Two arcs with ranges [U1, U2] and [V1, V2] correspondingly are
512     //considered to be parallel in the following case:
513     //  The range created by projection both points U1 and U2 of the
514     //  1st circle to the 2nd one intersects either the range [V1, V2] or
515     //  the range [V1-PI, V2-PI]. All ranges must be adjusted to correspond
516     //  periodic range before checking of intersection.
517
518     const gp_Circ aWorkCirc = ((Adaptor3d_Curve*) myC[1])->Circle();
519     const Standard_Real aPeriod = M_PI + M_PI;
520     gp_Vec aVTg1;
521     gp_Pnt aP11;
522     const gp_Pnt aP12 = Extrema_CurveTool::Value(*((Adaptor3d_Curve*) myC[0]), theUt12);
523     Extrema_CurveTool::D1(*((Adaptor3d_Curve*) myC[0]), theUt11, aP11, aVTg1);
524
525     const Bnd_Range aRange(theUt21, theUt22);
526     Bnd_Range aProjRng1;
527
528     // Project arc of the 1st circle between points theUt11 and theUt12 to the
529     // 2nd circle. It is necessary to chose correct arc from two possible ones.
530
531     Standard_Real aPar1 = ElCLib::InPeriod(ElCLib::Parameter(aWorkCirc, aP11),
532                                            theUt21, theUt21 + aPeriod);
533     const gp_Vec aVTg2 = Extrema_CurveTool::DN(*((Adaptor3d_Curve*) myC[1]), aPar1, 1);
534     
535     // Check if circles have same/opposite directions
536     const Standard_Boolean isOpposite(aVTg1.Dot(aVTg2) < 0.0);
537
538     Standard_Real aPar2 = ElCLib::InPeriod(ElCLib::Parameter(aWorkCirc, aP12),
539                                            theUt21, theUt21 + aPeriod);
540
541     if (isOpposite)
542     {
543       // Must be aPar2 < aPar1
544       if ((aRange.Delta() > Precision::Angular()) &&
545           ((aPar1 - aPar2) < Precision::Angular()))
546       {
547         aPar2 -= aPeriod;
548       }
549     }
550     else
551     {
552       // Must be aPar2 > aPar1
553       if ((aRange.Delta() > Precision::Angular()) &&
554           ((aPar2 - aPar1) < Precision::Angular()))
555       {
556         aPar1 -= aPeriod;
557       }
558     }
559
560     // Now the projection result is the range [aPar1, aPar2]
561     // if aPar1 < aPar2 or the range [aPar2, aPar1], otherwise.
562
563     Standard_Real aMinSquareDist = RealLast();
564
565     aProjRng1.Add(aPar1 - M_PI);
566     aProjRng1.Add(aPar2 - M_PI);
567     for (Standard_Integer i = 0; i < 2; i++)
568     {
569       // Repeat computation twice
570
571       Bnd_Range aRng = aProjRng1;
572       aRng.Common(aRange);
573
574       //Cases are possible and processed below:
575       //1. Extrema does not exist. In this case all common ranges are VOID.
576       //2. Arcs are parallel and distance between them is equal to sqrt(theSqDist).
577       //    In this case myIsPar = TRUE definitely.
578       //3. Arcs are parallel and distance between them is equal to (sqrt(theSqDist) + R),
579       //    where R is the least radius of the both circles. In this case myIsPar flag
580       //    will temporary be set to TRUE but check will be continued until less
581       //    distance will be found. At that, region with the least distance can be
582       //    either a local point or continuous range. In 1st case myIsPar = FALSE and
583       //    several (or single) extremas will be returned. In the 2nd one
584       //    myIsPar = TRUE and only the least distance will be returned.
585       //4. Arcs are not parallel. Then several (or single) extremas will be returned.
586
587       if (aRng.Delta() > Precision::Angular())
588       {
589         Standard_Real aPar = 0.0;
590         aRng.GetIntermediatePoint(0.5, aPar);
591         const gp_Pnt aPCirc2 = ElCLib::Value(aPar, aWorkCirc);
592         Extrema_ExtPElC ExtPCir(aPCirc2,
593                                 Extrema_CurveTool::Circle(*((Adaptor3d_Curve*) myC[0])),
594                                 Precision::Confusion(), theUt11, theUt12);
595
596         Standard_Real aMinSqD = ExtPCir.SquareDistance(1);
597         for (Standard_Integer anExtID = 2; anExtID <= ExtPCir.NbExt(); anExtID++)
598         {
599           aMinSqD = Min(aMinSqD, ExtPCir.SquareDistance(anExtID));
600         }
601
602         if (aMinSqD <= aMinSquareDist)
603         {
604           ClearSolutions();
605           mySqDist.Append(aMinSqD);
606           myIsPar = Standard_True;
607
608           const Standard_Real aDeltaSqDist = aMinSqD - theSqDist;
609           const Standard_Real aSqD = Max(aMinSqD, theSqDist);
610
611           //  0 <= Dist1-Dist2 <= Eps
612           //  0 <= Dist1^2 - Dist2^2 < Eps*(Dist1+Dist2)
613
614           //If Dist1 ~= Dist2 ==> Dist1+Dist2 ~= 2*Dist2.
615           //Consequently,
616           //  0 <= Dist1^2 - Dist2^2 <= 2*Dist2*Eps
617
618           //Or
619           //  (Dist1^2 - Dist2^2)^2 <= 4*Dist2^2*Eps^2
620
621           if (aDeltaSqDist*aDeltaSqDist < 4.0*aSqD*Precision::SquareConfusion())
622           {
623             // New solution is found
624             break;
625           }
626         }
627
628         //Nearer solution can be found
629       }
630       else if (!aRng.IsVoid())
631       {
632         //Check cases like this:
633
634         //  **************     aCirc1
635         //               o
636         //               o
637         //               ***************  aCirc2
638
639         Standard_Real aPar = 0.0;
640         aRng.GetIntermediatePoint(0.5, aPar);
641         const gp_Pnt aPCirc2 = ElCLib::Value(aPar, aWorkCirc);
642         const Extrema_POnCurv aP2(aPar, aPCirc2);
643
644         Extrema_ExtPElC ExtPCir(aPCirc2,
645                                 Extrema_CurveTool::Circle(*((Adaptor3d_Curve*) myC[0])),
646                                 Precision::Confusion(), theUt11, theUt12);
647
648         Standard_Boolean isFound = !myIsPar;
649
650         if (!isFound)
651         {
652           //If the flag myIsPar was set earlier then it does not mean that
653           //we have found the minimal distance. Here we check it. If there is
654           //a pair of points, which are in less distance then myIsPar flag
655           //was unset and the algorithm will return these nearest points.
656
657           for (Standard_Integer anExtID = 1; anExtID <= ExtPCir.NbExt(); anExtID++)
658           {
659             if (ExtPCir.SquareDistance(anExtID) < aMinSquareDist)
660             {
661               isFound = Standard_True;
662               break;
663             }
664           }
665         }
666
667         if (isFound)
668         {
669           ClearSolutions();
670           myIsPar = Standard_False;
671           for (Standard_Integer anExtID = 1; anExtID <= ExtPCir.NbExt(); anExtID++)
672           {
673             mypoints.Append(ExtPCir.Point(anExtID));
674             mypoints.Append(aP2);
675             mySqDist.Append(ExtPCir.SquareDistance(anExtID));
676             aMinSquareDist = Min(aMinSquareDist, ExtPCir.SquareDistance(anExtID));
677           }
678         }
679       }
680
681       aProjRng1.Shift(M_PI);
682     }
683   }
684 }
685
686 //=======================================================================
687 //function : Results
688 //purpose  : 
689 //=======================================================================
690
691 void Extrema_ExtCC::PrepareResults(const Extrema_ExtElC&  AlgExt,
692                                    const Standard_Boolean theIsInverse,
693                                    const Standard_Real    Ut11,
694                                    const Standard_Real    Ut12,
695                                    const Standard_Real    Ut21,
696                                    const Standard_Real    Ut22)
697 {
698   Standard_Integer i, NbExt;
699   Standard_Real Val, U, U2;
700   Extrema_POnCurv P1, P2;
701
702   myDone = AlgExt.IsDone();
703   if (myDone) {
704     myIsPar = AlgExt.IsParallel();
705     if (myIsPar) {
706       PrepareParallelResult(Ut11, Ut12, Ut21, Ut22, AlgExt.SquareDistance());
707     }
708     else {
709       NbExt = AlgExt.NbExt();
710       for (i = 1; i <= NbExt; i++) {
711         // Verification de la validite des parametres
712         AlgExt.Points(i, P1, P2);
713         if (!theIsInverse)
714         {
715           U = P1.Parameter();
716           U2 = P2.Parameter();
717         }
718         else {
719           U2 = P1.Parameter();
720           U = P2.Parameter();
721         }
722
723         if (Extrema_CurveTool::IsPeriodic(*((Adaptor3d_Curve*)myC[0]))) {
724           U = ElCLib::InPeriod(U, Ut11, Ut11+Extrema_CurveTool::Period(*((Adaptor3d_Curve*)myC[0])));
725         }
726         if (Extrema_CurveTool::IsPeriodic(*((Adaptor3d_Curve*)myC[1]))) {
727           U2 = ElCLib::InPeriod(U2, Ut21, Ut21+Extrema_CurveTool::Period(*((Adaptor3d_Curve*)myC[1])));
728         }
729
730         if ((U  >= Ut11 - RealEpsilon())  && 
731             (U  <= Ut12 + RealEpsilon())  &&
732             (U2 >= Ut21 - RealEpsilon())  &&
733             (U2 <= Ut22 + RealEpsilon())) {
734           Val = AlgExt.SquareDistance(i);
735           mySqDist.Append(Val);
736           if (!theIsInverse)
737           {
738             P1.SetValues(U, P1.Value());
739             P2.SetValues(U2, P2.Value());
740             mypoints.Append(P1);
741             mypoints.Append(P2);
742           }
743           else {
744             P1.SetValues(U2, P1.Value());
745             P2.SetValues(U, P2.Value());
746             mypoints.Append(P2);
747             mypoints.Append(P1);
748           }
749         }
750       }
751     }
752   }
753
754 }
755
756
757 //=======================================================================
758 //function : Results
759 //purpose  : 
760 //=======================================================================
761
762 void Extrema_ExtCC::PrepareResults(const Extrema_ECC&   AlgExt,
763                                    const Standard_Real  Ut11,
764                                    const Standard_Real  Ut12,
765                                    const Standard_Real  Ut21,
766                                    const Standard_Real  Ut22)
767 {
768   Standard_Integer i, NbExt;
769   Standard_Real Val, U, U2;
770   Extrema_POnCurv P1, P2;
771
772   myDone = AlgExt.IsDone();
773   if (myDone)
774   {
775     myIsPar = AlgExt.IsParallel();
776     if (myIsPar)
777     {
778       PrepareParallelResult(Ut11, Ut12, Ut21, Ut22, AlgExt.SquareDistance());
779     }
780     else
781     {
782       NbExt = AlgExt.NbExt();
783       for (i = 1; i <= NbExt; i++)
784       {
785         AlgExt.Points(i, P1, P2);
786         U = P1.Parameter();
787         U2 = P2.Parameter();
788
789         // Check points to be into param space.
790         if (Extrema_CurveTool::IsPeriodic(*((Adaptor3d_Curve*) myC[0])))
791         {
792           U = ElCLib::InPeriod(U, Ut11, Ut11 + Extrema_CurveTool::Period(*((Adaptor3d_Curve*) myC[0])));
793         }
794         if (Extrema_CurveTool::IsPeriodic(*((Adaptor3d_Curve*) myC[1])))
795         {
796           U2 = ElCLib::InPeriod(U2, Ut21, Ut21 + Extrema_CurveTool::Period(*((Adaptor3d_Curve*) myC[1])));
797         }
798
799         if ((U >= Ut11 - RealEpsilon()) &&
800             (U <= Ut12 + RealEpsilon()) &&
801             (U2 >= Ut21 - RealEpsilon()) &&
802             (U2 <= Ut22 + RealEpsilon()))
803         {
804           Val = AlgExt.SquareDistance(i);
805           mySqDist.Append(Val);
806           P1.SetValues(U, P1.Value());
807           P2.SetValues(U2, P2.Value());
808           mypoints.Append(P1);
809           mypoints.Append(P2);
810         }
811       }
812     }
813   }
814 }
815
816 //=======================================================================
817 //function : SetSingleSolutionFlag
818 //purpose  : 
819 //=======================================================================
820 void Extrema_ExtCC::SetSingleSolutionFlag(const Standard_Boolean theFlag)
821 {
822   myIsFindSingleSolution = theFlag;
823 }
824
825 //=======================================================================
826 //function : GetSingleSolutionFlag
827 //purpose  : 
828 //=======================================================================
829 Standard_Boolean Extrema_ExtCC::GetSingleSolutionFlag() const
830 {
831   return myIsFindSingleSolution;
832 }