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 #include <IntTools_EdgeEdge.ixx>
  16
  17 #include <gp_Dir.hxx>
  18 #include <gp_Lin.hxx>
  19
  20 #include <ElCLib.hxx>
  21
  22 #include <TopoDS_Iterator.hxx>
  23
  24 #include <Bnd_Box.hxx>
  25 #include <BndLib_Add3dCurve.hxx>
  26
  27 #include <GeomAPI_ProjectPointOnCurve.hxx>
  28
  29 #include <BRep_Tool.hxx>
  30 #include <BRepAdaptor_Curve.hxx>
  31
  32 #include <IntTools_Tools.hxx>
  33 #include <IntTools_CommonPrt.hxx>
  34
  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   void 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_Integer TypeToInteger(const GeomAbs_CurveType theCType);
  83
  84 //=======================================================================
  85 //function : Prepare
  86 //purpose  :
  87 //=======================================================================
  88 void IntTools_EdgeEdge::Prepare()
  89 {
  90   GeomAbs_CurveType aCT1, aCT2;
  91   Standard_Integer iCT1, iCT2;
  92   //
  93   myCurve1.Initialize(myEdge1);
  94   myCurve2.Initialize(myEdge2);
  95   //
  96   if (myRange1.First() == 0. && myRange1.Last() == 0.) {
  97     myRange1.SetFirst(myCurve1.FirstParameter());
  98     myRange1.SetLast (myCurve1.LastParameter());
  99   }
 100   //
 101   if (myRange2.First() == 0. && myRange2.Last() == 0.) {
 102     myRange2.SetFirst(myCurve2.FirstParameter());
 103     myRange2.SetLast (myCurve2.LastParameter());
 104   }
 105   //
 106   aCT1 = myCurve1.GetType();
 107   aCT2 = myCurve2.GetType();
 108   //
 109   iCT1 = TypeToInteger(aCT1);
 110   iCT2 = TypeToInteger(aCT2);
 111   //
 112   if (iCT1 == iCT2) {
 113     if (iCT1 != 0) {
 114       //compute deflection
 115       Standard_Integer i;
 116       Standard_Real aDt, aT, aT1, aT2;
 117       gp_Vec aV1, aV2;
 118       gp_Pnt aP;
 119       //
 120       Standard_Real aC1(10.), aC2(10.);
 121       for (i = 0; i < 2; ++i) {
 122         Standard_Real &aC = !i ? aC2 : aC1;
 123         aC = 0;
 124         IntTools_Range aR = !i ? myRange2 : myRange1;
 125         const BRepAdaptor_Curve& aBAC = !i ? myCurve2 : myCurve1;
 126         aR.Range(aT1, aT2);
 127         aDt = (aT2 - aT1) / 10.;
 128         aT = aT1;
 129         aBAC.D1(aT, aP, aV1);
 130         while (aT < aT2) {
 131           aT += aDt;
 132           aBAC.D1(aT, aP, aV2);
 133           if (aV1.Magnitude() > gp::Resolution() &&
 134               aV2.Magnitude() > gp::Resolution()) {
 135             gp_Dir aD1(aV1), aD2(aV2);
 136             aC += aD1.Angle(aD2);
 137           }
 138           aV1 = aV2;
 139         }
 140         //
 141         if (aC < Precision::Confusion()) {
 142           break;
 143         }
 144       }
 145       //
 146       if (aC1 < aC2) {
 147         --iCT1;
 148       }
 149     }
 150   }
 151   //
 152   if (iCT1 < iCT2) {
 153     TopoDS_Edge tmpE = myEdge1;
 154     myEdge1 = myEdge2;
 155     myEdge2 = tmpE;
 156     //
 157     BRepAdaptor_Curve tmpC = myCurve1;
 158     myCurve1 = myCurve2;
 159     myCurve2 = tmpC;
 160     //
 161     IntTools_Range tmpR = myRange1;
 162     myRange1 = myRange2;
 163     myRange2 = tmpR;
 164     //
 165     mySwap = Standard_True;
 166   }
 167   //
 168   myTol1 = myCurve1.Tolerance();
 169   myTol2 = myCurve2.Tolerance();
 170   myTol = myTol1 + myTol2;
 171   //
 172   myRes1 = myCurve1.Resolution(myTol);
 173   myRes2 = myCurve2.Resolution(myTol);
 174   //
 175   if (iCT1 != 0 || iCT2 != 0) {
 176     Standard_Real f, l, aTM;
 177     //
 178     myGeom1 = BRep_Tool::Curve(myEdge1, f, l);
 179     myGeom2 = BRep_Tool::Curve(myEdge2, f, l);
 180     //
 181     myPTol1 = 5.e-13;
 182     aTM = Max(fabs(myRange1.First()), fabs(myRange1.Last()));
 183     if (aTM > 999.) {
 184       myPTol1 = 5.e-16 * aTM;
 185     }
 186     //
 187     myPTol2 = 5.e-13;
 188     aTM = Max(fabs(myRange2.First()), fabs(myRange2.Last()));
 189     if (aTM > 999.) {
 190       myPTol2 = 5.e-16 * aTM;
 191     }
 192   }
 193 }
 194
 195 //=======================================================================
 196 //function : Perform
 197 //purpose  :
 198 //=======================================================================
 199 void IntTools_EdgeEdge::Perform()
 200 {
 201   //1. Check data
 202   CheckData();
 203   if (myErrorStatus) {
 204     return;
 205   }
 206   //
 207   //2. Prepare Data
 208   Prepare();
 209   //
 210   //3.1. Check Line/Line case
 211   if (myCurve1.GetType() == GeomAbs_Line &&
 212       myCurve2.GetType() == GeomAbs_Line) {
 213     ComputeLineLine();
 214     return;
 215   }
 216   //
 217   IntTools_SequenceOfRanges aRanges1, aRanges2;
 218   //
 219   //3.2. Find ranges containig solutions
 220   FindSolutions(myRange1, myRange2, aRanges1, aRanges2);
 221   //
 222   //4. Merge solutions and save common parts
 223   MergeSolutions(aRanges1, aRanges2);
 224 }
 225
 226 //=======================================================================
 227 //function : FindSolutions
 228 //purpose  :
 229 //=======================================================================
 230 void IntTools_EdgeEdge::FindSolutions(const IntTools_Range& theR1,
 231                                       const IntTools_Range& theR2,
 232                                       IntTools_SequenceOfRanges& theRanges1,
 233                                       IntTools_SequenceOfRanges& theRanges2)
 234 {
 235   Standard_Boolean bOut, bStop, bThin;
 236   Standard_Real aT11, aT12, aT21, aT22;
 237   Standard_Real aTB11, aTB12, aTB21, aTB22;
 238   Standard_Real aTol, aSmallStep1, aSmallStep2;
 239   Standard_Integer iCom;
 240   Bnd_Box aB1, aB2;
 241   //
 242   theR1.Range(aT11, aT12);
 243   theR2.Range(aT21, aT22);
 244   //
 245   BndBuildBox(myCurve1, aT11, aT12, myTol1, aB1);
 246   BndBuildBox(myCurve2, aT21, aT22, myTol2, aB2);
 247   if (aB1.IsOut(aB2)) {
 248     //no intersection;
 249     return;
 250   }
 251   //
 252   bOut  = Standard_False;
 253   bThin = Standard_False;
 254   bStop = Standard_False;
 255   aTol  = 2*myTol;
 256   iCom  = 1;
 257   //
 258   do {
 259     bOut = !FindParameters(myCurve2, aT21, aT22, myRes2, myPTol2, aB1, aTB21, aTB22);
 260     if (bOut) {
 261       break;
 262     }
 263     //
 264     bThin = (aTB22 - aTB21) < myRes2;
 265     if (bThin) {
 266       bOut = !FindParameters(myCurve1, aT11, aT12, myRes1, myPTol1, aB1, aTB11, aTB12);
 267       break;
 268     }
 269     //
 270     BndBuildBox(myCurve2, aTB21, aTB22, myTol2, aB2);
 271     //
 272     bOut = !FindParameters(myCurve1, aT11, aT12, myRes1, myPTol1, aB2, aTB11, aTB12);
 273     if (bOut) {
 274       break;
 275     }
 276     //
 277     bThin = ((aTB12 - aTB11) < myRes1) ||
 278       (aB2.IsXThin(aTol) && aB2.IsYThin(aTol) && aB2.IsZThin(aTol));
 279     //
 280     if (!bThin) {
 281       aSmallStep1 = (aT12 - aT11) / 250.;
 282       aSmallStep2 = (aT22 - aT21) / 250.;
 283       //
 284       if (aSmallStep1 < myRes1) {
 285         aSmallStep1 = myRes1;
 286       }
 287       if (aSmallStep2 < myRes2) {
 288         aSmallStep2 = myRes2;
 289       }
 290       //
 291       if (((aTB11 - aT11) < aSmallStep1) && ((aT12 - aTB12) < aSmallStep1) &&
 292           ((aTB21 - aT21) < aSmallStep2) && ((aT22 - aTB22) < aSmallStep2)) {
 293         bStop = Standard_True;
 294       } else {
 295         BndBuildBox(myCurve1, aTB11, aTB12, myTol1, aB1);
 296         bOut = aB1.IsOut(aB2);
 297         if (bOut) {
 298           break;
 299         }
 300       }
 301     }
 302     //
 303     aT11 = aTB11;
 304     aT12 = aTB12;
 305     aT21 = aTB21;
 306     aT22 = aTB22;
 307   } while (!bThin && !bStop);
 308   //
 309   if (bOut) {
 310     //no intersection;
 311     return;
 312   }
 313   //
 314   if (!bThin) {
 315     //check curves for coincidence on the ranges
 316     iCom = CheckCoincidence(aT11, aT12, aT21, aT22, myTol, myRes1);
 317     if (!iCom) {
 318       bThin = Standard_True;
 319     }
 320   }
 321   //
 322   if (bThin) {
 323     if (iCom != 0) {
 324       //check intermediate points
 325       Standard_Real aT1, aT2, aDist;
 326       gp_Pnt aP1, aP2;
 327       //
 328       aT1 = IntTools_Tools::IntermediatePoint(aT11, aT12);
 329       aT2 = IntTools_Tools::IntermediatePoint(aT21, aT22);
 330       //
 331       myGeom1->D0(aT1, aP1);
 332       myGeom2->D0(aT2, aP2);
 333       //
 334       aDist = aP1.Distance(aP2);
 335       if (aDist > myTol) {
 336         return;
 337       }
 338     }
 339     //add common part
 340     IntTools_Range aR1(aT11, aT12), aR2(aT21, aT22);
 341     //
 342     theRanges1.Append(aR1);
 343     theRanges2.Append(aR2);
 344     return;
 345   }
 346   //
 347   if (!IsIntersection(aT11, aT12, aT21, aT22)) {
 348     return;
 349   }
 350   //
 351   //split ranges on segments and repeat
 352   Standard_Integer i, aNb1;
 353   IntTools_SequenceOfRanges aSegments1;
 354   //
 355   IntTools_Range aR2(aT21, aT22);
 356   //
 357   SplitRangeOnSegments(aT11, aT12, myRes1, 3, aSegments1);
 358   aNb1 = aSegments1.Length();
 359   for (i = 1; i <= aNb1; ++i) {
 360     const IntTools_Range& aR1 = aSegments1(i);
 361     FindSolutions(aR1, aR2, theRanges1, theRanges2);
 362   }
 363 }
 364
 365 //=======================================================================
 366 //function : FindParameters
 367 //purpose  :
 368 //=======================================================================
 369 Standard_Boolean IntTools_EdgeEdge::FindParameters(const BRepAdaptor_Curve& theBAC,
 370                                                    const Standard_Real aT1,
 371                                                    const Standard_Real aT2,
 372                                                    const Standard_Real theRes,
 373                                                    const Standard_Real thePTol,
 374                                                    const Bnd_Box& theCBox,
 375                                                    Standard_Real& aTB1,
 376                                                    Standard_Real& aTB2)
 377 {
 378   Standard_Boolean bRet;
 379   Standard_Integer aC, i, k;
 380   Standard_Real aCf, aDiff, aDt, aT, aTB, aTOut, aTIn;
 381   Standard_Real aDist, aDistP, aDistTol, aTol;
 382   gp_Pnt aP;
 383   Bnd_Box aCBx;
 384   //
 385   bRet = Standard_False;
 386   aCf = 0.6180339887498948482045868343656;// =0.5*(1.+sqrt(5.))/2.;
 387   aDt = theRes;
 388   aTol = theBAC.Tolerance();
 389   aDistP = 0.;
 390   aDistTol = Precision::PConfusion();
 391   aCBx = theCBox;
 392   aCBx.Enlarge(aTol);
 393   //
 394   for (i = 0; i < 2; ++i) {
 395     aTB = !i ? aT1 : aT2;
 396     aT = !i ? aT2 : aTB1;
 397     aC = !i ? 1 : -1;
 398     bRet = Standard_False;
 399     k = 0;
 400     //looking for the point on the edge which is in the box;
 401     while (aC*(aT-aTB) >= 0) {
 402       theBAC.D0(aTB, aP);
 403       aDist = PointBoxDistance(theCBox, aP);
 404       if (aDist > aTol) {
 405         if (fabs(aDist - aDistP) < aDistTol) {
 406           aDt = theBAC.Resolution((++k)*aDist);
 407         } else {
 408           k = 0;
 409           aDt = theBAC.Resolution(aDist);
 410         }
 411         aTB += aC*aDt;
 412       } else {
 413         bRet = Standard_True;
 414         break;
 415       }
 416       aDistP = aDist;
 417     }
 418     //
 419     if (!bRet) {
 420       if (!i) {
 421         //edge is out of the box;
 422         return bRet;
 423       } else {
 424         bRet = !bRet;
 425         aTB = aTB1;
 426         aDt = aT2 - aTB1;
 427       }
 428     }
 429     //
 430     aT = !i ? aT1 : aT2;
 431     if (aTB != aT) {
 432       //one point IN, one point OUT; looking for the bounding point;
 433       aTIn = aTB;
 434       aTOut = aTB - aC*aDt;
 435       aDiff = aTIn - aTOut;
 436       while (fabs(aDiff) > thePTol) {
 437         aTB = aTOut + aDiff*aCf;
 438         theBAC.D0(aTB, aP);
 439         if (aCBx.IsOut(aP)) {
 440           aTOut = aTB;
 441         } else {
 442           aTIn = aTB;
 443         }
 444         aDiff = aTIn - aTOut;
 445       }
 446     }
 447     if (!i) {
 448       aTB1 = aTB;
 449     } else {
 450       aTB2 = aTB;
 451     }
 452   }
 453   return bRet;
 454 }
 455
 456 //=======================================================================
 457 //function : MergeSolutions
 458 //purpose  :
 459 //=======================================================================
 460 void IntTools_EdgeEdge::MergeSolutions(const IntTools_SequenceOfRanges& theRanges1,
 461                                        const IntTools_SequenceOfRanges& theRanges2)
 462 {
 463   IntTools_Range aRi1, aRi2, aRj1, aRj2;
 464   Standard_Integer aNbCP, i, j;
 465   TopAbs_ShapeEnum aType;
 466   Standard_Real aTi11, aTi12, aTi21, aTi22,
 467                 aTj11, aTj12, aTj21, aTj22;
 468   //
 469   aNbCP = theRanges1.Length();
 470   aType = TopAbs_VERTEX;
 471   //
 472   for (i = 1; i <= aNbCP; ) {
 473     aRi1 = theRanges1(i);
 474     aRi2 = theRanges2(i);
 475     //
 476     aRi1.Range(aTi11, aTi12);
 477     aRi2.Range(aTi21, aTi22);
 478     //
 479     for (j = i+1; j <= aNbCP; ++j) {
 480       aRj1 = theRanges1(j);
 481       aRj2 = theRanges2(j);
 482       //
 483       aRj1.Range(aTj11, aTj12);
 484       aRj2.Range(aTj21, aTj22);
 485       if (fabs(aTi12 - aTj11) < 10*myRes1 ||
 486           fabs(aTi22 - aTj21) < 10*myRes2) {
 487         aTi11 = Min(aTi11, aTj11);
 488         aTi12 = Max(aTi12, aTj12);
 489         aTi21 = Min(aTi21, aTj21);
 490         aTi22 = Max(aTi22, aTj22);
 491       } else {
 492         break;
 493       }
 494     }
 495     i = j;
 496     //
 497     if (aTi11 == myRange1.First() && aTi12 == myRange1.Last() &&
 498         aTi21 == myRange2.First() && aTi22 == myRange2.Last()) {
 499       aType = TopAbs_EDGE;
 500     }
 501     //
 502     AddSolution(aTi11, aTi12, aTi21, aTi22, aType);
 503   }
 504 }
 505
 506 //=======================================================================
 507 //function : AddSolution
 508 //purpose  :
 509 //=======================================================================
 510 void IntTools_EdgeEdge::AddSolution(const Standard_Real aT11,
 511                                     const Standard_Real aT12,
 512                                     const Standard_Real aT21,
 513                                     const Standard_Real aT22,
 514                                     const TopAbs_ShapeEnum theType)
 515 {
 516   IntTools_CommonPrt aCPart;
 517   //
 518   aCPart.SetType(theType);
 519   if (!mySwap) {
 520     aCPart.SetEdge1(myEdge1);
 521     aCPart.SetEdge2(myEdge2);
 522     aCPart.SetRange1(aT11, aT12);
 523     aCPart.AppendRange2(aT21, aT22);
 524   } else {
 525     aCPart.SetEdge1(myEdge2);
 526     aCPart.SetEdge2(myEdge1);
 527     aCPart.SetRange1(aT21, aT22);
 528     aCPart.AppendRange2(aT11, aT12);
 529   }
 530   //
 531   if (theType == TopAbs_VERTEX) {
 532     Standard_Real aT1, aT2;
 533     //
 534     FindBestSolution(aT11, aT12, aT21, aT22, aT1, aT2);
 535     //
 536     if (!mySwap) {
 537       aCPart.SetVertexParameter1(aT1);
 538       aCPart.SetVertexParameter2(aT2);
 539     } else {
 540       aCPart.SetVertexParameter1(aT2);
 541       aCPart.SetVertexParameter2(aT1);
 542     }
 543   }
 544   myCommonParts.Append(aCPart);
 545 }
 546
 547 //=======================================================================
 548 //function : FindBestSolution
 549 //purpose  :
 550 //=======================================================================
 551 void IntTools_EdgeEdge::FindBestSolution(const Standard_Real aT11,
 552                                          const Standard_Real aT12,
 553                                          const Standard_Real aT21,
 554                                          const Standard_Real aT22,
 555                                          Standard_Real& aT1,
 556                                          Standard_Real& aT2)
 557 {
 558   Standard_Integer i, aNbS, iErr;
 559   Standard_Real aDMin, aD, aCrit;
 560   Standard_Real aT1x, aT2x, aT1p, aT2p;
 561   GeomAPI_ProjectPointOnCurve aProj;
 562   IntTools_SequenceOfRanges aSeg1;
 563   //
 564   aT1 = IntTools_Tools::IntermediatePoint(aT11, aT12);
 565   aT2 = IntTools_Tools::IntermediatePoint(aT21, aT22);
 566   //
 567   aDMin = 100.;
 568   aD = 100.;
 569   aCrit = 0.;//1.e-16;
 570   //
 571   aNbS = 10;
 572   SplitRangeOnSegments(aT11, aT12, 3*myRes1, aNbS, aSeg1);
 573   aNbS = aSeg1.Length();
 574   //
 575   aProj.Init(myGeom2, aT21, aT22);
 576   for (i = 1; i <= aNbS; ++i) {
 577     const IntTools_Range& aR1 = aSeg1(i);
 578     aR1.Range(aT1x, aT2x);
 579     //
 580     iErr = FindDistPC(aT1x, aT2x, myGeom1, aCrit, myPTol1,
 581                       aProj, aD, aT1p, aT2p, Standard_False);
 582     if (iErr != 1 && aD < aDMin) {
 583       aT1 = aT1p;
 584       aT2 = aT2p;
 585       aDMin = aD;
 586       if (aDMin <= aCrit) {
 587         break;
 588       }
 589     }
 590   }
 591 }
 592
 593 //=======================================================================
 594 //function : ComputeLineLine
 595 //purpose  :
 596 //=======================================================================
 597 void IntTools_EdgeEdge::ComputeLineLine()
 598 {
 599   Standard_Boolean IsParallel, IsCoincide;
 600   Standard_Real aSin, aCos, aAng, aTol;
 601   Standard_Real aT1, aT2, aT11, aT12, aT21, aT22;
 602   gp_Pnt aP11, aP12;
 603   gp_Lin aL1, aL2;
 604   gp_Dir aD1, aD2;
 605   IntTools_CommonPrt aCommonPrt;
 606   //
 607   IsParallel = Standard_False;
 608   IsCoincide = Standard_False;
 609   aTol = myTol*myTol;
 610   aL1 = myCurve1.Line();
 611   aL2 = myCurve2.Line();
 612   aD1 = aL1.Position().Direction();
 613   aD2 = aL2.Position().Direction();
 614   myRange1.Range(aT11, aT12);
 615   myRange2.Range(aT21, aT22);
 616   //
 617   aCommonPrt.SetEdge1(myEdge1);
 618   aCommonPrt.SetEdge2(myEdge2);
 619   //
 620   aCos = aD1.Dot(aD2);
 621   aAng = (aCos >= 0.) ? 2.*(1. - aCos) : 2.*(1. + aCos);
 622   //
 623   if(aAng <= Precision::Angular()) {
 624     IsParallel = Standard_True;
 625     if(aL1.SquareDistance(aL2.Location()) <= aTol) {
 626       IsCoincide = Standard_True;
 627       aP11 = ElCLib::Value(aT11, aL1);
 628       aP12 = ElCLib::Value(aT12, aL1);
 629     }
 630   }
 631   else {
 632     aP11 = ElCLib::Value(aT11, aL1);
 633     aP12 = ElCLib::Value(aT12, aL1);
 634     if(aL2.SquareDistance(aP11) <= aTol && aL2.SquareDistance(aP12) <= aTol) {
 635       IsCoincide = Standard_True;
 636     }
 637   }
 638   //
 639   if (IsCoincide) {
 640     Standard_Real t21, t22;
 641     //
 642     t21 = ElCLib::Parameter(aL2, aP11);
 643     t22 = ElCLib::Parameter(aL2, aP12);
 644     if((t21 > aT22 && t22 > aT22) || (t21 < aT21 && t22 < aT21)) {
 645       return;
 646     }
 647     //
 648     Standard_Real temp;
 649     if(t21 > t22) {
 650       temp = t21;
 651       t21 = t22;
 652       t22 = temp;
 653     }
 654     //
 655     if(t21 >= aT21) {
 656       if(t22 <= aT22) {
 657         aCommonPrt.SetRange1(aT11, aT12);
 658         aCommonPrt.SetAllNullFlag(Standard_True);
 659         aCommonPrt.AppendRange2(t21, t22);
 660       }
 661       else {
 662         aCommonPrt.SetRange1(aT11, aT12 - (t22 - aT22));
 663         aCommonPrt.AppendRange2(t21, aT22);
 664       }
 665     }
 666     else {
 667       aCommonPrt.SetRange1(aT11 + (aT21 - t21), aT12);
 668       aCommonPrt.AppendRange2(aT21, t22);
 669     }
 670     aCommonPrt.SetType(TopAbs_EDGE);
 671     myCommonParts.Append(aCommonPrt);
 672     return;
 673   }
 674   //
 675   if (IsParallel) {
 676     return;
 677   }
 678   //
 679   {
 680     TopoDS_Iterator aIt1, aIt2;
 681     aIt1.Initialize(myEdge1);
 682     for (; aIt1.More(); aIt1.Next()) {
 683       const TopoDS_Shape& aV1 = aIt1.Value();
 684       aIt2.Initialize(myEdge2);
 685       for (; aIt2.More(); aIt2.Next()) {
 686         const TopoDS_Shape& aV2 = aIt2.Value();
 687         if (aV2.IsSame(aV1)) {
 688           return;
 689         }
 690       }
 691     }
 692   }
 693   //
 694   aSin = 1. - aCos*aCos;
 695   gp_Pnt O1 = aL1.Location();
 696   gp_Pnt O2 = aL2.Location();
 697   gp_Vec O1O2 (O1, O2);
 698   //
 699   aT2 = (aD1.XYZ()*(O1O2.Dot(aD1))-(O1O2.XYZ())).Dot(aD2.XYZ());
 700   aT2 /= aSin;
 701   //
 702   if(aT2 < aT21 || aT2 > aT22) {
 703     return;
 704   }
 705   //
 706   gp_Pnt aP2(ElCLib::Value(aT2, aL2));
 707   aT1 = (gp_Vec(O1, aP2)).Dot(aD1);
 708   //
 709   if(aT1 < aT11 || aT1 > aT12) {
 710     return;
 711   }
 712   //
 713   gp_Pnt aP1(ElCLib::Value(aT1, aL1));
 714   Standard_Real aDist = aP1.SquareDistance(aP2);
 715   //
 716   if (aDist > aTol) {
 717     return;
 718   }
 719   //
 720   aCommonPrt.SetRange1(aT1 - myTol, aT1 + myTol);
 721   aCommonPrt.AppendRange2(aT2 - myTol, aT2 + myTol);
 722   aCommonPrt.SetType(TopAbs_VERTEX);
 723   aCommonPrt.SetVertexParameter1(aT1);
 724   aCommonPrt.SetVertexParameter2(aT2);
 725   myCommonParts.Append(aCommonPrt);
 726 }
 727
 728 //=======================================================================
 729 //function : IsIntersection
 730 //purpose  :
 731 //=======================================================================
 732 Standard_Boolean IntTools_EdgeEdge::IsIntersection(const Standard_Real aT11,
 733                                                    const Standard_Real aT12,
 734                                                    const Standard_Real aT21,
 735                                                    const Standard_Real aT22)
 736 {
 737   Standard_Boolean bRet;
 738   gp_Pnt aP11, aP12, aP21, aP22;
 739   gp_Vec aV11, aV12, aV21, aV22;
 740   Standard_Real aD11_21, aD11_22, aD12_21, aD12_22, aCriteria, aCoef;
 741   Standard_Boolean bSmall_11_21, bSmall_11_22, bSmall_12_21, bSmall_12_22;
 742   //
 743   bRet = Standard_True;
 744   aCoef = 1.e+5;
 745   if (((aT12 - aT11) > aCoef*myRes1) && ((aT22 - aT21) > aCoef*myRes2)) {
 746     aCoef = 5000;
 747   } else {
 748     Standard_Real aTRMin = Min((aT12 - aT11)/myRes1, (aT22 - aT21)/myRes2);
 749     aCoef = aTRMin / 100.;
 750     if (aCoef < 1.) {
 751       aCoef = 1.;
 752     }
 753   }
 754   aCriteria = aCoef * myTol;
 755   aCriteria *= aCriteria;
 756   //
 757   myGeom1->D1(aT11, aP11, aV11);
 758   myGeom1->D1(aT12, aP12, aV12);
 759   myGeom2->D1(aT21, aP21, aV21);
 760   myGeom2->D1(aT22, aP22, aV22);
 761   //
 762   aD11_21 = aP11.SquareDistance(aP21);
 763   aD11_22 = aP11.SquareDistance(aP22);
 764   aD12_21 = aP12.SquareDistance(aP21);
 765   aD12_22 = aP12.SquareDistance(aP22);
 766   //
 767   bSmall_11_21 = aD11_21 < aCriteria;
 768   bSmall_11_22 = aD11_22 < aCriteria;
 769   bSmall_12_21 = aD12_21 < aCriteria;
 770   bSmall_12_22 = aD12_22 < aCriteria;
 771   //
 772   if ((bSmall_11_21 && bSmall_12_22) ||
 773       (bSmall_11_22 && bSmall_12_21)) {
 774     if (aCoef == 1.) {
 775       return bRet;
 776     }
 777     //
 778     Standard_Real anAngleCriteria;
 779     Standard_Real anAngle1, anAngle2;
 780     //
 781     anAngleCriteria = 5.e-3;
 782     if (bSmall_11_21 && bSmall_12_22) {
 783       anAngle1 = aV11.Angle(aV21);
 784       anAngle2 = aV12.Angle(aV22);
 785     } else {
 786       anAngle1 = aV11.Angle(aV22);
 787       anAngle2 = aV12.Angle(aV21);
 788     }
 789     //
 790     if (((anAngle1 < anAngleCriteria) || ((M_PI - anAngle1) < anAngleCriteria)) ||
 791         ((anAngle2 < anAngleCriteria) || ((M_PI - anAngle2) < anAngleCriteria))) {
 792       GeomAPI_ProjectPointOnCurve aProj;
 793       Standard_Integer iErr;
 794       Standard_Real aD, aT1p, aT2p;
 795       //
 796       aD = 100.;
 797       aProj.Init(myGeom2, aT21, aT22);
 798       iErr = FindDistPC(aT11, aT12, myGeom1, myTol, myRes1, aProj, aD, aT1p, aT2p, Standard_False);
 799       bRet = (iErr == 2);
 800     }
 801   }
 802   return bRet;
 803 }
 804
 805 //=======================================================================
 806 //function : CheckCoincidence
 807 //purpose  :
 808 //=======================================================================
 809 Standard_Integer IntTools_EdgeEdge::CheckCoincidence(const Standard_Real aT11,
 810                                                      const Standard_Real aT12,
 811                                                      const Standard_Real aT21,
 812                                                      const Standard_Real aT22,
 813                                                      const Standard_Real theCriteria,
 814                                                      const Standard_Real theCurveRes1)
 815 {
 816   Standard_Integer iErr, aNb, aNb1, i;
 817   Standard_Real aT1A, aT1B, aT1max, aT2max, aDmax;
 818   GeomAPI_ProjectPointOnCurve aProjPC;
 819   IntTools_SequenceOfRanges aSeg1;
 820   //
 821   iErr  = 0;
 822   aDmax = -1.;
 823   aProjPC.Init(myGeom2, aT21, aT22);
 824   //
 825   // 1. Express evaluation
 826   aNb = 10; // Number of intervals on the curve #1
 827   SplitRangeOnSegments(aT11, aT12, theCurveRes1, aNb, aSeg1);
 828   aNb1 = aSeg1.Length();
 829   for (i = 1; i < aNb1; ++i) {
 830     const IntTools_Range& aR1 = aSeg1(i);
 831     aR1.Range(aT1A, aT1B);
 832     //
 833     iErr = DistPC(aT1B, myGeom1, theCriteria, aProjPC, aDmax, aT2max);
 834     if (iErr) {
 835       return iErr;
 836     }
 837   }
 838   //
 839   // if the ranges in aSeg1 are less than theCurveRes1,
 840   // there is no need to do step 2 (deep evaluation)
 841   if (aNb1 < aNb) {
 842     return iErr;
 843   }
 844   //
 845   // 2. Deep evaluation
 846   for (i = 2; i < aNb1; ++i) {
 847     const IntTools_Range& aR1 = aSeg1(i);
 848     aR1.Range(aT1A, aT1B);
 849     //
 850     iErr = FindDistPC(aT1A, aT1B, myGeom1, theCriteria, theCurveRes1,
 851                       aProjPC, aDmax, aT1max, aT2max);
 852     if (iErr) {
 853       return iErr;
 854     }
 855   }
 856   // Possible values:
 857   // iErr == 0 - the patches are coincided
 858   // iErr == 1 - a point from aC1 can not be projected on aC2
 859   // iErr == 2 - the distance is too big
 860   return iErr;
 861 }
 862
 863 //=======================================================================
 864 //function : FindDistPC
 865 //purpose  :
 866 //=======================================================================
 867 Standard_Integer FindDistPC(const Standard_Real aT1A,
 868                             const Standard_Real aT1B,
 869                             const Handle(Geom_Curve)& theC1,
 870                             const Standard_Real theCriteria,
 871                             const Standard_Real theEps,
 872                             GeomAPI_ProjectPointOnCurve& theProjPC,
 873                             Standard_Real& aDmax,
 874                             Standard_Real& aT1max,
 875                             Standard_Real& aT2max,
 876                             const Standard_Boolean bMaxDist)
 877 {
 878   Standard_Integer iErr, iC;
 879   Standard_Real aGS, aXP, aA, aB, aXL, aYP, aYL, aT2P, aT2L;
 880   //
 881   iC = bMaxDist ? 1 : -1;
 882   iErr = 0;
 883   //
 884   aGS = 0.6180339887498948482045868343656;// =0.5*(1.+sqrt(5.))-1.;
 885   aA = aT1A;
 886   aB = aT1B;
 887   //
 888   // check bounds
 889   iErr = DistPC(aA, theC1, theCriteria, theProjPC, aYP, aT2P, aDmax, aT1max, aT2max, iC);
 890   if (iErr == 2) {
 891     return iErr;
 892   }
 893   //
 894   iErr = DistPC(aB, theC1, theCriteria, theProjPC, aYL, aT2L, aDmax, aT1max, aT2max, iC);
 895   if (iErr == 2) {
 896     return iErr;
 897   }
 898   //
 899   aXP = aA + (aB - aA)*aGS;
 900   aXL = aB - (aB - aA)*aGS;
 901   //
 902   iErr = DistPC(aXP, theC1, theCriteria, theProjPC, aYP, aT2P, aDmax, aT1max, aT2max, iC);
 903   if (iErr) {
 904     return iErr;
 905   }
 906   //
 907   iErr = DistPC(aXL, theC1, theCriteria, theProjPC, aYL, aT2L, aDmax, aT1max, aT2max, iC);
 908   if (iErr) {
 909     return iErr;
 910   }
 911   //
 912   for (;;) {
 913     if (iC*(aYP - aYL) > 0) {
 914       aA = aXL;
 915       aXL = aXP;
 916       aYL = aYP;
 917       aXP = aA + (aB - aA)*aGS;
 918       iErr = DistPC(aXP, theC1, theCriteria, theProjPC, aYP, aT2P, aDmax, aT1max, aT2max, iC);
 919       if (iErr) {
 920         return iErr;
 921       }
 922     }
 923     else {
 924       aB = aXP;
 925       aXP = aXL;
 926       aYP = aYL;
 927       aXL = aB - (aB - aA)*aGS;
 928       iErr = DistPC(aXL, theC1, theCriteria, theProjPC, aYL, aT2L, aDmax, aT1max, aT2max, iC);
 929       if (iErr) {
 930         return iErr;
 931       }
 932     }
 933     //
 934     if ((aB - aA) < theEps) {
 935       break;
 936     }
 937   }// for (;;) {
 938   //
 939   return iErr;
 940 }
 941 //=======================================================================
 942 //function : DistPC
 943 //purpose  :
 944 //=======================================================================
 945 Standard_Integer DistPC(const Standard_Real aT1,
 946                         const Handle(Geom_Curve)& theC1,
 947                         const Standard_Real theCriteria,
 948                         GeomAPI_ProjectPointOnCurve& theProjPC,
 949                         Standard_Real& aD,
 950                         Standard_Real& aT2,
 951                         Standard_Real& aDmax,
 952                         Standard_Real& aT1max,
 953                         Standard_Real& aT2max,
 954                         const Standard_Integer iC)
 955 {
 956   Standard_Integer iErr;
 957   //
 958   iErr = DistPC(aT1, theC1, theCriteria, theProjPC, aD, aT2, iC);
 959   if (iErr) {
 960     return iErr;
 961   }
 962   //
 963   if (iC*(aD - aDmax) > 0) {
 964     aDmax = aD;
 965     aT1max = aT1;
 966     aT2max = aT2;
 967   }
 968   //
 969   return iErr;
 970 }
 971 //=======================================================================
 972 //function : DistPC
 973 //purpose  :
 974 //=======================================================================
 975 Standard_Integer DistPC(const Standard_Real aT1,
 976                         const Handle(Geom_Curve)& theC1,
 977                         const Standard_Real theCriteria,
 978                         GeomAPI_ProjectPointOnCurve& theProjPC,
 979                         Standard_Real& aD,
 980                         Standard_Real& aT2,
 981                         const Standard_Integer iC)
 982 {
 983   Standard_Integer iErr, aNbP2;
 984   gp_Pnt aP1;
 985   //
 986   iErr = 0;
 987   theC1->D0(aT1, aP1);
 988   //
 989   theProjPC.Perform(aP1);
 990   aNbP2 = theProjPC.NbPoints();
 991   if (!aNbP2) {
 992     iErr = 1;// the point from aC1 can not be projected on aC2
 993     return iErr;
 994   }
 995   //
 996   aD  = theProjPC.LowerDistance();
 997   aT2 = theProjPC.LowerDistanceParameter();
 998   if (iC*(aD - theCriteria) > 0) {
 999     iErr = 2;// the distance is too big or small
1000   }
1001   //
1002   return iErr;
1003 }
1004
1005 //=======================================================================
1006 //function : SplitRangeOnSegments
1007 //purpose  :
1008 //=======================================================================
1009 void SplitRangeOnSegments(const Standard_Real aT1,
1010                           const Standard_Real aT2,
1011                           const Standard_Real theResolution,
1012                           const Standard_Integer theNbSeg,
1013                           IntTools_SequenceOfRanges& theSegments)
1014 {
1015   Standard_Real aDiff = aT2 - aT1;
1016   if (aDiff < theResolution) {
1017     theSegments.Append(IntTools_Range(aT1, aT2));
1018     return;
1019   }
1020   //
1021   Standard_Real aDt, aT1x, aT2x, aSeg;
1022   Standard_Integer aNbSegments, i;
1023   //
1024   aNbSegments = theNbSeg;
1025   aDt = aDiff / aNbSegments;
1026   if (aDt < theResolution) {
1027     aSeg = aDiff / theResolution;
1028     aNbSegments = Standard_Integer(aSeg) + 1;
1029     aDt = aDiff / aNbSegments;
1030   }
1031   //
1032   aT1x = aT1;
1033   for (i = 1; i < aNbSegments; ++i) {
1034     aT2x = aT1x + aDt;
1035     //
1036     IntTools_Range aR(aT1x, aT2x);
1037     theSegments.Append(aR);
1038     //
1039     aT1x = aT2x;
1040   }
1041   //
1042   IntTools_Range aR(aT1x, aT2);
1043   theSegments.Append(aR);
1044 }
1045
1046 //=======================================================================
1047 //function : BndBuildBox
1048 //purpose  :
1049 //=======================================================================
1050 void BndBuildBox(const BRepAdaptor_Curve& theBAC,
1051                  const Standard_Real aT1,
1052                  const Standard_Real aT2,
1053                  const Standard_Real theTol,
1054                  Bnd_Box& theBox)
1055 {
1056   Bnd_Box aB;
1057   BndLib_Add3dCurve::Add(theBAC, aT1, aT2, theTol, aB);
1058   theBox = aB;
1059 }
1060
1061 //=======================================================================
1062 //function : PointBoxDistance
1063 //purpose  :
1064 //=======================================================================
1065 Standard_Real PointBoxDistance(const Bnd_Box& aB,
1066                                const gp_Pnt& aP)
1067 {
1068   Standard_Real aPCoord[3];
1069   Standard_Real aBMinCoord[3], aBMaxCoord[3];
1070   Standard_Real aDist, aR1, aR2;
1071   Standard_Integer i;
1072   //
1073   aP.Coord(aPCoord[0], aPCoord[1], aPCoord[2]);
1074   aB.Get(aBMinCoord[0], aBMinCoord[1], aBMinCoord[2],
1075          aBMaxCoord[0], aBMaxCoord[1], aBMaxCoord[2]);
1076   //
1077   aDist = 0.;
1078   for (i = 0; i < 3; ++i) {
1079     aR1 = aBMinCoord[i] - aPCoord[i];
1080     if (aR1 > 0.) {
1081       aDist += aR1*aR1;
1082       continue;
1083     }
1084     //
1085     aR2 = aPCoord[i] - aBMaxCoord[i];
1086     if (aR2 > 0.) {
1087       aDist += aR2*aR2;
1088     }
1089   }
1090   //
1091   aDist = Sqrt(aDist);
1092   return aDist;
1093 }
1094
1095 //=======================================================================
1096 //function : TypeToInteger
1097 //purpose  :
1098 //=======================================================================
1099 Standard_Integer TypeToInteger(const GeomAbs_CurveType theCType)
1100 {
1101   Standard_Integer iRet;
1102   //
1103   switch(theCType) {
1104   case GeomAbs_Line:
1105     iRet=0;
1106     break;
1107   case GeomAbs_Circle:
1108     iRet=1;
1109     break;
1110   case GeomAbs_Ellipse:
1111   case GeomAbs_Hyperbola:
1112   case GeomAbs_Parabola:
1113     iRet=2;
1114     break;
1115   case GeomAbs_BezierCurve:
1116   case GeomAbs_BSplineCurve:
1117     iRet=3;
1118     break;
1119   default:
1120     iRet=4;
1121     break;
1122   }
1123   return iRet;
1124 }
1125