src/IntTools/IntTools_EdgeEdge.cxx

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