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