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