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