src/IntTools/IntTools_FaceFace.cxx

   1 // Created on: 2000-11-23
   2 // Created by: Michael KLOKOV
   3 // Copyright (c) 2000-2014 OPEN CASCADE SAS
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   7 // This library is free software; you can redistribute it and/or modify it under
   8 // the terms of the GNU Lesser General Public License version 2.1 as published
   9 // by the Free Software Foundation, with special exception defined in the file
  10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  11 // distribution for complete text of the license and disclaimer of any warranty.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement.
  15
  16 #include <IntTools_FaceFace.hxx>
  17
  18 #include <Precision.hxx>
  19
  20 #include <TColStd_HArray1OfReal.hxx>
  21 #include <TColStd_Array1OfReal.hxx>
  22 #include <TColStd_Array1OfInteger.hxx>
  23 #include <TColStd_SequenceOfReal.hxx>
  24 #include <TColStd_ListOfInteger.hxx>
  25 #include <TColStd_ListIteratorOfListOfInteger.hxx>
  26 #include <TColStd_Array1OfListOfInteger.hxx>
  27
  28 #include <gp_Lin2d.hxx>
  29 #include <gp_Ax22d.hxx>
  30 #include <gp_Circ2d.hxx>
  31 #include <gp_Torus.hxx>
  32 #include <gp_Cylinder.hxx>
  33
  34 #include <Bnd_Box.hxx>
  35
  36 #include <TColgp_HArray1OfPnt2d.hxx>
  37 #include <TColgp_SequenceOfPnt2d.hxx>
  38 #include <TColgp_Array1OfPnt.hxx>
  39 #include <TColgp_Array1OfPnt2d.hxx>
  40
  41 #include <IntAna_QuadQuadGeo.hxx>
  42
  43 #include <IntSurf_PntOn2S.hxx>
  44 #include <IntSurf_LineOn2S.hxx>
  45 #include <IntSurf_PntOn2S.hxx>
  46 #include <IntSurf_ListOfPntOn2S.hxx>
  47 #include <IntRes2d_Domain.hxx>
  48 #include <ProjLib_Plane.hxx>
  49
  50 #include <IntPatch_GLine.hxx>
  51 #include <IntPatch_RLine.hxx>
  52 #include <IntPatch_WLine.hxx>
  53 #include <IntPatch_ALine.hxx>
  54 #include <IntPatch_ALineToWLine.hxx>
  55
  56 #include <ElSLib.hxx>
  57 #include <ElCLib.hxx>
  58
  59 #include <Extrema_ExtCC.hxx>
  60 #include <Extrema_POnCurv.hxx>
  61 #include <BndLib_AddSurface.hxx>
  62
  63 #include <Adaptor3d_SurfacePtr.hxx>
  64 #include <Adaptor2d_HLine2d.hxx>
  65
  66 #include <GeomAbs_SurfaceType.hxx>
  67 #include <GeomAbs_CurveType.hxx>
  68
  69 #include <Geom_Surface.hxx>
  70 #include <Geom_Line.hxx>
  71 #include <Geom_Circle.hxx>
  72 #include <Geom_Ellipse.hxx>
  73 #include <Geom_Parabola.hxx>
  74 #include <Geom_Hyperbola.hxx>
  75 #include <Geom_TrimmedCurve.hxx>
  76 #include <Geom_BSplineCurve.hxx>
  77 #include <Geom_RectangularTrimmedSurface.hxx>
  78 #include <Geom_OffsetSurface.hxx>
  79 #include <Geom_Curve.hxx>
  80 #include <Geom_Conic.hxx>
  81
  82 #include <Geom2d_TrimmedCurve.hxx>
  83 #include <Geom2d_BSplineCurve.hxx>
  84 #include <Geom2d_Line.hxx>
  85 #include <Geom2d_Curve.hxx>
  86 #include <Geom2d_Circle.hxx>
  87
  88 #include <Geom2dAPI_InterCurveCurve.hxx>
  89 #include <Geom2dInt_GInter.hxx>
  90 #include <Geom2dAdaptor.hxx>
  91 #include <GeomAdaptor_Curve.hxx>
  92 #include <GeomAdaptor_HSurface.hxx>
  93 #include <GeomAdaptor_Surface.hxx>
  94 #include <GeomLib_CheckBSplineCurve.hxx>
  95 #include <GeomLib_Check2dBSplineCurve.hxx>
  96
  97 #include <GeomInt_WLApprox.hxx>
  98 #include <GeomProjLib.hxx>
  99 #include <GeomAPI_ProjectPointOnSurf.hxx>
 100 #include <Geom2dAdaptor_Curve.hxx>
 101 #include <TopoDS.hxx>
 102 #include <TopoDS_Edge.hxx>
 103 #include <TopExp_Explorer.hxx>
 104
 105 #include <BRep_Tool.hxx>
 106 #include <BRepTools.hxx>
 107 #include <BRepAdaptor_Surface.hxx>
 108
 109 #include <IntTools_Curve.hxx>
 110 #include <IntTools_Tools.hxx>
 111 #include <IntTools_Tools.hxx>
 112 #include <IntTools_TopolTool.hxx>
 113 #include <IntTools_PntOnFace.hxx>
 114 #include <IntTools_PntOn2Faces.hxx>
 115 #include <IntTools_Context.hxx>
 116 #include <IntSurf_ListIteratorOfListOfPntOn2S.hxx>
 117 #include <GeomInt.hxx>
 118
 119 #include <Approx_CurveOnSurface.hxx>
 120 #include <GeomAdaptor.hxx>
 121 #include <GeomInt_IntSS.hxx>
 122
 123 static
 124   void RefineVector(gp_Vec2d& aV2D);
 125 #ifdef OCCT_DEBUG_DUMPWLINE
 126 static
 127   void DumpWLine(const Handle(IntPatch_WLine)& aWLine);
 128 #endif
 129 //
 130 static
 131   void TolR3d(const TopoDS_Face& ,
 132               const TopoDS_Face& ,
 133               Standard_Real& );
 134 static
 135   Handle(Geom_Curve) MakeBSpline (const Handle(IntPatch_WLine)&,
 136                                   const Standard_Integer,
 137                                   const Standard_Integer);
 138
 139 static
 140   void Parameters(const Handle(GeomAdaptor_HSurface)&,
 141                   const Handle(GeomAdaptor_HSurface)&,
 142                   const gp_Pnt&,
 143                   Standard_Real&,
 144                   Standard_Real&,
 145                   Standard_Real&,
 146                   Standard_Real&);
 147
 148 static
 149   void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
 150                                 const Standard_Real theTolerance,
 151                                 Standard_Real&      theumin,
 152                                 Standard_Real&      theumax,
 153                                 Standard_Real&      thevmin,
 154                                 Standard_Real&      thevmax);
 155
 156 static
 157   Standard_Boolean NotUseSurfacesForApprox
 158           (const TopoDS_Face& aF1,
 159            const TopoDS_Face& aF2,
 160            const Handle(IntPatch_WLine)& WL,
 161            const Standard_Integer ifprm,
 162            const Standard_Integer ilprm);
 163
 164 static
 165   Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)       &theWLine,
 166                                             const Handle(GeomAdaptor_HSurface) &theS1,
 167                                             const Handle(GeomAdaptor_HSurface) &theS2);
 168
 169 static
 170   Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
 171                                             const Standard_Integer                         ideb,
 172                                             const Standard_Integer                         ifin,
 173                                             const Standard_Boolean                         onFirst);
 174
 175 static
 176   Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
 177                                         const Handle(GeomAdaptor_HSurface)&            theSurface1,
 178                                         const Handle(GeomAdaptor_HSurface)&            theSurface2,
 179                                         const TopoDS_Face&                             theFace1,
 180                                         const TopoDS_Face&                             theFace2,
 181                                         const GeomInt_LineConstructor&                 theLConstructor,
 182                                         const Standard_Boolean                         theAvoidLConstructor,
 183                                         IntPatch_SequenceOfLine&                       theNewLines,
 184                                         Standard_Real&                                 theReachedTol3d,
 185                                         const Handle(IntTools_Context)& );
 186
 187 static
 188   Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
 189                                           const Handle(Geom_Curve)& theCurve,
 190                                           const TopoDS_Face&        theFace1,
 191                                           const TopoDS_Face&        theFace2,
 192                                           const Standard_Real       theOtherParameter,
 193                                           const Standard_Boolean    bIncreasePar,
 194                                           Standard_Real&            theNewParameter,
 195                                           const Handle(IntTools_Context)& );
 196
 197 static
 198   Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve);
 199
 200 static
 201   Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
 202                                      const Standard_Real theFirstBoundary,
 203                                      const Standard_Real theSecondBoundary,
 204                                      const Standard_Real theResolution,
 205                                      Standard_Boolean&   IsOnFirstBoundary);
 206 static
 207   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
 208                              const gp_Pnt2d&     theLastPoint,
 209                              const Standard_Real theUmin,
 210                              const Standard_Real theUmax,
 211                              const Standard_Real theVmin,
 212                              const Standard_Real theVmax,
 213                              gp_Pnt2d&           theNewPoint);
 214
 215
 216 static
 217   Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
 218                                        const Handle(GeomAdaptor_HSurface)&  theSurface2,
 219                                        const TopoDS_Face&                   theFace1,
 220                                        const TopoDS_Face&                   theFace2,
 221                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
 222                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
 223                                        Handle(TColStd_HArray1OfReal)&       theResultRadius,
 224                                        const Handle(IntTools_Context)& );
 225
 226 static
 227   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
 228                              const gp_Pnt2d&     theLastPoint,
 229                              const Standard_Real theUmin,
 230                              const Standard_Real theUmax,
 231                              const Standard_Real theVmin,
 232                              const Standard_Real theVmax,
 233                              const gp_Pnt2d&     theTanZoneCenter,
 234                              const Standard_Real theZoneRadius,
 235                              Handle(GeomAdaptor_HSurface) theGASurface,
 236                              gp_Pnt2d&           theNewPoint);
 237
 238 static
 239   Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
 240                                    const gp_Pnt2d&     theTanZoneCenter,
 241                                    const Standard_Real theZoneRadius,
 242                                    Handle(GeomAdaptor_HSurface) theGASurface);
 243
 244 static
 245   gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
 246                              const gp_Pnt2d&     theOriginalPoint,
 247                              Handle(GeomAdaptor_HSurface) theGASurface);
 248 static
 249   Standard_Boolean  ApproxWithPCurves(const gp_Cylinder& theCyl,
 250                                       const gp_Sphere& theSph);
 251
 252 static void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
 253                            const Handle(GeomAdaptor_HSurface)& theS2,
 254                            const Standard_Real TolAng,
 255                            const Standard_Real TolTang,
 256                            const Standard_Boolean theApprox1,
 257                            const Standard_Boolean theApprox2,
 258                            IntTools_SequenceOfCurves& theSeqOfCurve,
 259                            Standard_Boolean& theTangentFaces);
 260
 261 static Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS,
 262                                       const gp_Lin2d& theLin2d,
 263                                       const Standard_Real theTol,
 264                                       Standard_Real& theP1,
 265                                       Standard_Real& theP2);
 266 //
 267 static
 268   void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
 269                         const Handle(GeomAdaptor_HSurface)& aHS2,
 270                         Standard_Integer& iDegMin,
 271                         Standard_Integer& iNbIter,
 272                         Standard_Integer& iDegMax);
 273
 274 static
 275   void Tolerances(const Handle(GeomAdaptor_HSurface)& aHS1,
 276                   const Handle(GeomAdaptor_HSurface)& aHS2,
 277                   Standard_Real& aTolTang);
 278
 279 static
 280   Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
 281                              const GeomAbs_SurfaceType aType2);
 282 static
 283   Standard_Integer IndexType(const GeomAbs_SurfaceType aType);
 284
 285 //
 286 static
 287   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
 288                                const TopoDS_Face& aFace);
 289
 290 static
 291   Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
 292                             const Standard_Real aT,
 293                             GeomAPI_ProjectPointOnSurf& theProjPS);
 294
 295 static
 296   Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
 297                                 const Standard_Real theFirst,
 298                                 const Standard_Real theLast,
 299                                 GeomAPI_ProjectPointOnSurf& theProjPS,
 300                                 const Standard_Real theEps);
 301
 302 static
 303   Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
 304                                 const Standard_Real theFirst,
 305                                 const Standard_Real theLast,
 306                                 const TopoDS_Face& theFace,
 307                                 const Handle(IntTools_Context)& theContext);
 308
 309 static
 310   void CorrectPlaneBoundaries(Standard_Real& aUmin,
 311                               Standard_Real& aUmax,
 312                               Standard_Real& aVmin,
 313                               Standard_Real& aVmax);
 314
 315 //=======================================================================
 316 //function :
 317 //purpose  :
 318 //=======================================================================
 319 IntTools_FaceFace::IntTools_FaceFace()
 320 {
 321   myIsDone=Standard_False;
 322   myTangentFaces=Standard_False;
 323   //
 324   myHS1 = new GeomAdaptor_HSurface ();
 325   myHS2 = new GeomAdaptor_HSurface ();
 326   myTolReached2d=0.;
 327   myTolReached3d=0.;
 328   SetParameters(Standard_True, Standard_True, Standard_True, 1.e-07);
 329
 330 }
 331 //=======================================================================
 332 //function : SetContext
 333 //purpose  :
 334 //=======================================================================
 335 void IntTools_FaceFace::SetContext(const Handle(IntTools_Context)& aContext)
 336 {
 337   myContext=aContext;
 338 }
 339 //=======================================================================
 340 //function : Context
 341 //purpose  :
 342 //=======================================================================
 343 const Handle(IntTools_Context)& IntTools_FaceFace::Context()const
 344 {
 345   return myContext;
 346 }
 347 //=======================================================================
 348 //function : Face1
 349 //purpose  :
 350 //=======================================================================
 351 const TopoDS_Face& IntTools_FaceFace::Face1() const
 352 {
 353   return myFace1;
 354 }
 355 //=======================================================================
 356 //function : Face2
 357 //purpose  :
 358 //=======================================================================
 359 const TopoDS_Face& IntTools_FaceFace::Face2() const
 360 {
 361   return myFace2;
 362 }
 363 //=======================================================================
 364 //function : TangentFaces
 365 //purpose  :
 366 //=======================================================================
 367 Standard_Boolean IntTools_FaceFace::TangentFaces() const
 368 {
 369   return myTangentFaces;
 370 }
 371 //=======================================================================
 372 //function : Points
 373 //purpose  :
 374 //=======================================================================
 375 const IntTools_SequenceOfPntOn2Faces& IntTools_FaceFace::Points() const
 376 {
 377   return myPnts;
 378 }
 379 //=======================================================================
 380 //function : IsDone
 381 //purpose  :
 382 //=======================================================================
 383 Standard_Boolean IntTools_FaceFace::IsDone() const
 384 {
 385   return myIsDone;
 386 }
 387 //=======================================================================
 388 //function : TolReached3d
 389 //purpose  :
 390 //=======================================================================
 391 Standard_Real IntTools_FaceFace::TolReached3d() const
 392 {
 393   return myTolReached3d;
 394 }
 395 //=======================================================================
 396 //function : Lines
 397 //purpose  : return lines of intersection
 398 //=======================================================================
 399 const IntTools_SequenceOfCurves& IntTools_FaceFace::Lines() const
 400 {
 401   StdFail_NotDone_Raise_if
 402     (!myIsDone,
 403      "IntTools_FaceFace::Lines() => myIntersector NOT DONE");
 404   return mySeqOfCurve;
 405 }
 406 //=======================================================================
 407 //function : TolReached2d
 408 //purpose  :
 409 //=======================================================================
 410 Standard_Real IntTools_FaceFace::TolReached2d() const
 411 {
 412   return myTolReached2d;
 413 }
 414 // =======================================================================
 415 // function: SetParameters
 416 //
 417 // =======================================================================
 418 void IntTools_FaceFace::SetParameters(const Standard_Boolean ToApproxC3d,
 419                                       const Standard_Boolean ToApproxC2dOnS1,
 420                                       const Standard_Boolean ToApproxC2dOnS2,
 421                                       const Standard_Real ApproximationTolerance)
 422 {
 423   myApprox = ToApproxC3d;
 424   myApprox1 = ToApproxC2dOnS1;
 425   myApprox2 = ToApproxC2dOnS2;
 426   myTolApprox = ApproximationTolerance;
 427 }
 428 //=======================================================================
 429 //function : SetList
 430 //purpose  :
 431 //=======================================================================
 432 void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
 433 {
 434   myListOfPnts = aListOfPnts;
 435 }
 436
 437
 438 static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
 439                                         const TopoDS_Face& theF2)
 440 {
 441   const Standard_Real Tolang = 1.e-8;
 442   const Standard_Real aTolF1=BRep_Tool::Tolerance(theF1);
 443   const Standard_Real aTolF2=BRep_Tool::Tolerance(theF2);
 444   const Standard_Real aTolSum = aTolF1 + aTolF2;
 445   Standard_Real aHigh = 0.0;
 446
 447   const BRepAdaptor_Surface aBAS1(theF1), aBAS2(theF2);
 448   const GeomAbs_SurfaceType aType1=aBAS1.GetType();
 449   const GeomAbs_SurfaceType aType2=aBAS2.GetType();
 450
 451   gp_Pln aS1;
 452   gp_Cylinder aS2;
 453   if(aType1 == GeomAbs_Plane)
 454   {
 455     aS1=aBAS1.Plane();
 456   }
 457   else if(aType2 == GeomAbs_Plane)
 458   {
 459     aS1=aBAS2.Plane();
 460   }
 461   else
 462   {
 463     return Standard_True;
 464   }
 465
 466   if(aType1 == GeomAbs_Cylinder)
 467   {
 468     aS2=aBAS1.Cylinder();
 469     const Standard_Real VMin = aBAS1.FirstVParameter();
 470     const Standard_Real VMax = aBAS1.LastVParameter();
 471
 472     if( Precision::IsNegativeInfinite(VMin) ||
 473         Precision::IsPositiveInfinite(VMax))
 474           return Standard_True;
 475     else
 476       aHigh = VMax - VMin;
 477   }
 478   else if(aType2 == GeomAbs_Cylinder)
 479   {
 480     aS2=aBAS2.Cylinder();
 481
 482     const Standard_Real VMin = aBAS2.FirstVParameter();
 483     const Standard_Real VMax = aBAS2.LastVParameter();
 484
 485     if( Precision::IsNegativeInfinite(VMin) ||
 486         Precision::IsPositiveInfinite(VMax))
 487           return Standard_True;
 488     else
 489       aHigh = VMax - VMin;
 490   }
 491   else
 492   {
 493     return Standard_True;
 494   }
 495
 496   IntAna_QuadQuadGeo inter;
 497   inter.Perform(aS1,aS2,Tolang,aTolSum, aHigh);
 498   if(inter.TypeInter() == IntAna_Ellipse)
 499   {
 500     const gp_Elips anEl = inter.Ellipse(1);
 501     const Standard_Real aMajorR = anEl.MajorRadius();
 502     const Standard_Real aMinorR = anEl.MinorRadius();
 503
 504     return (aMajorR < 100000.0 * aMinorR);
 505   }
 506   else
 507   {
 508     return inter.IsDone();
 509   }
 510 }
 511 //=======================================================================
 512 //function : Perform
 513 //purpose  : intersect surfaces of the faces
 514 //=======================================================================
 515 void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
 516                                 const TopoDS_Face& aF2)
 517 {
 518   Standard_Boolean RestrictLine = Standard_False, hasCone = Standard_False;
 519
 520   if (myContext.IsNull()) {
 521     myContext=new IntTools_Context;
 522   }
 523
 524   mySeqOfCurve.Clear();
 525   myTolReached2d=0.;
 526   myTolReached3d=0.;
 527   myIsDone = Standard_False;
 528   myNbrestr=0;//?
 529
 530   myFace1=aF1;
 531   myFace2=aF2;
 532
 533   const BRepAdaptor_Surface aBAS1(myFace1, Standard_False);
 534   const BRepAdaptor_Surface aBAS2(myFace2, Standard_False);
 535   GeomAbs_SurfaceType aType1=aBAS1.GetType();
 536   GeomAbs_SurfaceType aType2=aBAS2.GetType();
 537
 538   const Standard_Boolean bReverse=SortTypes(aType1, aType2);
 539   if (bReverse)
 540   {
 541     myFace1=aF2;
 542     myFace2=aF1;
 543     aType1=aBAS2.GetType();
 544     aType2=aBAS1.GetType();
 545
 546     if (myListOfPnts.Extent())
 547     {
 548       Standard_Real aU1,aV1,aU2,aV2;
 549       IntSurf_ListIteratorOfListOfPntOn2S aItP2S;
 550       //
 551       aItP2S.Initialize(myListOfPnts);
 552       for (; aItP2S.More(); aItP2S.Next())
 553       {
 554         IntSurf_PntOn2S& aP2S=aItP2S.Value();
 555         aP2S.Parameters(aU1,aV1,aU2,aV2);
 556         aP2S.SetValue(aU2,aV2,aU1,aV1);
 557       }
 558     }
 559     //
 560     Standard_Boolean anAproxTmp = myApprox1;
 561     myApprox1 = myApprox2;
 562     myApprox2 = anAproxTmp;
 563   }
 564
 565
 566   const Handle(Geom_Surface) S1=BRep_Tool::Surface(myFace1);
 567   const Handle(Geom_Surface) S2=BRep_Tool::Surface(myFace2);
 568
 569   const Standard_Real aTolF1=BRep_Tool::Tolerance(myFace1);
 570   const Standard_Real aTolF2=BRep_Tool::Tolerance(myFace2);
 571
 572   Standard_Real TolArc = aTolF1 + aTolF2;
 573   Standard_Real TolTang = TolArc;
 574
 575   const Standard_Boolean isFace1Quad = (aType1 == GeomAbs_Cylinder ||
 576                                         aType1 == GeomAbs_Cone ||
 577                                         aType1 == GeomAbs_Torus);
 578
 579   const Standard_Boolean isFace2Quad = (aType2 == GeomAbs_Cylinder ||
 580                                         aType2 == GeomAbs_Cone ||
 581                                         aType2 == GeomAbs_Torus);
 582
 583   if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane)  {
 584     Standard_Real umin, umax, vmin, vmax;
 585     //
 586     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 587     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 588     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 589     //
 590     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 591     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 592     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 593     //
 594     Standard_Real TolAng = 1.e-8;
 595     //
 596     PerformPlanes(myHS1, myHS2, TolAng, TolTang, myApprox1, myApprox2,
 597                   mySeqOfCurve, myTangentFaces);
 598     //
 599     myIsDone = Standard_True;
 600
 601     if(!myTangentFaces) {
 602       const Standard_Integer NbLinPP = mySeqOfCurve.Length();
 603       if(NbLinPP) {
 604         Standard_Real aTolFMax;
 605         myTolReached3d = 1.e-7;
 606         aTolFMax=Max(aTolF1, aTolF2);
 607         if (aTolFMax>myTolReached3d) {
 608           myTolReached3d=aTolFMax;
 609         }
 610         //
 611         myTolReached2d = myTolReached3d;
 612
 613         if (bReverse) {
 614           Handle(Geom2d_Curve) aC2D1, aC2D2;
 615           const Standard_Integer aNbLin = mySeqOfCurve.Length();
 616           for (Standard_Integer i = 1; i <= aNbLin; ++i) {
 617             IntTools_Curve& aIC=mySeqOfCurve(i);
 618             aC2D1=aIC.FirstCurve2d();
 619             aC2D2=aIC.SecondCurve2d();
 620             aIC.SetFirstCurve2d(aC2D2);
 621             aIC.SetSecondCurve2d(aC2D1);
 622           }
 623         }
 624       }
 625     }
 626     return;
 627   }//if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane){
 628
 629   if ((aType1==GeomAbs_Plane) && isFace2Quad)
 630   {
 631     Standard_Real umin, umax, vmin, vmax;
 632     // F1
 633     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 634     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 635     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 636     // F2
 637     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 638     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 639     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 640     //
 641     if( aType2==GeomAbs_Cone ) {
 642       TolArc = 0.0001;
 643       hasCone = Standard_True;
 644     }
 645   }
 646   else if ((aType2==GeomAbs_Plane) && isFace1Quad)
 647   {
 648     Standard_Real umin, umax, vmin, vmax;
 649     //F1
 650     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 651     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 652     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 653     // F2
 654     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 655     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 656     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 657     //
 658     if( aType1==GeomAbs_Cone ) {
 659       TolArc = 0.0001;
 660       hasCone = Standard_True;
 661     }
 662   }
 663   else
 664   {
 665     Standard_Real umin, umax, vmin, vmax;
 666     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 667     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 668     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 669     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 670     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 671     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 672   }
 673
 674   const Handle(IntTools_TopolTool) dom1 = new IntTools_TopolTool(myHS1);
 675   const Handle(IntTools_TopolTool) dom2 = new IntTools_TopolTool(myHS2);
 676
 677   myLConstruct.Load(dom1, dom2, myHS1, myHS2);
 678
 679
 680   Tolerances(myHS1, myHS2, TolTang);
 681
 682   {
 683     const Standard_Real UVMaxStep = 0.001;
 684     const Standard_Real Deflection = (hasCone) ? 0.085 : 0.1;
 685   myIntersector.SetTolerances(TolArc, TolTang, UVMaxStep, Deflection);
 686   }
 687
 688   if((myHS1->IsUClosed() && !myHS1->IsUPeriodic()) ||
 689      (myHS1->IsVClosed() && !myHS1->IsVPeriodic()) ||
 690      (myHS2->IsUClosed() && !myHS2->IsUPeriodic()) ||
 691      (myHS2->IsVClosed() && !myHS2->IsVPeriodic()))
 692   {
 693     RestrictLine = Standard_True;
 694   }
 695   //
 696   if((aType1 != GeomAbs_BSplineSurface) &&
 697       (aType1 != GeomAbs_BezierSurface)  &&
 698      (aType1 != GeomAbs_OtherSurface)  &&
 699      (aType2 != GeomAbs_BSplineSurface) &&
 700       (aType2 != GeomAbs_BezierSurface)  &&
 701      (aType2 != GeomAbs_OtherSurface))
 702   {
 703     RestrictLine = Standard_True;
 704
 705     if ((aType1 == GeomAbs_Torus) ||
 706         (aType2 == GeomAbs_Torus))
 707     {
 708       myListOfPnts.Clear();
 709     }
 710   }
 711
 712   //
 713   if(!RestrictLine)
 714   {
 715     TopExp_Explorer aExp;
 716     for(Standard_Integer i = 0; (!RestrictLine) && (i < 2); i++)
 717     {
 718       const TopoDS_Face& aF=(!i) ? myFace1 : myFace2;
 719       aExp.Init(aF, TopAbs_EDGE);
 720       for(; aExp.More(); aExp.Next())
 721       {
 722         const TopoDS_Edge& aE=TopoDS::Edge(aExp.Current());
 723
 724         if(BRep_Tool::Degenerated(aE))
 725         {
 726           RestrictLine = Standard_True;
 727           break;
 728         }
 729       }
 730     }
 731   }
 732
 733   const Standard_Boolean isGeomInt = isTreatAnalityc(aF1, aF2);
 734   myIntersector.Perform(myHS1, dom1, myHS2, dom2, TolArc, TolTang,
 735                                   myListOfPnts, RestrictLine, isGeomInt);
 736
 737   myIsDone = myIntersector.IsDone();
 738
 739   if (myIsDone)
 740   {
 741     myTangentFaces=myIntersector.TangentFaces();
 742     if (myTangentFaces) {
 743       return;
 744     }
 745     //
 746     if(RestrictLine) {
 747       myListOfPnts.Clear(); // to use LineConstructor
 748     }
 749     //
 750     const Standard_Integer aNbLin = myIntersector.NbLines();
 751     for (Standard_Integer i=1; i <= aNbLin; ++i) {
 752       MakeCurve(i, dom1, dom2);
 753     }
 754     //
 755     ComputeTolReached3d();
 756     //
 757     if (bReverse) {
 758       Handle(Geom2d_Curve) aC2D1, aC2D2;
 759       //
 760       const Standard_Integer aNbLin=mySeqOfCurve.Length();
 761       for (Standard_Integer i=1; i<=aNbLin; ++i)
 762       {
 763         IntTools_Curve& aIC=mySeqOfCurve(i);
 764         aC2D1=aIC.FirstCurve2d();
 765         aC2D2=aIC.SecondCurve2d();
 766         aIC.SetFirstCurve2d(aC2D2);
 767         aIC.SetSecondCurve2d(aC2D1);
 768       }
 769     }
 770
 771     // Points
 772     Standard_Boolean bValid2D1, bValid2D2;
 773     Standard_Real U1,V1,U2,V2;
 774     IntTools_PntOnFace aPntOnF1, aPntOnF2;
 775     IntTools_PntOn2Faces aPntOn2Faces;
 776     //
 777     const Standard_Integer aNbPnts = myIntersector.NbPnts();
 778     for (Standard_Integer i=1; i <= aNbPnts; ++i)
 779     {
 780       const IntSurf_PntOn2S& aISPnt=myIntersector.Point(i).PntOn2S();
 781       const gp_Pnt& aPnt=aISPnt.Value();
 782       aISPnt.Parameters(U1,V1,U2,V2);
 783       //
 784       // check the validity of the intersection point for the faces
 785       bValid2D1 = myContext->IsPointInOnFace(myFace1, gp_Pnt2d(U1, V1));
 786       if (!bValid2D1) {
 787         continue;
 788       }
 789       //
 790       bValid2D2 = myContext->IsPointInOnFace(myFace2, gp_Pnt2d(U2, V2));
 791       if (!bValid2D2) {
 792         continue;
 793       }
 794       //
 795       // add the intersection point
 796       aPntOnF1.Init(myFace1, aPnt, U1, V1);
 797       aPntOnF2.Init(myFace2, aPnt, U2, V2);
 798       //
 799       if (!bReverse)
 800       {
 801         aPntOn2Faces.SetP1(aPntOnF1);
 802         aPntOn2Faces.SetP2(aPntOnF2);
 803       }
 804       else
 805       {
 806         aPntOn2Faces.SetP2(aPntOnF1);
 807         aPntOn2Faces.SetP1(aPntOnF2);
 808       }
 809
 810       myPnts.Append(aPntOn2Faces);
 811     }
 812   }
 813 }
 814
 815 //=======================================================================
 816 //function : ComputeTolerance
 817 //purpose  :
 818 //=======================================================================
 819 Standard_Real IntTools_FaceFace::ComputeTolerance()
 820 {
 821   Standard_Integer i, j, aNbLin;
 822   Standard_Real aFirst, aLast, aD, aDMax, aT;
 823   Handle(Geom_Surface) aS1, aS2;
 824   //
 825   aDMax = 0;
 826   aNbLin = mySeqOfCurve.Length();
 827   //
 828   aS1 = myHS1->ChangeSurface().Surface();
 829   aS2 = myHS2->ChangeSurface().Surface();
 830   //
 831   for (i = 1; i <= aNbLin; ++i)
 832   {
 833     const IntTools_Curve& aIC = mySeqOfCurve(i);
 834     const Handle(Geom_Curve)& aC3D = aIC.Curve();
 835     if (aC3D.IsNull())
 836     {
 837       continue;
 838     }
 839     //
 840     aFirst = aC3D->FirstParameter();
 841     aLast  = aC3D->LastParameter();
 842     //
 843     const Handle(Geom2d_Curve)& aC2D1 = aIC.FirstCurve2d();
 844     const Handle(Geom2d_Curve)& aC2D2 = aIC.SecondCurve2d();
 845     //
 846     for (j = 0; j < 2; ++j)
 847     {
 848       const Handle(Geom2d_Curve)& aC2D = !j ? aC2D1 : aC2D2;
 849       const Handle(Geom_Surface)& aS = !j ? aS1 : aS2;
 850       //
 851       if (!aC2D.IsNull())
 852       {
 853         if (IntTools_Tools::ComputeTolerance
 854             (aC3D, aC2D, aS, aFirst, aLast, aD, aT))
 855         {
 856           if (aD > aDMax)
 857           {
 858             aDMax = aD;
 859           }
 860         }
 861       }
 862       else
 863       {
 864         const TopoDS_Face& aF = !j ? myFace1 : myFace2;
 865         aD = FindMaxDistance(aC3D, aFirst, aLast, aF, myContext);
 866         if (aD > aDMax)
 867         {
 868           aDMax = aD;
 869         }
 870       }
 871     }
 872   }
 873   //
 874   return aDMax;
 875 }
 876
 877 //=======================================================================
 878 //function :ComputeTolReached3d
 879 //purpose  :
 880 //=======================================================================
 881   void IntTools_FaceFace::ComputeTolReached3d()
 882 {
 883   Standard_Integer aNbLin;
 884   GeomAbs_SurfaceType aType1, aType2;
 885   //
 886   aNbLin=myIntersector.NbLines();
 887   if (!aNbLin) {
 888     return;
 889   }
 890   //
 891   aType1=myHS1->Surface().GetType();
 892   aType2=myHS2->Surface().GetType();
 893   //
 894   if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder)
 895   {
 896     if (aNbLin==2)
 897     {
 898       Handle(IntPatch_Line) aIL1, aIL2;
 899       IntPatch_IType aTL1, aTL2;
 900       //
 901       aIL1=myIntersector.Line(1);
 902       aIL2=myIntersector.Line(2);
 903       aTL1=aIL1->ArcType();
 904       aTL2=aIL2->ArcType();
 905       if (aTL1==IntPatch_Lin && aTL2==IntPatch_Lin) {
 906         Standard_Real aD, aDTresh, dTol;
 907         gp_Lin aL1, aL2;
 908         //
 909         dTol=1.e-8;
 910         aDTresh=1.5e-6;
 911         //
 912         aL1=Handle(IntPatch_GLine)::DownCast(aIL1)->Line();
 913         aL2=Handle(IntPatch_GLine)::DownCast(aIL2)->Line();
 914         aD=aL1.Distance(aL2);
 915         aD=0.5*aD;
 916         if (aD<aDTresh)
 917         {//In order to avoid creation too thin face
 918           myTolReached3d=aD+dTol;
 919         }
 920       }
 921     }
 922   }// if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
 923   //
 924
 925   Standard_Real aDMax = ComputeTolerance();
 926   if (aDMax > myTolReached3d)
 927   {
 928       myTolReached3d = aDMax;
 929     }
 930   }
 931
 932 //=======================================================================
 933 //function : MakeCurve
 934 //purpose  :
 935 //=======================================================================
 936   void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
 937                                     const Handle(Adaptor3d_TopolTool)& dom1,
 938                                     const Handle(Adaptor3d_TopolTool)& dom2)
 939 {
 940   Standard_Boolean bDone, rejectSurface, reApprox, bAvoidLineConstructor;
 941   Standard_Boolean ok, bPCurvesOk;
 942   Standard_Integer i, j, aNbParts;
 943   Standard_Real fprm, lprm;
 944   Standard_Real Tolpc;
 945   Handle(IntPatch_Line) L;
 946   IntPatch_IType typl;
 947   Handle(Geom_Curve) newc;
 948   //
 949   const Standard_Real TOLCHECK   =0.0000001;
 950   const Standard_Real TOLANGCHECK=0.1;
 951   //
 952   rejectSurface = Standard_False;
 953   reApprox = Standard_False;
 954   //
 955   bPCurvesOk = Standard_True;
 956
 957  reapprox:;
 958
 959   Tolpc = myTolApprox;
 960   bAvoidLineConstructor = Standard_False;
 961   L = myIntersector.Line(Index);
 962   typl = L->ArcType();
 963   //
 964   if(typl==IntPatch_Walking) {
 965     Handle(IntPatch_Line) anewL;
 966     //
 967     Handle(IntPatch_WLine) aWLine (Handle(IntPatch_WLine)::DownCast(L));
 968     anewL = ComputePurgedWLine(aWLine, myHS1, myHS2);
 969     if(anewL.IsNull()) {
 970       return;
 971     }
 972     L = anewL;
 973
 974     //Handle(IntPatch_WLine) aWLineX (Handle(IntPatch_WLine)::DownCast(L));
 975     //DumpWLine(aWLineX);
 976
 977     //
 978     if(!myListOfPnts.IsEmpty()) {
 979       bAvoidLineConstructor = Standard_True;
 980     }
 981
 982     Standard_Integer nbp = aWLine->NbPnts();
 983     const IntSurf_PntOn2S& p1 = aWLine->Point(1);
 984     const IntSurf_PntOn2S& p2 = aWLine->Point(nbp);
 985
 986     const gp_Pnt& P1 = p1.Value();
 987     const gp_Pnt& P2 = p2.Value();
 988
 989     if(P1.SquareDistance(P2) < 1.e-14) {
 990       bAvoidLineConstructor = Standard_False;
 991     }
 992   }
 993
 994   typl=L->ArcType();
 995
 996   //
 997   // Line Constructor
 998   if(!bAvoidLineConstructor) {
 999     myLConstruct.Perform(L);
1000     //
1001     bDone=myLConstruct.IsDone();
1002     if(!bDone)
1003     {
1004       return;
1005     }
1006
1007     if(typl != IntPatch_Restriction)
1008     {
1009       aNbParts=myLConstruct.NbParts();
1010       if (aNbParts <= 0)
1011       {
1012         return;
1013       }
1014     }
1015   }
1016   // Do the Curve
1017
1018
1019   switch (typl) {
1020   //########################################
1021   // Line, Parabola, Hyperbola
1022   //########################################
1023   case IntPatch_Lin:
1024   case IntPatch_Parabola:
1025   case IntPatch_Hyperbola: {
1026     if (typl == IntPatch_Lin) {
1027       newc =
1028         new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
1029     }
1030
1031     else if (typl == IntPatch_Parabola) {
1032       newc =
1033         new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
1034     }
1035
1036     else if (typl == IntPatch_Hyperbola) {
1037       newc =
1038         new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
1039     }
1040     //
1041     // myTolReached3d
1042     if (typl == IntPatch_Lin) {
1043       TolR3d (myFace1, myFace2, myTolReached3d);
1044     }
1045     //
1046     aNbParts=myLConstruct.NbParts();
1047     for (i=1; i<=aNbParts; i++) {
1048       Standard_Boolean bFNIt, bLPIt;
1049       //
1050       myLConstruct.Part(i, fprm, lprm);
1051         //
1052       bFNIt=Precision::IsNegativeInfinite(fprm);
1053       bLPIt=Precision::IsPositiveInfinite(lprm);
1054       //
1055       if (!bFNIt && !bLPIt) {
1056         //
1057         IntTools_Curve aCurve;
1058         //
1059         Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
1060         aCurve.SetCurve(aCT3D);
1061         if (typl == IntPatch_Parabola) {
1062           Standard_Real aTolF1, aTolF2, aTolBase;
1063
1064           aTolF1 = BRep_Tool::Tolerance(myFace1);
1065           aTolF2 = BRep_Tool::Tolerance(myFace2);
1066           aTolBase=aTolF1+aTolF2;
1067           myTolReached3d=IntTools_Tools::CurveTolerance(aCT3D, aTolBase);
1068         }
1069         //
1070         aCurve.SetCurve(new Geom_TrimmedCurve(newc, fprm, lprm));
1071         if(myApprox1) {
1072           Handle (Geom2d_Curve) C2d;
1073           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1074                 myHS1->ChangeSurface().Surface(), newc, C2d);
1075           if(Tolpc>myTolReached2d || myTolReached2d==0.) {
1076             myTolReached2d=Tolpc;
1077           }
1078             //
1079             aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
1080           }
1081           else {
1082             Handle(Geom2d_BSplineCurve) H1;
1083             //
1084             aCurve.SetFirstCurve2d(H1);
1085           }
1086         //
1087         if(myApprox2) {
1088           Handle (Geom2d_Curve) C2d;
1089           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1090                     myHS2->ChangeSurface().Surface(), newc, C2d);
1091           if(Tolpc>myTolReached2d || myTolReached2d==0.) {
1092             myTolReached2d=Tolpc;
1093           }
1094           //
1095           aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
1096           }
1097         else {
1098           Handle(Geom2d_BSplineCurve) H1;
1099           //
1100           aCurve.SetSecondCurve2d(H1);
1101         }
1102         mySeqOfCurve.Append(aCurve);
1103       } //if (!bFNIt && !bLPIt) {
1104       else {
1105         //  on regarde si on garde
1106         //
1107         Standard_Real aTestPrm, dT=100.;
1108         //
1109         aTestPrm=0.;
1110         if (bFNIt && !bLPIt) {
1111           aTestPrm=lprm-dT;
1112         }
1113         else if (!bFNIt && bLPIt) {
1114           aTestPrm=fprm+dT;
1115         }
1116         else {
1117           // i.e, if (bFNIt && bLPIt)
1118           aTestPrm=IntTools_Tools::IntermediatePoint(-dT, dT);
1119         }
1120         //
1121         gp_Pnt ptref(newc->Value(aTestPrm));
1122         //
1123         GeomAbs_SurfaceType typS1 = myHS1->GetType();
1124         GeomAbs_SurfaceType typS2 = myHS2->GetType();
1125         if( typS1 == GeomAbs_SurfaceOfExtrusion ||
1126             typS1 == GeomAbs_OffsetSurface ||
1127             typS1 == GeomAbs_SurfaceOfRevolution ||
1128             typS2 == GeomAbs_SurfaceOfExtrusion ||
1129             typS2 == GeomAbs_OffsetSurface ||
1130             typS2 == GeomAbs_SurfaceOfRevolution) {
1131           Handle(Geom2d_BSplineCurve) H1;
1132           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1133           continue;
1134         }
1135
1136         Standard_Real u1, v1, u2, v2, Tol;
1137
1138         Tol = Precision::Confusion();
1139         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
1140         ok = (dom1->Classify(gp_Pnt2d(u1, v1), Tol) != TopAbs_OUT);
1141         if(ok) {
1142           ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1143         }
1144         if (ok) {
1145           Handle(Geom2d_BSplineCurve) H1;
1146           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1147         }
1148       }
1149     }// for (i=1; i<=aNbParts; i++) {
1150   }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
1151     break;
1152
1153   //########################################
1154   // Circle and Ellipse
1155   //########################################
1156   case IntPatch_Circle:
1157   case IntPatch_Ellipse: {
1158
1159     if (typl == IntPatch_Circle) {
1160       newc = new Geom_Circle
1161         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
1162     }
1163     else { //IntPatch_Ellipse
1164       newc = new Geom_Ellipse
1165         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
1166     }
1167     //
1168     // myTolReached3d
1169     TolR3d (myFace1, myFace2, myTolReached3d);
1170     //
1171     aNbParts=myLConstruct.NbParts();
1172     //
1173     Standard_Real aPeriod, aNul;
1174     TColStd_SequenceOfReal aSeqFprm,  aSeqLprm;
1175
1176     aNul=0.;
1177     aPeriod=M_PI+M_PI;
1178
1179     for (i=1; i<=aNbParts; i++) {
1180       myLConstruct.Part(i, fprm, lprm);
1181
1182       if (fprm < aNul && lprm > aNul) {
1183         // interval that goes through 0. is divided on two intervals;
1184         while (fprm<aNul || fprm>aPeriod) fprm=fprm+aPeriod;
1185         while (lprm<aNul || lprm>aPeriod) lprm=lprm+aPeriod;
1186         //
1187         if((aPeriod - fprm) > Tolpc) {
1188           aSeqFprm.Append(fprm);
1189           aSeqLprm.Append(aPeriod);
1190         }
1191         else {
1192           gp_Pnt P1 = newc->Value(fprm);
1193           gp_Pnt P2 = newc->Value(aPeriod);
1194           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1195           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1196
1197           if(P1.Distance(P2) > aTolDist) {
1198             Standard_Real anewpar = fprm;
1199
1200             if(ParameterOutOfBoundary(fprm, newc, myFace1, myFace2,
1201                                       lprm, Standard_False, anewpar, myContext)) {
1202               fprm = anewpar;
1203             }
1204             aSeqFprm.Append(fprm);
1205             aSeqLprm.Append(aPeriod);
1206           }
1207         }
1208
1209         //
1210         if((lprm - aNul) > Tolpc) {
1211           aSeqFprm.Append(aNul);
1212           aSeqLprm.Append(lprm);
1213         }
1214         else {
1215           gp_Pnt P1 = newc->Value(aNul);
1216           gp_Pnt P2 = newc->Value(lprm);
1217           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1218           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1219
1220           if(P1.Distance(P2) > aTolDist) {
1221             Standard_Real anewpar = lprm;
1222
1223             if(ParameterOutOfBoundary(lprm, newc, myFace1, myFace2,
1224                                       fprm, Standard_True, anewpar, myContext)) {
1225               lprm = anewpar;
1226             }
1227             aSeqFprm.Append(aNul);
1228             aSeqLprm.Append(lprm);
1229           }
1230         }
1231       }
1232       else {
1233         // usual interval
1234         aSeqFprm.Append(fprm);
1235         aSeqLprm.Append(lprm);
1236       }
1237     }
1238
1239     //
1240     aNbParts=aSeqFprm.Length();
1241     for (i=1; i<=aNbParts; i++) {
1242       fprm=aSeqFprm(i);
1243       lprm=aSeqLprm(i);
1244       //
1245       Standard_Real aRealEpsilon=RealEpsilon();
1246       if (Abs(fprm) > aRealEpsilon || Abs(lprm-2.*M_PI) > aRealEpsilon) {
1247         //==============================================
1248         ////
1249         IntTools_Curve aCurve;
1250         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1251         aCurve.SetCurve(aTC3D);
1252         fprm=aTC3D->FirstParameter();
1253         lprm=aTC3D->LastParameter ();
1254         ////
1255         if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {////
1256           if(myApprox1) {
1257             Handle (Geom2d_Curve) C2d;
1258             GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1259                         myHS1->ChangeSurface().Surface(), newc, C2d);
1260             if(Tolpc>myTolReached2d || myTolReached2d==0) {
1261               myTolReached2d=Tolpc;
1262             }
1263             //
1264             aCurve.SetFirstCurve2d(C2d);
1265           }
1266           else { ////
1267             Handle(Geom2d_BSplineCurve) H1;
1268             aCurve.SetFirstCurve2d(H1);
1269           }
1270
1271
1272           if(myApprox2) {
1273             Handle (Geom2d_Curve) C2d;
1274             GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1275                     myHS2->ChangeSurface().Surface(),newc,C2d);
1276             if(Tolpc>myTolReached2d || myTolReached2d==0) {
1277               myTolReached2d=Tolpc;
1278             }
1279             //
1280             aCurve.SetSecondCurve2d(C2d);
1281           }
1282           else {
1283             Handle(Geom2d_BSplineCurve) H1;
1284             aCurve.SetSecondCurve2d(H1);
1285           }
1286         }
1287
1288         else {
1289           Handle(Geom2d_BSplineCurve) H1;
1290           aCurve.SetFirstCurve2d(H1);
1291           aCurve.SetSecondCurve2d(H1);
1292         }
1293         mySeqOfCurve.Append(aCurve);
1294           //==============================================
1295       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
1296
1297       else {
1298         //  on regarde si on garde
1299         //
1300         if (aNbParts==1) {
1301 //           if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
1302           if (Abs(fprm) <= aRealEpsilon && Abs(lprm-2.*M_PI) <= aRealEpsilon) {
1303             IntTools_Curve aCurve;
1304             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1305             aCurve.SetCurve(aTC3D);
1306             fprm=aTC3D->FirstParameter();
1307             lprm=aTC3D->LastParameter ();
1308
1309             if(myApprox1) {
1310               Handle (Geom2d_Curve) C2d;
1311               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1312                     myHS1->ChangeSurface().Surface(),newc,C2d);
1313               if(Tolpc>myTolReached2d || myTolReached2d==0) {
1314                 myTolReached2d=Tolpc;
1315               }
1316               //
1317               aCurve.SetFirstCurve2d(C2d);
1318             }
1319             else { ////
1320               Handle(Geom2d_BSplineCurve) H1;
1321               aCurve.SetFirstCurve2d(H1);
1322             }
1323
1324             if(myApprox2) {
1325               Handle (Geom2d_Curve) C2d;
1326               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1327                     myHS2->ChangeSurface().Surface(),newc,C2d);
1328               if(Tolpc>myTolReached2d || myTolReached2d==0) {
1329                 myTolReached2d=Tolpc;
1330               }
1331               //
1332               aCurve.SetSecondCurve2d(C2d);
1333             }
1334             else {
1335               Handle(Geom2d_BSplineCurve) H1;
1336               aCurve.SetSecondCurve2d(H1);
1337             }
1338             mySeqOfCurve.Append(aCurve);
1339             break;
1340           }
1341         }
1342         //
1343         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, Tol;
1344
1345         aTwoPIdiv17=2.*M_PI/17.;
1346
1347         for (j=0; j<=17; j++) {
1348           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
1349           Tol = Precision::Confusion();
1350
1351           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
1352           ok = (dom1->Classify(gp_Pnt2d(u1,v1),Tol) != TopAbs_OUT);
1353           if(ok) {
1354             ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1355           }
1356           if (ok) {
1357             IntTools_Curve aCurve;
1358             aCurve.SetCurve(newc);
1359             //==============================================
1360             if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {
1361
1362               if(myApprox1) {
1363                 Handle (Geom2d_Curve) C2d;
1364                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1365                         myHS1->ChangeSurface().Surface(), newc, C2d);
1366                 if(Tolpc>myTolReached2d || myTolReached2d==0) {
1367                   myTolReached2d=Tolpc;
1368                 }
1369                 //
1370                 aCurve.SetFirstCurve2d(C2d);
1371               }
1372               else {
1373                 Handle(Geom2d_BSplineCurve) H1;
1374                 aCurve.SetFirstCurve2d(H1);
1375               }
1376
1377               if(myApprox2) {
1378                 Handle (Geom2d_Curve) C2d;
1379                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1380                         myHS2->ChangeSurface().Surface(), newc, C2d);
1381                 if(Tolpc>myTolReached2d || myTolReached2d==0) {
1382                   myTolReached2d=Tolpc;
1383                 }
1384                 //
1385                 aCurve.SetSecondCurve2d(C2d);
1386               }
1387
1388               else {
1389                 Handle(Geom2d_BSplineCurve) H1;
1390                 aCurve.SetSecondCurve2d(H1);
1391               }
1392             }//  end of if (typl == IntPatch_Circle || typl == IntPatch_Ellipse)
1393
1394             else {
1395               Handle(Geom2d_BSplineCurve) H1;
1396               //
1397               aCurve.SetFirstCurve2d(H1);
1398               aCurve.SetSecondCurve2d(H1);
1399             }
1400             //==============================================
1401             //
1402             mySeqOfCurve.Append(aCurve);
1403             break;
1404
1405             }//  end of if (ok) {
1406           }//  end of for (Standard_Integer j=0; j<=17; j++)
1407         }//  end of else { on regarde si on garde
1408       }// for (i=1; i<=myLConstruct.NbParts(); i++)
1409     }// IntPatch_Circle: IntPatch_Ellipse:
1410     break;
1411
1412   case IntPatch_Analytic: {
1413     IntSurf_Quadric quad1,quad2;
1414     GeomAbs_SurfaceType typs = myHS1->Surface().GetType();
1415
1416     switch (typs) {
1417       case GeomAbs_Plane:
1418         quad1.SetValue(myHS1->Surface().Plane());
1419         break;
1420       case GeomAbs_Cylinder:
1421         quad1.SetValue(myHS1->Surface().Cylinder());
1422         break;
1423       case GeomAbs_Cone:
1424         quad1.SetValue(myHS1->Surface().Cone());
1425         break;
1426       case GeomAbs_Sphere:
1427         quad1.SetValue(myHS1->Surface().Sphere());
1428         break;
1429     case GeomAbs_Torus:
1430       quad1.SetValue(myHS1->Surface().Torus());
1431       break;
1432       default:
1433         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 1");
1434       }
1435
1436     typs = myHS2->Surface().GetType();
1437
1438     switch (typs) {
1439       case GeomAbs_Plane:
1440         quad2.SetValue(myHS2->Surface().Plane());
1441         break;
1442       case GeomAbs_Cylinder:
1443         quad2.SetValue(myHS2->Surface().Cylinder());
1444         break;
1445       case GeomAbs_Cone:
1446         quad2.SetValue(myHS2->Surface().Cone());
1447         break;
1448       case GeomAbs_Sphere:
1449         quad2.SetValue(myHS2->Surface().Sphere());
1450         break;
1451     case GeomAbs_Torus:
1452       quad2.SetValue(myHS2->Surface().Torus());
1453       break;
1454       default:
1455         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 2");
1456       }
1457     //
1458     //=========
1459     IntPatch_ALineToWLine convert (quad1, quad2);
1460
1461     if (!myApprox) {
1462       aNbParts=myLConstruct.NbParts();
1463       for (i=1; i<=aNbParts; i++) {
1464         myLConstruct.Part(i, fprm, lprm);
1465         Handle(IntPatch_WLine) WL =
1466           convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
1467         //
1468         Handle(Geom2d_BSplineCurve) H1;
1469         Handle(Geom2d_BSplineCurve) H2;
1470
1471         if(myApprox1) {
1472           H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1473         }
1474
1475         if(myApprox2) {
1476           H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1477         }
1478         //
1479         mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1480       }
1481     } // if (!myApprox)
1482
1483     else { // myApprox=TRUE
1484       GeomInt_WLApprox theapp3d;
1485       //
1486       Standard_Real tol2d = myTolApprox;
1487       //
1488       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
1489
1490       aNbParts=myLConstruct.NbParts();
1491       for (i=1; i<=aNbParts; i++) {
1492         myLConstruct.Part(i, fprm, lprm);
1493         Handle(IntPatch_WLine) WL =
1494           convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
1495
1496         theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
1497
1498         if (!theapp3d.IsDone()) {
1499           //
1500           Handle(Geom2d_BSplineCurve) H1;
1501           Handle(Geom2d_BSplineCurve) H2;
1502
1503           if(myApprox1) {
1504             H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1505           }
1506
1507           if(myApprox2) {
1508             H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1509           }
1510           //
1511           mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1512         }
1513
1514         else {
1515           if(myApprox1 || myApprox2) {
1516             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) {
1517               myTolReached2d = theapp3d.TolReached2d();
1518             }
1519           }
1520
1521           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) {
1522             myTolReached3d = theapp3d.TolReached3d();
1523           }
1524
1525           Standard_Integer aNbMultiCurves, nbpoles;
1526           aNbMultiCurves=theapp3d.NbMultiCurves();
1527           for (j=1; j<=aNbMultiCurves; j++) {
1528             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1529             nbpoles = mbspc.NbPoles();
1530
1531             TColgp_Array1OfPnt tpoles(1, nbpoles);
1532             mbspc.Curve(1, tpoles);
1533             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1534                                                                mbspc.Knots(),
1535                                                                mbspc.Multiplicities(),
1536                                                                mbspc.Degree());
1537
1538             GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1539             Check.FixTangent(Standard_True,Standard_True);
1540             //
1541             IntTools_Curve aCurve;
1542             aCurve.SetCurve(BS);
1543
1544             if(myApprox1) {
1545               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1546               mbspc.Curve(2,tpoles2d);
1547               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1548                                                                       mbspc.Knots(),
1549                                                                       mbspc.Multiplicities(),
1550                                                                       mbspc.Degree());
1551
1552               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1553               newCheck.FixTangent(Standard_True,Standard_True);
1554               //
1555               aCurve.SetFirstCurve2d(BS2);
1556             }
1557             else {
1558               Handle(Geom2d_BSplineCurve) H1;
1559               aCurve.SetFirstCurve2d(H1);
1560             }
1561
1562             if(myApprox2) {
1563               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
1564               Standard_Integer TwoOrThree;
1565               TwoOrThree=myApprox1 ? 3 : 2;
1566               mbspc.Curve(TwoOrThree, tpoles2d);
1567               Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
1568                                                                        mbspc.Knots(),
1569                                                                        mbspc.Multiplicities(),
1570                                                                        mbspc.Degree());
1571
1572               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1573               newCheck.FixTangent(Standard_True,Standard_True);
1574               //
1575               aCurve.SetSecondCurve2d(BS2);
1576             }
1577             else {
1578               Handle(Geom2d_BSplineCurve) H2;
1579               aCurve.SetSecondCurve2d(H2);
1580             }
1581             //
1582             mySeqOfCurve.Append(aCurve);
1583
1584           }// for (j=1; j<=aNbMultiCurves; j++) {
1585         }// else from if (!theapp3d.IsDone())
1586       }// for (i=1; i<=aNbParts; i++) {
1587     }// else { // myApprox=TRUE
1588   }// case IntPatch_Analytic:
1589     break;
1590
1591   case IntPatch_Walking:{
1592     Handle(IntPatch_WLine) WL =
1593       Handle(IntPatch_WLine)::DownCast(L);
1594
1595 #ifdef OCCT_DEBUG
1596     //WL->Dump();
1597 #endif
1598
1599     //
1600     Standard_Integer ifprm, ilprm;
1601     //
1602     if (!myApprox) {
1603       aNbParts = 1;
1604       if(!bAvoidLineConstructor){
1605         aNbParts=myLConstruct.NbParts();
1606       }
1607       for (i=1; i<=aNbParts; ++i) {
1608         Handle(Geom2d_BSplineCurve) H1, H2;
1609         Handle(Geom_Curve) aBSp;
1610         //
1611         if(bAvoidLineConstructor) {
1612           ifprm = 1;
1613           ilprm = WL->NbPnts();
1614         }
1615         else {
1616           myLConstruct.Part(i, fprm, lprm);
1617           ifprm=(Standard_Integer)fprm;
1618           ilprm=(Standard_Integer)lprm;
1619         }
1620         //
1621         if(myApprox1) {
1622           H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1623         }
1624         //
1625         if(myApprox2) {
1626           H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1627         }
1628         //
1629         aBSp=MakeBSpline(WL, ifprm, ilprm);
1630         IntTools_Curve aIC(aBSp, H1, H2);
1631         mySeqOfCurve.Append(aIC);
1632       }// for (i=1; i<=aNbParts; ++i) {
1633     }// if (!myApprox) {
1634     //
1635     else { // X
1636       Standard_Boolean bIsDecomposited;
1637       Standard_Integer nbiter, aNbSeqOfL;
1638       Standard_Real tol2d, aTolApproxImp;
1639       IntPatch_SequenceOfLine aSeqOfL;
1640       GeomInt_WLApprox theapp3d;
1641       Approx_ParametrizationType aParType = Approx_ChordLength;
1642       //
1643       Standard_Boolean anApprox1 = myApprox1;
1644       Standard_Boolean anApprox2 = myApprox2;
1645       //
1646       aTolApproxImp=1.e-5;
1647       tol2d = myTolApprox;
1648
1649       GeomAbs_SurfaceType typs1, typs2;
1650       typs1 = myHS1->Surface().GetType();
1651       typs2 = myHS2->Surface().GetType();
1652       Standard_Boolean anWithPC = Standard_True;
1653
1654       if(typs1 == GeomAbs_Cylinder && typs2 == GeomAbs_Sphere) {
1655         anWithPC =
1656           ApproxWithPCurves(myHS1->Surface().Cylinder(), myHS2->Surface().Sphere());
1657       }
1658       else if (typs1 == GeomAbs_Sphere && typs2 == GeomAbs_Cylinder) {
1659         anWithPC =
1660           ApproxWithPCurves(myHS2->Surface().Cylinder(), myHS1->Surface().Sphere());
1661       }
1662       //
1663       if(!anWithPC) {
1664         myTolApprox = aTolApproxImp;//1.e-5;
1665         anApprox1 = Standard_False;
1666         anApprox2 = Standard_False;
1667         //
1668         tol2d = myTolApprox;
1669       }
1670
1671       if(myHS1 == myHS2) {
1672         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1673         rejectSurface = Standard_True;
1674       }
1675       else {
1676         if(reApprox && !rejectSurface)
1677           theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1678         else {
1679           Standard_Integer iDegMax, iDegMin, iNbIter;
1680           //
1681           ApproxParameters(myHS1, myHS2, iDegMin, iDegMax, iNbIter);
1682           theapp3d.SetParameters(myTolApprox, tol2d, iDegMin, iDegMax, iNbIter, Standard_True, aParType);
1683         }
1684       }
1685       //
1686       Standard_Real aReachedTol = Precision::Confusion();
1687       bIsDecomposited=DecompositionOfWLine(WL,
1688                                            myHS1,
1689                                            myHS2,
1690                                            myFace1,
1691                                            myFace2,
1692                                            myLConstruct,
1693                                            bAvoidLineConstructor,
1694                                            aSeqOfL,
1695                                            aReachedTol,
1696                                            myContext);
1697       if ( bIsDecomposited && ( myTolReached3d < aReachedTol ) ) {
1698         myTolReached3d = aReachedTol;
1699       }
1700       //
1701       aNbSeqOfL=aSeqOfL.Length();
1702       //
1703       if (bIsDecomposited) {
1704         nbiter=aNbSeqOfL;
1705       }
1706       else {
1707         nbiter=1;
1708         aNbParts=1;
1709         if (!bAvoidLineConstructor) {
1710           aNbParts=myLConstruct.NbParts();
1711           nbiter=aNbParts;
1712         }
1713       }
1714       //
1715       for(i = 1; i <= nbiter; ++i) {
1716         if(bIsDecomposited) {
1717           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
1718           ifprm = 1;
1719           ilprm = WL->NbPnts();
1720         }
1721         else {
1722           if(bAvoidLineConstructor) {
1723             ifprm = 1;
1724             ilprm = WL->NbPnts();
1725           }
1726           else {
1727             myLConstruct.Part(i, fprm, lprm);
1728             ifprm = (Standard_Integer)fprm;
1729             ilprm = (Standard_Integer)lprm;
1730           }
1731         }
1732         //-- lbr :
1733         //-- Si une des surfaces est un plan , on approxime en 2d
1734         //-- sur cette surface et on remonte les points 2d en 3d.
1735         if(typs1 == GeomAbs_Plane) {
1736           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,Standard_True, myApprox2,ifprm,ilprm);
1737         }
1738         else if(typs2 == GeomAbs_Plane) {
1739           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,myApprox1,Standard_True,ifprm,ilprm);
1740         }
1741         else {
1742           //
1743           if (myHS1 != myHS2){
1744             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
1745                 (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
1746
1747               theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True, aParType);
1748
1749               Standard_Boolean bUseSurfaces;
1750               bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
1751               if (bUseSurfaces) {
1752                 // ######
1753                 rejectSurface = Standard_True;
1754                 // ######
1755                 theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1756               }
1757             }
1758           }
1759           //
1760           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,anApprox1,anApprox2,ifprm,ilprm);
1761         }
1762           //
1763         if (!theapp3d.IsDone()) {
1764           Handle(Geom2d_BSplineCurve) H1;
1765           Handle(Geom2d_BSplineCurve) H2;
1766           //
1767           Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
1768           //
1769           if(myApprox1) {
1770             H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1771           }
1772           //
1773           if(myApprox2) {
1774             H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1775           }
1776           //
1777           IntTools_Curve aIC(aBSp, H1, H2);
1778           mySeqOfCurve.Append(aIC);
1779         }
1780
1781         else {
1782           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) {
1783             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) {
1784               myTolReached2d = theapp3d.TolReached2d();
1785             }
1786           }
1787           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) {
1788             myTolReached3d = myTolReached2d;
1789             //
1790             if (typs1==GeomAbs_Torus || typs2==GeomAbs_Torus) {
1791               if (myTolReached3d<1.e-6) {
1792                 myTolReached3d = theapp3d.TolReached3d();
1793                 myTolReached3d=1.e-6;
1794               }
1795             }
1796           }
1797           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) {
1798             myTolReached3d = theapp3d.TolReached3d();
1799           }
1800
1801           Standard_Integer aNbMultiCurves, nbpoles;
1802           aNbMultiCurves=theapp3d.NbMultiCurves();
1803           for (j=1; j<=aNbMultiCurves; j++) {
1804             if(typs1 == GeomAbs_Plane) {
1805               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1806               nbpoles = mbspc.NbPoles();
1807
1808               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1809               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1810
1811               mbspc.Curve(1,tpoles2d);
1812               const gp_Pln&  Pln = myHS1->Surface().Plane();
1813               //
1814               Standard_Integer ik;
1815               for(ik = 1; ik<= nbpoles; ik++) {
1816                 tpoles.SetValue(ik,
1817                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1818                                               tpoles2d.Value(ik).Y(),
1819                                               Pln));
1820               }
1821               //
1822               Handle(Geom_BSplineCurve) BS =
1823                 new Geom_BSplineCurve(tpoles,
1824                                       mbspc.Knots(),
1825                                       mbspc.Multiplicities(),
1826                                       mbspc.Degree());
1827               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1828               Check.FixTangent(Standard_True, Standard_True);
1829               //
1830               IntTools_Curve aCurve;
1831               aCurve.SetCurve(BS);
1832
1833               if(myApprox1) {
1834                 Handle(Geom2d_BSplineCurve) BS1 =
1835                   new Geom2d_BSplineCurve(tpoles2d,
1836                                           mbspc.Knots(),
1837                                           mbspc.Multiplicities(),
1838                                           mbspc.Degree());
1839                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1840                 Check1.FixTangent(Standard_True,Standard_True);
1841                 //
1842                 // ############################################
1843                 if(!rejectSurface && !reApprox) {
1844                   Standard_Boolean isValid = IsCurveValid(BS1);
1845                   if(!isValid) {
1846                     reApprox = Standard_True;
1847                     goto reapprox;
1848                   }
1849                 }
1850                 // ############################################
1851                 aCurve.SetFirstCurve2d(BS1);
1852               }
1853               else {
1854                 Handle(Geom2d_BSplineCurve) H1;
1855                 aCurve.SetFirstCurve2d(H1);
1856               }
1857
1858               if(myApprox2) {
1859                 mbspc.Curve(2, tpoles2d);
1860
1861                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
1862                                                                           mbspc.Knots(),
1863                                                                           mbspc.Multiplicities(),
1864                                                                           mbspc.Degree());
1865                 GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1866                 newCheck.FixTangent(Standard_True,Standard_True);
1867
1868                 // ###########################################
1869                 if(!rejectSurface && !reApprox) {
1870                   Standard_Boolean isValid = IsCurveValid(BS2);
1871                   if(!isValid) {
1872                     reApprox = Standard_True;
1873                     goto reapprox;
1874                   }
1875                 }
1876                 // ###########################################
1877                 //
1878                 aCurve.SetSecondCurve2d(BS2);
1879               }
1880               else {
1881                 Handle(Geom2d_BSplineCurve) H2;
1882                 //
1883                   aCurve.SetSecondCurve2d(H2);
1884               }
1885               //
1886               mySeqOfCurve.Append(aCurve);
1887
1888             }//if(typs1 == GeomAbs_Plane) {
1889
1890             else if(typs2 == GeomAbs_Plane)
1891             {
1892               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1893               nbpoles = mbspc.NbPoles();
1894
1895               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1896               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1897               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
1898               const gp_Pln&  Pln = myHS2->Surface().Plane();
1899               //
1900               Standard_Integer ik;
1901               for(ik = 1; ik<= nbpoles; ik++) {
1902                 tpoles.SetValue(ik,
1903                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1904                                               tpoles2d.Value(ik).Y(),
1905                                               Pln));
1906
1907               }
1908               //
1909               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1910                                                                  mbspc.Knots(),
1911                                                                  mbspc.Multiplicities(),
1912                                                                  mbspc.Degree());
1913               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1914               Check.FixTangent(Standard_True,Standard_True);
1915               //
1916               IntTools_Curve aCurve;
1917               aCurve.SetCurve(BS);
1918
1919               if(myApprox2) {
1920                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
1921                                                                         mbspc.Knots(),
1922                                                                         mbspc.Multiplicities(),
1923                                                                         mbspc.Degree());
1924                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1925                 Check1.FixTangent(Standard_True,Standard_True);
1926                 //
1927                 // ###########################################
1928                 if(!rejectSurface && !reApprox) {
1929                   Standard_Boolean isValid = IsCurveValid(BS1);
1930                   if(!isValid) {
1931                     reApprox = Standard_True;
1932                     goto reapprox;
1933                   }
1934                 }
1935                 // ###########################################
1936                 bPCurvesOk = CheckPCurve(BS1, myFace2);
1937                 aCurve.SetSecondCurve2d(BS1);
1938               }
1939               else {
1940                 Handle(Geom2d_BSplineCurve) H2;
1941                 aCurve.SetSecondCurve2d(H2);
1942               }
1943
1944               if(myApprox1) {
1945                 mbspc.Curve(1,tpoles2d);
1946                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1947                                                                         mbspc.Knots(),
1948                                                                         mbspc.Multiplicities(),
1949                                                                         mbspc.Degree());
1950                 GeomLib_Check2dBSplineCurve Check2(BS2,TOLCHECK,TOLANGCHECK);
1951                 Check2.FixTangent(Standard_True,Standard_True);
1952                 //
1953                 // ###########################################
1954                 if(!rejectSurface && !reApprox) {
1955                   Standard_Boolean isValid = IsCurveValid(BS2);
1956                   if(!isValid) {
1957                     reApprox = Standard_True;
1958                     goto reapprox;
1959                   }
1960                 }
1961                 // ###########################################
1962                 bPCurvesOk = bPCurvesOk && CheckPCurve(BS2, myFace1);
1963                 aCurve.SetFirstCurve2d(BS2);
1964               }
1965               else {
1966                 Handle(Geom2d_BSplineCurve) H1;
1967                 //
1968                 aCurve.SetFirstCurve2d(H1);
1969               }
1970               //
1971               //if points of the pcurves are out of the faces bounds
1972               //create 3d and 2d curves without approximation
1973               if (!bPCurvesOk) {
1974                 Handle(Geom2d_BSplineCurve) H1, H2;
1975                 bPCurvesOk = Standard_True;
1976                 //
1977                 Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
1978
1979                 if(myApprox1) {
1980                   H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1981                   bPCurvesOk = CheckPCurve(H1, myFace1);
1982                 }
1983
1984                 if(myApprox2) {
1985                   H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1986                   bPCurvesOk = bPCurvesOk && CheckPCurve(H2, myFace2);
1987                 }
1988                 //
1989                 //if pcurves created without approximation are out of the
1990                 //faces bounds, use approximated 3d and 2d curves
1991                 if (bPCurvesOk) {
1992                   IntTools_Curve aIC(aBSp, H1, H2);
1993                   mySeqOfCurve.Append(aIC);
1994                 } else {
1995                   mySeqOfCurve.Append(aCurve);
1996                 }
1997               } else {
1998                 mySeqOfCurve.Append(aCurve);
1999               }
2000
2001             }// else if(typs2 == GeomAbs_Plane)
2002             //
2003             else { //typs2 != GeomAbs_Plane && typs1 != GeomAbs_Plane
2004               Standard_Boolean bIsValid1, bIsValid2;
2005               Handle(Geom_BSplineCurve) BS;
2006               Handle(Geom2d_BSplineCurve) aH2D;
2007               IntTools_Curve aCurve;
2008               //
2009               bIsValid1=Standard_True;
2010               bIsValid2=Standard_True;
2011               //
2012               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
2013               nbpoles = mbspc.NbPoles();
2014               TColgp_Array1OfPnt tpoles(1,nbpoles);
2015               mbspc.Curve(1,tpoles);
2016               BS=new Geom_BSplineCurve(tpoles,
2017                                                                  mbspc.Knots(),
2018                                                                  mbspc.Multiplicities(),
2019                                                                  mbspc.Degree());
2020               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
2021               Check.FixTangent(Standard_True,Standard_True);
2022               //
2023               aCurve.SetCurve(BS);
2024               aCurve.SetFirstCurve2d(aH2D);
2025               aCurve.SetSecondCurve2d(aH2D);
2026               //
2027               if(myApprox1) {
2028                 if(anApprox1) {
2029                   Handle(Geom2d_BSplineCurve) BS1;
2030                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2031                   mbspc.Curve(2,tpoles2d);
2032                   //
2033                   BS1=new Geom2d_BSplineCurve(tpoles2d,
2034                                                                         mbspc.Knots(),
2035                                                                         mbspc.Multiplicities(),
2036                                                                         mbspc.Degree());
2037                   GeomLib_Check2dBSplineCurve newCheck(BS1,TOLCHECK,TOLANGCHECK);
2038                   newCheck.FixTangent(Standard_True,Standard_True);
2039                   //
2040                   if (!reApprox) {
2041                     bIsValid1=CheckPCurve(BS1, myFace1);
2042                   }
2043                   //
2044                   aCurve.SetFirstCurve2d(BS1);
2045                 }
2046                 else {
2047                   Handle(Geom2d_BSplineCurve) BS1;
2048                   fprm = BS->FirstParameter();
2049                   lprm = BS->LastParameter();
2050
2051                   Handle(Geom2d_Curve) C2d;
2052                   Standard_Real aTol = myTolApprox;
2053                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
2054                             myHS1->ChangeSurface().Surface(), BS, C2d);
2055                   BS1 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2056                   aCurve.SetFirstCurve2d(BS1);
2057                 }
2058               } // if(myApprox1) {
2059                 //
2060               if(myApprox2) {
2061                 if(anApprox2) {
2062                   Handle(Geom2d_BSplineCurve) BS2;
2063                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2064                   mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
2065                   BS2=new Geom2d_BSplineCurve(tpoles2d,
2066                                                                         mbspc.Knots(),
2067                                                                         mbspc.Multiplicities(),
2068                                                                         mbspc.Degree());
2069                   GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
2070                   newCheck.FixTangent(Standard_True,Standard_True);
2071                 //
2072                   if (!reApprox) {
2073                     bIsValid2=CheckPCurve(BS2, myFace2);
2074                   }
2075                   aCurve.SetSecondCurve2d(BS2);
2076                 }
2077                 else {
2078                   Handle(Geom2d_BSplineCurve) BS2;
2079                   fprm = BS->FirstParameter();
2080                   lprm = BS->LastParameter();
2081
2082                   Handle(Geom2d_Curve) C2d;
2083                   Standard_Real aTol = myTolApprox;
2084                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
2085                             myHS2->ChangeSurface().Surface(), BS, C2d);
2086                   BS2 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2087                   aCurve.SetSecondCurve2d(BS2);
2088                 }
2089               } //if(myApprox2) {
2090               if (!bIsValid1 || !bIsValid2) {
2091                 myTolApprox=aTolApproxImp;//1.e-5;
2092                 tol2d = myTolApprox;
2093                 reApprox = Standard_True;
2094                 goto reapprox;
2095               }
2096                 //
2097               mySeqOfCurve.Append(aCurve);
2098             }
2099           }
2100         }
2101       }
2102     }// else { // X
2103   }// case IntPatch_Walking:{
2104     break;
2105
2106   case IntPatch_Restriction:
2107     {
2108       Handle(IntPatch_RLine) RL =
2109         Handle(IntPatch_RLine)::DownCast(L);
2110       Handle(Geom_Curve) aC3d;
2111       Handle(Geom2d_Curve) aC2d1, aC2d2;
2112       Standard_Real aTolReached;
2113       GeomInt_IntSS::TreatRLine(RL, myHS1, myHS2, aC3d,
2114                                   aC2d1, aC2d2, aTolReached);
2115
2116       if(aC3d.IsNull())
2117         break;
2118
2119       Bnd_Box2d aBox1, aBox2;
2120
2121       const Standard_Real aU1f = myHS1->FirstUParameter(),
2122                           aV1f = myHS1->FirstVParameter(),
2123                           aU1l = myHS1->LastUParameter(),
2124                           aV1l = myHS1->LastVParameter();
2125       const Standard_Real aU2f = myHS2->FirstUParameter(),
2126                           aV2f = myHS2->FirstVParameter(),
2127                           aU2l = myHS2->LastUParameter(),
2128                           aV2l = myHS2->LastVParameter();
2129
2130       aBox1.Add(gp_Pnt2d(aU1f, aV1f));
2131       aBox1.Add(gp_Pnt2d(aU1l, aV1l));
2132       aBox2.Add(gp_Pnt2d(aU2f, aV2f));
2133       aBox2.Add(gp_Pnt2d(aU2l, aV2l));
2134
2135       GeomInt_VectorOfReal anArrayOfParameters;
2136
2137       //We consider here that the intersection line is same-parameter-line
2138       anArrayOfParameters.Append(aC3d->FirstParameter());
2139       anArrayOfParameters.Append(aC3d->LastParameter());
2140
2141       GeomInt_IntSS::
2142         TrimILineOnSurfBoundaries(aC2d1, aC2d2, aBox1, aBox2, anArrayOfParameters);
2143
2144       const Standard_Integer aNbIntersSolutionsm1 = anArrayOfParameters.Length() - 1;
2145
2146       //Trim RLine found.
2147       for(Standard_Integer anInd = 0; anInd < aNbIntersSolutionsm1; anInd++)
2148       {
2149         const Standard_Real aParF = anArrayOfParameters(anInd),
2150                             aParL = anArrayOfParameters(anInd+1);
2151
2152         if((aParL - aParF) <= Precision::PConfusion())
2153           continue;
2154
2155         const Standard_Real aPar = 0.5*(aParF + aParL);
2156         gp_Pnt2d aPt;
2157
2158         Handle(Geom2d_Curve) aCurv2d1, aCurv2d2;
2159         if(!aC2d1.IsNull())
2160         {
2161           aC2d1->D0(aPar, aPt);
2162
2163           if(aBox1.IsOut(aPt))
2164             continue;
2165
2166           if(myApprox1)
2167             aCurv2d1 = new Geom2d_TrimmedCurve(aC2d1, aParF, aParL);
2168         }
2169
2170         if(!aC2d2.IsNull())
2171         {
2172           aC2d2->D0(aPar, aPt);
2173
2174           if(aBox2.IsOut(aPt))
2175             continue;
2176
2177           if(myApprox2)
2178             aCurv2d2 = new Geom2d_TrimmedCurve(aC2d2, aParF, aParL);
2179         }
2180
2181         Handle(Geom_Curve) aCurv3d = new Geom_TrimmedCurve(aC3d, aParF, aParL);
2182
2183         IntTools_Curve aIC(aCurv3d, aCurv2d1, aCurv2d2);
2184         mySeqOfCurve.Append(aIC);
2185       }
2186     }
2187     break;
2188   default:
2189     break;
2190
2191   }
2192 }
2193
2194 //=======================================================================
2195 //function : Parameters
2196 //purpose  :
2197 //=======================================================================
2198  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
2199                  const Handle(GeomAdaptor_HSurface)& HS2,
2200                  const gp_Pnt& Ptref,
2201                  Standard_Real& U1,
2202                  Standard_Real& V1,
2203                  Standard_Real& U2,
2204                  Standard_Real& V2)
2205 {
2206
2207   IntSurf_Quadric quad1,quad2;
2208   GeomAbs_SurfaceType typs = HS1->Surface().GetType();
2209
2210   switch (typs) {
2211   case GeomAbs_Plane:
2212     quad1.SetValue(HS1->Surface().Plane());
2213     break;
2214   case GeomAbs_Cylinder:
2215     quad1.SetValue(HS1->Surface().Cylinder());
2216     break;
2217   case GeomAbs_Cone:
2218     quad1.SetValue(HS1->Surface().Cone());
2219     break;
2220   case GeomAbs_Sphere:
2221     quad1.SetValue(HS1->Surface().Sphere());
2222     break;
2223   case GeomAbs_Torus:
2224     quad1.SetValue(HS1->Surface().Torus());
2225     break;
2226   default:
2227     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2228   }
2229
2230   typs = HS2->Surface().GetType();
2231   switch (typs) {
2232   case GeomAbs_Plane:
2233     quad2.SetValue(HS2->Surface().Plane());
2234     break;
2235   case GeomAbs_Cylinder:
2236     quad2.SetValue(HS2->Surface().Cylinder());
2237     break;
2238   case GeomAbs_Cone:
2239     quad2.SetValue(HS2->Surface().Cone());
2240     break;
2241   case GeomAbs_Sphere:
2242     quad2.SetValue(HS2->Surface().Sphere());
2243     break;
2244   case GeomAbs_Torus:
2245     quad2.SetValue(HS2->Surface().Torus());
2246     break;
2247   default:
2248     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2249   }
2250
2251   quad1.Parameters(Ptref,U1,V1);
2252   quad2.Parameters(Ptref,U2,V2);
2253 }
2254
2255 //=======================================================================
2256 //function : MakeBSpline
2257 //purpose  :
2258 //=======================================================================
2259 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
2260                                  const Standard_Integer ideb,
2261                                  const Standard_Integer ifin)
2262 {
2263   Standard_Integer i,nbpnt = ifin-ideb+1;
2264   TColgp_Array1OfPnt poles(1,nbpnt);
2265   TColStd_Array1OfReal knots(1,nbpnt);
2266   TColStd_Array1OfInteger mults(1,nbpnt);
2267   Standard_Integer ipidebm1;
2268   for(i=1,ipidebm1=i+ideb-1; i<=nbpnt;ipidebm1++, i++) {
2269     poles(i) = WL->Point(ipidebm1).Value();
2270     mults(i) = 1;
2271     knots(i) = i-1;
2272   }
2273   mults(1) = mults(nbpnt) = 2;
2274   return
2275     new Geom_BSplineCurve(poles,knots,mults,1);
2276 }
2277 //
2278
2279 //=======================================================================
2280 //function : MakeBSpline2d
2281 //purpose  :
2282 //=======================================================================
2283 Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
2284                                           const Standard_Integer ideb,
2285                                           const Standard_Integer ifin,
2286                                           const Standard_Boolean onFirst)
2287 {
2288   Standard_Integer i, nbpnt = ifin-ideb+1;
2289   TColgp_Array1OfPnt2d poles(1,nbpnt);
2290   TColStd_Array1OfReal knots(1,nbpnt);
2291   TColStd_Array1OfInteger mults(1,nbpnt);
2292   Standard_Integer ipidebm1;
2293
2294   for(i = 1, ipidebm1 = i+ideb-1; i <= nbpnt; ipidebm1++, i++) {
2295       Standard_Real U, V;
2296       if(onFirst)
2297         theWLine->Point(ipidebm1).ParametersOnS1(U, V);
2298       else
2299         theWLine->Point(ipidebm1).ParametersOnS2(U, V);
2300       poles(i).SetCoord(U, V);
2301       mults(i) = 1;
2302       knots(i) = i-1;
2303     }
2304     mults(1) = mults(nbpnt) = 2;
2305
2306   return new Geom2d_BSplineCurve(poles,knots,mults,1);
2307 }
2308 //=======================================================================
2309 //function : PrepareLines3D
2310 //purpose  :
2311 //=======================================================================
2312   void IntTools_FaceFace::PrepareLines3D(const Standard_Boolean bToSplit)
2313 {
2314   Standard_Integer i, aNbCurves;
2315   GeomAbs_SurfaceType aType1, aType2;
2316   IntTools_SequenceOfCurves aNewCvs;
2317   //
2318   // 1. Treatment closed  curves
2319   aNbCurves=mySeqOfCurve.Length();
2320   for (i=1; i<=aNbCurves; ++i) {
2321     const IntTools_Curve& aIC=mySeqOfCurve(i);
2322     //
2323     if (bToSplit) {
2324       Standard_Integer j, aNbC;
2325       IntTools_SequenceOfCurves aSeqCvs;
2326       //
2327       aNbC=IntTools_Tools::SplitCurve(aIC, aSeqCvs);
2328       if (aNbC) {
2329         for (j=1; j<=aNbC; ++j) {
2330           const IntTools_Curve& aICNew=aSeqCvs(j);
2331           aNewCvs.Append(aICNew);
2332         }
2333       }
2334       else {
2335         aNewCvs.Append(aIC);
2336       }
2337     }
2338     else {
2339       aNewCvs.Append(aIC);
2340     }
2341   }
2342   //
2343   // 2. Plane\Cone intersection when we had 4 curves
2344   aType1=myHS1->GetType();
2345   aType2=myHS2->GetType();
2346   aNbCurves=aNewCvs.Length();
2347   //
2348   if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone) ||
2349       (aType2==GeomAbs_Plane && aType1==GeomAbs_Cone)) {
2350     if (aNbCurves==4) {
2351       GeomAbs_CurveType aCType1;
2352       //
2353       aCType1=aNewCvs(1).Type();
2354       if (aCType1==GeomAbs_Line) {
2355         IntTools_SequenceOfCurves aSeqIn, aSeqOut;
2356         //
2357         for (i=1; i<=aNbCurves; ++i) {
2358           const IntTools_Curve& aIC=aNewCvs(i);
2359           aSeqIn.Append(aIC);
2360         }
2361         //
2362         IntTools_Tools::RejectLines(aSeqIn, aSeqOut);
2363         //
2364         aNewCvs.Clear();
2365         aNbCurves=aSeqOut.Length();
2366         for (i=1; i<=aNbCurves; ++i) {
2367           const IntTools_Curve& aIC=aSeqOut(i);
2368           aNewCvs.Append(aIC);
2369         }
2370       }
2371     }
2372   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone)...
2373   //
2374   // 3. Fill  mySeqOfCurve
2375   mySeqOfCurve.Clear();
2376   aNbCurves=aNewCvs.Length();
2377   for (i=1; i<=aNbCurves; ++i) {
2378     const IntTools_Curve& aIC=aNewCvs(i);
2379     mySeqOfCurve.Append(aIC);
2380   }
2381 }
2382 //=======================================================================
2383 //function : CorrectSurfaceBoundaries
2384 //purpose  :
2385 //=======================================================================
2386  void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
2387                               const Standard_Real theTolerance,
2388                               Standard_Real&      theumin,
2389                               Standard_Real&      theumax,
2390                               Standard_Real&      thevmin,
2391                               Standard_Real&      thevmax)
2392 {
2393   Standard_Boolean enlarge, isuperiodic, isvperiodic;
2394   Standard_Real uinf, usup, vinf, vsup, delta;
2395   GeomAbs_SurfaceType aType;
2396   Handle(Geom_Surface) aSurface;
2397   //
2398   aSurface = BRep_Tool::Surface(theFace);
2399   aSurface->Bounds(uinf, usup, vinf, vsup);
2400   delta = theTolerance;
2401   enlarge = Standard_False;
2402   //
2403   GeomAdaptor_Surface anAdaptorSurface(aSurface);
2404   //
2405   if(aSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
2406     Handle(Geom_Surface) aBasisSurface =
2407       (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface))->BasisSurface();
2408
2409     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2410        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2411       return;
2412     }
2413   }
2414   //
2415   if(aSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2416     Handle(Geom_Surface) aBasisSurface =
2417       (Handle(Geom_OffsetSurface)::DownCast(aSurface))->BasisSurface();
2418
2419     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2420        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2421       return;
2422     }
2423   }
2424   //
2425   isuperiodic = anAdaptorSurface.IsUPeriodic();
2426   isvperiodic = anAdaptorSurface.IsVPeriodic();
2427   //
2428   aType=anAdaptorSurface.GetType();
2429   if((aType==GeomAbs_BezierSurface) ||
2430      (aType==GeomAbs_BSplineSurface) ||
2431      (aType==GeomAbs_SurfaceOfExtrusion) ||
2432      (aType==GeomAbs_SurfaceOfRevolution) ||
2433      (aType==GeomAbs_Cylinder)) {
2434     enlarge=Standard_True;
2435   }
2436   //
2437   if(!isuperiodic && enlarge) {
2438
2439     if(!Precision::IsInfinite(theumin) &&
2440         ((theumin - uinf) > delta))
2441       theumin -= delta;
2442     else {
2443       theumin = uinf;
2444     }
2445
2446     if(!Precision::IsInfinite(theumax) &&
2447         ((usup - theumax) > delta))
2448       theumax += delta;
2449     else
2450       theumax = usup;
2451   }
2452   //
2453   if(!isvperiodic && enlarge) {
2454     if(!Precision::IsInfinite(thevmin) &&
2455         ((thevmin - vinf) > delta)) {
2456       thevmin -= delta;
2457     }
2458     else {
2459       thevmin = vinf;
2460     }
2461     if(!Precision::IsInfinite(thevmax) &&
2462         ((vsup - thevmax) > delta)) {
2463       thevmax += delta;
2464     }
2465     else {
2466       thevmax = vsup;
2467     }
2468   }
2469   //
2470   {
2471     Standard_Integer aNbP;
2472     Standard_Real aXP, dXfact, aXmid, aX1, aX2, aTolPA;
2473     //
2474     aTolPA=Precision::Angular();
2475     // U
2476     if (isuperiodic) {
2477       aXP=anAdaptorSurface.UPeriod();
2478       dXfact=theumax-theumin;
2479       if (dXfact-aTolPA>aXP) {
2480         aXmid=0.5*(theumax+theumin);
2481         aNbP=RealToInt(aXmid/aXP);
2482         if (aXmid<0.) {
2483           aNbP=aNbP-1;
2484         }
2485         aX1=aNbP*aXP;
2486         if (theumin>aTolPA) {
2487           aX1=theumin+aNbP*aXP;
2488         }
2489         aX2=aX1+aXP;
2490         if (theumin<aX1) {
2491           theumin=aX1;
2492         }
2493         if (theumax>aX2) {
2494           theumax=aX2;
2495         }
2496       }
2497     }
2498     // V
2499     if (isvperiodic) {
2500       aXP=anAdaptorSurface.VPeriod();
2501       dXfact=thevmax-thevmin;
2502       if (dXfact-aTolPA>aXP) {
2503         aXmid=0.5*(thevmax+thevmin);
2504         aNbP=RealToInt(aXmid/aXP);
2505         if (aXmid<0.) {
2506           aNbP=aNbP-1;
2507         }
2508         aX1=aNbP*aXP;
2509         if (thevmin>aTolPA) {
2510           aX1=thevmin+aNbP*aXP;
2511         }
2512         aX2=aX1+aXP;
2513         if (thevmin<aX1) {
2514           thevmin=aX1;
2515         }
2516         if (thevmax>aX2) {
2517           thevmax=aX2;
2518         }
2519       }
2520     }
2521   }
2522   //
2523   if(isuperiodic || isvperiodic) {
2524     Standard_Boolean correct = Standard_False;
2525     Standard_Boolean correctU = Standard_False;
2526     Standard_Boolean correctV = Standard_False;
2527     Bnd_Box2d aBox;
2528     TopExp_Explorer anExp;
2529
2530     for(anExp.Init(theFace, TopAbs_EDGE); anExp.More(); anExp.Next()) {
2531       if(BRep_Tool::IsClosed(TopoDS::Edge(anExp.Current()), theFace)) {
2532         correct = Standard_True;
2533         Standard_Real f, l;
2534         TopoDS_Edge anEdge = TopoDS::Edge(anExp.Current());
2535
2536         for(Standard_Integer i = 0; i < 2; i++) {
2537           if(i==0) {
2538             anEdge.Orientation(TopAbs_FORWARD);
2539           }
2540           else {
2541             anEdge.Orientation(TopAbs_REVERSED);
2542           }
2543           Handle(Geom2d_Curve) aCurve = BRep_Tool::CurveOnSurface(anEdge, theFace, f, l);
2544
2545           if(aCurve.IsNull()) {
2546             correct = Standard_False;
2547             break;
2548           }
2549           Handle(Geom2d_Line) aLine = Handle(Geom2d_Line)::DownCast(aCurve);
2550
2551           if(aLine.IsNull()) {
2552             correct = Standard_False;
2553             break;
2554           }
2555           gp_Dir2d anUDir(1., 0.);
2556           gp_Dir2d aVDir(0., 1.);
2557           Standard_Real anAngularTolerance = Precision::Angular();
2558
2559           correctU = correctU || aLine->Position().Direction().IsParallel(aVDir, anAngularTolerance);
2560           correctV = correctV || aLine->Position().Direction().IsParallel(anUDir, anAngularTolerance);
2561
2562           gp_Pnt2d pp1 = aCurve->Value(f);
2563           aBox.Add(pp1);
2564           gp_Pnt2d pp2 = aCurve->Value(l);
2565           aBox.Add(pp2);
2566         }
2567         if(!correct)
2568           break;
2569       }
2570     }
2571
2572     if(correct) {
2573       Standard_Real umin, vmin, umax, vmax;
2574       aBox.Get(umin, vmin, umax, vmax);
2575
2576       if(isuperiodic && correctU) {
2577
2578         if(theumin < umin)
2579           theumin = umin;
2580
2581         if(theumax > umax) {
2582           theumax = umax;
2583         }
2584       }
2585       if(isvperiodic && correctV) {
2586
2587         if(thevmin < vmin)
2588           thevmin = vmin;
2589         if(thevmax > vmax)
2590           thevmax = vmax;
2591       }
2592     }
2593   }
2594 }
2595 //
2596 //
2597 // The block is dedicated to determine whether WLine [ifprm, ilprm]
2598 // crosses the degenerated zone on each given surface or not.
2599 // If Yes -> We will not use info about surfaces during approximation
2600 // because inside degenerated zone of the surface the approx. algo.
2601 // uses wrong values of normal, etc., and resulting curve will have
2602 // oscillations that we would not like to have.
2603
2604
2605
2606 static
2607   Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2608                                      const Handle(Geom_Surface)& aS,
2609                                      const Standard_Integer iDir);
2610 static
2611   Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2612                                             const TopoDS_Face& aF1,
2613                                             const TopoDS_Face& aF2);
2614 //=======================================================================
2615 //function :  NotUseSurfacesForApprox
2616 //purpose  :
2617 //=======================================================================
2618 Standard_Boolean NotUseSurfacesForApprox(const TopoDS_Face& aF1,
2619                                          const TopoDS_Face& aF2,
2620                                          const Handle(IntPatch_WLine)& WL,
2621                                          const Standard_Integer ifprm,
2622                                          const Standard_Integer ilprm)
2623 {
2624   Standard_Boolean bPInDZ;
2625
2626   Handle(IntSurf_LineOn2S) aLineOn2S=WL->Curve();
2627
2628   const IntSurf_PntOn2S& aP2Sfprm=aLineOn2S->Value(ifprm);
2629   bPInDZ=IsPointInDegeneratedZone(aP2Sfprm, aF1, aF2);
2630   if (bPInDZ) {
2631     return bPInDZ;
2632   }
2633
2634   const IntSurf_PntOn2S& aP2Slprm=aLineOn2S->Value(ilprm);
2635   bPInDZ=IsPointInDegeneratedZone(aP2Slprm, aF1, aF2);
2636
2637   return bPInDZ;
2638 }
2639 //=======================================================================
2640 //function : IsPointInDegeneratedZone
2641 //purpose  :
2642 //=======================================================================
2643 Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2644                                           const TopoDS_Face& aF1,
2645                                           const TopoDS_Face& aF2)
2646
2647 {
2648   Standard_Boolean bFlag=Standard_True;
2649   Standard_Real US11, US12, VS11, VS12, US21, US22, VS21, VS22;
2650   Standard_Real U1, V1, U2, V2, aDelta, aD;
2651   gp_Pnt2d aP2d;
2652
2653   Handle(Geom_Surface)aS1 = BRep_Tool::Surface(aF1);
2654   aS1->Bounds(US11, US12, VS11, VS12);
2655   GeomAdaptor_Surface aGAS1(aS1);
2656
2657   Handle(Geom_Surface)aS2 = BRep_Tool::Surface(aF2);
2658   aS1->Bounds(US21, US22, VS21, VS22);
2659   GeomAdaptor_Surface aGAS2(aS2);
2660   //
2661   //const gp_Pnt& aP=aP2S.Value();
2662   aP2S.Parameters(U1, V1, U2, V2);
2663   //
2664   aDelta=1.e-7;
2665   // Check on Surf 1
2666   aD=aGAS1.UResolution(aDelta);
2667   aP2d.SetCoord(U1, V1);
2668   if (fabs(U1-US11) < aD) {
2669     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2670     if (bFlag) {
2671       return bFlag;
2672     }
2673   }
2674   if (fabs(U1-US12) < aD) {
2675     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2676     if (bFlag) {
2677       return bFlag;
2678     }
2679   }
2680   aD=aGAS1.VResolution(aDelta);
2681   if (fabs(V1-VS11) < aDelta) {
2682     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2683     if (bFlag) {
2684       return bFlag;
2685     }
2686   }
2687   if (fabs(V1-VS12) < aDelta) {
2688     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2689     if (bFlag) {
2690       return bFlag;
2691     }
2692   }
2693   // Check on Surf 2
2694   aD=aGAS2.UResolution(aDelta);
2695   aP2d.SetCoord(U2, V2);
2696   if (fabs(U2-US21) < aDelta) {
2697     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2698     if (bFlag) {
2699       return bFlag;
2700     }
2701   }
2702   if (fabs(U2-US22) < aDelta) {
2703     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2704     if (bFlag) {
2705       return bFlag;
2706     }
2707   }
2708   aD=aGAS2.VResolution(aDelta);
2709   if (fabs(V2-VS21) < aDelta) {
2710     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2711     if (bFlag) {
2712       return bFlag;
2713     }
2714   }
2715   if (fabs(V2-VS22) < aDelta) {
2716     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2717     if (bFlag) {
2718       return bFlag;
2719     }
2720   }
2721   return !bFlag;
2722 }
2723
2724 //=======================================================================
2725 //function : IsDegeneratedZone
2726 //purpose  :
2727 //=======================================================================
2728 Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2729                                    const Handle(Geom_Surface)& aS,
2730                                    const Standard_Integer iDir)
2731 {
2732   Standard_Boolean bFlag=Standard_True;
2733   Standard_Real US1, US2, VS1, VS2, dY, dX, d1, d2, dD;
2734   Standard_Real aXm, aYm, aXb, aYb, aXe, aYe;
2735   aS->Bounds(US1, US2, VS1, VS2);
2736
2737   gp_Pnt aPm, aPb, aPe;
2738
2739   aXm=aP2d.X();
2740   aYm=aP2d.Y();
2741
2742   aS->D0(aXm, aYm, aPm);
2743
2744   dX=1.e-5;
2745   dY=1.e-5;
2746   dD=1.e-12;
2747
2748   if (iDir==1) {
2749     aXb=aXm;
2750     aXe=aXm;
2751     aYb=aYm-dY;
2752     if (aYb < VS1) {
2753       aYb=VS1;
2754     }
2755     aYe=aYm+dY;
2756     if (aYe > VS2) {
2757       aYe=VS2;
2758     }
2759     aS->D0(aXb, aYb, aPb);
2760     aS->D0(aXe, aYe, aPe);
2761
2762     d1=aPm.Distance(aPb);
2763     d2=aPm.Distance(aPe);
2764     if (d1 < dD && d2 < dD) {
2765       return bFlag;
2766     }
2767     return !bFlag;
2768   }
2769   //
2770   else if (iDir==2) {
2771     aYb=aYm;
2772     aYe=aYm;
2773     aXb=aXm-dX;
2774     if (aXb < US1) {
2775       aXb=US1;
2776     }
2777     aXe=aXm+dX;
2778     if (aXe > US2) {
2779       aXe=US2;
2780     }
2781     aS->D0(aXb, aYb, aPb);
2782     aS->D0(aXe, aYe, aPe);
2783
2784     d1=aPm.Distance(aPb);
2785     d2=aPm.Distance(aPe);
2786     if (d1 < dD && d2 < dD) {
2787       return bFlag;
2788     }
2789     return !bFlag;
2790   }
2791   return !bFlag;
2792 }
2793
2794 // Check if aNextPnt lies inside of tube build on aBasePnt and aBaseVec.
2795 // In 2d space.
2796 static Standard_Boolean IsInsideIn2d(const gp_Pnt2d& aBasePnt,
2797                                      const gp_Vec2d& aBaseVec,
2798                                      const gp_Pnt2d& aNextPnt,
2799                                      const Standard_Real aSquareMaxDist)
2800 {
2801   gp_Vec2d aVec2d(aBasePnt, aNextPnt);
2802
2803   //d*d = (basevec^(nextpnt-basepnt))**2 / basevec**2
2804   Standard_Real aCross = aVec2d.Crossed(aBaseVec);
2805   Standard_Real aSquareDist = aCross * aCross
2806                             / aBaseVec.SquareMagnitude();
2807
2808   return (aSquareDist <= aSquareMaxDist);
2809 }
2810
2811 // Check if aNextPnt lies inside of tube build on aBasePnt and aBaseVec.
2812 // In 3d space.
2813 static Standard_Boolean IsInsideIn3d(const gp_Pnt& aBasePnt,
2814                                      const gp_Vec& aBaseVec,
2815                                      const gp_Pnt& aNextPnt,
2816                                      const Standard_Real aSquareMaxDist)
2817 {
2818   gp_Vec aVec(aBasePnt, aNextPnt);
2819
2820   //d*d = (basevec^(nextpnt-basepnt))**2 / basevec**2
2821   Standard_Real aSquareDist = aVec.CrossSquareMagnitude(aBaseVec)
2822                             / aBaseVec.SquareMagnitude();
2823
2824   return (aSquareDist <= aSquareMaxDist);
2825 }
2826
2827 //=========================================================================
2828 // static function : ComputePurgedWLine
2829 // purpose : Removes equal points (leave one of equal points) from theWLine
2830 //           and recompute vertex parameters.
2831 //           Removes exceed points using tube criteria:
2832 //           delete 7D point if it lies near to expected lines in 2d and 3d.
2833 //           Each task (2d, 2d, 3d) have its own tolerance and checked separately.
2834 //           Returns new WLine or null WLine if the number
2835 //           of the points is less than 2.
2836 //=========================================================================
2837 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)       &theWLine,
2838                                           const Handle(GeomAdaptor_HSurface) &theS1,
2839                                           const Handle(GeomAdaptor_HSurface) &theS2)
2840 {
2841   Standard_Integer i, k, v, nb, nbvtx;
2842   Handle(IntPatch_WLine) aResult;
2843   nbvtx = theWLine->NbVertex();
2844   nb = theWLine->NbPnts();
2845   if (nb==2) {
2846     const IntSurf_PntOn2S& p1 = theWLine->Point(1);
2847     const IntSurf_PntOn2S& p2 = theWLine->Point(2);
2848     if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
2849       return aResult;
2850     }
2851   }
2852   //
2853   Handle(IntPatch_WLine) aLocalWLine;
2854   Handle(IntPatch_WLine) aTmpWLine = theWLine;
2855   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
2856   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
2857   for(i = 1; i <= nb; i++) {
2858     aLineOn2S->Add(theWLine->Point(i));
2859   }
2860
2861   for(v = 1; v <= nbvtx; v++) {
2862     aLocalWLine->AddVertex(theWLine->Vertex(v));
2863   }
2864   for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
2865     Standard_Integer aStartIndex = i + 1;
2866     Standard_Integer anEndIndex = i + 5;
2867     nb = aLineOn2S->NbPoints();
2868     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
2869
2870     if((aStartIndex > nb) || (anEndIndex <= 1)) {
2871       continue;
2872     }
2873     k = aStartIndex;
2874
2875     while(k <= anEndIndex) {
2876
2877       if(i != k) {
2878         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
2879         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
2880
2881         Standard_Real UV[8];
2882         p1.Parameters(UV[0], UV[1], UV[2], UV[3]);
2883         p2.Parameters(UV[4], UV[5], UV[6], UV[7]);
2884
2885         Standard_Real aMax = Abs(UV[0]);
2886         for(Standard_Integer anIdx = 1; anIdx < 8; anIdx++)
2887         {
2888           if (aMax < Abs(UV[anIdx]))
2889             aMax = Abs(UV[anIdx]);
2890         }
2891
2892         if(p1.Value().IsEqual(p2.Value(), gp::Resolution()) ||
2893            Abs(UV[0] - UV[4]) + Abs(UV[1] - UV[5]) < 1.0e-16 * aMax ||
2894            Abs(UV[2] - UV[6]) + Abs(UV[3] - UV[7]) < 1.0e-16 * aMax )
2895         {
2896           aTmpWLine = aLocalWLine;
2897           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
2898
2899           for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
2900             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
2901             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
2902
2903             if(avertexindex >= k) {
2904               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
2905             }
2906             aLocalWLine->AddVertex(aVertex);
2907           }
2908           aLineOn2S->RemovePoint(k);
2909           anEndIndex--;
2910           continue;
2911         }
2912       }
2913       k++;
2914     }
2915   }
2916
2917   if (aLineOn2S->NbPoints() <= 2)
2918   {
2919     if (aLineOn2S->NbPoints() == 2)
2920       return aLocalWLine;
2921     else
2922       return aResult;
2923   }
2924
2925   const Standard_Integer aMinNbBadDistr = 15;
2926   const Standard_Integer aNbSingleBezier = 30;
2927   // Avoid purge in case of C0 continuity:
2928   // Intersection approximator may produce invalid curve after purge, example:
2929   // bugs modalg_5 bug24731,
2930   // Do not run purger when base number of points is too small.
2931   if (theS1->UContinuity() == GeomAbs_C0 ||
2932       theS1->VContinuity() == GeomAbs_C0 ||
2933       theS2->UContinuity() == GeomAbs_C0 ||
2934       theS2->VContinuity() == GeomAbs_C0 ||
2935       nb < aNbSingleBezier)
2936   {
2937     return aLocalWLine;
2938   }
2939
2940   // 1 - Delete point.
2941   // 0 - Store point.
2942   // -1 - Vertex point (not delete).
2943   NCollection_Array1<Standard_Integer> aNewPointsHash(1, aLineOn2S->NbPoints());
2944   for(i = 1; i <= aLineOn2S->NbPoints(); i++)
2945     aNewPointsHash.SetValue(i, 0);
2946
2947   for(v = 1; v <= aLocalWLine->NbVertex(); v++)
2948   {
2949     IntPatch_Point aVertex = aLocalWLine->Vertex(v);
2950     Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
2951     aNewPointsHash.SetValue(avertexindex, -1);
2952   }
2953
2954   // Workaround to handle case of small amount points after purge.
2955   // Test "boolean boptuc_complex B5" and similar.
2956   Standard_Integer aNbPnt = 0;
2957
2958   // Inital computations.
2959   Standard_Real UonS1[3], VonS1[3], UonS2[3], VonS2[3];
2960   aLineOn2S->Value(1).ParametersOnS1(UonS1[0], VonS1[0]);
2961   aLineOn2S->Value(2).ParametersOnS1(UonS1[1], VonS1[1]);
2962   aLineOn2S->Value(1).ParametersOnS2(UonS2[0], VonS2[0]);
2963   aLineOn2S->Value(2).ParametersOnS2(UonS2[1], VonS2[1]);
2964
2965   gp_Pnt2d aBase2dPnt1(UonS1[0], VonS1[0]);
2966   gp_Pnt2d aBase2dPnt2(UonS2[0], VonS2[0]);
2967   gp_Vec2d aBase2dVec1(UonS1[1] - UonS1[0], VonS1[1] - VonS1[0]);
2968   gp_Vec2d aBase2dVec2(UonS2[1] - UonS2[0], VonS2[1] - VonS2[0]);
2969   gp_Pnt   aBase3dPnt = aLineOn2S->Value(1).Value();
2970   gp_Vec   aBase3dVec(aLineOn2S->Value(1).Value(), aLineOn2S->Value(2).Value());
2971
2972   // Choose base tolerance and scale it to pipe algorithm.
2973   const Standard_Real aBaseTolerance = Precision::Approximation();
2974   Standard_Real aResS1Tol = Min(theS1->UResolution(aBaseTolerance),
2975                                 theS1->VResolution(aBaseTolerance));
2976   Standard_Real aResS2Tol = Min(theS2->UResolution(aBaseTolerance),
2977                                 theS2->VResolution(aBaseTolerance));
2978   Standard_Real aTol1 = aResS1Tol * aResS1Tol;
2979   Standard_Real aTol2 = aResS2Tol * aResS2Tol;
2980   Standard_Real aTol3d = aBaseTolerance * aBaseTolerance;
2981
2982   const Standard_Real aLimitCoeff = 0.99 * 0.99;
2983   for(i = 3; i <= aLineOn2S->NbPoints(); i++)
2984   {
2985     Standard_Boolean isDeleteState = Standard_False;
2986
2987     aLineOn2S->Value(i).ParametersOnS1(UonS1[2], VonS1[2]);
2988     aLineOn2S->Value(i).ParametersOnS2(UonS2[2], VonS2[2]);
2989     gp_Pnt2d aPnt2dOnS1(UonS1[2], VonS1[2]);
2990     gp_Pnt2d aPnt2dOnS2(UonS2[2], VonS2[2]);
2991     const gp_Pnt& aPnt3d = aLineOn2S->Value(i).Value();
2992
2993     if (aNewPointsHash(i - 1) != - 1 &&
2994         IsInsideIn2d(aBase2dPnt1, aBase2dVec1, aPnt2dOnS1, aTol1) &&
2995         IsInsideIn2d(aBase2dPnt2, aBase2dVec2, aPnt2dOnS2, aTol2) &&
2996         IsInsideIn3d(aBase3dPnt, aBase3dVec, aPnt3d, aTol3d) )
2997     {
2998       // Handle possible uneven parametrization on one of 2d subspaces.
2999       // Delete point only when expected lengths are close to each other (aLimitCoeff).
3000       // Example:
3001       // c2d1 - line
3002       // c3d - line
3003       // c2d2 - geometrically line, but have uneven parametrization -> c2d2 is bspline.
3004       gp_XY aPntOnS1[2]= { gp_XY(UonS1[1] - UonS1[0], VonS1[1] - VonS1[0])
3005                          , gp_XY(UonS1[2] - UonS1[1], VonS1[2] - VonS1[1])};
3006       gp_XY aPntOnS2[2]= { gp_XY(UonS2[1] - UonS2[0], VonS2[1] - VonS2[0])
3007                          , gp_XY(UonS2[2] - UonS2[1], VonS2[2] - VonS2[1])};
3008
3009       Standard_Real aStepOnS1 = aPntOnS1[0].SquareModulus() / aPntOnS1[1].SquareModulus();
3010       Standard_Real aStepOnS2 = aPntOnS2[0].SquareModulus() / aPntOnS2[1].SquareModulus();
3011
3012       Standard_Real aStepCoeff = Min(aStepOnS1, aStepOnS2) / Max(aStepOnS1, aStepOnS2);
3013
3014       if (aStepCoeff > aLimitCoeff)
3015       {
3016         // Set hash flag to "Delete" state.
3017         isDeleteState = Standard_True;
3018         aNewPointsHash.SetValue(i - 1, 1);
3019
3020         // Change middle point.
3021         UonS1[1] = UonS1[2];
3022         UonS2[1] = UonS2[2];
3023         VonS1[1] = VonS1[2];
3024         VonS2[1] = VonS2[2];
3025       }
3026     }
3027
3028     if (!isDeleteState)
3029     {
3030       // Compute new pipe parameters.
3031       UonS1[0] = UonS1[1];
3032       VonS1[0] = VonS1[1];
3033       UonS2[0] = UonS2[1];
3034       VonS2[0] = VonS2[1];
3035
3036       UonS1[1] = UonS1[2];
3037       VonS1[1] = VonS1[2];
3038       UonS2[1] = UonS2[2];
3039       VonS2[1] = VonS2[2];
3040
3041       aBase2dPnt1.SetCoord(UonS1[0], VonS1[0]);
3042       aBase2dPnt2.SetCoord(UonS2[0], VonS2[0]);
3043       aBase2dVec1.SetCoord(UonS1[1] - UonS1[0], VonS1[1] - VonS1[0]);
3044       aBase2dVec2.SetCoord(UonS2[1] - UonS2[0], VonS2[1] - VonS2[0]);
3045       aBase3dPnt = aLineOn2S->Value(i - 1).Value();
3046       aBase3dVec = gp_Vec(aLineOn2S->Value(i - 1).Value(), aLineOn2S->Value(i).Value());
3047
3048       aNbPnt++;
3049     }
3050   }
3051
3052   // Workaround to handle case of small amount of points after purge.
3053   // Test "boolean boptuc_complex B5" and similar.
3054   // This is possible since there are at least two points.
3055   if (aNewPointsHash(1) == -1 &&
3056       aNewPointsHash(2) == -1 &&
3057       aNbPnt <= 3)
3058   {
3059     // Delete first.
3060     aNewPointsHash(1) = 1;
3061   }
3062   if (aNewPointsHash(aLineOn2S->NbPoints() - 1) == -1 &&
3063       aNewPointsHash(aLineOn2S->NbPoints()    ) == -1 &&
3064       aNbPnt <= 3)
3065   {
3066     // Delete last.
3067     aNewPointsHash(aLineOn2S->NbPoints() ) = 1;
3068   }
3069
3070   // Purgre when too small amount of points left.
3071   if (aNbPnt <= 2)
3072   {
3073     for(i = aNewPointsHash.Lower(); i <= aNewPointsHash.Upper(); i++)
3074     {
3075       if (aNewPointsHash(i) != -1)
3076       {
3077         aNewPointsHash(i) = 1;
3078       }
3079     }
3080   }
3081
3082   // Handle possible bad distribution of points,
3083   // which are will converted into one single bezier curve (less than 30 points).
3084   // Make distribution more even:
3085   // max step will be nearly to 0.1 of param distance.
3086   if (aNbPnt + 2 > aMinNbBadDistr &&
3087       aNbPnt + 2 < aNbSingleBezier )
3088   {
3089     for(Standard_Integer anIdx = 1; anIdx <= 8; anIdx++)
3090     {
3091       Standard_Integer aHashIdx =
3092         Standard_Integer(anIdx * aLineOn2S->NbPoints() / 9);
3093
3094       //Store this point.
3095       aNewPointsHash(aHashIdx) = 0;
3096     }
3097   }
3098
3099   aTmpWLine = aLocalWLine;
3100   Handle(IntSurf_LineOn2S) aPurgedLineOn2S = new IntSurf_LineOn2S();
3101   aLocalWLine = new IntPatch_WLine(aPurgedLineOn2S, Standard_False);
3102   Standard_Integer anOldLineIdx = 1, aVertexIdx = 1;
3103   for(i = 1; i <= aNewPointsHash.Upper(); i++)
3104   {
3105     if (aNewPointsHash(i) == 0)
3106     {
3107       // Store this point.
3108       aPurgedLineOn2S->Add(aLineOn2S->Value(i));
3109       anOldLineIdx++;
3110     }
3111     else if (aNewPointsHash(i) == -1)
3112     {
3113       // Add vertex.
3114       IntPatch_Point aVertex = aTmpWLine->Vertex(aVertexIdx++);
3115       aVertex.SetParameter(anOldLineIdx++);
3116       aLocalWLine->AddVertex(aVertex);
3117       aPurgedLineOn2S->Add(aLineOn2S->Value(i));
3118     }
3119   }
3120
3121   if(aPurgedLineOn2S->NbPoints() > 1) {
3122     aResult = aLocalWLine;
3123   }
3124   return aResult;
3125 }
3126
3127 //=======================================================================
3128 //function : TolR3d
3129 //purpose  :
3130 //=======================================================================
3131 void TolR3d(const TopoDS_Face& aF1,
3132             const TopoDS_Face& aF2,
3133             Standard_Real& myTolReached3d)
3134 {
3135   Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
3136
3137   aTolTresh=2.999999e-3;
3138   aTolF1 = BRep_Tool::Tolerance(aF1);
3139   aTolF2 = BRep_Tool::Tolerance(aF2);
3140   aTolFMax=Max(aTolF1, aTolF2);
3141
3142   if (aTolFMax>aTolTresh) {
3143     myTolReached3d=aTolFMax;
3144   }
3145 }
3146 //=======================================================================
3147 //function : IsPointOnBoundary
3148 //purpose  :
3149 //=======================================================================
3150 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
3151                                    const Standard_Real theFirstBoundary,
3152                                    const Standard_Real theSecondBoundary,
3153                                    const Standard_Real theResolution,
3154                                    Standard_Boolean&   IsOnFirstBoundary)
3155 {
3156   Standard_Boolean bRet;
3157   Standard_Integer i;
3158   Standard_Real adist;
3159   //
3160   bRet=Standard_False;
3161   for(i = 0; i < 2; ++i) {
3162     IsOnFirstBoundary = (i == 0);
3163     if (IsOnFirstBoundary) {
3164       adist = fabs(theParameter - theFirstBoundary);
3165     }
3166     else {
3167       adist = fabs(theParameter - theSecondBoundary);
3168     }
3169     if(adist < theResolution) {
3170       return !bRet;
3171     }
3172   }
3173   return bRet;
3174 }
3175 // ------------------------------------------------------------------------------------------------
3176 // static function: FindPoint
3177 // purpose:
3178 // ------------------------------------------------------------------------------------------------
3179 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
3180                            const gp_Pnt2d&     theLastPoint,
3181                            const Standard_Real theUmin,
3182                            const Standard_Real theUmax,
3183                            const Standard_Real theVmin,
3184                            const Standard_Real theVmax,
3185                            gp_Pnt2d&           theNewPoint) {
3186
3187   gp_Vec2d aVec(theFirstPoint, theLastPoint);
3188   Standard_Integer i = 0, j = 0;
3189
3190   for(i = 0; i < 4; i++) {
3191     gp_Vec2d anOtherVec;
3192     gp_Vec2d anOtherVecNormal;
3193     gp_Pnt2d aprojpoint = theLastPoint;
3194
3195     if((i % 2) == 0) {
3196       anOtherVec.SetX(0.);
3197       anOtherVec.SetY(1.);
3198       anOtherVecNormal.SetX(1.);
3199       anOtherVecNormal.SetY(0.);
3200
3201       if(i < 2)
3202         aprojpoint.SetX(theUmin);
3203       else
3204         aprojpoint.SetX(theUmax);
3205     }
3206     else {
3207       anOtherVec.SetX(1.);
3208       anOtherVec.SetY(0.);
3209       anOtherVecNormal.SetX(0.);
3210       anOtherVecNormal.SetY(1.);
3211
3212       if(i < 2)
3213         aprojpoint.SetY(theVmin);
3214       else
3215         aprojpoint.SetY(theVmax);
3216     }
3217     gp_Vec2d anormvec = aVec;
3218     anormvec.Normalize();
3219     RefineVector(anormvec);
3220     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
3221
3222     if(fabs(adot1) < Precision::Angular())
3223       continue;
3224     Standard_Real adist = 0.;
3225     Standard_Boolean bIsOut = Standard_False;
3226
3227     if((i % 2) == 0) {
3228       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
3229       bIsOut = (i < 2) ? (theLastPoint.X() < theUmin) : (theLastPoint.X() > theUmax);
3230     }
3231     else {
3232       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
3233       bIsOut = (i < 2) ? (theLastPoint.Y() < theVmin) : (theLastPoint.Y() > theVmax);
3234     }
3235     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
3236
3237     for(j = 0; j < 2; j++) {
3238       anoffset = (j == 0) ? anoffset : -anoffset;
3239       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
3240       gp_Vec2d acurvec(theLastPoint, acurpoint);
3241       if ( bIsOut )
3242         acurvec.Reverse();
3243
3244       Standard_Real aDotX, anAngleX;
3245       //
3246       aDotX = aVec.Dot(acurvec);
3247       anAngleX = aVec.Angle(acurvec);
3248       //
3249       if(aDotX > 0. && fabs(anAngleX) < Precision::PConfusion()) {
3250         if((i % 2) == 0) {
3251           if((acurpoint.Y() >= theVmin) &&
3252              (acurpoint.Y() <= theVmax)) {
3253             theNewPoint = acurpoint;
3254             return Standard_True;
3255           }
3256         }
3257         else {
3258           if((acurpoint.X() >= theUmin) &&
3259              (acurpoint.X() <= theUmax)) {
3260             theNewPoint = acurpoint;
3261             return Standard_True;
3262           }
3263         }
3264       }
3265     }
3266   }
3267   return Standard_False;
3268 }
3269
3270
3271 // ------------------------------------------------------------------------------------------------
3272 // static function: FindPoint
3273 // purpose: Find point on the boundary of radial tangent zone
3274 // ------------------------------------------------------------------------------------------------
3275 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
3276                            const gp_Pnt2d&     theLastPoint,
3277                            const Standard_Real theUmin,
3278                            const Standard_Real theUmax,
3279                            const Standard_Real theVmin,
3280                            const Standard_Real theVmax,
3281                            const gp_Pnt2d&     theTanZoneCenter,
3282                            const Standard_Real theZoneRadius,
3283                            Handle(GeomAdaptor_HSurface) theGASurface,
3284                            gp_Pnt2d&           theNewPoint) {
3285   theNewPoint = theLastPoint;
3286
3287   if ( !IsInsideTanZone( theLastPoint, theTanZoneCenter, theZoneRadius, theGASurface) )
3288     return Standard_False;
3289
3290   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3291   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3292
3293   Standard_Real aRadius = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3294   gp_Ax22d anAxis( theTanZoneCenter, gp_Dir2d(1, 0), gp_Dir2d(0, 1) );
3295   gp_Circ2d aCircle( anAxis, aRadius );
3296
3297   //
3298   gp_Vec2d aDir( theLastPoint.XY() - theFirstPoint.XY() );
3299   Standard_Real aLength = aDir.Magnitude();
3300   if ( aLength <= gp::Resolution() )
3301     return Standard_False;
3302   gp_Lin2d aLine( theFirstPoint, aDir );
3303
3304   //
3305   Handle(Geom2d_Line) aCLine = new Geom2d_Line( aLine );
3306   Handle(Geom2d_TrimmedCurve) aC1 = new Geom2d_TrimmedCurve( aCLine, 0, aLength );
3307   Handle(Geom2d_Circle) aC2 = new Geom2d_Circle( aCircle );
3308
3309   Standard_Real aTol = aRadius * 0.001;
3310   aTol = ( aTol < Precision::PConfusion() ) ? Precision::PConfusion() : aTol;
3311
3312   Geom2dAPI_InterCurveCurve anIntersector;
3313   anIntersector.Init( aC1, aC2, aTol );
3314
3315   if ( anIntersector.NbPoints() == 0 )
3316     return Standard_False;
3317
3318   Standard_Boolean aFound = Standard_False;
3319   Standard_Real aMinDist = aLength * aLength;
3320   Standard_Integer i = 0;
3321   for ( i = 1; i <= anIntersector.NbPoints(); i++ ) {
3322     gp_Pnt2d aPInt = anIntersector.Point( i );
3323     if ( aPInt.SquareDistance( theFirstPoint ) < aMinDist ) {
3324       if ( ( aPInt.X() >= theUmin ) && ( aPInt.X() <= theUmax ) &&
3325            ( aPInt.Y() >= theVmin ) && ( aPInt.Y() <= theVmax ) ) {
3326         theNewPoint = aPInt;
3327         aFound = Standard_True;
3328       }
3329     }
3330   }
3331
3332   return aFound;
3333 }
3334
3335 // ------------------------------------------------------------------------------------------------
3336 // static function: IsInsideTanZone
3337 // purpose: Check if point is inside a radial tangent zone
3338 // ------------------------------------------------------------------------------------------------
3339 Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
3340                                  const gp_Pnt2d&     theTanZoneCenter,
3341                                  const Standard_Real theZoneRadius,
3342                                  Handle(GeomAdaptor_HSurface) theGASurface) {
3343
3344   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3345   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3346   Standard_Real aRadiusSQR = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3347   aRadiusSQR *= aRadiusSQR;
3348   if ( thePoint.SquareDistance( theTanZoneCenter ) <= aRadiusSQR )
3349     return Standard_True;
3350   return Standard_False;
3351 }
3352
3353 // ------------------------------------------------------------------------------------------------
3354 // static function: CheckTangentZonesExist
3355 // purpose: Check if tangent zone exists
3356 // ------------------------------------------------------------------------------------------------
3357 Standard_Boolean CheckTangentZonesExist( const Handle(GeomAdaptor_HSurface)& theSurface1,
3358                                         const Handle(GeomAdaptor_HSurface)&  theSurface2 )
3359 {
3360   if ( ( theSurface1->GetType() != GeomAbs_Torus ) ||
3361       ( theSurface2->GetType() != GeomAbs_Torus ) )
3362     return Standard_False;
3363
3364   gp_Torus aTor1 = theSurface1->Torus();
3365   gp_Torus aTor2 = theSurface2->Torus();
3366
3367   if ( aTor1.Location().Distance( aTor2.Location() ) > Precision::Confusion() )
3368     return Standard_False;
3369
3370   if ( ( fabs( aTor1.MajorRadius() - aTor2.MajorRadius() ) > Precision::Confusion() ) ||
3371        ( fabs( aTor1.MinorRadius() - aTor2.MinorRadius() ) > Precision::Confusion() ) )
3372     return Standard_False;
3373
3374   if ( ( aTor1.MajorRadius() < aTor1.MinorRadius() ) ||
3375        ( aTor2.MajorRadius() < aTor2.MinorRadius() ) )
3376     return Standard_False;
3377   return Standard_True;
3378 }
3379
3380 // ------------------------------------------------------------------------------------------------
3381 // static function: ComputeTangentZones
3382 // purpose:
3383 // ------------------------------------------------------------------------------------------------
3384 Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
3385                                      const Handle(GeomAdaptor_HSurface)&  theSurface2,
3386                                      const TopoDS_Face&                   theFace1,
3387                                      const TopoDS_Face&                   theFace2,
3388                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
3389                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
3390                                      Handle(TColStd_HArray1OfReal)&       theResultRadius,
3391                                      const Handle(IntTools_Context)& aContext)
3392 {
3393   Standard_Integer aResult = 0;
3394   if ( !CheckTangentZonesExist( theSurface1, theSurface2 ) )
3395     return aResult;
3396
3397
3398   TColgp_SequenceOfPnt2d aSeqResultS1, aSeqResultS2;
3399   TColStd_SequenceOfReal aSeqResultRad;
3400
3401   gp_Torus aTor1 = theSurface1->Torus();
3402   gp_Torus aTor2 = theSurface2->Torus();
3403
3404   gp_Ax2 anax1( aTor1.Location(), aTor1.Axis().Direction() );
3405   gp_Ax2 anax2( aTor2.Location(), aTor2.Axis().Direction() );
3406   Standard_Integer j = 0;
3407
3408   for ( j = 0; j < 2; j++ ) {
3409     Standard_Real aCoef = ( j == 0 ) ? -1 : 1;
3410     Standard_Real aRadius1 = fabs(aTor1.MajorRadius() + aCoef * aTor1.MinorRadius());
3411     Standard_Real aRadius2 = fabs(aTor2.MajorRadius() + aCoef * aTor2.MinorRadius());
3412
3413     gp_Circ aCircle1( anax1, aRadius1 );
3414     gp_Circ aCircle2( anax2, aRadius2 );
3415
3416     // roughly compute radius of tangent zone for perpendicular case
3417     Standard_Real aCriteria = Precision::Confusion() * 0.5;
3418
3419     Standard_Real aT1 = aCriteria;
3420     Standard_Real aT2 = aCriteria;
3421     if ( j == 0 ) {
3422       // internal tangency
3423       Standard_Real aR = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3424       //aT1 = aCriteria * aCriteria + aR * aR - ( aR - aCriteria ) * ( aR - aCriteria );
3425       aT1 = 2. * aR * aCriteria;
3426       aT2 = aT1;
3427     }
3428     else {
3429       // external tangency
3430       Standard_Real aRb = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3431       Standard_Real aRm = ( aRadius1 < aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3432       Standard_Real aDelta = aRb - aCriteria;
3433       aDelta *= aDelta;
3434       aDelta -= aRm * aRm;
3435       aDelta /= 2. * (aRb - aRm);
3436       aDelta -= 0.5 * (aRb - aRm);
3437
3438       aT1 = 2. * aRm * (aRm - aDelta);
3439       aT2 = aT1;
3440     }
3441     aCriteria = ( aT1 > aT2) ? aT1 : aT2;
3442     if ( aCriteria > 0 )
3443       aCriteria = sqrt( aCriteria );
3444
3445     if ( aCriteria > 0.5 * aTor1.MinorRadius() ) {
3446       // too big zone -> drop to minimum
3447       aCriteria = Precision::Confusion();
3448     }
3449
3450     GeomAdaptor_Curve aC1( new Geom_Circle(aCircle1) );
3451     GeomAdaptor_Curve aC2( new Geom_Circle(aCircle2) );
3452     Extrema_ExtCC anExtrema(aC1, aC2, 0, 2. * M_PI, 0, 2. * M_PI,
3453                             Precision::PConfusion(), Precision::PConfusion());
3454
3455     if ( anExtrema.IsDone() ) {
3456
3457       Standard_Integer i = 0;
3458       for ( i = 1; i <= anExtrema.NbExt(); i++ ) {
3459         if ( anExtrema.SquareDistance(i) > aCriteria * aCriteria )
3460           continue;
3461
3462         Extrema_POnCurv P1, P2;
3463         anExtrema.Points( i, P1, P2 );
3464
3465         Standard_Boolean bFoundResult = Standard_True;
3466         gp_Pnt2d pr1, pr2;
3467
3468         Standard_Integer surfit = 0;
3469         for ( surfit = 0; surfit < 2; surfit++ ) {
3470           GeomAPI_ProjectPointOnSurf& aProjector =
3471             (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
3472
3473           gp_Pnt aP3d = (surfit == 0) ? P1.Value() : P2.Value();
3474           aProjector.Perform(aP3d);
3475
3476           if(!aProjector.IsDone())
3477             bFoundResult = Standard_False;
3478           else {
3479             if(aProjector.LowerDistance() > aCriteria) {
3480               bFoundResult = Standard_False;
3481             }
3482             else {
3483               Standard_Real foundU = 0, foundV = 0;
3484               aProjector.LowerDistanceParameters(foundU, foundV);
3485               if ( surfit == 0 )
3486                 pr1 = gp_Pnt2d( foundU, foundV );
3487               else
3488                 pr2 = gp_Pnt2d( foundU, foundV );
3489             }
3490           }
3491         }
3492         if ( bFoundResult ) {
3493           aSeqResultS1.Append( pr1 );
3494           aSeqResultS2.Append( pr2 );
3495           aSeqResultRad.Append( aCriteria );
3496
3497           // torus is u and v periodic
3498           const Standard_Real twoPI = M_PI + M_PI;
3499           Standard_Real arr1tmp[2] = {pr1.X(), pr1.Y()};
3500           Standard_Real arr2tmp[2] = {pr2.X(), pr2.Y()};
3501
3502           // iteration on period bounds
3503           for ( Standard_Integer k1 = 0; k1 < 2; k1++ ) {
3504             Standard_Real aBound = ( k1 == 0 ) ? 0 : twoPI;
3505             Standard_Real aShift = ( k1 == 0 ) ? twoPI : -twoPI;
3506
3507             // iteration on surfaces
3508             for ( Standard_Integer k2 = 0; k2 < 2; k2++ ) {
3509               Standard_Real* arr1 = ( k2 == 0 ) ? arr1tmp : arr2tmp;
3510               Standard_Real* arr2 = ( k2 != 0 ) ? arr1tmp : arr2tmp;
3511               TColgp_SequenceOfPnt2d& aSeqS1 = ( k2 == 0 ) ? aSeqResultS1 : aSeqResultS2;
3512               TColgp_SequenceOfPnt2d& aSeqS2 = ( k2 != 0 ) ? aSeqResultS1 : aSeqResultS2;
3513
3514               if (fabs(arr1[0] - aBound) < Precision::PConfusion()) {
3515                 aSeqS1.Append( gp_Pnt2d( arr1[0] + aShift, arr1[1] ) );
3516                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
3517                 aSeqResultRad.Append( aCriteria );
3518               }
3519               if (fabs(arr1[1] - aBound) < Precision::PConfusion()) {
3520                 aSeqS1.Append( gp_Pnt2d( arr1[0], arr1[1] + aShift) );
3521                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
3522                 aSeqResultRad.Append( aCriteria );
3523               }
3524             }
3525           } //
3526         }
3527       }
3528     }
3529   }
3530   aResult = aSeqResultRad.Length();
3531
3532   if ( aResult > 0 ) {
3533     theResultOnS1 = new TColgp_HArray1OfPnt2d( 1, aResult );
3534     theResultOnS2 = new TColgp_HArray1OfPnt2d( 1, aResult );
3535     theResultRadius = new TColStd_HArray1OfReal( 1, aResult );
3536
3537     for ( Standard_Integer i = 1 ; i <= aResult; i++ ) {
3538       theResultOnS1->SetValue( i, aSeqResultS1.Value(i) );
3539       theResultOnS2->SetValue( i, aSeqResultS2.Value(i) );
3540       theResultRadius->SetValue( i, aSeqResultRad.Value(i) );
3541     }
3542   }
3543   return aResult;
3544 }
3545
3546 // ------------------------------------------------------------------------------------------------
3547 // static function: AdjustByNeighbour
3548 // purpose:
3549 // ------------------------------------------------------------------------------------------------
3550 gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
3551                            const gp_Pnt2d&     theOriginalPoint,
3552                            Handle(GeomAdaptor_HSurface) theGASurface) {
3553
3554   gp_Pnt2d ap1 = theaNeighbourPoint;
3555   gp_Pnt2d ap2 = theOriginalPoint;
3556
3557   if ( theGASurface->IsUPeriodic() ) {
3558     Standard_Real aPeriod     = theGASurface->UPeriod();
3559     gp_Pnt2d aPTest = ap2;
3560     Standard_Real aSqDistMin = 1.e+100;
3561
3562     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
3563       aPTest.SetX( theOriginalPoint.X() + aPeriod * pIt );
3564       Standard_Real dd = ap1.SquareDistance( aPTest );
3565
3566       if ( dd < aSqDistMin ) {
3567         ap2 = aPTest;
3568         aSqDistMin = dd;
3569       }
3570     }
3571   }
3572   if ( theGASurface->IsVPeriodic() ) {
3573     Standard_Real aPeriod     = theGASurface->VPeriod();
3574     gp_Pnt2d aPTest = ap2;
3575     Standard_Real aSqDistMin = 1.e+100;
3576
3577     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
3578       aPTest.SetY( theOriginalPoint.Y() + aPeriod * pIt );
3579       Standard_Real dd = ap1.SquareDistance( aPTest );
3580
3581       if ( dd < aSqDistMin ) {
3582         ap2 = aPTest;
3583         aSqDistMin = dd;
3584       }
3585     }
3586   }
3587   return ap2;
3588 }
3589
3590 // ------------------------------------------------------------------------------------------------
3591 //function: DecompositionOfWLine
3592 // purpose:
3593 // ------------------------------------------------------------------------------------------------
3594 Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
3595                                       const Handle(GeomAdaptor_HSurface)&            theSurface1,
3596                                       const Handle(GeomAdaptor_HSurface)&            theSurface2,
3597                                       const TopoDS_Face&                             theFace1,
3598                                       const TopoDS_Face&                             theFace2,
3599                                       const GeomInt_LineConstructor&                 theLConstructor,
3600                                       const Standard_Boolean                         theAvoidLConstructor,
3601                                       IntPatch_SequenceOfLine&                       theNewLines,
3602                                       Standard_Real&                                 theReachedTol3d,
3603                                       const Handle(IntTools_Context)& aContext)
3604 {
3605
3606   Standard_Boolean bRet, bAvoidLineConstructor;
3607   Standard_Integer aNbPnts, aNbParts;
3608   //
3609   bRet=Standard_False;
3610   aNbPnts=theWLine->NbPnts();
3611   bAvoidLineConstructor=theAvoidLConstructor;
3612   //
3613   if(!aNbPnts){
3614     return bRet;
3615   }
3616   if (!bAvoidLineConstructor) {
3617     aNbParts=theLConstructor.NbParts();
3618     if (!aNbParts) {
3619       return bRet;
3620     }
3621   }
3622   //
3623   Standard_Boolean bIsPrevPointOnBoundary, bIsPointOnBoundary, bIsCurrentPointOnBoundary;
3624   Standard_Integer nblines, pit, i, j;
3625   Standard_Real aTol;
3626   TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts);
3627   TColStd_Array1OfInteger       anArrayOfLineType(1, aNbPnts);
3628   TColStd_ListOfInteger aListOfPointIndex;
3629
3630   Handle(TColgp_HArray1OfPnt2d) aTanZoneS1;
3631   Handle(TColgp_HArray1OfPnt2d) aTanZoneS2;
3632   Handle(TColStd_HArray1OfReal) aTanZoneRadius;
3633   Standard_Integer aNbZone = ComputeTangentZones( theSurface1, theSurface2, theFace1, theFace2,
3634                                                  aTanZoneS1, aTanZoneS2, aTanZoneRadius, aContext);
3635
3636   //
3637   nblines=0;
3638   aTol=Precision::Confusion();
3639   aTol=0.5*aTol;
3640   bIsPrevPointOnBoundary=Standard_False;
3641   bIsPointOnBoundary=Standard_False;
3642   //
3643   // 1. ...
3644   //
3645   // Points
3646   for(pit = 1; pit <= aNbPnts; ++pit) {
3647     Standard_Boolean bIsOnFirstBoundary, isperiodic;
3648     Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
3649     Standard_Real aParameter, anoffset, anAdjustPar;
3650     Standard_Real umin, umax, vmin, vmax;
3651     //
3652     bIsCurrentPointOnBoundary = Standard_False;
3653     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
3654     //
3655     // Surface
3656     for(i = 0; i < 2; ++i) {
3657       Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
3658       aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
3659       if(!i) {
3660         aPoint.ParametersOnS1(U, V);
3661       }
3662       else {
3663         aPoint.ParametersOnS2(U, V);
3664       }
3665       // U, V
3666       for(j = 0; j < 2; j++) {
3667         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3668         if(!isperiodic){
3669           continue;
3670         }
3671         //
3672         if (!j) {
3673           aResolution=aGASurface->UResolution(aTol);
3674           aPeriod=aGASurface->UPeriod();
3675           alowerboundary=umin;
3676           aupperboundary=umax;
3677           aParameter=U;
3678         }
3679         else {
3680           aResolution=aGASurface->VResolution(aTol);
3681           aPeriod=aGASurface->VPeriod();
3682           alowerboundary=vmin;
3683           aupperboundary=vmax;
3684           aParameter=V;
3685         }
3686
3687         GeomInt::AdjustPeriodic(aParameter,
3688                                        alowerboundary,
3689                                        aupperboundary,
3690                                        aPeriod,
3691                                        anAdjustPar,
3692                                        anoffset);
3693         //
3694         bIsOnFirstBoundary = Standard_True;// ?
3695         bIsPointOnBoundary=
3696           IsPointOnBoundary(anAdjustPar,
3697                             alowerboundary,
3698                             aupperboundary,
3699                             aResolution,
3700                             bIsOnFirstBoundary);
3701         //
3702         if(bIsPointOnBoundary) {
3703           bIsCurrentPointOnBoundary = Standard_True;
3704           break;
3705         }
3706         else {
3707           // check if a point belong to a tangent zone. Begin
3708           Standard_Integer zIt = 0;
3709           for ( zIt = 1; zIt <= aNbZone; zIt++ ) {
3710             gp_Pnt2d aPZone = (i == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
3711             Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
3712
3713             if ( IsInsideTanZone(gp_Pnt2d( U, V ), aPZone, aZoneRadius, aGASurface ) ) {
3714               // set boundary flag to split the curve by a tangent zone
3715               bIsPointOnBoundary = Standard_True;
3716               bIsCurrentPointOnBoundary = Standard_True;
3717               if ( theReachedTol3d < aZoneRadius ) {
3718                 theReachedTol3d = aZoneRadius;
3719               }
3720               break;
3721             }
3722           }
3723         }
3724       }//for(j = 0; j < 2; j++) {
3725
3726       if(bIsCurrentPointOnBoundary){
3727         break;
3728       }
3729     }//for(i = 0; i < 2; ++i) {
3730     //
3731     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
3732       if(!aListOfPointIndex.IsEmpty()) {
3733         nblines++;
3734         anArrayOfLines.SetValue(nblines, aListOfPointIndex);
3735         anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
3736         aListOfPointIndex.Clear();
3737       }
3738       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
3739     }
3740     aListOfPointIndex.Append(pit);
3741   } //for(pit = 1; pit <= aNbPnts; ++pit) {
3742   //
3743   if(!aListOfPointIndex.IsEmpty()) {
3744     nblines++;
3745     anArrayOfLines.SetValue(nblines, aListOfPointIndex);
3746     anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
3747     aListOfPointIndex.Clear();
3748   }
3749   //
3750   if(nblines<=1) {
3751     return bRet; //Standard_False;
3752   }
3753   //
3754   //
3755   // 2. Correct wlines.begin
3756   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
3757   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S();
3758   //
3759   for(i = 1; i <= nblines; i++) {
3760     if(anArrayOfLineType.Value(i) != 0) {
3761       continue;
3762     }
3763     const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
3764     TColStd_ListOfInteger aListOfFLIndex;
3765
3766     for(j = 0; j < 2; j++) {
3767       Standard_Integer aneighbourindex = (j == 0) ? (i - 1) : (i + 1);
3768
3769       if((aneighbourindex < 1) || (aneighbourindex > nblines))
3770         continue;
3771
3772       if(anArrayOfLineType.Value(aneighbourindex) == 0)
3773         continue;
3774       const TColStd_ListOfInteger& aNeighbour = anArrayOfLines.Value(aneighbourindex);
3775       Standard_Integer anIndex = (j == 0) ? aNeighbour.Last() : aNeighbour.First();
3776       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
3777
3778       IntSurf_PntOn2S aNewP = aPoint;
3779       if(aListOfIndex.Extent() < 2) {
3780         aSeqOfPntOn2S->Add(aNewP);
3781         aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
3782         continue;
3783       }
3784       //
3785       Standard_Integer iFirst = aListOfIndex.First();
3786       Standard_Integer iLast  = aListOfIndex.Last();
3787       //
3788       for(Standard_Integer surfit = 0; surfit < 2; surfit++) {
3789
3790         Handle(GeomAdaptor_HSurface) aGASurface = (surfit == 0) ? theSurface1 : theSurface2;
3791         Standard_Real umin=0., umax=0., vmin=0., vmax=0.;
3792         aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
3793         Standard_Real U=0., V=0.;
3794
3795         if(surfit == 0)
3796           aNewP.ParametersOnS1(U, V);
3797         else
3798           aNewP.ParametersOnS2(U, V);
3799         Standard_Integer nbboundaries = 0;
3800
3801         Standard_Boolean bIsNearBoundary = Standard_False;
3802         Standard_Integer aZoneIndex = 0;
3803         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
3804         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
3805
3806
3807         for(Standard_Integer parit = 0; parit < 2; parit++) {
3808           Standard_Boolean isperiodic = (parit == 0) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3809
3810           Standard_Real aResolution = (parit == 0) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
3811           Standard_Real alowerboundary = (parit == 0) ? umin : vmin;
3812           Standard_Real aupperboundary = (parit == 0) ? umax : vmax;
3813
3814           Standard_Real aParameter = (parit == 0) ? U : V;
3815           Standard_Boolean bIsOnFirstBoundary = Standard_True;
3816
3817           if(!isperiodic) {
3818             bIsPointOnBoundary=
3819               IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
3820             if(bIsPointOnBoundary) {
3821               bIsUBoundary = (parit == 0);
3822               bIsFirstBoundary = bIsOnFirstBoundary;
3823               nbboundaries++;
3824             }
3825           }
3826           else {
3827             Standard_Real aPeriod     = (parit == 0) ? aGASurface->UPeriod() : aGASurface->VPeriod();
3828             Standard_Real anoffset, anAdjustPar;
3829             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
3830                                            aPeriod, anAdjustPar, anoffset);
3831
3832             bIsPointOnBoundary=
3833               IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
3834             if(bIsPointOnBoundary) {
3835               bIsUBoundary = (parit == 0);
3836               bIsFirstBoundary = bIsOnFirstBoundary;
3837               nbboundaries++;
3838             }
3839             else {
3840             //check neighbourhood of boundary
3841             Standard_Real anEpsilon = aResolution * 100.;
3842             Standard_Real aPart = ( aupperboundary - alowerboundary ) * 0.1;
3843             anEpsilon = ( anEpsilon > aPart ) ? aPart : anEpsilon;
3844
3845               bIsNearBoundary = IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary,
3846                                                   anEpsilon, bIsOnFirstBoundary);
3847
3848             }
3849           }
3850         }
3851
3852         // check if a point belong to a tangent zone. Begin
3853         for ( Standard_Integer zIt = 1; zIt <= aNbZone; zIt++ ) {
3854           gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
3855           Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
3856
3857           Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
3858           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3859           Standard_Real nU1, nV1;
3860
3861           if(surfit == 0)
3862             aNeighbourPoint.ParametersOnS1(nU1, nV1);
3863           else
3864             aNeighbourPoint.ParametersOnS2(nU1, nV1);
3865           gp_Pnt2d ap1(nU1, nV1);
3866           gp_Pnt2d ap2 = AdjustByNeighbour( ap1, gp_Pnt2d( U, V ), aGASurface );
3867
3868
3869           if ( IsInsideTanZone( ap2, aPZone, aZoneRadius, aGASurface ) ) {
3870             aZoneIndex = zIt;
3871             bIsNearBoundary = Standard_True;
3872             if ( theReachedTol3d < aZoneRadius ) {
3873               theReachedTol3d = aZoneRadius;
3874             }
3875           }
3876         }
3877         // check if a point belong to a tangent zone. End
3878         Standard_Boolean bComputeLineEnd = Standard_False;
3879
3880         if(nbboundaries == 2) {
3881           //xf
3882           bComputeLineEnd = Standard_True;
3883           //xt
3884         }
3885         else if(nbboundaries == 1) {
3886           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3887
3888           if(isperiodic) {
3889             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
3890             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
3891             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
3892             Standard_Real aParameter = (bIsUBoundary) ? U : V;
3893             Standard_Real anoffset, anAdjustPar;
3894             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
3895                                            aPeriod, anAdjustPar, anoffset);
3896
3897             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
3898             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
3899             anotherPar += anoffset;
3900             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
3901             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
3902             Standard_Real nU1, nV1;
3903
3904             if(surfit == 0)
3905               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3906             else
3907               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3908
3909             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
3910             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
3911             bComputeLineEnd = Standard_True;
3912             Standard_Boolean bCheckAngle1 = Standard_False;
3913             Standard_Boolean bCheckAngle2 = Standard_False;
3914             gp_Vec2d aNewVec;
3915             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
3916             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
3917
3918             if(((adist1 - adist2) > Precision::PConfusion()) &&
3919                (adist2 < (aPeriod / 4.))) {
3920               bCheckAngle1 = Standard_True;
3921               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
3922
3923               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
3924                 aNewP.SetValue((surfit == 0), anewU, anewV);
3925                 bCheckAngle1 = Standard_False;
3926               }
3927             }
3928             else if(adist1 < (aPeriod / 4.)) {
3929               bCheckAngle2 = Standard_True;
3930               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
3931
3932               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
3933                 bCheckAngle2 = Standard_False;
3934               }
3935             }
3936
3937             if(bCheckAngle1 || bCheckAngle2) {
3938               // assume there are at least two points in line (see "if" above)
3939               Standard_Integer anindexother = aneighbourpointindex;
3940
3941               while((anindexother <= iLast) && (anindexother >= iFirst)) {
3942                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
3943                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
3944                 Standard_Real nU2, nV2;
3945
3946                 if(surfit == 0)
3947                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
3948                 else
3949                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
3950                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
3951
3952                 if(aVecOld.SquareMagnitude() <= gp::Resolution()) {
3953                   continue;
3954                 }
3955                 else {
3956                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
3957
3958                   if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
3959
3960                     if(bCheckAngle1) {
3961                       Standard_Real U1, U2, V1, V2;
3962                       IntSurf_PntOn2S atmppoint = aNewP;
3963                       atmppoint.SetValue((surfit == 0), anewU, anewV);
3964                       atmppoint.Parameters(U1, V1, U2, V2);
3965                       gp_Pnt P1 = theSurface1->Value(U1, V1);
3966                       gp_Pnt P2 = theSurface2->Value(U2, V2);
3967                       gp_Pnt P0 = aPoint.Value();
3968
3969                       if(P0.IsEqual(P1, aTol) &&
3970                          P0.IsEqual(P2, aTol) &&
3971                          P1.IsEqual(P2, aTol)) {
3972                         bComputeLineEnd = Standard_False;
3973                         aNewP.SetValue((surfit == 0), anewU, anewV);
3974                       }
3975                     }
3976
3977                     if(bCheckAngle2) {
3978                       bComputeLineEnd = Standard_False;
3979                     }
3980                   }
3981                   break;
3982                 }
3983               } // end while(anindexother...)
3984             }
3985           }
3986         }
3987         else if ( bIsNearBoundary ) {
3988           bComputeLineEnd = Standard_True;
3989         }
3990
3991         if(bComputeLineEnd) {
3992
3993           gp_Pnt2d anewpoint;
3994           Standard_Boolean found = Standard_False;
3995
3996           if ( bIsNearBoundary ) {
3997             // re-compute point near natural boundary or near tangent zone
3998             Standard_Real u1, v1, u2, v2;
3999             aNewP.Parameters( u1, v1, u2, v2 );
4000             if(surfit == 0)
4001               anewpoint = gp_Pnt2d( u1, v1 );
4002             else
4003               anewpoint = gp_Pnt2d( u2, v2 );
4004
4005             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
4006             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
4007             Standard_Real nU1, nV1;
4008
4009             if(surfit == 0)
4010               aNeighbourPoint.ParametersOnS1(nU1, nV1);
4011             else
4012               aNeighbourPoint.ParametersOnS2(nU1, nV1);
4013             gp_Pnt2d ap1(nU1, nV1);
4014             gp_Pnt2d ap2;
4015
4016
4017             if ( aZoneIndex ) {
4018               // exclude point from a tangent zone
4019               anewpoint = AdjustByNeighbour( ap1, anewpoint, aGASurface );
4020               gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(aZoneIndex) : aTanZoneS2->Value(aZoneIndex);
4021               Standard_Real aZoneRadius = aTanZoneRadius->Value(aZoneIndex);
4022
4023               if ( FindPoint(ap1, anewpoint, umin, umax, vmin, vmax,
4024                              aPZone, aZoneRadius, aGASurface, ap2) ) {
4025                 anewpoint = ap2;
4026                 found = Standard_True;
4027               }
4028             }
4029             else if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) {
4030               // re-compute point near boundary if shifted on a period
4031               ap2 = AdjustByNeighbour( ap1, anewpoint, aGASurface );
4032
4033               if ( ( ap2.X() < umin ) || ( ap2.X() > umax ) ||
4034                   ( ap2.Y() < vmin ) || ( ap2.Y() > vmax ) ) {
4035                 found = FindPoint(ap1, ap2, umin, umax, vmin, vmax, anewpoint);
4036               }
4037               else {
4038                 anewpoint = ap2;
4039                 aNewP.SetValue( (surfit == 0), anewpoint.X(), anewpoint.Y() );
4040               }
4041             }
4042           }
4043           else {
4044
4045             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
4046             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
4047             Standard_Real nU1, nV1;
4048
4049             if(surfit == 0)
4050               aNeighbourPoint.ParametersOnS1(nU1, nV1);
4051             else
4052               aNeighbourPoint.ParametersOnS2(nU1, nV1);
4053             gp_Pnt2d ap1(nU1, nV1);
4054             gp_Pnt2d ap2(nU1, nV1);
4055             Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
4056
4057             while((aneighbourpointindex2 <= iLast) && (aneighbourpointindex2 >= iFirst)) {
4058               aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
4059               const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
4060               Standard_Real nU2, nV2;
4061
4062               if(surfit == 0)
4063                 aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
4064               else
4065                 aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
4066               ap2.SetX(nU2);
4067               ap2.SetY(nV2);
4068
4069               if(ap1.SquareDistance(ap2) > gp::Resolution()) {
4070                 break;
4071               }
4072             }
4073             found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
4074           }
4075
4076           if(found) {
4077             // check point
4078             Standard_Real aCriteria = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
4079             GeomAPI_ProjectPointOnSurf& aProjector =
4080               (surfit == 0) ? aContext->ProjPS(theFace2) : aContext->ProjPS(theFace1);
4081             Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
4082
4083             Handle(GeomAdaptor_HSurface) aSurfaceOther = (surfit == 0) ? theSurface2 : theSurface1;
4084
4085             gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
4086             aProjector.Perform(aP3d);
4087
4088             if(aProjector.IsDone()) {
4089               if(aProjector.LowerDistance() < aCriteria) {
4090                 Standard_Real foundU = U, foundV = V;
4091                 aProjector.LowerDistanceParameters(foundU, foundV);
4092
4093                 //Correction of projected coordinates. Begin
4094                 //Note, it may be shifted on a period
4095                 Standard_Integer aneindex1 = (j == 0) ? iFirst : iLast;
4096                 const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneindex1);
4097                 Standard_Real nUn, nVn;
4098
4099                 if(surfit == 0)
4100                   aNeighbourPoint.ParametersOnS2(nUn, nVn);
4101                 else
4102                   aNeighbourPoint.ParametersOnS1(nUn, nVn);
4103                 gp_Pnt2d aNeighbour2d(nUn, nVn);
4104                 gp_Pnt2d anAdjustedPoint = AdjustByNeighbour( aNeighbour2d, gp_Pnt2d(foundU, foundV), aSurfaceOther );
4105                 foundU = anAdjustedPoint.X();
4106                 foundV = anAdjustedPoint.Y();
4107
4108                 if ( ( anAdjustedPoint.X() < umin ) && ( anAdjustedPoint.X() > umax ) &&
4109                     ( anAdjustedPoint.Y() < vmin ) && ( anAdjustedPoint.Y() > vmax ) ) {
4110                   // attempt to roughly re-compute point
4111                   foundU = ( foundU < umin ) ? umin : foundU;
4112                   foundU = ( foundU > umax ) ? umax : foundU;
4113                   foundV = ( foundV < vmin ) ? vmin : foundV;
4114                   foundV = ( foundV > vmax ) ? vmax : foundV;
4115
4116                   GeomAPI_ProjectPointOnSurf& aProjector2 =
4117                     (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
4118
4119                   aP3d = aSurfaceOther->Value(foundU, foundV);
4120                   aProjector2.Perform(aP3d);
4121
4122                   if(aProjector2.IsDone()) {
4123                     if(aProjector2.LowerDistance() < aCriteria) {
4124                       Standard_Real foundU2 = anewpoint.X(), foundV2 = anewpoint.Y();
4125                       aProjector2.LowerDistanceParameters(foundU2, foundV2);
4126                       anewpoint.SetX(foundU2);
4127                       anewpoint.SetY(foundV2);
4128                     }
4129                   }
4130                 }
4131                 //Correction of projected coordinates. End
4132
4133                 if(surfit == 0)
4134                   aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
4135                 else
4136                   aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
4137               }
4138             }
4139           }
4140         }
4141       }
4142       aSeqOfPntOn2S->Add(aNewP);
4143       aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
4144     }
4145     anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
4146   }
4147   // Correct wlines.end
4148
4149   // Split wlines.begin
4150   Standard_Integer nbiter;
4151   //
4152   nbiter=1;
4153   if (!bAvoidLineConstructor) {
4154     nbiter=theLConstructor.NbParts();
4155   }
4156   //
4157   for(j = 1; j <= nbiter; ++j) {
4158     Standard_Real fprm, lprm;
4159     Standard_Integer ifprm, ilprm;
4160     //
4161     if(bAvoidLineConstructor) {
4162       ifprm = 1;
4163       ilprm = theWLine->NbPnts();
4164     }
4165     else {
4166       theLConstructor.Part(j, fprm, lprm);
4167       ifprm = (Standard_Integer)fprm;
4168       ilprm = (Standard_Integer)lprm;
4169     }
4170
4171     Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
4172     //
4173     for(i = 1; i <= nblines; i++) {
4174       if(anArrayOfLineType.Value(i) != 0) {
4175         continue;
4176       }
4177       const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
4178       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
4179       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
4180       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
4181
4182       if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
4183         bhasfirstpoint = (i != 1);
4184       }
4185
4186       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
4187         bhaslastpoint = (i != nblines);
4188       }
4189
4190       Standard_Integer iFirst = aListOfIndex.First();
4191       Standard_Integer iLast  = aListOfIndex.Last();
4192       Standard_Boolean bIsFirstInside = ((ifprm >= iFirst) && (ifprm <= iLast));
4193       Standard_Boolean bIsLastInside =  ((ilprm >= iFirst) && (ilprm <= iLast));
4194
4195       if(!bIsFirstInside && !bIsLastInside) {
4196         if((ifprm < iFirst) && (ilprm > iLast)) {
4197           // append whole line, and boundaries if neccesary
4198           if(bhasfirstpoint) {
4199             pit = aListOfFLIndex.First();
4200             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4201             aLineOn2S->Add(aP);
4202           }
4203           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4204
4205           for(; anIt.More(); anIt.Next()) {
4206             pit = anIt.Value();
4207             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4208             aLineOn2S->Add(aP);
4209           }
4210
4211           if(bhaslastpoint) {
4212             pit = aListOfFLIndex.Last();
4213             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4214             aLineOn2S->Add(aP);
4215           }
4216
4217           // check end of split line (end is almost always)
4218           Standard_Integer aneighbour = i + 1;
4219           Standard_Boolean bIsEndOfLine = Standard_True;
4220
4221           if(aneighbour <= nblines) {
4222             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
4223
4224             if((anArrayOfLineType.Value(aneighbour) != 0) &&
4225                (aListOfNeighbourIndex.IsEmpty())) {
4226               bIsEndOfLine = Standard_False;
4227             }
4228           }
4229
4230           if(bIsEndOfLine) {
4231             if(aLineOn2S->NbPoints() > 1) {
4232               Handle(IntPatch_WLine) aNewWLine =
4233                 new IntPatch_WLine(aLineOn2S, Standard_False);
4234               theNewLines.Append(aNewWLine);
4235             }
4236             aLineOn2S = new IntSurf_LineOn2S();
4237           }
4238         }
4239         continue;
4240       }
4241       // end if(!bIsFirstInside && !bIsLastInside)
4242
4243       if(bIsFirstInside && bIsLastInside) {
4244         // append inside points between ifprm and ilprm
4245         TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4246
4247         for(; anIt.More(); anIt.Next()) {
4248           pit = anIt.Value();
4249           if((pit < ifprm) || (pit > ilprm))
4250             continue;
4251           const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4252           aLineOn2S->Add(aP);
4253         }
4254       }
4255       else {
4256
4257         if(bIsFirstInside) {
4258           // append points from ifprm to last point + boundary point
4259           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4260
4261           for(; anIt.More(); anIt.Next()) {
4262             pit = anIt.Value();
4263             if(pit < ifprm)
4264               continue;
4265             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4266             aLineOn2S->Add(aP);
4267           }
4268
4269           if(bhaslastpoint) {
4270             pit = aListOfFLIndex.Last();
4271             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4272             aLineOn2S->Add(aP);
4273           }
4274           // check end of split line (end is almost always)
4275           Standard_Integer aneighbour = i + 1;
4276           Standard_Boolean bIsEndOfLine = Standard_True;
4277
4278           if(aneighbour <= nblines) {
4279             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
4280
4281             if((anArrayOfLineType.Value(aneighbour) != 0) &&
4282                (aListOfNeighbourIndex.IsEmpty())) {
4283               bIsEndOfLine = Standard_False;
4284             }
4285           }
4286
4287           if(bIsEndOfLine) {
4288             if(aLineOn2S->NbPoints() > 1) {
4289               Handle(IntPatch_WLine) aNewWLine =
4290                 new IntPatch_WLine(aLineOn2S, Standard_False);
4291               theNewLines.Append(aNewWLine);
4292             }
4293             aLineOn2S = new IntSurf_LineOn2S();
4294           }
4295         }
4296         // end if(bIsFirstInside)
4297
4298         if(bIsLastInside) {
4299           // append points from first boundary point to ilprm
4300           if(bhasfirstpoint) {
4301             pit = aListOfFLIndex.First();
4302             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4303             aLineOn2S->Add(aP);
4304           }
4305           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4306
4307           for(; anIt.More(); anIt.Next()) {
4308             pit = anIt.Value();
4309             if(pit > ilprm)
4310               continue;
4311             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4312             aLineOn2S->Add(aP);
4313           }
4314         }
4315         //end if(bIsLastInside)
4316       }
4317     }
4318
4319     if(aLineOn2S->NbPoints() > 1) {
4320       Handle(IntPatch_WLine) aNewWLine =
4321         new IntPatch_WLine(aLineOn2S, Standard_False);
4322       theNewLines.Append(aNewWLine);
4323     }
4324   }
4325   // Split wlines.end
4326
4327   return Standard_True;
4328 }
4329
4330 // ------------------------------------------------------------------------------------------------
4331 // static function: ParameterOutOfBoundary
4332 // purpose:         Computes a new parameter for given curve. The corresponding 2d points
4333 //                  does not lay on any boundary of given faces
4334 // ------------------------------------------------------------------------------------------------
4335 Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
4336                                         const Handle(Geom_Curve)& theCurve,
4337                                         const TopoDS_Face&        theFace1,
4338                                         const TopoDS_Face&        theFace2,
4339                                         const Standard_Real       theOtherParameter,
4340                                         const Standard_Boolean    bIncreasePar,
4341                                         Standard_Real&            theNewParameter,
4342                                         const Handle(IntTools_Context)& aContext)
4343 {
4344   Standard_Boolean bIsComputed = Standard_False;
4345   theNewParameter = theParameter;
4346
4347   Standard_Real acurpar = theParameter;
4348   TopAbs_State aState = TopAbs_ON;
4349   Standard_Integer iter = 0;
4350   Standard_Real asumtol = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
4351   Standard_Real adelta = asumtol * 0.1;
4352   adelta = (adelta < Precision::Confusion()) ? Precision::Confusion() : adelta;
4353   Handle(Geom_Surface) aSurf1 = BRep_Tool::Surface(theFace1);
4354   Handle(Geom_Surface) aSurf2 = BRep_Tool::Surface(theFace2);
4355
4356   Standard_Real u1, u2, v1, v2;
4357
4358   GeomAPI_ProjectPointOnSurf aPrj1;
4359   aSurf1->Bounds(u1, u2, v1, v2);
4360   aPrj1.Init(aSurf1, u1, u2, v1, v2);
4361
4362   GeomAPI_ProjectPointOnSurf aPrj2;
4363   aSurf2->Bounds(u1, u2, v1, v2);
4364   aPrj2.Init(aSurf2, u1, u2, v1, v2);
4365
4366   while(aState == TopAbs_ON) {
4367     if(bIncreasePar)
4368       acurpar += adelta;
4369     else
4370       acurpar -= adelta;
4371     gp_Pnt aPCurrent = theCurve->Value(acurpar);
4372     aPrj1.Perform(aPCurrent);
4373     Standard_Real U=0., V=0.;
4374
4375     if(aPrj1.IsDone()) {
4376       aPrj1.LowerDistanceParameters(U, V);
4377       aState = aContext->StatePointFace(theFace1, gp_Pnt2d(U, V));
4378     }
4379
4380     if(aState != TopAbs_ON) {
4381       aPrj2.Perform(aPCurrent);
4382
4383       if(aPrj2.IsDone()) {
4384         aPrj2.LowerDistanceParameters(U, V);
4385         aState = aContext->StatePointFace(theFace2, gp_Pnt2d(U, V));
4386       }
4387     }
4388
4389     if(iter > 11) {
4390       break;
4391     }
4392     iter++;
4393   }
4394
4395   if(iter <= 11) {
4396     theNewParameter = acurpar;
4397     bIsComputed = Standard_True;
4398
4399     if(bIncreasePar) {
4400       if(acurpar >= theOtherParameter)
4401         theNewParameter = theOtherParameter;
4402     }
4403     else {
4404       if(acurpar <= theOtherParameter)
4405         theNewParameter = theOtherParameter;
4406     }
4407   }
4408   return bIsComputed;
4409 }
4410
4411 //=======================================================================
4412 //function : IsCurveValid
4413 //purpose  :
4414 //=======================================================================
4415 Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve)
4416 {
4417   if(thePCurve.IsNull())
4418     return Standard_False;
4419
4420   Standard_Real tolint = 1.e-10;
4421   Geom2dAdaptor_Curve PCA;
4422   IntRes2d_Domain PCD;
4423   Geom2dInt_GInter PCI;
4424
4425   Standard_Real pf = 0., pl = 0.;
4426   gp_Pnt2d pntf, pntl;
4427
4428   if(!thePCurve->IsClosed() && !thePCurve->IsPeriodic()) {
4429     pf = thePCurve->FirstParameter();
4430     pl = thePCurve->LastParameter();
4431     pntf = thePCurve->Value(pf);
4432     pntl = thePCurve->Value(pl);
4433     PCA.Load(thePCurve);
4434     if(!PCA.IsPeriodic()) {
4435       if(PCA.FirstParameter() > pf) pf = PCA.FirstParameter();
4436       if(PCA.LastParameter()  < pl) pl = PCA.LastParameter();
4437     }
4438     PCD.SetValues(pntf,pf,tolint,pntl,pl,tolint);
4439     PCI.Perform(PCA,PCD,tolint,tolint);
4440     if(PCI.IsDone())
4441       if(PCI.NbPoints() > 0) {
4442         return Standard_False;
4443       }
4444   }
4445
4446   return Standard_True;
4447 }
4448
4449 //=======================================================================
4450 //static function : ApproxWithPCurves
4451 //purpose  : for bug 20964 only
4452 //=======================================================================
4453 Standard_Boolean ApproxWithPCurves(const gp_Cylinder& theCyl,
4454                                    const gp_Sphere& theSph)
4455 {
4456   Standard_Boolean bRes = Standard_True;
4457   Standard_Real R1 = theCyl.Radius(), R2 = theSph.Radius();
4458   //
4459   {
4460     Standard_Real aD2, aRc2, aEps;
4461     gp_Pnt aApexSph;
4462     //
4463     aEps=1.E-7;
4464     aRc2=R1*R1;
4465     //
4466     const gp_Ax3& aAx3Sph=theSph.Position();
4467     const gp_Pnt& aLocSph=aAx3Sph.Location();
4468     const gp_Dir& aDirSph=aAx3Sph.Direction();
4469     //
4470     const gp_Ax1& aAx1Cyl=theCyl.Axis();
4471     gp_Lin aLinCyl(aAx1Cyl);
4472     //
4473     aApexSph.SetXYZ(aLocSph.XYZ()+R2*aDirSph.XYZ());
4474     aD2=aLinCyl.SquareDistance(aApexSph);
4475     if (fabs(aD2-aRc2)<aEps) {
4476       return !bRes;
4477     }
4478     //
4479     aApexSph.SetXYZ(aLocSph.XYZ()-R2*aDirSph.XYZ());
4480     aD2=aLinCyl.SquareDistance(aApexSph);
4481     if (fabs(aD2-aRc2)<aEps) {
4482       return !bRes;
4483     }
4484   }
4485   //
4486
4487   if(R1 < 2.*R2) {
4488     return bRes;
4489   }
4490   gp_Lin anCylAx(theCyl.Axis());
4491
4492   Standard_Real aDist = anCylAx.Distance(theSph.Location());
4493   Standard_Real aDRel = Abs(aDist - R1)/R2;
4494
4495   if(aDRel > .2) return bRes;
4496
4497   Standard_Real par = ElCLib::Parameter(anCylAx, theSph.Location());
4498   gp_Pnt aP = ElCLib::Value(par, anCylAx);
4499   gp_Vec aV(aP, theSph.Location());
4500
4501   Standard_Real dd = aV.Dot(theSph.Position().XDirection());
4502
4503   if(aDist < R1 && dd > 0.) return Standard_False;
4504   if(aDist > R1 && dd < 0.) return Standard_False;
4505
4506
4507   return bRes;
4508 }
4509 //=======================================================================
4510 //function : PerformPlanes
4511 //purpose  :
4512 //=======================================================================
4513 void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
4514                     const Handle(GeomAdaptor_HSurface)& theS2,
4515                     const Standard_Real TolAng,
4516                     const Standard_Real TolTang,
4517                     const Standard_Boolean theApprox1,
4518                     const Standard_Boolean theApprox2,
4519                     IntTools_SequenceOfCurves& theSeqOfCurve,
4520                     Standard_Boolean& theTangentFaces)
4521 {
4522
4523   gp_Pln aPln1 = theS1->Surface().Plane();
4524   gp_Pln aPln2 = theS2->Surface().Plane();
4525
4526   IntAna_QuadQuadGeo aPlnInter(aPln1, aPln2, TolAng, TolTang);
4527
4528   if(!aPlnInter.IsDone()) {
4529     theTangentFaces = Standard_False;
4530     return;
4531   }
4532
4533   IntAna_ResultType aResType = aPlnInter.TypeInter();
4534
4535   if(aResType == IntAna_Same) {
4536     theTangentFaces = Standard_True;
4537     return;
4538   }
4539
4540   theTangentFaces = Standard_False;
4541
4542   if(aResType == IntAna_Empty) {
4543     return;
4544   }
4545
4546   gp_Lin aLin = aPlnInter.Line(1);
4547
4548   ProjLib_Plane aProj;
4549
4550   aProj.Init(aPln1);
4551   aProj.Project(aLin);
4552   gp_Lin2d aLin2d1 = aProj.Line();
4553   //
4554   aProj.Init(aPln2);
4555   aProj.Project(aLin);
4556   gp_Lin2d aLin2d2 = aProj.Line();
4557   //
4558   //classify line2d1 relatively first plane
4559   Standard_Real P11, P12;
4560   Standard_Boolean IsCrossed = ClassifyLin2d(theS1, aLin2d1, TolTang, P11, P12);
4561   if(!IsCrossed) return;
4562   //classify line2d2 relatively second plane
4563   Standard_Real P21, P22;
4564   IsCrossed = ClassifyLin2d(theS2, aLin2d2, TolTang, P21, P22);
4565   if(!IsCrossed) return;
4566
4567   //Analysis of parametric intervals: must have common part
4568
4569   if(P21 >= P12) return;
4570   if(P22 <= P11) return;
4571
4572   Standard_Real pmin, pmax;
4573   pmin = Max(P11, P21);
4574   pmax = Min(P12, P22);
4575
4576   if(pmax - pmin <= TolTang) return;
4577
4578   Handle(Geom_Line) aGLin = new Geom_Line(aLin);
4579
4580   IntTools_Curve aCurve;
4581   Handle(Geom_TrimmedCurve) aGTLin = new Geom_TrimmedCurve(aGLin, pmin, pmax);
4582
4583   aCurve.SetCurve(aGTLin);
4584
4585   if(theApprox1) {
4586     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d1);
4587     aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
4588   }
4589   else {
4590     Handle(Geom2d_Curve) H1;
4591     aCurve.SetFirstCurve2d(H1);
4592   }
4593   if(theApprox2) {
4594     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d2);
4595     aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
4596   }
4597   else {
4598     Handle(Geom2d_Curve) H1;
4599     aCurve.SetFirstCurve2d(H1);
4600   }
4601
4602   theSeqOfCurve.Append(aCurve);
4603
4604 }
4605
4606 //=======================================================================
4607 //function : ClassifyLin2d
4608 //purpose  :
4609 //=======================================================================
4610 static inline Standard_Boolean INTER(const Standard_Real d1,
4611                                      const Standard_Real d2,
4612                                      const Standard_Real tol)
4613 {
4614   return (d1 > tol && d2 < -tol) ||
4615          (d1 < -tol && d2 > tol) ||
4616          ((d1 <= tol && d1 >= -tol) && (d2 > tol || d2 < -tol)) ||
4617          ((d2 <= tol && d2 >= -tol) && (d1 > tol || d1 < -tol));
4618 }
4619 static inline Standard_Boolean COINC(const Standard_Real d1,
4620                                      const Standard_Real d2,
4621                                      const Standard_Real tol)
4622 {
4623   return (d1 <= tol && d1 >= -tol) && (d2 <= tol && d2 >= -tol);
4624 }
4625 Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS,
4626                                const gp_Lin2d& theLin2d,
4627                                const Standard_Real theTol,
4628                                Standard_Real& theP1,
4629                                Standard_Real& theP2)
4630
4631 {
4632   Standard_Real xmin, xmax, ymin, ymax, d1, d2, A, B, C;
4633   Standard_Real par[2];
4634   Standard_Integer nbi = 0;
4635
4636   xmin = theS->Surface().FirstUParameter();
4637   xmax = theS->Surface().LastUParameter();
4638   ymin = theS->Surface().FirstVParameter();
4639   ymax = theS->Surface().LastVParameter();
4640
4641   theLin2d.Coefficients(A, B, C);
4642
4643   //xmin, ymin <-> xmin, ymax
4644   d1 = A*xmin + B*ymin + C;
4645   d2 = A*xmin + B*ymax + C;
4646
4647   if(INTER(d1, d2, theTol)) {
4648     //Intersection with boundary
4649     Standard_Real y = -(C + A*xmin)/B;
4650     par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, y));
4651     nbi++;
4652   }
4653   else if (COINC(d1, d2, theTol)) {
4654     //Coincidence with boundary
4655     par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
4656     par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
4657     nbi = 2;
4658   }
4659
4660   if(nbi == 2) {
4661
4662     if(fabs(par[0]-par[1]) > theTol) {
4663       theP1 = Min(par[0], par[1]);
4664       theP2 = Max(par[0], par[1]);
4665       return Standard_True;
4666     }
4667     else return Standard_False;
4668
4669   }
4670
4671   //xmin, ymax <-> xmax, ymax
4672   d1 = d2;
4673   d2 = A*xmax + B*ymax + C;
4674
4675   if(d1 > theTol || d1 < -theTol) {//to avoid checking of
4676                                    //coincidence with the same point
4677     if(INTER(d1, d2, theTol)) {
4678       Standard_Real x = -(C + B*ymax)/A;
4679       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymax));
4680       nbi++;
4681     }
4682     else if (COINC(d1, d2, theTol)) {
4683       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
4684       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
4685       nbi = 2;
4686     }
4687   }
4688
4689   if(nbi == 2) {
4690
4691     if(fabs(par[0]-par[1]) > theTol) {
4692       theP1 = Min(par[0], par[1]);
4693       theP2 = Max(par[0], par[1]);
4694       return Standard_True;
4695     }
4696     else return Standard_False;
4697
4698   }
4699
4700   //xmax, ymax <-> xmax, ymin
4701   d1 = d2;
4702   d2 = A*xmax + B*ymin + C;
4703
4704   if(d1 > theTol || d1 < -theTol) {
4705     if(INTER(d1, d2, theTol)) {
4706       Standard_Real y = -(C + A*xmax)/B;
4707       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, y));
4708       nbi++;
4709     }
4710     else if (COINC(d1, d2, theTol)) {
4711       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
4712       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
4713       nbi = 2;
4714     }
4715   }
4716
4717   if(nbi == 2) {
4718     if(fabs(par[0]-par[1]) > theTol) {
4719       theP1 = Min(par[0], par[1]);
4720       theP2 = Max(par[0], par[1]);
4721       return Standard_True;
4722     }
4723     else return Standard_False;
4724   }
4725
4726   //xmax, ymin <-> xmin, ymin
4727   d1 = d2;
4728   d2 = A*xmin + B*ymin + C;
4729
4730   if(d1 > theTol || d1 < -theTol) {
4731     if(INTER(d1, d2, theTol)) {
4732       Standard_Real x = -(C + B*ymin)/A;
4733       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymin));
4734       nbi++;
4735     }
4736     else if (COINC(d1, d2, theTol)) {
4737       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
4738       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
4739       nbi = 2;
4740     }
4741   }
4742
4743   if(nbi == 2) {
4744     if(fabs(par[0]-par[1]) > theTol) {
4745       theP1 = Min(par[0], par[1]);
4746       theP2 = Max(par[0], par[1]);
4747       return Standard_True;
4748     }
4749     else return Standard_False;
4750   }
4751
4752   return Standard_False;
4753
4754 }
4755 //
4756 //=======================================================================
4757 //function : ApproxParameters
4758 //purpose  :
4759 //=======================================================================
4760 void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
4761                       const Handle(GeomAdaptor_HSurface)& aHS2,
4762                       Standard_Integer& iDegMin,
4763                       Standard_Integer& iDegMax,
4764                       Standard_Integer& iNbIter)
4765
4766 {
4767   GeomAbs_SurfaceType aTS1, aTS2;
4768
4769   //
4770   iNbIter=0;
4771   iDegMin=4;
4772   iDegMax=8;
4773   //
4774   aTS1=aHS1->Surface().GetType();
4775   aTS2=aHS2->Surface().GetType();
4776   //
4777   // Cylinder/Torus
4778   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
4779       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
4780     Standard_Real aRC, aRT, dR, aPC;
4781     gp_Cylinder aCylinder;
4782     gp_Torus aTorus;
4783     //
4784     aPC=Precision::Confusion();
4785     //
4786     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Surface().Cylinder() : aHS2->Surface().Cylinder();
4787     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Surface().Torus() : aHS2->Surface().Torus();
4788     //
4789     aRC=aCylinder.Radius();
4790     aRT=aTorus.MinorRadius();
4791     dR=aRC-aRT;
4792     if (dR<0.) {
4793       dR=-dR;
4794     }
4795     //
4796     if (dR<aPC) {
4797       iDegMax=6;
4798     }
4799   }
4800   if (aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Cylinder) {
4801     iNbIter=1;
4802   }
4803 }
4804 //=======================================================================
4805 //function : Tolerances
4806 //purpose  :
4807 //=======================================================================
4808 void Tolerances(const Handle(GeomAdaptor_HSurface)& aHS1,
4809                 const Handle(GeomAdaptor_HSurface)& aHS2,
4810                 Standard_Real& aTolTang)
4811 {
4812   GeomAbs_SurfaceType aTS1, aTS2;
4813   //
4814   aTS1=aHS1->Surface().GetType();
4815   aTS2=aHS2->Surface().GetType();
4816   //
4817   // Cylinder/Torus
4818   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
4819       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
4820     Standard_Real aRC, aRT, dR, aPC;
4821     gp_Cylinder aCylinder;
4822     gp_Torus aTorus;
4823     //
4824     aPC=Precision::Confusion();
4825     //
4826     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Surface().Cylinder() : aHS2->Surface().Cylinder();
4827     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Surface().Torus() : aHS2->Surface().Torus();
4828     //
4829     aRC=aCylinder.Radius();
4830     aRT=aTorus.MinorRadius();
4831     dR=aRC-aRT;
4832     if (dR<0.) {
4833       dR=-dR;
4834     }
4835     //
4836     if (dR<aPC) {
4837       aTolTang=0.1*aTolTang;
4838     }
4839   }
4840 }
4841 //=======================================================================
4842 //function : SortTypes
4843 //purpose  :
4844 //=======================================================================
4845 Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
4846                            const GeomAbs_SurfaceType aType2)
4847 {
4848   Standard_Boolean bRet;
4849   Standard_Integer aI1, aI2;
4850   //
4851   bRet=Standard_False;
4852   //
4853   aI1=IndexType(aType1);
4854   aI2=IndexType(aType2);
4855   if (aI1<aI2){
4856     bRet=!bRet;
4857   }
4858   return bRet;
4859 }
4860 //=======================================================================
4861 //function : IndexType
4862 //purpose  :
4863 //=======================================================================
4864 Standard_Integer IndexType(const GeomAbs_SurfaceType aType)
4865 {
4866   Standard_Integer aIndex;
4867   //
4868   aIndex=11;
4869   //
4870   if (aType==GeomAbs_Plane) {
4871     aIndex=0;
4872   }
4873   else if (aType==GeomAbs_Cylinder) {
4874     aIndex=1;
4875   }
4876   else if (aType==GeomAbs_Cone) {
4877     aIndex=2;
4878   }
4879   else if (aType==GeomAbs_Sphere) {
4880     aIndex=3;
4881   }
4882   else if (aType==GeomAbs_Torus) {
4883     aIndex=4;
4884   }
4885   else if (aType==GeomAbs_BezierSurface) {
4886     aIndex=5;
4887   }
4888   else if (aType==GeomAbs_BSplineSurface) {
4889     aIndex=6;
4890   }
4891   else if (aType==GeomAbs_SurfaceOfRevolution) {
4892     aIndex=7;
4893   }
4894   else if (aType==GeomAbs_SurfaceOfExtrusion) {
4895     aIndex=8;
4896   }
4897   else if (aType==GeomAbs_OffsetSurface) {
4898     aIndex=9;
4899   }
4900   else if (aType==GeomAbs_OtherSurface) {
4901     aIndex=10;
4902   }
4903   return aIndex;
4904 }
4905 #ifdef OCCT_DEBUG_DUMPWLINE
4906 //=======================================================================
4907 //function : DumpWLine
4908 //purpose  :
4909 //=======================================================================
4910 void DumpWLine(const Handle(IntPatch_WLine)& aWLine)
4911 {
4912   Standard_Integer i, aNbPnts;
4913   Standard_Real aX, aY, aZ, aU1, aV1, aU2, aV2;
4914   //
4915   printf(" *WLine\n");
4916   aNbPnts=aWLine->NbPnts();
4917   for (i=1; i<=aNbPnts; ++i) {
4918     const IntSurf_PntOn2S aPntOn2S=aWLine->Point(i);
4919     const gp_Pnt& aP3D=aPntOn2S.Value();
4920     aP3D.Coord(aX, aY, aZ);
4921     aPntOn2S.Parameters(aU1, aV1, aU2, aV2);
4922     //
4923     printf("point p_%d %lf %lf %lf\n", i, aX, aY, aZ);
4924     //printf("point p_%d %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf\n",
4925         //   i, aX, aY, aZ, aU1, aV1, aU2, aV2);
4926   }
4927 }
4928 #endif
4929 //=======================================================================
4930 //function : RefineVector
4931 //purpose  :
4932 //=======================================================================
4933 void RefineVector(gp_Vec2d& aV2D)
4934 {
4935   Standard_Integer k,m;
4936   Standard_Real aC[2], aEps, aR1, aR2, aNum;
4937   //
4938   aEps=RealEpsilon();
4939   aR1=1.-aEps;
4940   aR2=1.+aEps;
4941   //
4942   aV2D.Coord(aC[0], aC[1]);
4943   //
4944   for (k=0; k<2; ++k) {
4945     m=(k+1)%2;
4946     aNum=fabs(aC[k]);
4947     if (aNum>aR1 && aNum<aR2) {
4948       if (aC[k]<0.) {
4949         aC[k]=-1.;
4950       }
4951       else {
4952         aC[k]=1.;
4953       }
4954       aC[m]=0.;
4955       break;
4956     }
4957   }
4958   aV2D.SetCoord(aC[0], aC[1]);
4959 }
4960
4961 //=======================================================================
4962 // Function : FindMaxDistance
4963 // purpose :
4964 //=======================================================================
4965 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
4966                               const Standard_Real theFirst,
4967                               const Standard_Real theLast,
4968                               const TopoDS_Face& theFace,
4969                               const Handle(IntTools_Context)& theContext)
4970 {
4971   Standard_Integer aNbS;
4972   Standard_Real aT1, aT2, aDt, aD, aDMax, anEps;
4973   //
4974   aNbS = 11;
4975   aDt = (theLast - theFirst) / aNbS;
4976   aDMax = 0.;
4977   anEps = 1.e-4 * aDt;
4978   //
4979   GeomAPI_ProjectPointOnSurf& aProjPS = theContext->ProjPS(theFace);
4980   aT2 = theFirst;
4981   for (;;) {
4982     aT1 = aT2;
4983     aT2 += aDt;
4984     //
4985     if (aT2 > theLast) {
4986       break;
4987     }
4988     //
4989     aD = FindMaxDistance(theCurve, aT1, aT2, aProjPS, anEps);
4990     if (aD > aDMax) {
4991       aDMax = aD;
4992     }
4993   }
4994   //
4995   return aDMax;
4996 }
4997
4998 //=======================================================================
4999 // Function : FindMaxDistance
5000 // purpose :
5001 //=======================================================================
5002 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
5003                               const Standard_Real theFirst,
5004                               const Standard_Real theLast,
5005                               GeomAPI_ProjectPointOnSurf& theProjPS,
5006                               const Standard_Real theEps)
5007 {
5008   Standard_Real aA, aB, aCf, aX, aX1, aX2, aF1, aF2, aF;
5009   //
5010   aCf = 0.61803398874989484820458683436564;//(sqrt(5.)-1)/2.;
5011   aA = theFirst;
5012   aB = theLast;
5013   //
5014   aX1 = aB - aCf * (aB - aA);
5015   aF1 = MaxDistance(theC, aX1, theProjPS);
5016   aX2 = aA + aCf * (aB - aA);
5017   aF2 = MaxDistance(theC, aX2, theProjPS);
5018   //
5019   for (;;) {
5020     if ((aB - aA) < theEps) {
5021       break;
5022     }
5023     //
5024     if (aF1 > aF2) {
5025       aB = aX2;
5026       aX2 = aX1;
5027       aF2 = aF1;
5028       aX1 = aB - aCf * (aB - aA);
5029       aF1 = MaxDistance(theC, aX1, theProjPS);
5030     }
5031     else {
5032       aA = aX1;
5033       aX1 = aX2;
5034       aF1 = aF2;
5035       aX2 = aA + aCf * (aB - aA);
5036       aF2 = MaxDistance(theC, aX2, theProjPS);
5037     }
5038   }
5039   //
5040   aX = 0.5 * (aA + aB);
5041   aF = MaxDistance(theC, aX, theProjPS);
5042   //
5043   if (aF1 > aF) {
5044     aF = aF1;
5045   }
5046   //
5047   if (aF2 > aF) {
5048     aF = aF2;
5049   }
5050   //
5051   return aF;
5052 }
5053
5054 //=======================================================================
5055 // Function : MaxDistance
5056 // purpose :
5057 //=======================================================================
5058 Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
5059                           const Standard_Real aT,
5060                           GeomAPI_ProjectPointOnSurf& theProjPS)
5061 {
5062   Standard_Real aD;
5063   gp_Pnt aP;
5064   //
5065   theC->D0(aT, aP);
5066   theProjPS.Perform(aP);
5067   aD = theProjPS.NbPoints() ? theProjPS.LowerDistance() : 0.;
5068   //
5069   return aD;
5070 }
5071
5072 //=======================================================================
5073 //function : CheckPCurve
5074 //purpose  : Checks if points of the pcurve are out of the face bounds.
5075 //=======================================================================
5076   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
5077                                const TopoDS_Face& aFace)
5078 {
5079   const Standard_Integer NPoints = 23;
5080   Standard_Integer i;
5081   Standard_Real umin,umax,vmin,vmax;
5082
5083   BRepTools::UVBounds(aFace, umin, umax, vmin, vmax);
5084   Standard_Real tolU = Max ((umax-umin)*0.01, Precision::Confusion());
5085   Standard_Real tolV = Max ((vmax-vmin)*0.01, Precision::Confusion());
5086   Standard_Real fp = aPC->FirstParameter();
5087   Standard_Real lp = aPC->LastParameter();
5088
5089
5090   // adjust domain for periodic surfaces
5091   TopLoc_Location aLoc;
5092   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace, aLoc);
5093   if (aSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
5094     aSurf = (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurf))->BasisSurface();
5095   }
5096   gp_Pnt2d pnt = aPC->Value((fp+lp)/2);
5097   Standard_Real u,v;
5098   pnt.Coord(u,v);
5099   //
5100   if (aSurf->IsUPeriodic()) {
5101     Standard_Real aPer = aSurf->UPeriod();
5102     Standard_Integer nshift = (Standard_Integer) ((u-umin)/aPer);
5103     if (u < umin+aPer*nshift) nshift--;
5104     umin += aPer*nshift;
5105     umax += aPer*nshift;
5106   }
5107   if (aSurf->IsVPeriodic()) {
5108     Standard_Real aPer = aSurf->VPeriod();
5109     Standard_Integer nshift = (Standard_Integer) ((v-vmin)/aPer);
5110     if (v < vmin+aPer*nshift) nshift--;
5111     vmin += aPer*nshift;
5112     vmax += aPer*nshift;
5113   }
5114   //
5115   //--------------------------------------------------------
5116   Standard_Boolean bRet;
5117   Standard_Integer j, aNbIntervals;
5118   Standard_Real aT, dT;
5119   gp_Pnt2d aP2D;
5120   //
5121   Geom2dAdaptor_Curve aGAC(aPC);
5122   aNbIntervals=aGAC.NbIntervals(GeomAbs_CN);
5123   //
5124   TColStd_Array1OfReal aTI(1, aNbIntervals+1);
5125   aGAC.Intervals(aTI,GeomAbs_CN);
5126   //
5127   bRet=Standard_False;
5128   //
5129   aT=aGAC.FirstParameter();
5130   for (j=1; j<=aNbIntervals; ++j) {
5131     dT=(aTI(j+1)-aTI(j))/NPoints;
5132     //
5133     for (i=1; i<NPoints; i++) {
5134       aT=aT+dT;
5135       aGAC.D0(aT, aP2D);
5136       aP2D.Coord(u,v);
5137     if (umin-u > tolU || u-umax > tolU ||
5138           vmin-v > tolV || v-vmax > tolV) {
5139         return bRet;
5140   }
5141 }
5142   }
5143   return !bRet;
5144 }
5145 //=======================================================================
5146 //function : CorrectPlaneBoundaries
5147 //purpose  :
5148 //=======================================================================
5149  void CorrectPlaneBoundaries(Standard_Real& aUmin,
5150                              Standard_Real& aUmax,
5151                              Standard_Real& aVmin,
5152                              Standard_Real& aVmax)
5153 {
5154   if (!(Precision::IsInfinite(aUmin) ||
5155         Precision::IsInfinite(aUmax))) {
5156     Standard_Real dU;
5157     //
5158     dU=0.1*(aUmax-aUmin);
5159     aUmin=aUmin-dU;
5160     aUmax=aUmax+dU;
5161   }
5162   if (!(Precision::IsInfinite(aVmin) ||
5163         Precision::IsInfinite(aVmax))) {
5164     Standard_Real dV;
5165     //
5166     dV=0.1*(aVmax-aVmin);
5167     aVmin=aVmin-dV;
5168     aVmax=aVmax+dV;
5169   }
5170 }