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