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
 167 static
 168   Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
 169                                             const Standard_Integer                         ideb,
 170                                             const Standard_Integer                         ifin,
 171                                             const Standard_Boolean                         onFirst);
 172
 173 static
 174   Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
 175                                         const Handle(GeomAdaptor_HSurface)&            theSurface1,
 176                                         const Handle(GeomAdaptor_HSurface)&            theSurface2,
 177                                         const TopoDS_Face&                             theFace1,
 178                                         const TopoDS_Face&                             theFace2,
 179                                         const GeomInt_LineConstructor&                 theLConstructor,
 180                                         const Standard_Boolean                         theAvoidLConstructor,
 181                                         IntPatch_SequenceOfLine&                       theNewLines,
 182                                         Standard_Real&                                 theReachedTol3d,
 183                                         const Handle(IntTools_Context)& );
 184
 185 static
 186   Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
 187                                           const Handle(Geom_Curve)& theCurve,
 188                                           const TopoDS_Face&        theFace1,
 189                                           const TopoDS_Face&        theFace2,
 190                                           const Standard_Real       theOtherParameter,
 191                                           const Standard_Boolean    bIncreasePar,
 192                                           Standard_Real&            theNewParameter,
 193                                           const Handle(IntTools_Context)& );
 194
 195 static
 196   Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve);
 197
 198 static
 199   Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
 200                                      const Standard_Real theFirstBoundary,
 201                                      const Standard_Real theSecondBoundary,
 202                                      const Standard_Real theResolution,
 203                                      Standard_Boolean&   IsOnFirstBoundary);
 204 static
 205   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
 206                              const gp_Pnt2d&     theLastPoint,
 207                              const Standard_Real theUmin,
 208                              const Standard_Real theUmax,
 209                              const Standard_Real theVmin,
 210                              const Standard_Real theVmax,
 211                              gp_Pnt2d&           theNewPoint);
 212
 213
 214 static
 215   Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
 216                                        const Handle(GeomAdaptor_HSurface)&  theSurface2,
 217                                        const TopoDS_Face&                   theFace1,
 218                                        const TopoDS_Face&                   theFace2,
 219                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
 220                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
 221                                        Handle(TColStd_HArray1OfReal)&       theResultRadius,
 222                                        const Handle(IntTools_Context)& );
 223
 224 static
 225   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
 226                              const gp_Pnt2d&     theLastPoint,
 227                              const Standard_Real theUmin,
 228                              const Standard_Real theUmax,
 229                              const Standard_Real theVmin,
 230                              const Standard_Real theVmax,
 231                              const gp_Pnt2d&     theTanZoneCenter,
 232                              const Standard_Real theZoneRadius,
 233                              Handle(GeomAdaptor_HSurface) theGASurface,
 234                              gp_Pnt2d&           theNewPoint);
 235
 236 static
 237   Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
 238                                    const gp_Pnt2d&     theTanZoneCenter,
 239                                    const Standard_Real theZoneRadius,
 240                                    Handle(GeomAdaptor_HSurface) theGASurface);
 241
 242 static
 243   gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
 244                              const gp_Pnt2d&     theOriginalPoint,
 245                              Handle(GeomAdaptor_HSurface) theGASurface);
 246 static
 247   Standard_Boolean  ApproxWithPCurves(const gp_Cylinder& theCyl,
 248                                       const gp_Sphere& theSph);
 249
 250 static void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
 251                            const Handle(GeomAdaptor_HSurface)& theS2,
 252                            const Standard_Real TolAng,
 253                            const Standard_Real TolTang,
 254                            const Standard_Boolean theApprox1,
 255                            const Standard_Boolean theApprox2,
 256                            IntTools_SequenceOfCurves& theSeqOfCurve,
 257                            Standard_Boolean& theTangentFaces);
 258
 259 static Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS,
 260                                       const gp_Lin2d& theLin2d,
 261                                       const Standard_Real theTol,
 262                                       Standard_Real& theP1,
 263                                       Standard_Real& theP2);
 264 //
 265 static
 266   void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
 267                         const Handle(GeomAdaptor_HSurface)& aHS2,
 268                         Standard_Integer& iDegMin,
 269                         Standard_Integer& iNbIter,
 270                         Standard_Integer& iDegMax);
 271
 272 static
 273   void Tolerances(const Handle(GeomAdaptor_HSurface)& aHS1,
 274                   const Handle(GeomAdaptor_HSurface)& aHS2,
 275                   Standard_Real& aTolTang);
 276
 277 static
 278   Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
 279                              const GeomAbs_SurfaceType aType2);
 280 static
 281   Standard_Integer IndexType(const GeomAbs_SurfaceType aType);
 282
 283 //
 284 static
 285   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
 286                                const TopoDS_Face& aFace);
 287
 288 static
 289   Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
 290                             const Standard_Real aT,
 291                             GeomAPI_ProjectPointOnSurf& theProjPS);
 292
 293 static
 294   Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
 295                                 const Standard_Real theFirst,
 296                                 const Standard_Real theLast,
 297                                 GeomAPI_ProjectPointOnSurf& theProjPS,
 298                                 const Standard_Real theEps);
 299
 300 static
 301   Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
 302                                 const Standard_Real theFirst,
 303                                 const Standard_Real theLast,
 304                                 const TopoDS_Face& theFace,
 305                                 const Handle(IntTools_Context)& theContext);
 306
 307 static
 308   void CorrectPlaneBoundaries(Standard_Real& aUmin,
 309                               Standard_Real& aUmax,
 310                               Standard_Real& aVmin,
 311                               Standard_Real& aVmax);
 312
 313 //=======================================================================
 314 //function :
 315 //purpose  :
 316 //=======================================================================
 317 IntTools_FaceFace::IntTools_FaceFace()
 318 {
 319   myIsDone=Standard_False;
 320   myTangentFaces=Standard_False;
 321   //
 322   myHS1 = new GeomAdaptor_HSurface ();
 323   myHS2 = new GeomAdaptor_HSurface ();
 324   myTolReached2d=0.;
 325   myTolReached3d=0.;
 326   SetParameters(Standard_True, Standard_True, Standard_True, 1.e-07);
 327
 328 }
 329 //=======================================================================
 330 //function : SetContext
 331 //purpose  :
 332 //=======================================================================
 333 void IntTools_FaceFace::SetContext(const Handle(IntTools_Context)& aContext)
 334 {
 335   myContext=aContext;
 336 }
 337 //=======================================================================
 338 //function : Context
 339 //purpose  :
 340 //=======================================================================
 341 const Handle(IntTools_Context)& IntTools_FaceFace::Context()const
 342 {
 343   return myContext;
 344 }
 345 //=======================================================================
 346 //function : Face1
 347 //purpose  :
 348 //=======================================================================
 349 const TopoDS_Face& IntTools_FaceFace::Face1() const
 350 {
 351   return myFace1;
 352 }
 353 //=======================================================================
 354 //function : Face2
 355 //purpose  :
 356 //=======================================================================
 357 const TopoDS_Face& IntTools_FaceFace::Face2() const
 358 {
 359   return myFace2;
 360 }
 361 //=======================================================================
 362 //function : TangentFaces
 363 //purpose  :
 364 //=======================================================================
 365 Standard_Boolean IntTools_FaceFace::TangentFaces() const
 366 {
 367   return myTangentFaces;
 368 }
 369 //=======================================================================
 370 //function : Points
 371 //purpose  :
 372 //=======================================================================
 373 const IntTools_SequenceOfPntOn2Faces& IntTools_FaceFace::Points() const
 374 {
 375   return myPnts;
 376 }
 377 //=======================================================================
 378 //function : IsDone
 379 //purpose  :
 380 //=======================================================================
 381 Standard_Boolean IntTools_FaceFace::IsDone() const
 382 {
 383   return myIsDone;
 384 }
 385 //=======================================================================
 386 //function : TolReached3d
 387 //purpose  :
 388 //=======================================================================
 389 Standard_Real IntTools_FaceFace::TolReached3d() const
 390 {
 391   return myTolReached3d;
 392 }
 393 //=======================================================================
 394 //function : Lines
 395 //purpose  : return lines of intersection
 396 //=======================================================================
 397 const IntTools_SequenceOfCurves& IntTools_FaceFace::Lines() const
 398 {
 399   StdFail_NotDone_Raise_if
 400     (!myIsDone,
 401      "IntTools_FaceFace::Lines() => myIntersector NOT DONE");
 402   return mySeqOfCurve;
 403 }
 404 //=======================================================================
 405 //function : TolReached2d
 406 //purpose  :
 407 //=======================================================================
 408 Standard_Real IntTools_FaceFace::TolReached2d() const
 409 {
 410   return myTolReached2d;
 411 }
 412 // =======================================================================
 413 // function: SetParameters
 414 //
 415 // =======================================================================
 416 void IntTools_FaceFace::SetParameters(const Standard_Boolean ToApproxC3d,
 417                                       const Standard_Boolean ToApproxC2dOnS1,
 418                                       const Standard_Boolean ToApproxC2dOnS2,
 419                                       const Standard_Real ApproximationTolerance)
 420 {
 421   myApprox = ToApproxC3d;
 422   myApprox1 = ToApproxC2dOnS1;
 423   myApprox2 = ToApproxC2dOnS2;
 424   myTolApprox = ApproximationTolerance;
 425 }
 426 //=======================================================================
 427 //function : SetList
 428 //purpose  :
 429 //=======================================================================
 430 void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
 431 {
 432   myListOfPnts = aListOfPnts;
 433 }
 434
 435
 436 static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
 437                                         const TopoDS_Face& theF2)
 438 {
 439   const Standard_Real Tolang = 1.e-8;
 440   const Standard_Real aTolF1=BRep_Tool::Tolerance(theF1);
 441   const Standard_Real aTolF2=BRep_Tool::Tolerance(theF2);
 442   const Standard_Real aTolSum = aTolF1 + aTolF2;
 443   Standard_Real aHigh = 0.0;
 444
 445   const BRepAdaptor_Surface aBAS1(theF1), aBAS2(theF2);
 446   const GeomAbs_SurfaceType aType1=aBAS1.GetType();
 447   const GeomAbs_SurfaceType aType2=aBAS2.GetType();
 448
 449   gp_Pln aS1;
 450   gp_Cylinder aS2;
 451   if(aType1 == GeomAbs_Plane)
 452   {
 453     aS1=aBAS1.Plane();
 454   }
 455   else if(aType2 == GeomAbs_Plane)
 456   {
 457     aS1=aBAS2.Plane();
 458   }
 459   else
 460   {
 461     return Standard_True;
 462   }
 463
 464   if(aType1 == GeomAbs_Cylinder)
 465   {
 466     aS2=aBAS1.Cylinder();
 467     const Standard_Real VMin = aBAS1.FirstVParameter();
 468     const Standard_Real VMax = aBAS1.LastVParameter();
 469
 470     if( Precision::IsNegativeInfinite(VMin) ||
 471         Precision::IsPositiveInfinite(VMax))
 472           return Standard_True;
 473     else
 474       aHigh = VMax - VMin;
 475   }
 476   else if(aType2 == GeomAbs_Cylinder)
 477   {
 478     aS2=aBAS2.Cylinder();
 479
 480     const Standard_Real VMin = aBAS2.FirstVParameter();
 481     const Standard_Real VMax = aBAS2.LastVParameter();
 482
 483     if( Precision::IsNegativeInfinite(VMin) ||
 484         Precision::IsPositiveInfinite(VMax))
 485           return Standard_True;
 486     else
 487       aHigh = VMax - VMin;
 488   }
 489   else
 490   {
 491     return Standard_True;
 492   }
 493
 494   IntAna_QuadQuadGeo inter;
 495   inter.Perform(aS1,aS2,Tolang,aTolSum, aHigh);
 496   if(inter.TypeInter() == IntAna_Ellipse)
 497   {
 498     const gp_Elips anEl = inter.Ellipse(1);
 499     const Standard_Real aMajorR = anEl.MajorRadius();
 500     const Standard_Real aMinorR = anEl.MinorRadius();
 501
 502     return (aMajorR < 100000.0 * aMinorR);
 503   }
 504   else
 505   {
 506     return inter.IsDone();
 507   }
 508 }
 509 //=======================================================================
 510 //function : Perform
 511 //purpose  : intersect surfaces of the faces
 512 //=======================================================================
 513 void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
 514                                 const TopoDS_Face& aF2)
 515 {
 516   Standard_Boolean RestrictLine = Standard_False, hasCone = Standard_False;
 517
 518   if (myContext.IsNull()) {
 519     myContext=new IntTools_Context;
 520   }
 521
 522   mySeqOfCurve.Clear();
 523   myTolReached2d=0.;
 524   myTolReached3d=0.;
 525   myIsDone = Standard_False;
 526   myNbrestr=0;//?
 527
 528   myFace1=aF1;
 529   myFace2=aF2;
 530
 531   const BRepAdaptor_Surface aBAS1(myFace1, Standard_False);
 532   const BRepAdaptor_Surface aBAS2(myFace2, Standard_False);
 533   GeomAbs_SurfaceType aType1=aBAS1.GetType();
 534   GeomAbs_SurfaceType aType2=aBAS2.GetType();
 535
 536   const Standard_Boolean bReverse=SortTypes(aType1, aType2);
 537   if (bReverse)
 538   {
 539     myFace1=aF2;
 540     myFace2=aF1;
 541     aType1=aBAS2.GetType();
 542     aType2=aBAS1.GetType();
 543
 544     if (myListOfPnts.Extent())
 545     {
 546       Standard_Real aU1,aV1,aU2,aV2;
 547       IntSurf_ListIteratorOfListOfPntOn2S aItP2S;
 548       //
 549       aItP2S.Initialize(myListOfPnts);
 550       for (; aItP2S.More(); aItP2S.Next())
 551       {
 552         IntSurf_PntOn2S& aP2S=aItP2S.Value();
 553         aP2S.Parameters(aU1,aV1,aU2,aV2);
 554         aP2S.SetValue(aU2,aV2,aU1,aV1);
 555       }
 556     }
 557     //
 558     Standard_Boolean anAproxTmp = myApprox1;
 559     myApprox1 = myApprox2;
 560     myApprox2 = anAproxTmp;
 561   }
 562
 563
 564   const Handle(Geom_Surface) S1=BRep_Tool::Surface(myFace1);
 565   const Handle(Geom_Surface) S2=BRep_Tool::Surface(myFace2);
 566
 567   const Standard_Real aTolF1=BRep_Tool::Tolerance(myFace1);
 568   const Standard_Real aTolF2=BRep_Tool::Tolerance(myFace2);
 569
 570   Standard_Real TolArc = aTolF1 + aTolF2;
 571   Standard_Real TolTang = TolArc;
 572
 573   const Standard_Boolean isFace1Quad = (aType1 == GeomAbs_Cylinder ||
 574                                         aType1 == GeomAbs_Cone ||
 575                                         aType1 == GeomAbs_Torus);
 576
 577   const Standard_Boolean isFace2Quad = (aType2 == GeomAbs_Cylinder ||
 578                                         aType2 == GeomAbs_Cone ||
 579                                         aType2 == GeomAbs_Torus);
 580
 581   if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane)  {
 582     Standard_Real umin, umax, vmin, vmax;
 583     //
 584     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 585     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 586     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 587     //
 588     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 589     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 590     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 591     //
 592     Standard_Real TolAng = 1.e-8;
 593     //
 594     PerformPlanes(myHS1, myHS2, TolAng, TolTang, myApprox1, myApprox2,
 595                   mySeqOfCurve, myTangentFaces);
 596     //
 597     myIsDone = Standard_True;
 598
 599     if(!myTangentFaces) {
 600       const Standard_Integer NbLinPP = mySeqOfCurve.Length();
 601       if(NbLinPP) {
 602         Standard_Real aTolFMax;
 603         myTolReached3d = 1.e-7;
 604         aTolFMax=Max(aTolF1, aTolF2);
 605         if (aTolFMax>myTolReached3d) {
 606           myTolReached3d=aTolFMax;
 607         }
 608         //
 609         myTolReached2d = myTolReached3d;
 610
 611         if (bReverse) {
 612           Handle(Geom2d_Curve) aC2D1, aC2D2;
 613           const Standard_Integer aNbLin = mySeqOfCurve.Length();
 614           for (Standard_Integer i = 1; i <= aNbLin; ++i) {
 615             IntTools_Curve& aIC=mySeqOfCurve(i);
 616             aC2D1=aIC.FirstCurve2d();
 617             aC2D2=aIC.SecondCurve2d();
 618             aIC.SetFirstCurve2d(aC2D2);
 619             aIC.SetSecondCurve2d(aC2D1);
 620           }
 621         }
 622       }
 623     }
 624     return;
 625   }//if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane){
 626
 627   if ((aType1==GeomAbs_Plane) && isFace2Quad)
 628   {
 629     Standard_Real umin, umax, vmin, vmax;
 630     // F1
 631     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 632     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 633     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 634     // F2
 635     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 636     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 637     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 638     //
 639     if( aType2==GeomAbs_Cone ) {
 640       TolArc = 0.0001;
 641       hasCone = Standard_True;
 642     }
 643   }
 644   else if ((aType2==GeomAbs_Plane) && isFace1Quad)
 645   {
 646     Standard_Real umin, umax, vmin, vmax;
 647     //F1
 648     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 649     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 650     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 651     // F2
 652     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 653     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 654     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 655     //
 656     if( aType1==GeomAbs_Cone ) {
 657       TolArc = 0.0001;
 658       hasCone = Standard_True;
 659     }
 660   }
 661   else
 662   {
 663     Standard_Real umin, umax, vmin, vmax;
 664     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 665     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 666     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 667     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 668     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 669     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 670   }
 671
 672   const Handle(IntTools_TopolTool) dom1 = new IntTools_TopolTool(myHS1);
 673   const Handle(IntTools_TopolTool) dom2 = new IntTools_TopolTool(myHS2);
 674
 675   myLConstruct.Load(dom1, dom2, myHS1, myHS2);
 676
 677
 678   Tolerances(myHS1, myHS2, TolTang);
 679
 680   {
 681     const Standard_Real UVMaxStep = 0.001;
 682     const Standard_Real Deflection = (hasCone) ? 0.085 : 0.1;
 683     myIntersector.SetTolerances(TolArc, TolTang, UVMaxStep, Deflection);
 684   }
 685
 686   if((myHS1->IsUClosed() && !myHS1->IsUPeriodic()) ||
 687      (myHS1->IsVClosed() && !myHS1->IsVPeriodic()) ||
 688      (myHS2->IsUClosed() && !myHS2->IsUPeriodic()) ||
 689      (myHS2->IsVClosed() && !myHS2->IsVPeriodic()))
 690   {
 691     RestrictLine = Standard_True;
 692   }
 693   //
 694   if((aType1 != GeomAbs_BSplineSurface) &&
 695       (aType1 != GeomAbs_BezierSurface)  &&
 696      (aType1 != GeomAbs_OtherSurface)  &&
 697      (aType2 != GeomAbs_BSplineSurface) &&
 698       (aType2 != GeomAbs_BezierSurface)  &&
 699      (aType2 != GeomAbs_OtherSurface))
 700   {
 701     RestrictLine = Standard_True;
 702
 703     if ((aType1 == GeomAbs_Torus) ||
 704         (aType2 == GeomAbs_Torus))
 705     {
 706       myListOfPnts.Clear();
 707     }
 708   }
 709
 710   //
 711   if(!RestrictLine)
 712   {
 713     TopExp_Explorer aExp;
 714     for(Standard_Integer i = 0; (!RestrictLine) && (i < 2); i++)
 715     {
 716       const TopoDS_Face& aF=(!i) ? myFace1 : myFace2;
 717       aExp.Init(aF, TopAbs_EDGE);
 718       for(; aExp.More(); aExp.Next())
 719       {
 720         const TopoDS_Edge& aE=TopoDS::Edge(aExp.Current());
 721
 722         if(BRep_Tool::Degenerated(aE))
 723         {
 724           RestrictLine = Standard_True;
 725           break;
 726         }
 727       }
 728     }
 729   }
 730
 731   const Standard_Boolean isGeomInt = isTreatAnalityc(aF1, aF2);
 732   myIntersector.Perform(myHS1, dom1, myHS2, dom2, TolArc, TolTang,
 733                                   myListOfPnts, RestrictLine, isGeomInt);
 734
 735   myIsDone = myIntersector.IsDone();
 736
 737   if (myIsDone)
 738   {
 739     myTangentFaces=myIntersector.TangentFaces();
 740     if (myTangentFaces) {
 741       return;
 742     }
 743     //
 744     if(RestrictLine) {
 745       myListOfPnts.Clear(); // to use LineConstructor
 746     }
 747     //
 748     const Standard_Integer aNbLin = myIntersector.NbLines();
 749     for (Standard_Integer i=1; i <= aNbLin; ++i) {
 750       MakeCurve(i, dom1, dom2);
 751     }
 752     //
 753     ComputeTolReached3d();
 754     //
 755     if (bReverse) {
 756       Handle(Geom2d_Curve) aC2D1, aC2D2;
 757       //
 758       const Standard_Integer aNbLin=mySeqOfCurve.Length();
 759       for (Standard_Integer i=1; i<=aNbLin; ++i)
 760       {
 761         IntTools_Curve& aIC=mySeqOfCurve(i);
 762         aC2D1=aIC.FirstCurve2d();
 763         aC2D2=aIC.SecondCurve2d();
 764         aIC.SetFirstCurve2d(aC2D2);
 765         aIC.SetSecondCurve2d(aC2D1);
 766       }
 767     }
 768
 769     // Points
 770     Standard_Boolean bValid2D1, bValid2D2;
 771     Standard_Real U1,V1,U2,V2;
 772     IntTools_PntOnFace aPntOnF1, aPntOnF2;
 773     IntTools_PntOn2Faces aPntOn2Faces;
 774     //
 775     const Standard_Integer aNbPnts = myIntersector.NbPnts();
 776     for (Standard_Integer i=1; i <= aNbPnts; ++i)
 777     {
 778       const IntSurf_PntOn2S& aISPnt=myIntersector.Point(i).PntOn2S();
 779       const gp_Pnt& aPnt=aISPnt.Value();
 780       aISPnt.Parameters(U1,V1,U2,V2);
 781       //
 782       // check the validity of the intersection point for the faces
 783       bValid2D1 = myContext->IsPointInOnFace(myFace1, gp_Pnt2d(U1, V1));
 784       if (!bValid2D1) {
 785         continue;
 786       }
 787       //
 788       bValid2D2 = myContext->IsPointInOnFace(myFace2, gp_Pnt2d(U2, V2));
 789       if (!bValid2D2) {
 790         continue;
 791       }
 792       //
 793       // add the intersection point
 794       aPntOnF1.Init(myFace1, aPnt, U1, V1);
 795       aPntOnF2.Init(myFace2, aPnt, U2, V2);
 796       //
 797       if (!bReverse)
 798       {
 799         aPntOn2Faces.SetP1(aPntOnF1);
 800         aPntOn2Faces.SetP2(aPntOnF2);
 801       }
 802       else
 803       {
 804         aPntOn2Faces.SetP2(aPntOnF1);
 805         aPntOn2Faces.SetP1(aPntOnF2);
 806       }
 807
 808       myPnts.Append(aPntOn2Faces);
 809     }
 810   }
 811 }
 812
 813 //=======================================================================
 814 //function : ComputeTolerance
 815 //purpose  :
 816 //=======================================================================
 817 Standard_Real IntTools_FaceFace::ComputeTolerance()
 818 {
 819   Standard_Integer i, j, aNbLin;
 820   Standard_Real aFirst, aLast, aD, aDMax, aT;
 821   Handle(Geom_Surface) aS1, aS2;
 822   //
 823   aDMax = 0;
 824   aNbLin = mySeqOfCurve.Length();
 825   //
 826   aS1 = myHS1->ChangeSurface().Surface();
 827   aS2 = myHS2->ChangeSurface().Surface();
 828   //
 829   for (i = 1; i <= aNbLin; ++i)
 830   {
 831     const IntTools_Curve& aIC = mySeqOfCurve(i);
 832     const Handle(Geom_Curve)& aC3D = aIC.Curve();
 833     if (aC3D.IsNull())
 834     {
 835       continue;
 836     }
 837     //
 838     aFirst = aC3D->FirstParameter();
 839     aLast  = aC3D->LastParameter();
 840     //
 841     const Handle(Geom2d_Curve)& aC2D1 = aIC.FirstCurve2d();
 842     const Handle(Geom2d_Curve)& aC2D2 = aIC.SecondCurve2d();
 843     //
 844     for (j = 0; j < 2; ++j)
 845     {
 846       const Handle(Geom2d_Curve)& aC2D = !j ? aC2D1 : aC2D2;
 847       const Handle(Geom_Surface)& aS = !j ? aS1 : aS2;
 848       //
 849       if (!aC2D.IsNull())
 850       {
 851         if (IntTools_Tools::ComputeTolerance
 852             (aC3D, aC2D, aS, aFirst, aLast, aD, aT))
 853         {
 854           if (aD > aDMax)
 855           {
 856             aDMax = aD;
 857           }
 858         }
 859       }
 860       else
 861       {
 862         const TopoDS_Face& aF = !j ? myFace1 : myFace2;
 863         aD = FindMaxDistance(aC3D, aFirst, aLast, aF, myContext);
 864         if (aD > aDMax)
 865         {
 866           aDMax = aD;
 867         }
 868       }
 869     }
 870   }
 871   //
 872   return aDMax;
 873 }
 874
 875 //=======================================================================
 876 //function :ComputeTolReached3d
 877 //purpose  :
 878 //=======================================================================
 879 void IntTools_FaceFace::ComputeTolReached3d()
 880 {
 881   Standard_Integer aNbLin;
 882   GeomAbs_SurfaceType aType1, aType2;
 883   //
 884   aNbLin=myIntersector.NbLines();
 885   if (!aNbLin) {
 886     return;
 887   }
 888   //
 889   aType1=myHS1->Surface().GetType();
 890   aType2=myHS2->Surface().GetType();
 891   //
 892   if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder)
 893   {
 894     if (aNbLin==2)
 895     {
 896       Handle(IntPatch_Line) aIL1, aIL2;
 897       IntPatch_IType aTL1, aTL2;
 898       //
 899       aIL1=myIntersector.Line(1);
 900       aIL2=myIntersector.Line(2);
 901       aTL1=aIL1->ArcType();
 902       aTL2=aIL2->ArcType();
 903       if (aTL1==IntPatch_Lin && aTL2==IntPatch_Lin) {
 904         Standard_Real aD, aDTresh, dTol;
 905         gp_Lin aL1, aL2;
 906         //
 907         dTol=1.e-8;
 908         aDTresh=1.5e-6;
 909         //
 910         aL1=Handle(IntPatch_GLine)::DownCast(aIL1)->Line();
 911         aL2=Handle(IntPatch_GLine)::DownCast(aIL2)->Line();
 912         aD=aL1.Distance(aL2);
 913         aD=0.5*aD;
 914         if (aD<aDTresh)
 915         {//In order to avoid creation too thin face
 916           myTolReached3d=aD+dTol;
 917         }
 918       }
 919     }
 920   }// if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
 921   //
 922
 923   Standard_Real aDMax = ComputeTolerance();
 924   if (aDMax > myTolReached3d)
 925   {
 926     myTolReached3d = aDMax;
 927   }
 928 }
 929
 930 //=======================================================================
 931 //function : MakeCurve
 932 //purpose  :
 933 //=======================================================================
 934   void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
 935                                     const Handle(Adaptor3d_TopolTool)& dom1,
 936                                     const Handle(Adaptor3d_TopolTool)& dom2)
 937 {
 938   Standard_Boolean bDone, rejectSurface, reApprox, bAvoidLineConstructor;
 939   Standard_Boolean ok, bPCurvesOk;
 940   Standard_Integer i, j, aNbParts;
 941   Standard_Real fprm, lprm;
 942   Standard_Real Tolpc;
 943   Handle(IntPatch_Line) L;
 944   IntPatch_IType typl;
 945   Handle(Geom_Curve) newc;
 946   //
 947   const Standard_Real TOLCHECK   =0.0000001;
 948   const Standard_Real TOLANGCHECK=0.1;
 949   //
 950   rejectSurface = Standard_False;
 951   reApprox = Standard_False;
 952   //
 953   bPCurvesOk = Standard_True;
 954
 955  reapprox:;
 956
 957   Tolpc = myTolApprox;
 958   bAvoidLineConstructor = Standard_False;
 959   L = myIntersector.Line(Index);
 960   typl = L->ArcType();
 961   //
 962   if(typl==IntPatch_Walking) {
 963     Handle(IntPatch_Line) anewL;
 964     //
 965     Handle(IntPatch_WLine) aWLine (Handle(IntPatch_WLine)::DownCast(L));
 966     //DumpWLine(aWLine);
 967
 968     anewL = ComputePurgedWLine(aWLine);
 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             }//if(typs1 == GeomAbs_Plane) {
1888
1889             else if(typs2 == GeomAbs_Plane) {
1890               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1891               nbpoles = mbspc.NbPoles();
1892
1893               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1894               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1895               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
1896               const gp_Pln&  Pln = myHS2->Surface().Plane();
1897               //
1898               Standard_Integer ik;
1899               for(ik = 1; ik<= nbpoles; ik++) {
1900                 tpoles.SetValue(ik,
1901                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1902                                               tpoles2d.Value(ik).Y(),
1903                                               Pln));
1904
1905               }
1906               //
1907               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1908                                                                  mbspc.Knots(),
1909                                                                  mbspc.Multiplicities(),
1910                                                                  mbspc.Degree());
1911               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1912               Check.FixTangent(Standard_True,Standard_True);
1913               //
1914               IntTools_Curve aCurve;
1915               aCurve.SetCurve(BS);
1916
1917               if(myApprox2) {
1918                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
1919                                                                         mbspc.Knots(),
1920                                                                         mbspc.Multiplicities(),
1921                                                                         mbspc.Degree());
1922                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1923                 Check1.FixTangent(Standard_True,Standard_True);
1924                 //
1925                 // ###########################################
1926                 if(!rejectSurface && !reApprox) {
1927                   Standard_Boolean isValid = IsCurveValid(BS1);
1928                   if(!isValid) {
1929                     reApprox = Standard_True;
1930                     goto reapprox;
1931                   }
1932                 }
1933                 // ###########################################
1934                 bPCurvesOk = CheckPCurve(BS1, myFace2);
1935                 aCurve.SetSecondCurve2d(BS1);
1936               }
1937               else {
1938                 Handle(Geom2d_BSplineCurve) H2;
1939                 aCurve.SetSecondCurve2d(H2);
1940               }
1941
1942               if(myApprox1) {
1943                 mbspc.Curve(1,tpoles2d);
1944                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1945                                                                         mbspc.Knots(),
1946                                                                         mbspc.Multiplicities(),
1947                                                                         mbspc.Degree());
1948                 GeomLib_Check2dBSplineCurve Check2(BS2,TOLCHECK,TOLANGCHECK);
1949                 Check2.FixTangent(Standard_True,Standard_True);
1950                 //
1951                 // ###########################################
1952                 if(!rejectSurface && !reApprox) {
1953                   Standard_Boolean isValid = IsCurveValid(BS2);
1954                   if(!isValid) {
1955                     reApprox = Standard_True;
1956                     goto reapprox;
1957                   }
1958                 }
1959                 // ###########################################
1960                 bPCurvesOk = bPCurvesOk && CheckPCurve(BS2, myFace1);
1961                 aCurve.SetFirstCurve2d(BS2);
1962               }
1963               else {
1964                 Handle(Geom2d_BSplineCurve) H1;
1965                 //
1966                 aCurve.SetFirstCurve2d(H1);
1967               }
1968               //
1969               //if points of the pcurves are out of the faces bounds
1970               //create 3d and 2d curves without approximation
1971               if (!bPCurvesOk) {
1972                 Handle(Geom2d_BSplineCurve) H1, H2;
1973                 bPCurvesOk = Standard_True;
1974                 //
1975                 Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
1976
1977                 if(myApprox1) {
1978                   H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1979                   bPCurvesOk = CheckPCurve(H1, myFace1);
1980                 }
1981
1982                 if(myApprox2) {
1983                   H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1984                   bPCurvesOk = bPCurvesOk && CheckPCurve(H2, myFace2);
1985                 }
1986                 //
1987                 //if pcurves created without approximation are out of the
1988                 //faces bounds, use approximated 3d and 2d curves
1989                 if (bPCurvesOk) {
1990                   IntTools_Curve aIC(aBSp, H1, H2);
1991                   mySeqOfCurve.Append(aIC);
1992                 } else {
1993                   mySeqOfCurve.Append(aCurve);
1994                 }
1995               } else {
1996                 mySeqOfCurve.Append(aCurve);
1997               }
1998             }// else if(typs2 == GeomAbs_Plane)
1999             //
2000             else { //typs2 != GeomAbs_Plane && typs1 != GeomAbs_Plane
2001               Standard_Boolean bIsValid1, bIsValid2;
2002               Handle(Geom_BSplineCurve) BS;
2003               Handle(Geom2d_BSplineCurve) aH2D;
2004               IntTools_Curve aCurve;
2005               //
2006               bIsValid1=Standard_True;
2007               bIsValid2=Standard_True;
2008               //
2009               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
2010               nbpoles = mbspc.NbPoles();
2011               TColgp_Array1OfPnt tpoles(1,nbpoles);
2012               mbspc.Curve(1,tpoles);
2013               BS=new Geom_BSplineCurve(tpoles,
2014                                                                  mbspc.Knots(),
2015                                                                  mbspc.Multiplicities(),
2016                                                                  mbspc.Degree());
2017               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
2018               Check.FixTangent(Standard_True,Standard_True);
2019               //
2020               aCurve.SetCurve(BS);
2021               aCurve.SetFirstCurve2d(aH2D);
2022               aCurve.SetSecondCurve2d(aH2D);
2023               //
2024               if(myApprox1) {
2025                 if(anApprox1) {
2026                   Handle(Geom2d_BSplineCurve) BS1;
2027                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2028                   mbspc.Curve(2,tpoles2d);
2029                   //
2030                   BS1=new Geom2d_BSplineCurve(tpoles2d,
2031                                                                         mbspc.Knots(),
2032                                                                         mbspc.Multiplicities(),
2033                                                                         mbspc.Degree());
2034                   GeomLib_Check2dBSplineCurve newCheck(BS1,TOLCHECK,TOLANGCHECK);
2035                   newCheck.FixTangent(Standard_True,Standard_True);
2036                   //
2037                   if (!reApprox) {
2038                     bIsValid1=CheckPCurve(BS1, myFace1);
2039                   }
2040                   //
2041                   aCurve.SetFirstCurve2d(BS1);
2042                 }
2043                 else {
2044                   Handle(Geom2d_BSplineCurve) BS1;
2045                   fprm = BS->FirstParameter();
2046                   lprm = BS->LastParameter();
2047
2048                   Handle(Geom2d_Curve) C2d;
2049                   Standard_Real aTol = myTolApprox;
2050                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
2051                             myHS1->ChangeSurface().Surface(), BS, C2d);
2052                   BS1 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2053                   aCurve.SetFirstCurve2d(BS1);
2054                 }
2055               } // if(myApprox1) {
2056                 //
2057               if(myApprox2) {
2058                 if(anApprox2) {
2059                   Handle(Geom2d_BSplineCurve) BS2;
2060                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2061                   mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
2062                   BS2=new Geom2d_BSplineCurve(tpoles2d,
2063                                                                         mbspc.Knots(),
2064                                                                         mbspc.Multiplicities(),
2065                                                                         mbspc.Degree());
2066                   GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
2067                   newCheck.FixTangent(Standard_True,Standard_True);
2068                 //
2069                   if (!reApprox) {
2070                     bIsValid2=CheckPCurve(BS2, myFace2);
2071                   }
2072                   aCurve.SetSecondCurve2d(BS2);
2073                 }
2074                 else {
2075                   Handle(Geom2d_BSplineCurve) BS2;
2076                   fprm = BS->FirstParameter();
2077                   lprm = BS->LastParameter();
2078
2079                   Handle(Geom2d_Curve) C2d;
2080                   Standard_Real aTol = myTolApprox;
2081                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
2082                             myHS2->ChangeSurface().Surface(), BS, C2d);
2083                   BS2 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2084                   aCurve.SetSecondCurve2d(BS2);
2085                 }
2086               } //if(myApprox2) {
2087               if (!bIsValid1 || !bIsValid2) {
2088                 myTolApprox=aTolApproxImp;//1.e-5;
2089                 tol2d = myTolApprox;
2090                 reApprox = Standard_True;
2091                 goto reapprox;
2092               }
2093                 //
2094               mySeqOfCurve.Append(aCurve);
2095             }
2096           }
2097         }
2098       }
2099     }// else { // X
2100   }// case IntPatch_Walking:{
2101     break;
2102
2103   case IntPatch_Restriction:
2104     {
2105       Handle(IntPatch_RLine) RL =
2106         Handle(IntPatch_RLine)::DownCast(L);
2107       Handle(Geom_Curve) aC3d;
2108       Handle(Geom2d_Curve) aC2d1, aC2d2;
2109       Standard_Real aTolReached;
2110       GeomInt_IntSS::TreatRLine(RL, myHS1, myHS2, aC3d,
2111                                   aC2d1, aC2d2, aTolReached);
2112
2113       if(aC3d.IsNull())
2114         break;
2115
2116       Bnd_Box2d aBox1, aBox2;
2117
2118       const Standard_Real aU1f = myHS1->FirstUParameter(),
2119                           aV1f = myHS1->FirstVParameter(),
2120                           aU1l = myHS1->LastUParameter(),
2121                           aV1l = myHS1->LastVParameter();
2122       const Standard_Real aU2f = myHS2->FirstUParameter(),
2123                           aV2f = myHS2->FirstVParameter(),
2124                           aU2l = myHS2->LastUParameter(),
2125                           aV2l = myHS2->LastVParameter();
2126
2127       aBox1.Add(gp_Pnt2d(aU1f, aV1f));
2128       aBox1.Add(gp_Pnt2d(aU1l, aV1l));
2129       aBox2.Add(gp_Pnt2d(aU2f, aV2f));
2130       aBox2.Add(gp_Pnt2d(aU2l, aV2l));
2131
2132       GeomInt_VectorOfReal anArrayOfParameters;
2133
2134       //We consider here that the intersection line is same-parameter-line
2135       anArrayOfParameters.Append(aC3d->FirstParameter());
2136       anArrayOfParameters.Append(aC3d->LastParameter());
2137
2138       GeomInt_IntSS::
2139         TrimILineOnSurfBoundaries(aC2d1, aC2d2, aBox1, aBox2, anArrayOfParameters);
2140
2141       const Standard_Integer aNbIntersSolutionsm1 = anArrayOfParameters.Length() - 1;
2142
2143       //Trim RLine found.
2144       for(Standard_Integer anInd = 0; anInd < aNbIntersSolutionsm1; anInd++)
2145       {
2146         const Standard_Real aParF = anArrayOfParameters(anInd),
2147                             aParL = anArrayOfParameters(anInd+1);
2148
2149         if((aParL - aParF) <= Precision::PConfusion())
2150           continue;
2151
2152         const Standard_Real aPar = 0.5*(aParF + aParL);
2153         gp_Pnt2d aPt;
2154
2155         Handle(Geom2d_Curve) aCurv2d1, aCurv2d2;
2156         if(!aC2d1.IsNull())
2157         {
2158           aC2d1->D0(aPar, aPt);
2159
2160           if(aBox1.IsOut(aPt))
2161             continue;
2162
2163           if(myApprox1)
2164             aCurv2d1 = new Geom2d_TrimmedCurve(aC2d1, aParF, aParL);
2165         }
2166
2167         if(!aC2d2.IsNull())
2168         {
2169           aC2d2->D0(aPar, aPt);
2170
2171           if(aBox2.IsOut(aPt))
2172             continue;
2173
2174           if(myApprox2)
2175             aCurv2d2 = new Geom2d_TrimmedCurve(aC2d2, aParF, aParL);
2176         }
2177
2178         Handle(Geom_Curve) aCurv3d = new Geom_TrimmedCurve(aC3d, aParF, aParL);
2179
2180         IntTools_Curve aIC(aCurv3d, aCurv2d1, aCurv2d2);
2181         mySeqOfCurve.Append(aIC);
2182       }
2183     }
2184     break;
2185   default:
2186     break;
2187
2188   }
2189 }
2190
2191 //=======================================================================
2192 //function : Parameters
2193 //purpose  :
2194 //=======================================================================
2195  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
2196                  const Handle(GeomAdaptor_HSurface)& HS2,
2197                  const gp_Pnt& Ptref,
2198                  Standard_Real& U1,
2199                  Standard_Real& V1,
2200                  Standard_Real& U2,
2201                  Standard_Real& V2)
2202 {
2203
2204   IntSurf_Quadric quad1,quad2;
2205   GeomAbs_SurfaceType typs = HS1->Surface().GetType();
2206
2207   switch (typs) {
2208   case GeomAbs_Plane:
2209     quad1.SetValue(HS1->Surface().Plane());
2210     break;
2211   case GeomAbs_Cylinder:
2212     quad1.SetValue(HS1->Surface().Cylinder());
2213     break;
2214   case GeomAbs_Cone:
2215     quad1.SetValue(HS1->Surface().Cone());
2216     break;
2217   case GeomAbs_Sphere:
2218     quad1.SetValue(HS1->Surface().Sphere());
2219     break;
2220   case GeomAbs_Torus:
2221     quad1.SetValue(HS1->Surface().Torus());
2222     break;
2223   default:
2224     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2225   }
2226
2227   typs = HS2->Surface().GetType();
2228   switch (typs) {
2229   case GeomAbs_Plane:
2230     quad2.SetValue(HS2->Surface().Plane());
2231     break;
2232   case GeomAbs_Cylinder:
2233     quad2.SetValue(HS2->Surface().Cylinder());
2234     break;
2235   case GeomAbs_Cone:
2236     quad2.SetValue(HS2->Surface().Cone());
2237     break;
2238   case GeomAbs_Sphere:
2239     quad2.SetValue(HS2->Surface().Sphere());
2240     break;
2241   case GeomAbs_Torus:
2242     quad2.SetValue(HS2->Surface().Torus());
2243     break;
2244   default:
2245     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2246   }
2247
2248   quad1.Parameters(Ptref,U1,V1);
2249   quad2.Parameters(Ptref,U2,V2);
2250 }
2251
2252 //=======================================================================
2253 //function : MakeBSpline
2254 //purpose  :
2255 //=======================================================================
2256 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
2257                                  const Standard_Integer ideb,
2258                                  const Standard_Integer ifin)
2259 {
2260   Standard_Integer i,nbpnt = ifin-ideb+1;
2261   TColgp_Array1OfPnt poles(1,nbpnt);
2262   TColStd_Array1OfReal knots(1,nbpnt);
2263   TColStd_Array1OfInteger mults(1,nbpnt);
2264   Standard_Integer ipidebm1;
2265   for(i=1,ipidebm1=i+ideb-1; i<=nbpnt;ipidebm1++, i++) {
2266     poles(i) = WL->Point(ipidebm1).Value();
2267     mults(i) = 1;
2268     knots(i) = i-1;
2269   }
2270   mults(1) = mults(nbpnt) = 2;
2271   return
2272     new Geom_BSplineCurve(poles,knots,mults,1);
2273 }
2274 //
2275
2276 //=======================================================================
2277 //function : MakeBSpline2d
2278 //purpose  :
2279 //=======================================================================
2280 Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
2281                                           const Standard_Integer ideb,
2282                                           const Standard_Integer ifin,
2283                                           const Standard_Boolean onFirst)
2284 {
2285   Standard_Integer i, nbpnt = ifin-ideb+1;
2286   TColgp_Array1OfPnt2d poles(1,nbpnt);
2287   TColStd_Array1OfReal knots(1,nbpnt);
2288   TColStd_Array1OfInteger mults(1,nbpnt);
2289   Standard_Integer ipidebm1;
2290
2291   for(i = 1, ipidebm1 = i+ideb-1; i <= nbpnt; ipidebm1++, i++) {
2292       Standard_Real U, V;
2293       if(onFirst)
2294         theWLine->Point(ipidebm1).ParametersOnS1(U, V);
2295       else
2296         theWLine->Point(ipidebm1).ParametersOnS2(U, V);
2297       poles(i).SetCoord(U, V);
2298       mults(i) = 1;
2299       knots(i) = i-1;
2300     }
2301     mults(1) = mults(nbpnt) = 2;
2302
2303   return new Geom2d_BSplineCurve(poles,knots,mults,1);
2304 }
2305 //=======================================================================
2306 //function : PrepareLines3D
2307 //purpose  :
2308 //=======================================================================
2309   void IntTools_FaceFace::PrepareLines3D(const Standard_Boolean bToSplit)
2310 {
2311   Standard_Integer i, aNbCurves;
2312   GeomAbs_SurfaceType aType1, aType2;
2313   IntTools_SequenceOfCurves aNewCvs;
2314   //
2315   // 1. Treatment closed  curves
2316   aNbCurves=mySeqOfCurve.Length();
2317   for (i=1; i<=aNbCurves; ++i) {
2318     const IntTools_Curve& aIC=mySeqOfCurve(i);
2319     //
2320     if (bToSplit) {
2321       Standard_Integer j, aNbC;
2322       IntTools_SequenceOfCurves aSeqCvs;
2323       //
2324       aNbC=IntTools_Tools::SplitCurve(aIC, aSeqCvs);
2325       if (aNbC) {
2326         for (j=1; j<=aNbC; ++j) {
2327           const IntTools_Curve& aICNew=aSeqCvs(j);
2328           aNewCvs.Append(aICNew);
2329         }
2330       }
2331       else {
2332         aNewCvs.Append(aIC);
2333       }
2334     }
2335     else {
2336       aNewCvs.Append(aIC);
2337     }
2338   }
2339   //
2340   // 2. Plane\Cone intersection when we had 4 curves
2341   aType1=myHS1->GetType();
2342   aType2=myHS2->GetType();
2343   aNbCurves=aNewCvs.Length();
2344   //
2345   if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone) ||
2346       (aType2==GeomAbs_Plane && aType1==GeomAbs_Cone)) {
2347     if (aNbCurves==4) {
2348       GeomAbs_CurveType aCType1;
2349       //
2350       aCType1=aNewCvs(1).Type();
2351       if (aCType1==GeomAbs_Line) {
2352         IntTools_SequenceOfCurves aSeqIn, aSeqOut;
2353         //
2354         for (i=1; i<=aNbCurves; ++i) {
2355           const IntTools_Curve& aIC=aNewCvs(i);
2356           aSeqIn.Append(aIC);
2357         }
2358         //
2359         IntTools_Tools::RejectLines(aSeqIn, aSeqOut);
2360         //
2361         aNewCvs.Clear();
2362         aNbCurves=aSeqOut.Length();
2363         for (i=1; i<=aNbCurves; ++i) {
2364           const IntTools_Curve& aIC=aSeqOut(i);
2365           aNewCvs.Append(aIC);
2366         }
2367       }
2368     }
2369   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone)...
2370   //
2371   // 3. Fill  mySeqOfCurve
2372   mySeqOfCurve.Clear();
2373   aNbCurves=aNewCvs.Length();
2374   for (i=1; i<=aNbCurves; ++i) {
2375     const IntTools_Curve& aIC=aNewCvs(i);
2376     mySeqOfCurve.Append(aIC);
2377   }
2378 }
2379 //=======================================================================
2380 //function : CorrectSurfaceBoundaries
2381 //purpose  :
2382 //=======================================================================
2383  void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
2384                               const Standard_Real theTolerance,
2385                               Standard_Real&      theumin,
2386                               Standard_Real&      theumax,
2387                               Standard_Real&      thevmin,
2388                               Standard_Real&      thevmax)
2389 {
2390   Standard_Boolean enlarge, isuperiodic, isvperiodic;
2391   Standard_Real uinf, usup, vinf, vsup, delta;
2392   GeomAbs_SurfaceType aType;
2393   Handle(Geom_Surface) aSurface;
2394   //
2395   aSurface = BRep_Tool::Surface(theFace);
2396   aSurface->Bounds(uinf, usup, vinf, vsup);
2397   delta = theTolerance;
2398   enlarge = Standard_False;
2399   //
2400   GeomAdaptor_Surface anAdaptorSurface(aSurface);
2401   //
2402   if(aSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
2403     Handle(Geom_Surface) aBasisSurface =
2404       (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface))->BasisSurface();
2405
2406     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2407        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2408       return;
2409     }
2410   }
2411   //
2412   if(aSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2413     Handle(Geom_Surface) aBasisSurface =
2414       (Handle(Geom_OffsetSurface)::DownCast(aSurface))->BasisSurface();
2415
2416     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2417        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2418       return;
2419     }
2420   }
2421   //
2422   isuperiodic = anAdaptorSurface.IsUPeriodic();
2423   isvperiodic = anAdaptorSurface.IsVPeriodic();
2424   //
2425   aType=anAdaptorSurface.GetType();
2426   if((aType==GeomAbs_BezierSurface) ||
2427      (aType==GeomAbs_BSplineSurface) ||
2428      (aType==GeomAbs_SurfaceOfExtrusion) ||
2429      (aType==GeomAbs_SurfaceOfRevolution) ||
2430      (aType==GeomAbs_Cylinder)) {
2431     enlarge=Standard_True;
2432   }
2433   //
2434   if(!isuperiodic && enlarge) {
2435
2436     if(!Precision::IsInfinite(theumin) &&
2437         ((theumin - uinf) > delta))
2438       theumin -= delta;
2439     else {
2440       theumin = uinf;
2441     }
2442
2443     if(!Precision::IsInfinite(theumax) &&
2444         ((usup - theumax) > delta))
2445       theumax += delta;
2446     else
2447       theumax = usup;
2448   }
2449   //
2450   if(!isvperiodic && enlarge) {
2451     if(!Precision::IsInfinite(thevmin) &&
2452         ((thevmin - vinf) > delta)) {
2453       thevmin -= delta;
2454     }
2455     else {
2456       thevmin = vinf;
2457     }
2458     if(!Precision::IsInfinite(thevmax) &&
2459         ((vsup - thevmax) > delta)) {
2460       thevmax += delta;
2461     }
2462     else {
2463       thevmax = vsup;
2464     }
2465   }
2466   //
2467   {
2468     Standard_Integer aNbP;
2469     Standard_Real aXP, dXfact, aXmid, aX1, aX2, aTolPA;
2470     //
2471     aTolPA=Precision::Angular();
2472     // U
2473     if (isuperiodic) {
2474       aXP=anAdaptorSurface.UPeriod();
2475       dXfact=theumax-theumin;
2476       if (dXfact-aTolPA>aXP) {
2477         aXmid=0.5*(theumax+theumin);
2478         aNbP=RealToInt(aXmid/aXP);
2479         if (aXmid<0.) {
2480           aNbP=aNbP-1;
2481         }
2482         aX1=aNbP*aXP;
2483         if (theumin>aTolPA) {
2484           aX1=theumin+aNbP*aXP;
2485         }
2486         aX2=aX1+aXP;
2487         if (theumin<aX1) {
2488           theumin=aX1;
2489         }
2490         if (theumax>aX2) {
2491           theumax=aX2;
2492         }
2493       }
2494     }
2495     // V
2496     if (isvperiodic) {
2497       aXP=anAdaptorSurface.VPeriod();
2498       dXfact=thevmax-thevmin;
2499       if (dXfact-aTolPA>aXP) {
2500         aXmid=0.5*(thevmax+thevmin);
2501         aNbP=RealToInt(aXmid/aXP);
2502         if (aXmid<0.) {
2503           aNbP=aNbP-1;
2504         }
2505         aX1=aNbP*aXP;
2506         if (thevmin>aTolPA) {
2507           aX1=thevmin+aNbP*aXP;
2508         }
2509         aX2=aX1+aXP;
2510         if (thevmin<aX1) {
2511           thevmin=aX1;
2512         }
2513         if (thevmax>aX2) {
2514           thevmax=aX2;
2515         }
2516       }
2517     }
2518   }
2519   //
2520   if(isuperiodic || isvperiodic) {
2521     Standard_Boolean correct = Standard_False;
2522     Standard_Boolean correctU = Standard_False;
2523     Standard_Boolean correctV = Standard_False;
2524     Bnd_Box2d aBox;
2525     TopExp_Explorer anExp;
2526
2527     for(anExp.Init(theFace, TopAbs_EDGE); anExp.More(); anExp.Next()) {
2528       if(BRep_Tool::IsClosed(TopoDS::Edge(anExp.Current()), theFace)) {
2529         correct = Standard_True;
2530         Standard_Real f, l;
2531         TopoDS_Edge anEdge = TopoDS::Edge(anExp.Current());
2532
2533         for(Standard_Integer i = 0; i < 2; i++) {
2534           if(i==0) {
2535             anEdge.Orientation(TopAbs_FORWARD);
2536           }
2537           else {
2538             anEdge.Orientation(TopAbs_REVERSED);
2539           }
2540           Handle(Geom2d_Curve) aCurve = BRep_Tool::CurveOnSurface(anEdge, theFace, f, l);
2541
2542           if(aCurve.IsNull()) {
2543             correct = Standard_False;
2544             break;
2545           }
2546           Handle(Geom2d_Line) aLine = Handle(Geom2d_Line)::DownCast(aCurve);
2547
2548           if(aLine.IsNull()) {
2549             correct = Standard_False;
2550             break;
2551           }
2552           gp_Dir2d anUDir(1., 0.);
2553           gp_Dir2d aVDir(0., 1.);
2554           Standard_Real anAngularTolerance = Precision::Angular();
2555
2556           correctU = correctU || aLine->Position().Direction().IsParallel(aVDir, anAngularTolerance);
2557           correctV = correctV || aLine->Position().Direction().IsParallel(anUDir, anAngularTolerance);
2558
2559           gp_Pnt2d pp1 = aCurve->Value(f);
2560           aBox.Add(pp1);
2561           gp_Pnt2d pp2 = aCurve->Value(l);
2562           aBox.Add(pp2);
2563         }
2564         if(!correct)
2565           break;
2566       }
2567     }
2568
2569     if(correct) {
2570       Standard_Real umin, vmin, umax, vmax;
2571       aBox.Get(umin, vmin, umax, vmax);
2572
2573       if(isuperiodic && correctU) {
2574
2575         if(theumin < umin)
2576           theumin = umin;
2577
2578         if(theumax > umax) {
2579           theumax = umax;
2580         }
2581       }
2582       if(isvperiodic && correctV) {
2583
2584         if(thevmin < vmin)
2585           thevmin = vmin;
2586         if(thevmax > vmax)
2587           thevmax = vmax;
2588       }
2589     }
2590   }
2591 }
2592 //
2593 //
2594 // The block is dedicated to determine whether WLine [ifprm, ilprm]
2595 // crosses the degenerated zone on each given surface or not.
2596 // If Yes -> We will not use info about surfaces during approximation
2597 // because inside degenerated zone of the surface the approx. algo.
2598 // uses wrong values of normal, etc., and resulting curve will have
2599 // oscillations that we would not like to have.
2600
2601
2602
2603 static
2604   Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2605                                      const Handle(Geom_Surface)& aS,
2606                                      const Standard_Integer iDir);
2607 static
2608   Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2609                                             const TopoDS_Face& aF1,
2610                                             const TopoDS_Face& aF2);
2611 //=======================================================================
2612 //function :  NotUseSurfacesForApprox
2613 //purpose  :
2614 //=======================================================================
2615 Standard_Boolean NotUseSurfacesForApprox(const TopoDS_Face& aF1,
2616                                          const TopoDS_Face& aF2,
2617                                          const Handle(IntPatch_WLine)& WL,
2618                                          const Standard_Integer ifprm,
2619                                          const Standard_Integer ilprm)
2620 {
2621   Standard_Boolean bPInDZ;
2622
2623   Handle(IntSurf_LineOn2S) aLineOn2S=WL->Curve();
2624
2625   const IntSurf_PntOn2S& aP2Sfprm=aLineOn2S->Value(ifprm);
2626   bPInDZ=IsPointInDegeneratedZone(aP2Sfprm, aF1, aF2);
2627   if (bPInDZ) {
2628     return bPInDZ;
2629   }
2630
2631   const IntSurf_PntOn2S& aP2Slprm=aLineOn2S->Value(ilprm);
2632   bPInDZ=IsPointInDegeneratedZone(aP2Slprm, aF1, aF2);
2633
2634   return bPInDZ;
2635 }
2636 //=======================================================================
2637 //function : IsPointInDegeneratedZone
2638 //purpose  :
2639 //=======================================================================
2640 Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2641                                           const TopoDS_Face& aF1,
2642                                           const TopoDS_Face& aF2)
2643
2644 {
2645   Standard_Boolean bFlag=Standard_True;
2646   Standard_Real US11, US12, VS11, VS12, US21, US22, VS21, VS22;
2647   Standard_Real U1, V1, U2, V2, aDelta, aD;
2648   gp_Pnt2d aP2d;
2649
2650   Handle(Geom_Surface)aS1 = BRep_Tool::Surface(aF1);
2651   aS1->Bounds(US11, US12, VS11, VS12);
2652   GeomAdaptor_Surface aGAS1(aS1);
2653
2654   Handle(Geom_Surface)aS2 = BRep_Tool::Surface(aF2);
2655   aS1->Bounds(US21, US22, VS21, VS22);
2656   GeomAdaptor_Surface aGAS2(aS2);
2657   //
2658   //const gp_Pnt& aP=aP2S.Value();
2659   aP2S.Parameters(U1, V1, U2, V2);
2660   //
2661   aDelta=1.e-7;
2662   // Check on Surf 1
2663   aD=aGAS1.UResolution(aDelta);
2664   aP2d.SetCoord(U1, V1);
2665   if (fabs(U1-US11) < aD) {
2666     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2667     if (bFlag) {
2668       return bFlag;
2669     }
2670   }
2671   if (fabs(U1-US12) < aD) {
2672     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2673     if (bFlag) {
2674       return bFlag;
2675     }
2676   }
2677   aD=aGAS1.VResolution(aDelta);
2678   if (fabs(V1-VS11) < aDelta) {
2679     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2680     if (bFlag) {
2681       return bFlag;
2682     }
2683   }
2684   if (fabs(V1-VS12) < aDelta) {
2685     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2686     if (bFlag) {
2687       return bFlag;
2688     }
2689   }
2690   // Check on Surf 2
2691   aD=aGAS2.UResolution(aDelta);
2692   aP2d.SetCoord(U2, V2);
2693   if (fabs(U2-US21) < aDelta) {
2694     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2695     if (bFlag) {
2696       return bFlag;
2697     }
2698   }
2699   if (fabs(U2-US22) < aDelta) {
2700     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2701     if (bFlag) {
2702       return bFlag;
2703     }
2704   }
2705   aD=aGAS2.VResolution(aDelta);
2706   if (fabs(V2-VS21) < aDelta) {
2707     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2708     if (bFlag) {
2709       return bFlag;
2710     }
2711   }
2712   if (fabs(V2-VS22) < aDelta) {
2713     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2714     if (bFlag) {
2715       return bFlag;
2716     }
2717   }
2718   return !bFlag;
2719 }
2720
2721 //=======================================================================
2722 //function : IsDegeneratedZone
2723 //purpose  :
2724 //=======================================================================
2725 Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2726                                    const Handle(Geom_Surface)& aS,
2727                                    const Standard_Integer iDir)
2728 {
2729   Standard_Boolean bFlag=Standard_True;
2730   Standard_Real US1, US2, VS1, VS2, dY, dX, d1, d2, dD;
2731   Standard_Real aXm, aYm, aXb, aYb, aXe, aYe;
2732   aS->Bounds(US1, US2, VS1, VS2);
2733
2734   gp_Pnt aPm, aPb, aPe;
2735
2736   aXm=aP2d.X();
2737   aYm=aP2d.Y();
2738
2739   aS->D0(aXm, aYm, aPm);
2740
2741   dX=1.e-5;
2742   dY=1.e-5;
2743   dD=1.e-12;
2744
2745   if (iDir==1) {
2746     aXb=aXm;
2747     aXe=aXm;
2748     aYb=aYm-dY;
2749     if (aYb < VS1) {
2750       aYb=VS1;
2751     }
2752     aYe=aYm+dY;
2753     if (aYe > VS2) {
2754       aYe=VS2;
2755     }
2756     aS->D0(aXb, aYb, aPb);
2757     aS->D0(aXe, aYe, aPe);
2758
2759     d1=aPm.Distance(aPb);
2760     d2=aPm.Distance(aPe);
2761     if (d1 < dD && d2 < dD) {
2762       return bFlag;
2763     }
2764     return !bFlag;
2765   }
2766   //
2767   else if (iDir==2) {
2768     aYb=aYm;
2769     aYe=aYm;
2770     aXb=aXm-dX;
2771     if (aXb < US1) {
2772       aXb=US1;
2773     }
2774     aXe=aXm+dX;
2775     if (aXe > US2) {
2776       aXe=US2;
2777     }
2778     aS->D0(aXb, aYb, aPb);
2779     aS->D0(aXe, aYe, aPe);
2780
2781     d1=aPm.Distance(aPb);
2782     d2=aPm.Distance(aPe);
2783     if (d1 < dD && d2 < dD) {
2784       return bFlag;
2785     }
2786     return !bFlag;
2787   }
2788   return !bFlag;
2789 }
2790
2791 //=========================================================================
2792 // static function : ComputePurgedWLine
2793 // purpose : Removes equal points (leave one of equal points) from theWLine
2794 //           and recompute vertex parameters.
2795 //           Returns new WLine or null WLine if the number
2796 //           of the points is less than 2.
2797 //=========================================================================
2798 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine) {
2799
2800   Standard_Integer i, k, v, nb, nbvtx;
2801   Handle(IntPatch_WLine) aResult;
2802   nbvtx = theWLine->NbVertex();
2803   nb = theWLine->NbPnts();
2804   if (nb==2) {
2805     const IntSurf_PntOn2S& p1 = theWLine->Point(1);
2806     const IntSurf_PntOn2S& p2 = theWLine->Point(2);
2807     if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
2808       return aResult;
2809     }
2810   }
2811   //
2812   Handle(IntPatch_WLine) aLocalWLine;
2813   Handle(IntPatch_WLine) aTmpWLine = theWLine;
2814   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
2815   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
2816   for(i = 1; i <= nb; i++) {
2817     aLineOn2S->Add(theWLine->Point(i));
2818   }
2819
2820   for(v = 1; v <= nbvtx; v++) {
2821     aLocalWLine->AddVertex(theWLine->Vertex(v));
2822   }
2823
2824   for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
2825     Standard_Integer aStartIndex = i + 1;
2826     Standard_Integer anEndIndex = i + 5;
2827     nb = aLineOn2S->NbPoints();
2828     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
2829
2830     if((aStartIndex > nb) || (anEndIndex <= 1)) {
2831       continue;
2832     }
2833     k = aStartIndex;
2834
2835     while(k <= anEndIndex) {
2836
2837       if(i != k) {
2838         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
2839         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
2840
2841         if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
2842           aTmpWLine = aLocalWLine;
2843           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
2844
2845           for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
2846             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
2847             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
2848
2849             if(avertexindex >= k) {
2850               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
2851             }
2852             aLocalWLine->AddVertex(aVertex);
2853           }
2854           aLineOn2S->RemovePoint(k);
2855           anEndIndex--;
2856           continue;
2857         }
2858       }
2859       k++;
2860     }
2861   }
2862
2863   if(aLineOn2S->NbPoints() > 1) {
2864     aResult = aLocalWLine;
2865   }
2866   return aResult;
2867 }
2868
2869 //=======================================================================
2870 //function : TolR3d
2871 //purpose  :
2872 //=======================================================================
2873 void TolR3d(const TopoDS_Face& aF1,
2874             const TopoDS_Face& aF2,
2875             Standard_Real& myTolReached3d)
2876 {
2877   Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
2878
2879   aTolTresh=2.999999e-3;
2880   aTolF1 = BRep_Tool::Tolerance(aF1);
2881   aTolF2 = BRep_Tool::Tolerance(aF2);
2882   aTolFMax=Max(aTolF1, aTolF2);
2883
2884   if (aTolFMax>aTolTresh) {
2885     myTolReached3d=aTolFMax;
2886   }
2887 }
2888 //=======================================================================
2889 //function : IsPointOnBoundary
2890 //purpose  :
2891 //=======================================================================
2892 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
2893                                    const Standard_Real theFirstBoundary,
2894                                    const Standard_Real theSecondBoundary,
2895                                    const Standard_Real theResolution,
2896                                    Standard_Boolean&   IsOnFirstBoundary)
2897 {
2898   Standard_Boolean bRet;
2899   Standard_Integer i;
2900   Standard_Real adist;
2901   //
2902   bRet=Standard_False;
2903   for(i = 0; i < 2; ++i) {
2904     IsOnFirstBoundary = (i == 0);
2905     if (IsOnFirstBoundary) {
2906       adist = fabs(theParameter - theFirstBoundary);
2907     }
2908     else {
2909       adist = fabs(theParameter - theSecondBoundary);
2910     }
2911     if(adist < theResolution) {
2912       return !bRet;
2913     }
2914   }
2915   return bRet;
2916 }
2917 // ------------------------------------------------------------------------------------------------
2918 // static function: FindPoint
2919 // purpose:
2920 // ------------------------------------------------------------------------------------------------
2921 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
2922                            const gp_Pnt2d&     theLastPoint,
2923                            const Standard_Real theUmin,
2924                            const Standard_Real theUmax,
2925                            const Standard_Real theVmin,
2926                            const Standard_Real theVmax,
2927                            gp_Pnt2d&           theNewPoint) {
2928
2929   gp_Vec2d aVec(theFirstPoint, theLastPoint);
2930   Standard_Integer i = 0, j = 0;
2931
2932   for(i = 0; i < 4; i++) {
2933     gp_Vec2d anOtherVec;
2934     gp_Vec2d anOtherVecNormal;
2935     gp_Pnt2d aprojpoint = theLastPoint;
2936
2937     if((i % 2) == 0) {
2938       anOtherVec.SetX(0.);
2939       anOtherVec.SetY(1.);
2940       anOtherVecNormal.SetX(1.);
2941       anOtherVecNormal.SetY(0.);
2942
2943       if(i < 2)
2944         aprojpoint.SetX(theUmin);
2945       else
2946         aprojpoint.SetX(theUmax);
2947     }
2948     else {
2949       anOtherVec.SetX(1.);
2950       anOtherVec.SetY(0.);
2951       anOtherVecNormal.SetX(0.);
2952       anOtherVecNormal.SetY(1.);
2953
2954       if(i < 2)
2955         aprojpoint.SetY(theVmin);
2956       else
2957         aprojpoint.SetY(theVmax);
2958     }
2959     gp_Vec2d anormvec = aVec;
2960     anormvec.Normalize();
2961     RefineVector(anormvec);
2962     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
2963
2964     if(fabs(adot1) < Precision::Angular())
2965       continue;
2966     Standard_Real adist = 0.;
2967     Standard_Boolean bIsOut = Standard_False;
2968
2969     if((i % 2) == 0) {
2970       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
2971       bIsOut = (i < 2) ? (theLastPoint.X() < theUmin) : (theLastPoint.X() > theUmax);
2972     }
2973     else {
2974       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
2975       bIsOut = (i < 2) ? (theLastPoint.Y() < theVmin) : (theLastPoint.Y() > theVmax);
2976     }
2977     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
2978
2979     for(j = 0; j < 2; j++) {
2980       anoffset = (j == 0) ? anoffset : -anoffset;
2981       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
2982       gp_Vec2d acurvec(theLastPoint, acurpoint);
2983       if ( bIsOut )
2984         acurvec.Reverse();
2985
2986       Standard_Real aDotX, anAngleX;
2987       //
2988       aDotX = aVec.Dot(acurvec);
2989       anAngleX = aVec.Angle(acurvec);
2990       //
2991       if(aDotX > 0. && fabs(anAngleX) < Precision::PConfusion()) {
2992         if((i % 2) == 0) {
2993           if((acurpoint.Y() >= theVmin) &&
2994              (acurpoint.Y() <= theVmax)) {
2995             theNewPoint = acurpoint;
2996             return Standard_True;
2997           }
2998         }
2999         else {
3000           if((acurpoint.X() >= theUmin) &&
3001              (acurpoint.X() <= theUmax)) {
3002             theNewPoint = acurpoint;
3003             return Standard_True;
3004           }
3005         }
3006       }
3007     }
3008   }
3009   return Standard_False;
3010 }
3011
3012
3013 // ------------------------------------------------------------------------------------------------
3014 // static function: FindPoint
3015 // purpose: Find point on the boundary of radial tangent zone
3016 // ------------------------------------------------------------------------------------------------
3017 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
3018                            const gp_Pnt2d&     theLastPoint,
3019                            const Standard_Real theUmin,
3020                            const Standard_Real theUmax,
3021                            const Standard_Real theVmin,
3022                            const Standard_Real theVmax,
3023                            const gp_Pnt2d&     theTanZoneCenter,
3024                            const Standard_Real theZoneRadius,
3025                            Handle(GeomAdaptor_HSurface) theGASurface,
3026                            gp_Pnt2d&           theNewPoint) {
3027   theNewPoint = theLastPoint;
3028
3029   if ( !IsInsideTanZone( theLastPoint, theTanZoneCenter, theZoneRadius, theGASurface) )
3030     return Standard_False;
3031
3032   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3033   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3034
3035   Standard_Real aRadius = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3036   gp_Ax22d anAxis( theTanZoneCenter, gp_Dir2d(1, 0), gp_Dir2d(0, 1) );
3037   gp_Circ2d aCircle( anAxis, aRadius );
3038
3039   //
3040   gp_Vec2d aDir( theLastPoint.XY() - theFirstPoint.XY() );
3041   Standard_Real aLength = aDir.Magnitude();
3042   if ( aLength <= gp::Resolution() )
3043     return Standard_False;
3044   gp_Lin2d aLine( theFirstPoint, aDir );
3045
3046   //
3047   Handle(Geom2d_Line) aCLine = new Geom2d_Line( aLine );
3048   Handle(Geom2d_TrimmedCurve) aC1 = new Geom2d_TrimmedCurve( aCLine, 0, aLength );
3049   Handle(Geom2d_Circle) aC2 = new Geom2d_Circle( aCircle );
3050
3051   Standard_Real aTol = aRadius * 0.001;
3052   aTol = ( aTol < Precision::PConfusion() ) ? Precision::PConfusion() : aTol;
3053
3054   Geom2dAPI_InterCurveCurve anIntersector;
3055   anIntersector.Init( aC1, aC2, aTol );
3056
3057   if ( anIntersector.NbPoints() == 0 )
3058     return Standard_False;
3059
3060   Standard_Boolean aFound = Standard_False;
3061   Standard_Real aMinDist = aLength * aLength;
3062   Standard_Integer i = 0;
3063   for ( i = 1; i <= anIntersector.NbPoints(); i++ ) {
3064     gp_Pnt2d aPInt = anIntersector.Point( i );
3065     if ( aPInt.SquareDistance( theFirstPoint ) < aMinDist ) {
3066       if ( ( aPInt.X() >= theUmin ) && ( aPInt.X() <= theUmax ) &&
3067            ( aPInt.Y() >= theVmin ) && ( aPInt.Y() <= theVmax ) ) {
3068         theNewPoint = aPInt;
3069         aFound = Standard_True;
3070       }
3071     }
3072   }
3073
3074   return aFound;
3075 }
3076
3077 // ------------------------------------------------------------------------------------------------
3078 // static function: IsInsideTanZone
3079 // purpose: Check if point is inside a radial tangent zone
3080 // ------------------------------------------------------------------------------------------------
3081 Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
3082                                  const gp_Pnt2d&     theTanZoneCenter,
3083                                  const Standard_Real theZoneRadius,
3084                                  Handle(GeomAdaptor_HSurface) theGASurface) {
3085
3086   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3087   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3088   Standard_Real aRadiusSQR = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3089   aRadiusSQR *= aRadiusSQR;
3090   if ( thePoint.SquareDistance( theTanZoneCenter ) <= aRadiusSQR )
3091     return Standard_True;
3092   return Standard_False;
3093 }
3094
3095 // ------------------------------------------------------------------------------------------------
3096 // static function: CheckTangentZonesExist
3097 // purpose: Check if tangent zone exists
3098 // ------------------------------------------------------------------------------------------------
3099 Standard_Boolean CheckTangentZonesExist( const Handle(GeomAdaptor_HSurface)& theSurface1,
3100                                         const Handle(GeomAdaptor_HSurface)&  theSurface2 )
3101 {
3102   if ( ( theSurface1->GetType() != GeomAbs_Torus ) ||
3103       ( theSurface2->GetType() != GeomAbs_Torus ) )
3104     return Standard_False;
3105
3106   gp_Torus aTor1 = theSurface1->Torus();
3107   gp_Torus aTor2 = theSurface2->Torus();
3108
3109   if ( aTor1.Location().Distance( aTor2.Location() ) > Precision::Confusion() )
3110     return Standard_False;
3111
3112   if ( ( fabs( aTor1.MajorRadius() - aTor2.MajorRadius() ) > Precision::Confusion() ) ||
3113        ( fabs( aTor1.MinorRadius() - aTor2.MinorRadius() ) > Precision::Confusion() ) )
3114     return Standard_False;
3115
3116   if ( ( aTor1.MajorRadius() < aTor1.MinorRadius() ) ||
3117        ( aTor2.MajorRadius() < aTor2.MinorRadius() ) )
3118     return Standard_False;
3119   return Standard_True;
3120 }
3121
3122 // ------------------------------------------------------------------------------------------------
3123 // static function: ComputeTangentZones
3124 // purpose:
3125 // ------------------------------------------------------------------------------------------------
3126 Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
3127                                      const Handle(GeomAdaptor_HSurface)&  theSurface2,
3128                                      const TopoDS_Face&                   theFace1,
3129                                      const TopoDS_Face&                   theFace2,
3130                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
3131                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
3132                                      Handle(TColStd_HArray1OfReal)&       theResultRadius,
3133                                      const Handle(IntTools_Context)& aContext)
3134 {
3135   Standard_Integer aResult = 0;
3136   if ( !CheckTangentZonesExist( theSurface1, theSurface2 ) )
3137     return aResult;
3138
3139
3140   TColgp_SequenceOfPnt2d aSeqResultS1, aSeqResultS2;
3141   TColStd_SequenceOfReal aSeqResultRad;
3142
3143   gp_Torus aTor1 = theSurface1->Torus();
3144   gp_Torus aTor2 = theSurface2->Torus();
3145
3146   gp_Ax2 anax1( aTor1.Location(), aTor1.Axis().Direction() );
3147   gp_Ax2 anax2( aTor2.Location(), aTor2.Axis().Direction() );
3148   Standard_Integer j = 0;
3149
3150   for ( j = 0; j < 2; j++ ) {
3151     Standard_Real aCoef = ( j == 0 ) ? -1 : 1;
3152     Standard_Real aRadius1 = fabs(aTor1.MajorRadius() + aCoef * aTor1.MinorRadius());
3153     Standard_Real aRadius2 = fabs(aTor2.MajorRadius() + aCoef * aTor2.MinorRadius());
3154
3155     gp_Circ aCircle1( anax1, aRadius1 );
3156     gp_Circ aCircle2( anax2, aRadius2 );
3157
3158     // roughly compute radius of tangent zone for perpendicular case
3159     Standard_Real aCriteria = Precision::Confusion() * 0.5;
3160
3161     Standard_Real aT1 = aCriteria;
3162     Standard_Real aT2 = aCriteria;
3163     if ( j == 0 ) {
3164       // internal tangency
3165       Standard_Real aR = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3166       //aT1 = aCriteria * aCriteria + aR * aR - ( aR - aCriteria ) * ( aR - aCriteria );
3167       aT1 = 2. * aR * aCriteria;
3168       aT2 = aT1;
3169     }
3170     else {
3171       // external tangency
3172       Standard_Real aRb = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3173       Standard_Real aRm = ( aRadius1 < aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3174       Standard_Real aDelta = aRb - aCriteria;
3175       aDelta *= aDelta;
3176       aDelta -= aRm * aRm;
3177       aDelta /= 2. * (aRb - aRm);
3178       aDelta -= 0.5 * (aRb - aRm);
3179
3180       aT1 = 2. * aRm * (aRm - aDelta);
3181       aT2 = aT1;
3182     }
3183     aCriteria = ( aT1 > aT2) ? aT1 : aT2;
3184     if ( aCriteria > 0 )
3185       aCriteria = sqrt( aCriteria );
3186
3187     if ( aCriteria > 0.5 * aTor1.MinorRadius() ) {
3188       // too big zone -> drop to minimum
3189       aCriteria = Precision::Confusion();
3190     }
3191
3192     GeomAdaptor_Curve aC1( new Geom_Circle(aCircle1) );
3193     GeomAdaptor_Curve aC2( new Geom_Circle(aCircle2) );
3194     Extrema_ExtCC anExtrema(aC1, aC2, 0, 2. * M_PI, 0, 2. * M_PI,
3195                             Precision::PConfusion(), Precision::PConfusion());
3196
3197     if ( anExtrema.IsDone() ) {
3198
3199       Standard_Integer i = 0;
3200       for ( i = 1; i <= anExtrema.NbExt(); i++ ) {
3201         if ( anExtrema.SquareDistance(i) > aCriteria * aCriteria )
3202           continue;
3203
3204         Extrema_POnCurv P1, P2;
3205         anExtrema.Points( i, P1, P2 );
3206
3207         Standard_Boolean bFoundResult = Standard_True;
3208         gp_Pnt2d pr1, pr2;
3209
3210         Standard_Integer surfit = 0;
3211         for ( surfit = 0; surfit < 2; surfit++ ) {
3212           GeomAPI_ProjectPointOnSurf& aProjector =
3213             (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
3214
3215           gp_Pnt aP3d = (surfit == 0) ? P1.Value() : P2.Value();
3216           aProjector.Perform(aP3d);
3217
3218           if(!aProjector.IsDone())
3219             bFoundResult = Standard_False;
3220           else {
3221             if(aProjector.LowerDistance() > aCriteria) {
3222               bFoundResult = Standard_False;
3223             }
3224             else {
3225               Standard_Real foundU = 0, foundV = 0;
3226               aProjector.LowerDistanceParameters(foundU, foundV);
3227               if ( surfit == 0 )
3228                 pr1 = gp_Pnt2d( foundU, foundV );
3229               else
3230                 pr2 = gp_Pnt2d( foundU, foundV );
3231             }
3232           }
3233         }
3234         if ( bFoundResult ) {
3235           aSeqResultS1.Append( pr1 );
3236           aSeqResultS2.Append( pr2 );
3237           aSeqResultRad.Append( aCriteria );
3238
3239           // torus is u and v periodic
3240           const Standard_Real twoPI = M_PI + M_PI;
3241           Standard_Real arr1tmp[2] = {pr1.X(), pr1.Y()};
3242           Standard_Real arr2tmp[2] = {pr2.X(), pr2.Y()};
3243
3244           // iteration on period bounds
3245           for ( Standard_Integer k1 = 0; k1 < 2; k1++ ) {
3246             Standard_Real aBound = ( k1 == 0 ) ? 0 : twoPI;
3247             Standard_Real aShift = ( k1 == 0 ) ? twoPI : -twoPI;
3248
3249             // iteration on surfaces
3250             for ( Standard_Integer k2 = 0; k2 < 2; k2++ ) {
3251               Standard_Real* arr1 = ( k2 == 0 ) ? arr1tmp : arr2tmp;
3252               Standard_Real* arr2 = ( k2 != 0 ) ? arr1tmp : arr2tmp;
3253               TColgp_SequenceOfPnt2d& aSeqS1 = ( k2 == 0 ) ? aSeqResultS1 : aSeqResultS2;
3254               TColgp_SequenceOfPnt2d& aSeqS2 = ( k2 != 0 ) ? aSeqResultS1 : aSeqResultS2;
3255
3256               if (fabs(arr1[0] - aBound) < Precision::PConfusion()) {
3257                 aSeqS1.Append( gp_Pnt2d( arr1[0] + aShift, arr1[1] ) );
3258                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
3259                 aSeqResultRad.Append( aCriteria );
3260               }
3261               if (fabs(arr1[1] - aBound) < Precision::PConfusion()) {
3262                 aSeqS1.Append( gp_Pnt2d( arr1[0], arr1[1] + aShift) );
3263                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
3264                 aSeqResultRad.Append( aCriteria );
3265               }
3266             }
3267           } //
3268         }
3269       }
3270     }
3271   }
3272   aResult = aSeqResultRad.Length();
3273
3274   if ( aResult > 0 ) {
3275     theResultOnS1 = new TColgp_HArray1OfPnt2d( 1, aResult );
3276     theResultOnS2 = new TColgp_HArray1OfPnt2d( 1, aResult );
3277     theResultRadius = new TColStd_HArray1OfReal( 1, aResult );
3278
3279     for ( Standard_Integer i = 1 ; i <= aResult; i++ ) {
3280       theResultOnS1->SetValue( i, aSeqResultS1.Value(i) );
3281       theResultOnS2->SetValue( i, aSeqResultS2.Value(i) );
3282       theResultRadius->SetValue( i, aSeqResultRad.Value(i) );
3283     }
3284   }
3285   return aResult;
3286 }
3287
3288 // ------------------------------------------------------------------------------------------------
3289 // static function: AdjustByNeighbour
3290 // purpose:
3291 // ------------------------------------------------------------------------------------------------
3292 gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
3293                            const gp_Pnt2d&     theOriginalPoint,
3294                            Handle(GeomAdaptor_HSurface) theGASurface) {
3295
3296   gp_Pnt2d ap1 = theaNeighbourPoint;
3297   gp_Pnt2d ap2 = theOriginalPoint;
3298
3299   if ( theGASurface->IsUPeriodic() ) {
3300     Standard_Real aPeriod     = theGASurface->UPeriod();
3301     gp_Pnt2d aPTest = ap2;
3302     Standard_Real aSqDistMin = 1.e+100;
3303
3304     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
3305       aPTest.SetX( theOriginalPoint.X() + aPeriod * pIt );
3306       Standard_Real dd = ap1.SquareDistance( aPTest );
3307
3308       if ( dd < aSqDistMin ) {
3309         ap2 = aPTest;
3310         aSqDistMin = dd;
3311       }
3312     }
3313   }
3314   if ( theGASurface->IsVPeriodic() ) {
3315     Standard_Real aPeriod     = theGASurface->VPeriod();
3316     gp_Pnt2d aPTest = ap2;
3317     Standard_Real aSqDistMin = 1.e+100;
3318
3319     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
3320       aPTest.SetY( theOriginalPoint.Y() + aPeriod * pIt );
3321       Standard_Real dd = ap1.SquareDistance( aPTest );
3322
3323       if ( dd < aSqDistMin ) {
3324         ap2 = aPTest;
3325         aSqDistMin = dd;
3326       }
3327     }
3328   }
3329   return ap2;
3330 }
3331
3332 // ------------------------------------------------------------------------------------------------
3333 //function: DecompositionOfWLine
3334 // purpose:
3335 // ------------------------------------------------------------------------------------------------
3336 Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
3337                                       const Handle(GeomAdaptor_HSurface)&            theSurface1,
3338                                       const Handle(GeomAdaptor_HSurface)&            theSurface2,
3339                                       const TopoDS_Face&                             theFace1,
3340                                       const TopoDS_Face&                             theFace2,
3341                                       const GeomInt_LineConstructor&                 theLConstructor,
3342                                       const Standard_Boolean                         theAvoidLConstructor,
3343                                       IntPatch_SequenceOfLine&                       theNewLines,
3344                                       Standard_Real&                                 theReachedTol3d,
3345                                       const Handle(IntTools_Context)& aContext)
3346 {
3347
3348   Standard_Boolean bRet, bAvoidLineConstructor;
3349   Standard_Integer aNbPnts, aNbParts;
3350   //
3351   bRet=Standard_False;
3352   aNbPnts=theWLine->NbPnts();
3353   bAvoidLineConstructor=theAvoidLConstructor;
3354   //
3355   if(!aNbPnts){
3356     return bRet;
3357   }
3358   if (!bAvoidLineConstructor) {
3359     aNbParts=theLConstructor.NbParts();
3360     if (!aNbParts) {
3361       return bRet;
3362     }
3363   }
3364   //
3365   Standard_Boolean bIsPrevPointOnBoundary, bIsPointOnBoundary, bIsCurrentPointOnBoundary;
3366   Standard_Integer nblines, pit, i, j;
3367   Standard_Real aTol;
3368   TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts);
3369   TColStd_Array1OfInteger       anArrayOfLineType(1, aNbPnts);
3370   TColStd_ListOfInteger aListOfPointIndex;
3371
3372   Handle(TColgp_HArray1OfPnt2d) aTanZoneS1;
3373   Handle(TColgp_HArray1OfPnt2d) aTanZoneS2;
3374   Handle(TColStd_HArray1OfReal) aTanZoneRadius;
3375   Standard_Integer aNbZone = ComputeTangentZones( theSurface1, theSurface2, theFace1, theFace2,
3376                                                  aTanZoneS1, aTanZoneS2, aTanZoneRadius, aContext);
3377
3378   //
3379   nblines=0;
3380   aTol=Precision::Confusion();
3381   aTol=0.5*aTol;
3382   bIsPrevPointOnBoundary=Standard_False;
3383   bIsPointOnBoundary=Standard_False;
3384   //
3385   // 1. ...
3386   //
3387   // Points
3388   for(pit = 1; pit <= aNbPnts; ++pit) {
3389     Standard_Boolean bIsOnFirstBoundary, isperiodic;
3390     Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
3391     Standard_Real aParameter, anoffset, anAdjustPar;
3392     Standard_Real umin, umax, vmin, vmax;
3393     //
3394     bIsCurrentPointOnBoundary = Standard_False;
3395     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
3396     //
3397     // Surface
3398     for(i = 0; i < 2; ++i) {
3399       Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
3400       aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
3401       if(!i) {
3402         aPoint.ParametersOnS1(U, V);
3403       }
3404       else {
3405         aPoint.ParametersOnS2(U, V);
3406       }
3407       // U, V
3408       for(j = 0; j < 2; j++) {
3409         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3410         if(!isperiodic){
3411           continue;
3412         }
3413         //
3414         if (!j) {
3415           aResolution=aGASurface->UResolution(aTol);
3416           aPeriod=aGASurface->UPeriod();
3417           alowerboundary=umin;
3418           aupperboundary=umax;
3419           aParameter=U;
3420         }
3421         else {
3422           aResolution=aGASurface->VResolution(aTol);
3423           aPeriod=aGASurface->VPeriod();
3424           alowerboundary=vmin;
3425           aupperboundary=vmax;
3426           aParameter=V;
3427         }
3428
3429         GeomInt::AdjustPeriodic(aParameter,
3430                                        alowerboundary,
3431                                        aupperboundary,
3432                                        aPeriod,
3433                                        anAdjustPar,
3434                                        anoffset);
3435         //
3436         bIsOnFirstBoundary = Standard_True;// ?
3437         bIsPointOnBoundary=
3438           IsPointOnBoundary(anAdjustPar,
3439                             alowerboundary,
3440                             aupperboundary,
3441                             aResolution,
3442                             bIsOnFirstBoundary);
3443         //
3444         if(bIsPointOnBoundary) {
3445           bIsCurrentPointOnBoundary = Standard_True;
3446           break;
3447         }
3448         else {
3449           // check if a point belong to a tangent zone. Begin
3450           Standard_Integer zIt = 0;
3451           for ( zIt = 1; zIt <= aNbZone; zIt++ ) {
3452             gp_Pnt2d aPZone = (i == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
3453             Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
3454
3455             if ( IsInsideTanZone(gp_Pnt2d( U, V ), aPZone, aZoneRadius, aGASurface ) ) {
3456               // set boundary flag to split the curve by a tangent zone
3457               bIsPointOnBoundary = Standard_True;
3458               bIsCurrentPointOnBoundary = Standard_True;
3459               if ( theReachedTol3d < aZoneRadius ) {
3460                 theReachedTol3d = aZoneRadius;
3461               }
3462               break;
3463             }
3464           }
3465         }
3466       }//for(j = 0; j < 2; j++) {
3467
3468       if(bIsCurrentPointOnBoundary){
3469         break;
3470       }
3471     }//for(i = 0; i < 2; ++i) {
3472     //
3473     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
3474       if(!aListOfPointIndex.IsEmpty()) {
3475         nblines++;
3476         anArrayOfLines.SetValue(nblines, aListOfPointIndex);
3477         anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
3478         aListOfPointIndex.Clear();
3479       }
3480       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
3481     }
3482     aListOfPointIndex.Append(pit);
3483   } //for(pit = 1; pit <= aNbPnts; ++pit) {
3484   //
3485   if(!aListOfPointIndex.IsEmpty()) {
3486     nblines++;
3487     anArrayOfLines.SetValue(nblines, aListOfPointIndex);
3488     anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
3489     aListOfPointIndex.Clear();
3490   }
3491   //
3492   if(nblines<=1) {
3493     return bRet; //Standard_False;
3494   }
3495   //
3496   //
3497   // 2. Correct wlines.begin
3498   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
3499   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S();
3500   //
3501   for(i = 1; i <= nblines; i++) {
3502     if(anArrayOfLineType.Value(i) != 0) {
3503       continue;
3504     }
3505     const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
3506     TColStd_ListOfInteger aListOfFLIndex;
3507
3508     for(j = 0; j < 2; j++) {
3509       Standard_Integer aneighbourindex = (j == 0) ? (i - 1) : (i + 1);
3510
3511       if((aneighbourindex < 1) || (aneighbourindex > nblines))
3512         continue;
3513
3514       if(anArrayOfLineType.Value(aneighbourindex) == 0)
3515         continue;
3516       const TColStd_ListOfInteger& aNeighbour = anArrayOfLines.Value(aneighbourindex);
3517       Standard_Integer anIndex = (j == 0) ? aNeighbour.Last() : aNeighbour.First();
3518       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
3519
3520       IntSurf_PntOn2S aNewP = aPoint;
3521       if(aListOfIndex.Extent() < 2) {
3522         aSeqOfPntOn2S->Add(aNewP);
3523         aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
3524         continue;
3525       }
3526       //
3527       Standard_Integer iFirst = aListOfIndex.First();
3528       Standard_Integer iLast  = aListOfIndex.Last();
3529       //
3530       for(Standard_Integer surfit = 0; surfit < 2; surfit++) {
3531
3532         Handle(GeomAdaptor_HSurface) aGASurface = (surfit == 0) ? theSurface1 : theSurface2;
3533         Standard_Real umin=0., umax=0., vmin=0., vmax=0.;
3534         aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
3535         Standard_Real U=0., V=0.;
3536
3537         if(surfit == 0)
3538           aNewP.ParametersOnS1(U, V);
3539         else
3540           aNewP.ParametersOnS2(U, V);
3541         Standard_Integer nbboundaries = 0;
3542
3543         Standard_Boolean bIsNearBoundary = Standard_False;
3544         Standard_Integer aZoneIndex = 0;
3545         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
3546         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
3547
3548
3549         for(Standard_Integer parit = 0; parit < 2; parit++) {
3550           Standard_Boolean isperiodic = (parit == 0) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3551
3552           Standard_Real aResolution = (parit == 0) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
3553           Standard_Real alowerboundary = (parit == 0) ? umin : vmin;
3554           Standard_Real aupperboundary = (parit == 0) ? umax : vmax;
3555
3556           Standard_Real aParameter = (parit == 0) ? U : V;
3557           Standard_Boolean bIsOnFirstBoundary = Standard_True;
3558
3559           if(!isperiodic) {
3560             bIsPointOnBoundary=
3561               IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
3562             if(bIsPointOnBoundary) {
3563               bIsUBoundary = (parit == 0);
3564               bIsFirstBoundary = bIsOnFirstBoundary;
3565               nbboundaries++;
3566             }
3567           }
3568           else {
3569             Standard_Real aPeriod     = (parit == 0) ? aGASurface->UPeriod() : aGASurface->VPeriod();
3570             Standard_Real anoffset, anAdjustPar;
3571             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
3572                                            aPeriod, anAdjustPar, anoffset);
3573
3574             bIsPointOnBoundary=
3575               IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
3576             if(bIsPointOnBoundary) {
3577               bIsUBoundary = (parit == 0);
3578               bIsFirstBoundary = bIsOnFirstBoundary;
3579               nbboundaries++;
3580             }
3581             else {
3582             //check neighbourhood of boundary
3583             Standard_Real anEpsilon = aResolution * 100.;
3584             Standard_Real aPart = ( aupperboundary - alowerboundary ) * 0.1;
3585             anEpsilon = ( anEpsilon > aPart ) ? aPart : anEpsilon;
3586
3587               bIsNearBoundary = IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary,
3588                                                   anEpsilon, bIsOnFirstBoundary);
3589
3590             }
3591           }
3592         }
3593
3594         // check if a point belong to a tangent zone. Begin
3595         for ( Standard_Integer zIt = 1; zIt <= aNbZone; zIt++ ) {
3596           gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
3597           Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
3598
3599           Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
3600           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3601           Standard_Real nU1, nV1;
3602
3603           if(surfit == 0)
3604             aNeighbourPoint.ParametersOnS1(nU1, nV1);
3605           else
3606             aNeighbourPoint.ParametersOnS2(nU1, nV1);
3607           gp_Pnt2d ap1(nU1, nV1);
3608           gp_Pnt2d ap2 = AdjustByNeighbour( ap1, gp_Pnt2d( U, V ), aGASurface );
3609
3610
3611           if ( IsInsideTanZone( ap2, aPZone, aZoneRadius, aGASurface ) ) {
3612             aZoneIndex = zIt;
3613             bIsNearBoundary = Standard_True;
3614             if ( theReachedTol3d < aZoneRadius ) {
3615               theReachedTol3d = aZoneRadius;
3616             }
3617           }
3618         }
3619         // check if a point belong to a tangent zone. End
3620         Standard_Boolean bComputeLineEnd = Standard_False;
3621
3622         if(nbboundaries == 2) {
3623           //xf
3624           bComputeLineEnd = Standard_True;
3625           //xt
3626         }
3627         else if(nbboundaries == 1) {
3628           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3629
3630           if(isperiodic) {
3631             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
3632             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
3633             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
3634             Standard_Real aParameter = (bIsUBoundary) ? U : V;
3635             Standard_Real anoffset, anAdjustPar;
3636             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
3637                                            aPeriod, anAdjustPar, anoffset);
3638
3639             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
3640             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
3641             anotherPar += anoffset;
3642             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
3643             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
3644             Standard_Real nU1, nV1;
3645
3646             if(surfit == 0)
3647               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3648             else
3649               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3650
3651             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
3652             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
3653             bComputeLineEnd = Standard_True;
3654             Standard_Boolean bCheckAngle1 = Standard_False;
3655             Standard_Boolean bCheckAngle2 = Standard_False;
3656             gp_Vec2d aNewVec;
3657             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
3658             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
3659
3660             if(((adist1 - adist2) > Precision::PConfusion()) &&
3661                (adist2 < (aPeriod / 4.))) {
3662               bCheckAngle1 = Standard_True;
3663               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
3664
3665               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
3666                 aNewP.SetValue((surfit == 0), anewU, anewV);
3667                 bCheckAngle1 = Standard_False;
3668               }
3669             }
3670             else if(adist1 < (aPeriod / 4.)) {
3671               bCheckAngle2 = Standard_True;
3672               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
3673
3674               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
3675                 bCheckAngle2 = Standard_False;
3676               }
3677             }
3678
3679             if(bCheckAngle1 || bCheckAngle2) {
3680               // assume there are at least two points in line (see "if" above)
3681               Standard_Integer anindexother = aneighbourpointindex;
3682
3683               while((anindexother <= iLast) && (anindexother >= iFirst)) {
3684                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
3685                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
3686                 Standard_Real nU2, nV2;
3687
3688                 if(surfit == 0)
3689                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
3690                 else
3691                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
3692                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
3693
3694                 if(aVecOld.SquareMagnitude() <= gp::Resolution()) {
3695                   continue;
3696                 }
3697                 else {
3698                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
3699
3700                   if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
3701
3702                     if(bCheckAngle1) {
3703                       Standard_Real U1, U2, V1, V2;
3704                       IntSurf_PntOn2S atmppoint = aNewP;
3705                       atmppoint.SetValue((surfit == 0), anewU, anewV);
3706                       atmppoint.Parameters(U1, V1, U2, V2);
3707                       gp_Pnt P1 = theSurface1->Value(U1, V1);
3708                       gp_Pnt P2 = theSurface2->Value(U2, V2);
3709                       gp_Pnt P0 = aPoint.Value();
3710
3711                       if(P0.IsEqual(P1, aTol) &&
3712                          P0.IsEqual(P2, aTol) &&
3713                          P1.IsEqual(P2, aTol)) {
3714                         bComputeLineEnd = Standard_False;
3715                         aNewP.SetValue((surfit == 0), anewU, anewV);
3716                       }
3717                     }
3718
3719                     if(bCheckAngle2) {
3720                       bComputeLineEnd = Standard_False;
3721                     }
3722                   }
3723                   break;
3724                 }
3725               } // end while(anindexother...)
3726             }
3727           }
3728         }
3729         else if ( bIsNearBoundary ) {
3730           bComputeLineEnd = Standard_True;
3731         }
3732
3733         if(bComputeLineEnd) {
3734
3735           gp_Pnt2d anewpoint;
3736           Standard_Boolean found = Standard_False;
3737
3738           if ( bIsNearBoundary ) {
3739             // re-compute point near natural boundary or near tangent zone
3740             Standard_Real u1, v1, u2, v2;
3741             aNewP.Parameters( u1, v1, u2, v2 );
3742             if(surfit == 0)
3743               anewpoint = gp_Pnt2d( u1, v1 );
3744             else
3745               anewpoint = gp_Pnt2d( u2, v2 );
3746
3747             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
3748             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3749             Standard_Real nU1, nV1;
3750
3751             if(surfit == 0)
3752               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3753             else
3754               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3755             gp_Pnt2d ap1(nU1, nV1);
3756             gp_Pnt2d ap2;
3757
3758
3759             if ( aZoneIndex ) {
3760               // exclude point from a tangent zone
3761               anewpoint = AdjustByNeighbour( ap1, anewpoint, aGASurface );
3762               gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(aZoneIndex) : aTanZoneS2->Value(aZoneIndex);
3763               Standard_Real aZoneRadius = aTanZoneRadius->Value(aZoneIndex);
3764
3765               if ( FindPoint(ap1, anewpoint, umin, umax, vmin, vmax,
3766                              aPZone, aZoneRadius, aGASurface, ap2) ) {
3767                 anewpoint = ap2;
3768                 found = Standard_True;
3769               }
3770             }
3771             else if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) {
3772               // re-compute point near boundary if shifted on a period
3773               ap2 = AdjustByNeighbour( ap1, anewpoint, aGASurface );
3774
3775               if ( ( ap2.X() < umin ) || ( ap2.X() > umax ) ||
3776                   ( ap2.Y() < vmin ) || ( ap2.Y() > vmax ) ) {
3777                 found = FindPoint(ap1, ap2, umin, umax, vmin, vmax, anewpoint);
3778               }
3779               else {
3780                 anewpoint = ap2;
3781                 aNewP.SetValue( (surfit == 0), anewpoint.X(), anewpoint.Y() );
3782               }
3783             }
3784           }
3785           else {
3786
3787             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
3788             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3789             Standard_Real nU1, nV1;
3790
3791             if(surfit == 0)
3792               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3793             else
3794               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3795             gp_Pnt2d ap1(nU1, nV1);
3796             gp_Pnt2d ap2(nU1, nV1);
3797             Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
3798
3799             while((aneighbourpointindex2 <= iLast) && (aneighbourpointindex2 >= iFirst)) {
3800               aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
3801               const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
3802               Standard_Real nU2, nV2;
3803
3804               if(surfit == 0)
3805                 aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
3806               else
3807                 aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
3808               ap2.SetX(nU2);
3809               ap2.SetY(nV2);
3810
3811               if(ap1.SquareDistance(ap2) > gp::Resolution()) {
3812                 break;
3813               }
3814             }
3815             found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
3816           }
3817
3818           if(found) {
3819             // check point
3820             Standard_Real aCriteria = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
3821             GeomAPI_ProjectPointOnSurf& aProjector =
3822               (surfit == 0) ? aContext->ProjPS(theFace2) : aContext->ProjPS(theFace1);
3823             Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
3824
3825             Handle(GeomAdaptor_HSurface) aSurfaceOther = (surfit == 0) ? theSurface2 : theSurface1;
3826
3827             gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
3828             aProjector.Perform(aP3d);
3829
3830             if(aProjector.IsDone()) {
3831               if(aProjector.LowerDistance() < aCriteria) {
3832                 Standard_Real foundU = U, foundV = V;
3833                 aProjector.LowerDistanceParameters(foundU, foundV);
3834
3835                 //Correction of projected coordinates. Begin
3836                 //Note, it may be shifted on a period
3837                 Standard_Integer aneindex1 = (j == 0) ? iFirst : iLast;
3838                 const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneindex1);
3839                 Standard_Real nUn, nVn;
3840
3841                 if(surfit == 0)
3842                   aNeighbourPoint.ParametersOnS2(nUn, nVn);
3843                 else
3844                   aNeighbourPoint.ParametersOnS1(nUn, nVn);
3845                 gp_Pnt2d aNeighbour2d(nUn, nVn);
3846                 gp_Pnt2d anAdjustedPoint = AdjustByNeighbour( aNeighbour2d, gp_Pnt2d(foundU, foundV), aSurfaceOther );
3847                 foundU = anAdjustedPoint.X();
3848                 foundV = anAdjustedPoint.Y();
3849
3850                 if ( ( anAdjustedPoint.X() < umin ) && ( anAdjustedPoint.X() > umax ) &&
3851                     ( anAdjustedPoint.Y() < vmin ) && ( anAdjustedPoint.Y() > vmax ) ) {
3852                   // attempt to roughly re-compute point
3853                   foundU = ( foundU < umin ) ? umin : foundU;
3854                   foundU = ( foundU > umax ) ? umax : foundU;
3855                   foundV = ( foundV < vmin ) ? vmin : foundV;
3856                   foundV = ( foundV > vmax ) ? vmax : foundV;
3857
3858                   GeomAPI_ProjectPointOnSurf& aProjector2 =
3859                     (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
3860
3861                   aP3d = aSurfaceOther->Value(foundU, foundV);
3862                   aProjector2.Perform(aP3d);
3863
3864                   if(aProjector2.IsDone()) {
3865                     if(aProjector2.LowerDistance() < aCriteria) {
3866                       Standard_Real foundU2 = anewpoint.X(), foundV2 = anewpoint.Y();
3867                       aProjector2.LowerDistanceParameters(foundU2, foundV2);
3868                       anewpoint.SetX(foundU2);
3869                       anewpoint.SetY(foundV2);
3870                     }
3871                   }
3872                 }
3873                 //Correction of projected coordinates. End
3874
3875                 if(surfit == 0)
3876                   aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
3877                 else
3878                   aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
3879               }
3880             }
3881           }
3882         }
3883       }
3884       aSeqOfPntOn2S->Add(aNewP);
3885       aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
3886     }
3887     anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
3888   }
3889   // Correct wlines.end
3890
3891   // Split wlines.begin
3892   Standard_Integer nbiter;
3893   //
3894   nbiter=1;
3895   if (!bAvoidLineConstructor) {
3896     nbiter=theLConstructor.NbParts();
3897   }
3898   //
3899   for(j = 1; j <= nbiter; ++j) {
3900     Standard_Real fprm, lprm;
3901     Standard_Integer ifprm, ilprm;
3902     //
3903     if(bAvoidLineConstructor) {
3904       ifprm = 1;
3905       ilprm = theWLine->NbPnts();
3906     }
3907     else {
3908       theLConstructor.Part(j, fprm, lprm);
3909       ifprm = (Standard_Integer)fprm;
3910       ilprm = (Standard_Integer)lprm;
3911     }
3912
3913     Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
3914     //
3915     for(i = 1; i <= nblines; i++) {
3916       if(anArrayOfLineType.Value(i) != 0) {
3917         continue;
3918       }
3919       const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
3920       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
3921       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
3922       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
3923
3924       if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
3925         bhasfirstpoint = (i != 1);
3926       }
3927
3928       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
3929         bhaslastpoint = (i != nblines);
3930       }
3931
3932       Standard_Integer iFirst = aListOfIndex.First();
3933       Standard_Integer iLast  = aListOfIndex.Last();
3934       Standard_Boolean bIsFirstInside = ((ifprm >= iFirst) && (ifprm <= iLast));
3935       Standard_Boolean bIsLastInside =  ((ilprm >= iFirst) && (ilprm <= iLast));
3936
3937       if(!bIsFirstInside && !bIsLastInside) {
3938         if((ifprm < iFirst) && (ilprm > iLast)) {
3939           // append whole line, and boundaries if neccesary
3940           if(bhasfirstpoint) {
3941             pit = aListOfFLIndex.First();
3942             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
3943             aLineOn2S->Add(aP);
3944           }
3945           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
3946
3947           for(; anIt.More(); anIt.Next()) {
3948             pit = anIt.Value();
3949             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
3950             aLineOn2S->Add(aP);
3951           }
3952
3953           if(bhaslastpoint) {
3954             pit = aListOfFLIndex.Last();
3955             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
3956             aLineOn2S->Add(aP);
3957           }
3958
3959           // check end of split line (end is almost always)
3960           Standard_Integer aneighbour = i + 1;
3961           Standard_Boolean bIsEndOfLine = Standard_True;
3962
3963           if(aneighbour <= nblines) {
3964             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
3965
3966             if((anArrayOfLineType.Value(aneighbour) != 0) &&
3967                (aListOfNeighbourIndex.IsEmpty())) {
3968               bIsEndOfLine = Standard_False;
3969             }
3970           }
3971
3972           if(bIsEndOfLine) {
3973             if(aLineOn2S->NbPoints() > 1) {
3974               Handle(IntPatch_WLine) aNewWLine =
3975                 new IntPatch_WLine(aLineOn2S, Standard_False);
3976               theNewLines.Append(aNewWLine);
3977             }
3978             aLineOn2S = new IntSurf_LineOn2S();
3979           }
3980         }
3981         continue;
3982       }
3983       // end if(!bIsFirstInside && !bIsLastInside)
3984
3985       if(bIsFirstInside && bIsLastInside) {
3986         // append inside points between ifprm and ilprm
3987         TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
3988
3989         for(; anIt.More(); anIt.Next()) {
3990           pit = anIt.Value();
3991           if((pit < ifprm) || (pit > ilprm))
3992             continue;
3993           const IntSurf_PntOn2S& aP = theWLine->Point(pit);
3994           aLineOn2S->Add(aP);
3995         }
3996       }
3997       else {
3998
3999         if(bIsFirstInside) {
4000           // append points from ifprm to last point + boundary point
4001           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4002
4003           for(; anIt.More(); anIt.Next()) {
4004             pit = anIt.Value();
4005             if(pit < ifprm)
4006               continue;
4007             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4008             aLineOn2S->Add(aP);
4009           }
4010
4011           if(bhaslastpoint) {
4012             pit = aListOfFLIndex.Last();
4013             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4014             aLineOn2S->Add(aP);
4015           }
4016           // check end of split line (end is almost always)
4017           Standard_Integer aneighbour = i + 1;
4018           Standard_Boolean bIsEndOfLine = Standard_True;
4019
4020           if(aneighbour <= nblines) {
4021             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
4022
4023             if((anArrayOfLineType.Value(aneighbour) != 0) &&
4024                (aListOfNeighbourIndex.IsEmpty())) {
4025               bIsEndOfLine = Standard_False;
4026             }
4027           }
4028
4029           if(bIsEndOfLine) {
4030             if(aLineOn2S->NbPoints() > 1) {
4031               Handle(IntPatch_WLine) aNewWLine =
4032                 new IntPatch_WLine(aLineOn2S, Standard_False);
4033               theNewLines.Append(aNewWLine);
4034             }
4035             aLineOn2S = new IntSurf_LineOn2S();
4036           }
4037         }
4038         // end if(bIsFirstInside)
4039
4040         if(bIsLastInside) {
4041           // append points from first boundary point to ilprm
4042           if(bhasfirstpoint) {
4043             pit = aListOfFLIndex.First();
4044             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4045             aLineOn2S->Add(aP);
4046           }
4047           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4048
4049           for(; anIt.More(); anIt.Next()) {
4050             pit = anIt.Value();
4051             if(pit > ilprm)
4052               continue;
4053             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4054             aLineOn2S->Add(aP);
4055           }
4056         }
4057         //end if(bIsLastInside)
4058       }
4059     }
4060
4061     if(aLineOn2S->NbPoints() > 1) {
4062       Handle(IntPatch_WLine) aNewWLine =
4063         new IntPatch_WLine(aLineOn2S, Standard_False);
4064       theNewLines.Append(aNewWLine);
4065     }
4066   }
4067   // Split wlines.end
4068
4069   return Standard_True;
4070 }
4071
4072 // ------------------------------------------------------------------------------------------------
4073 // static function: ParameterOutOfBoundary
4074 // purpose:         Computes a new parameter for given curve. The corresponding 2d points
4075 //                  does not lay on any boundary of given faces
4076 // ------------------------------------------------------------------------------------------------
4077 Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
4078                                         const Handle(Geom_Curve)& theCurve,
4079                                         const TopoDS_Face&        theFace1,
4080                                         const TopoDS_Face&        theFace2,
4081                                         const Standard_Real       theOtherParameter,
4082                                         const Standard_Boolean    bIncreasePar,
4083                                         Standard_Real&            theNewParameter,
4084                                         const Handle(IntTools_Context)& aContext)
4085 {
4086   Standard_Boolean bIsComputed = Standard_False;
4087   theNewParameter = theParameter;
4088
4089   Standard_Real acurpar = theParameter;
4090   TopAbs_State aState = TopAbs_ON;
4091   Standard_Integer iter = 0;
4092   Standard_Real asumtol = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
4093   Standard_Real adelta = asumtol * 0.1;
4094   adelta = (adelta < Precision::Confusion()) ? Precision::Confusion() : adelta;
4095   Handle(Geom_Surface) aSurf1 = BRep_Tool::Surface(theFace1);
4096   Handle(Geom_Surface) aSurf2 = BRep_Tool::Surface(theFace2);
4097
4098   Standard_Real u1, u2, v1, v2;
4099
4100   GeomAPI_ProjectPointOnSurf aPrj1;
4101   aSurf1->Bounds(u1, u2, v1, v2);
4102   aPrj1.Init(aSurf1, u1, u2, v1, v2);
4103
4104   GeomAPI_ProjectPointOnSurf aPrj2;
4105   aSurf2->Bounds(u1, u2, v1, v2);
4106   aPrj2.Init(aSurf2, u1, u2, v1, v2);
4107
4108   while(aState == TopAbs_ON) {
4109     if(bIncreasePar)
4110       acurpar += adelta;
4111     else
4112       acurpar -= adelta;
4113     gp_Pnt aPCurrent = theCurve->Value(acurpar);
4114     aPrj1.Perform(aPCurrent);
4115     Standard_Real U=0., V=0.;
4116
4117     if(aPrj1.IsDone()) {
4118       aPrj1.LowerDistanceParameters(U, V);
4119       aState = aContext->StatePointFace(theFace1, gp_Pnt2d(U, V));
4120     }
4121
4122     if(aState != TopAbs_ON) {
4123       aPrj2.Perform(aPCurrent);
4124
4125       if(aPrj2.IsDone()) {
4126         aPrj2.LowerDistanceParameters(U, V);
4127         aState = aContext->StatePointFace(theFace2, gp_Pnt2d(U, V));
4128       }
4129     }
4130
4131     if(iter > 11) {
4132       break;
4133     }
4134     iter++;
4135   }
4136
4137   if(iter <= 11) {
4138     theNewParameter = acurpar;
4139     bIsComputed = Standard_True;
4140
4141     if(bIncreasePar) {
4142       if(acurpar >= theOtherParameter)
4143         theNewParameter = theOtherParameter;
4144     }
4145     else {
4146       if(acurpar <= theOtherParameter)
4147         theNewParameter = theOtherParameter;
4148     }
4149   }
4150   return bIsComputed;
4151 }
4152
4153 //=======================================================================
4154 //function : IsCurveValid
4155 //purpose  :
4156 //=======================================================================
4157 Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve)
4158 {
4159   if(thePCurve.IsNull())
4160     return Standard_False;
4161
4162   Standard_Real tolint = 1.e-10;
4163   Geom2dAdaptor_Curve PCA;
4164   IntRes2d_Domain PCD;
4165   Geom2dInt_GInter PCI;
4166
4167   Standard_Real pf = 0., pl = 0.;
4168   gp_Pnt2d pntf, pntl;
4169
4170   if(!thePCurve->IsClosed() && !thePCurve->IsPeriodic()) {
4171     pf = thePCurve->FirstParameter();
4172     pl = thePCurve->LastParameter();
4173     pntf = thePCurve->Value(pf);
4174     pntl = thePCurve->Value(pl);
4175     PCA.Load(thePCurve);
4176     if(!PCA.IsPeriodic()) {
4177       if(PCA.FirstParameter() > pf) pf = PCA.FirstParameter();
4178       if(PCA.LastParameter()  < pl) pl = PCA.LastParameter();
4179     }
4180     PCD.SetValues(pntf,pf,tolint,pntl,pl,tolint);
4181     PCI.Perform(PCA,PCD,tolint,tolint);
4182     if(PCI.IsDone())
4183       if(PCI.NbPoints() > 0) {
4184         return Standard_False;
4185       }
4186   }
4187
4188   return Standard_True;
4189 }
4190
4191 //=======================================================================
4192 //static function : ApproxWithPCurves
4193 //purpose  : for bug 20964 only
4194 //=======================================================================
4195 Standard_Boolean ApproxWithPCurves(const gp_Cylinder& theCyl,
4196                                    const gp_Sphere& theSph)
4197 {
4198   Standard_Boolean bRes = Standard_True;
4199   Standard_Real R1 = theCyl.Radius(), R2 = theSph.Radius();
4200   //
4201   {
4202     Standard_Real aD2, aRc2, aEps;
4203     gp_Pnt aApexSph;
4204     //
4205     aEps=1.E-7;
4206     aRc2=R1*R1;
4207     //
4208     const gp_Ax3& aAx3Sph=theSph.Position();
4209     const gp_Pnt& aLocSph=aAx3Sph.Location();
4210     const gp_Dir& aDirSph=aAx3Sph.Direction();
4211     //
4212     const gp_Ax1& aAx1Cyl=theCyl.Axis();
4213     gp_Lin aLinCyl(aAx1Cyl);
4214     //
4215     aApexSph.SetXYZ(aLocSph.XYZ()+R2*aDirSph.XYZ());
4216     aD2=aLinCyl.SquareDistance(aApexSph);
4217     if (fabs(aD2-aRc2)<aEps) {
4218       return !bRes;
4219     }
4220     //
4221     aApexSph.SetXYZ(aLocSph.XYZ()-R2*aDirSph.XYZ());
4222     aD2=aLinCyl.SquareDistance(aApexSph);
4223     if (fabs(aD2-aRc2)<aEps) {
4224       return !bRes;
4225     }
4226   }
4227   //
4228
4229   if(R1 < 2.*R2) {
4230     return bRes;
4231   }
4232   gp_Lin anCylAx(theCyl.Axis());
4233
4234   Standard_Real aDist = anCylAx.Distance(theSph.Location());
4235   Standard_Real aDRel = Abs(aDist - R1)/R2;
4236
4237   if(aDRel > .2) return bRes;
4238
4239   Standard_Real par = ElCLib::Parameter(anCylAx, theSph.Location());
4240   gp_Pnt aP = ElCLib::Value(par, anCylAx);
4241   gp_Vec aV(aP, theSph.Location());
4242
4243   Standard_Real dd = aV.Dot(theSph.Position().XDirection());
4244
4245   if(aDist < R1 && dd > 0.) return Standard_False;
4246   if(aDist > R1 && dd < 0.) return Standard_False;
4247
4248
4249   return bRes;
4250 }
4251 //=======================================================================
4252 //function : PerformPlanes
4253 //purpose  :
4254 //=======================================================================
4255 void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
4256                     const Handle(GeomAdaptor_HSurface)& theS2,
4257                     const Standard_Real TolAng,
4258                     const Standard_Real TolTang,
4259                     const Standard_Boolean theApprox1,
4260                     const Standard_Boolean theApprox2,
4261                     IntTools_SequenceOfCurves& theSeqOfCurve,
4262                     Standard_Boolean& theTangentFaces)
4263 {
4264
4265   gp_Pln aPln1 = theS1->Surface().Plane();
4266   gp_Pln aPln2 = theS2->Surface().Plane();
4267
4268   IntAna_QuadQuadGeo aPlnInter(aPln1, aPln2, TolAng, TolTang);
4269
4270   if(!aPlnInter.IsDone()) {
4271     theTangentFaces = Standard_False;
4272     return;
4273   }
4274
4275   IntAna_ResultType aResType = aPlnInter.TypeInter();
4276
4277   if(aResType == IntAna_Same) {
4278     theTangentFaces = Standard_True;
4279     return;
4280   }
4281
4282   theTangentFaces = Standard_False;
4283
4284   if(aResType == IntAna_Empty) {
4285     return;
4286   }
4287
4288   gp_Lin aLin = aPlnInter.Line(1);
4289
4290   ProjLib_Plane aProj;
4291
4292   aProj.Init(aPln1);
4293   aProj.Project(aLin);
4294   gp_Lin2d aLin2d1 = aProj.Line();
4295   //
4296   aProj.Init(aPln2);
4297   aProj.Project(aLin);
4298   gp_Lin2d aLin2d2 = aProj.Line();
4299   //
4300   //classify line2d1 relatively first plane
4301   Standard_Real P11, P12;
4302   Standard_Boolean IsCrossed = ClassifyLin2d(theS1, aLin2d1, TolTang, P11, P12);
4303   if(!IsCrossed) return;
4304   //classify line2d2 relatively second plane
4305   Standard_Real P21, P22;
4306   IsCrossed = ClassifyLin2d(theS2, aLin2d2, TolTang, P21, P22);
4307   if(!IsCrossed) return;
4308
4309   //Analysis of parametric intervals: must have common part
4310
4311   if(P21 >= P12) return;
4312   if(P22 <= P11) return;
4313
4314   Standard_Real pmin, pmax;
4315   pmin = Max(P11, P21);
4316   pmax = Min(P12, P22);
4317
4318   if(pmax - pmin <= TolTang) return;
4319
4320   Handle(Geom_Line) aGLin = new Geom_Line(aLin);
4321
4322   IntTools_Curve aCurve;
4323   Handle(Geom_TrimmedCurve) aGTLin = new Geom_TrimmedCurve(aGLin, pmin, pmax);
4324
4325   aCurve.SetCurve(aGTLin);
4326
4327   if(theApprox1) {
4328     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d1);
4329     aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
4330   }
4331   else {
4332     Handle(Geom2d_Curve) H1;
4333     aCurve.SetFirstCurve2d(H1);
4334   }
4335   if(theApprox2) {
4336     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d2);
4337     aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
4338   }
4339   else {
4340     Handle(Geom2d_Curve) H1;
4341     aCurve.SetFirstCurve2d(H1);
4342   }
4343
4344   theSeqOfCurve.Append(aCurve);
4345
4346 }
4347
4348 //=======================================================================
4349 //function : ClassifyLin2d
4350 //purpose  :
4351 //=======================================================================
4352 static inline Standard_Boolean INTER(const Standard_Real d1,
4353                                      const Standard_Real d2,
4354                                      const Standard_Real tol)
4355 {
4356   return (d1 > tol && d2 < -tol) ||
4357          (d1 < -tol && d2 > tol) ||
4358          ((d1 <= tol && d1 >= -tol) && (d2 > tol || d2 < -tol)) ||
4359          ((d2 <= tol && d2 >= -tol) && (d1 > tol || d1 < -tol));
4360 }
4361 static inline Standard_Boolean COINC(const Standard_Real d1,
4362                                      const Standard_Real d2,
4363                                      const Standard_Real tol)
4364 {
4365   return (d1 <= tol && d1 >= -tol) && (d2 <= tol && d2 >= -tol);
4366 }
4367 Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS,
4368                                const gp_Lin2d& theLin2d,
4369                                const Standard_Real theTol,
4370                                Standard_Real& theP1,
4371                                Standard_Real& theP2)
4372
4373 {
4374   Standard_Real xmin, xmax, ymin, ymax, d1, d2, A, B, C;
4375   Standard_Real par[2];
4376   Standard_Integer nbi = 0;
4377
4378   xmin = theS->Surface().FirstUParameter();
4379   xmax = theS->Surface().LastUParameter();
4380   ymin = theS->Surface().FirstVParameter();
4381   ymax = theS->Surface().LastVParameter();
4382
4383   theLin2d.Coefficients(A, B, C);
4384
4385   //xmin, ymin <-> xmin, ymax
4386   d1 = A*xmin + B*ymin + C;
4387   d2 = A*xmin + B*ymax + C;
4388
4389   if(INTER(d1, d2, theTol)) {
4390     //Intersection with boundary
4391     Standard_Real y = -(C + A*xmin)/B;
4392     par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, y));
4393     nbi++;
4394   }
4395   else if (COINC(d1, d2, theTol)) {
4396     //Coincidence with boundary
4397     par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
4398     par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
4399     nbi = 2;
4400   }
4401
4402   if(nbi == 2) {
4403
4404     if(fabs(par[0]-par[1]) > theTol) {
4405       theP1 = Min(par[0], par[1]);
4406       theP2 = Max(par[0], par[1]);
4407       return Standard_True;
4408     }
4409     else return Standard_False;
4410
4411   }
4412
4413   //xmin, ymax <-> xmax, ymax
4414   d1 = d2;
4415   d2 = A*xmax + B*ymax + C;
4416
4417   if(d1 > theTol || d1 < -theTol) {//to avoid checking of
4418                                    //coincidence with the same point
4419     if(INTER(d1, d2, theTol)) {
4420       Standard_Real x = -(C + B*ymax)/A;
4421       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymax));
4422       nbi++;
4423     }
4424     else if (COINC(d1, d2, theTol)) {
4425       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
4426       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
4427       nbi = 2;
4428     }
4429   }
4430
4431   if(nbi == 2) {
4432
4433     if(fabs(par[0]-par[1]) > theTol) {
4434       theP1 = Min(par[0], par[1]);
4435       theP2 = Max(par[0], par[1]);
4436       return Standard_True;
4437     }
4438     else return Standard_False;
4439
4440   }
4441
4442   //xmax, ymax <-> xmax, ymin
4443   d1 = d2;
4444   d2 = A*xmax + B*ymin + C;
4445
4446   if(d1 > theTol || d1 < -theTol) {
4447     if(INTER(d1, d2, theTol)) {
4448       Standard_Real y = -(C + A*xmax)/B;
4449       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, y));
4450       nbi++;
4451     }
4452     else if (COINC(d1, d2, theTol)) {
4453       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
4454       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
4455       nbi = 2;
4456     }
4457   }
4458
4459   if(nbi == 2) {
4460     if(fabs(par[0]-par[1]) > theTol) {
4461       theP1 = Min(par[0], par[1]);
4462       theP2 = Max(par[0], par[1]);
4463       return Standard_True;
4464     }
4465     else return Standard_False;
4466   }
4467
4468   //xmax, ymin <-> xmin, ymin
4469   d1 = d2;
4470   d2 = A*xmin + B*ymin + C;
4471
4472   if(d1 > theTol || d1 < -theTol) {
4473     if(INTER(d1, d2, theTol)) {
4474       Standard_Real x = -(C + B*ymin)/A;
4475       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymin));
4476       nbi++;
4477     }
4478     else if (COINC(d1, d2, theTol)) {
4479       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
4480       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
4481       nbi = 2;
4482     }
4483   }
4484
4485   if(nbi == 2) {
4486     if(fabs(par[0]-par[1]) > theTol) {
4487       theP1 = Min(par[0], par[1]);
4488       theP2 = Max(par[0], par[1]);
4489       return Standard_True;
4490     }
4491     else return Standard_False;
4492   }
4493
4494   return Standard_False;
4495
4496 }
4497 //
4498 //=======================================================================
4499 //function : ApproxParameters
4500 //purpose  :
4501 //=======================================================================
4502 void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
4503                       const Handle(GeomAdaptor_HSurface)& aHS2,
4504                       Standard_Integer& iDegMin,
4505                       Standard_Integer& iDegMax,
4506                       Standard_Integer& iNbIter)
4507
4508 {
4509   GeomAbs_SurfaceType aTS1, aTS2;
4510
4511   //
4512   iNbIter=0;
4513   iDegMin=4;
4514   iDegMax=8;
4515   //
4516   aTS1=aHS1->Surface().GetType();
4517   aTS2=aHS2->Surface().GetType();
4518   //
4519   // Cylinder/Torus
4520   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
4521       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
4522     Standard_Real aRC, aRT, dR, aPC;
4523     gp_Cylinder aCylinder;
4524     gp_Torus aTorus;
4525     //
4526     aPC=Precision::Confusion();
4527     //
4528     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Surface().Cylinder() : aHS2->Surface().Cylinder();
4529     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Surface().Torus() : aHS2->Surface().Torus();
4530     //
4531     aRC=aCylinder.Radius();
4532     aRT=aTorus.MinorRadius();
4533     dR=aRC-aRT;
4534     if (dR<0.) {
4535       dR=-dR;
4536     }
4537     //
4538     if (dR<aPC) {
4539       iDegMax=6;
4540     }
4541   }
4542   if (aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Cylinder) {
4543     iNbIter=1;
4544   }
4545 }
4546 //=======================================================================
4547 //function : Tolerances
4548 //purpose  :
4549 //=======================================================================
4550 void Tolerances(const Handle(GeomAdaptor_HSurface)& aHS1,
4551                 const Handle(GeomAdaptor_HSurface)& aHS2,
4552                 Standard_Real& aTolTang)
4553 {
4554   GeomAbs_SurfaceType aTS1, aTS2;
4555   //
4556   aTS1=aHS1->Surface().GetType();
4557   aTS2=aHS2->Surface().GetType();
4558   //
4559   // Cylinder/Torus
4560   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
4561       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
4562     Standard_Real aRC, aRT, dR, aPC;
4563     gp_Cylinder aCylinder;
4564     gp_Torus aTorus;
4565     //
4566     aPC=Precision::Confusion();
4567     //
4568     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Surface().Cylinder() : aHS2->Surface().Cylinder();
4569     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Surface().Torus() : aHS2->Surface().Torus();
4570     //
4571     aRC=aCylinder.Radius();
4572     aRT=aTorus.MinorRadius();
4573     dR=aRC-aRT;
4574     if (dR<0.) {
4575       dR=-dR;
4576     }
4577     //
4578     if (dR<aPC) {
4579       aTolTang=0.1*aTolTang;
4580     }
4581   }
4582 }
4583 //=======================================================================
4584 //function : SortTypes
4585 //purpose  :
4586 //=======================================================================
4587 Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
4588                            const GeomAbs_SurfaceType aType2)
4589 {
4590   Standard_Boolean bRet;
4591   Standard_Integer aI1, aI2;
4592   //
4593   bRet=Standard_False;
4594   //
4595   aI1=IndexType(aType1);
4596   aI2=IndexType(aType2);
4597   if (aI1<aI2){
4598     bRet=!bRet;
4599   }
4600   return bRet;
4601 }
4602 //=======================================================================
4603 //function : IndexType
4604 //purpose  :
4605 //=======================================================================
4606 Standard_Integer IndexType(const GeomAbs_SurfaceType aType)
4607 {
4608   Standard_Integer aIndex;
4609   //
4610   aIndex=11;
4611   //
4612   if (aType==GeomAbs_Plane) {
4613     aIndex=0;
4614   }
4615   else if (aType==GeomAbs_Cylinder) {
4616     aIndex=1;
4617   }
4618   else if (aType==GeomAbs_Cone) {
4619     aIndex=2;
4620   }
4621   else if (aType==GeomAbs_Sphere) {
4622     aIndex=3;
4623   }
4624   else if (aType==GeomAbs_Torus) {
4625     aIndex=4;
4626   }
4627   else if (aType==GeomAbs_BezierSurface) {
4628     aIndex=5;
4629   }
4630   else if (aType==GeomAbs_BSplineSurface) {
4631     aIndex=6;
4632   }
4633   else if (aType==GeomAbs_SurfaceOfRevolution) {
4634     aIndex=7;
4635   }
4636   else if (aType==GeomAbs_SurfaceOfExtrusion) {
4637     aIndex=8;
4638   }
4639   else if (aType==GeomAbs_OffsetSurface) {
4640     aIndex=9;
4641   }
4642   else if (aType==GeomAbs_OtherSurface) {
4643     aIndex=10;
4644   }
4645   return aIndex;
4646 }
4647 #ifdef OCCT_DEBUG_DUMPWLINE
4648 //=======================================================================
4649 //function : DumpWLine
4650 //purpose  :
4651 //=======================================================================
4652 void DumpWLine(const Handle(IntPatch_WLine)& aWLine)
4653 {
4654   Standard_Integer i, aNbPnts;
4655   Standard_Real aX, aY, aZ, aU1, aV1, aU2, aV2;
4656   //
4657   printf(" *WLine\n");
4658   aNbPnts=aWLine->NbPnts();
4659   for (i=1; i<=aNbPnts; ++i) {
4660     const IntSurf_PntOn2S aPntOn2S=aWLine->Point(i);
4661     const gp_Pnt& aP3D=aPntOn2S.Value();
4662     aP3D.Coord(aX, aY, aZ);
4663     aPntOn2S.Parameters(aU1, aV1, aU2, aV2);
4664     //
4665     printf("point p_%d %lf %lf %lf\n", i, aX, aY, aZ);
4666     //printf("point p_%d %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf\n",
4667         //   i, aX, aY, aZ, aU1, aV1, aU2, aV2);
4668   }
4669 }
4670 #endif
4671 //=======================================================================
4672 //function : RefineVector
4673 //purpose  :
4674 //=======================================================================
4675 void RefineVector(gp_Vec2d& aV2D)
4676 {
4677   Standard_Integer k,m;
4678   Standard_Real aC[2], aEps, aR1, aR2, aNum;
4679   //
4680   aEps=RealEpsilon();
4681   aR1=1.-aEps;
4682   aR2=1.+aEps;
4683   //
4684   aV2D.Coord(aC[0], aC[1]);
4685   //
4686   for (k=0; k<2; ++k) {
4687     m=(k+1)%2;
4688     aNum=fabs(aC[k]);
4689     if (aNum>aR1 && aNum<aR2) {
4690       if (aC[k]<0.) {
4691         aC[k]=-1.;
4692       }
4693       else {
4694         aC[k]=1.;
4695       }
4696       aC[m]=0.;
4697       break;
4698     }
4699   }
4700   aV2D.SetCoord(aC[0], aC[1]);
4701 }
4702
4703 //=======================================================================
4704 // Function : FindMaxDistance
4705 // purpose :
4706 //=======================================================================
4707 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
4708                               const Standard_Real theFirst,
4709                               const Standard_Real theLast,
4710                               const TopoDS_Face& theFace,
4711                               const Handle(IntTools_Context)& theContext)
4712 {
4713   Standard_Integer aNbS;
4714   Standard_Real aT1, aT2, aDt, aD, aDMax, anEps;
4715   //
4716   aNbS = 11;
4717   aDt = (theLast - theFirst) / aNbS;
4718   aDMax = 0.;
4719   anEps = 1.e-4 * aDt;
4720   //
4721   GeomAPI_ProjectPointOnSurf& aProjPS = theContext->ProjPS(theFace);
4722   aT2 = theFirst;
4723   for (;;) {
4724     aT1 = aT2;
4725     aT2 += aDt;
4726     //
4727     if (aT2 > theLast) {
4728       break;
4729     }
4730     //
4731     aD = FindMaxDistance(theCurve, aT1, aT2, aProjPS, anEps);
4732     if (aD > aDMax) {
4733       aDMax = aD;
4734     }
4735   }
4736   //
4737   return aDMax;
4738 }
4739
4740 //=======================================================================
4741 // Function : FindMaxDistance
4742 // purpose :
4743 //=======================================================================
4744 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
4745                               const Standard_Real theFirst,
4746                               const Standard_Real theLast,
4747                               GeomAPI_ProjectPointOnSurf& theProjPS,
4748                               const Standard_Real theEps)
4749 {
4750   Standard_Real aA, aB, aCf, aX, aX1, aX2, aF1, aF2, aF;
4751   //
4752   aCf = 0.61803398874989484820458683436564;//(sqrt(5.)-1)/2.;
4753   aA = theFirst;
4754   aB = theLast;
4755   //
4756   aX1 = aB - aCf * (aB - aA);
4757   aF1 = MaxDistance(theC, aX1, theProjPS);
4758   aX2 = aA + aCf * (aB - aA);
4759   aF2 = MaxDistance(theC, aX2, theProjPS);
4760   //
4761   for (;;) {
4762     if ((aB - aA) < theEps) {
4763       break;
4764     }
4765     //
4766     if (aF1 > aF2) {
4767       aB = aX2;
4768       aX2 = aX1;
4769       aF2 = aF1;
4770       aX1 = aB - aCf * (aB - aA);
4771       aF1 = MaxDistance(theC, aX1, theProjPS);
4772     }
4773     else {
4774       aA = aX1;
4775       aX1 = aX2;
4776       aF1 = aF2;
4777       aX2 = aA + aCf * (aB - aA);
4778       aF2 = MaxDistance(theC, aX2, theProjPS);
4779     }
4780   }
4781   //
4782   aX = 0.5 * (aA + aB);
4783   aF = MaxDistance(theC, aX, theProjPS);
4784   //
4785   if (aF1 > aF) {
4786     aF = aF1;
4787   }
4788   //
4789   if (aF2 > aF) {
4790     aF = aF2;
4791   }
4792   //
4793   return aF;
4794 }
4795
4796 //=======================================================================
4797 // Function : MaxDistance
4798 // purpose :
4799 //=======================================================================
4800 Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
4801                           const Standard_Real aT,
4802                           GeomAPI_ProjectPointOnSurf& theProjPS)
4803 {
4804   Standard_Real aD;
4805   gp_Pnt aP;
4806   //
4807   theC->D0(aT, aP);
4808   theProjPS.Perform(aP);
4809   aD = theProjPS.NbPoints() ? theProjPS.LowerDistance() : 0.;
4810   //
4811   return aD;
4812 }
4813
4814 //=======================================================================
4815 //function : CheckPCurve
4816 //purpose  : Checks if points of the pcurve are out of the face bounds.
4817 //=======================================================================
4818   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
4819                                const TopoDS_Face& aFace)
4820 {
4821   const Standard_Integer NPoints = 23;
4822   Standard_Integer i;
4823   Standard_Real umin,umax,vmin,vmax;
4824
4825   BRepTools::UVBounds(aFace, umin, umax, vmin, vmax);
4826   Standard_Real tolU = Max ((umax-umin)*0.01, Precision::Confusion());
4827   Standard_Real tolV = Max ((vmax-vmin)*0.01, Precision::Confusion());
4828   Standard_Real fp = aPC->FirstParameter();
4829   Standard_Real lp = aPC->LastParameter();
4830
4831
4832   // adjust domain for periodic surfaces
4833   TopLoc_Location aLoc;
4834   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace, aLoc);
4835   if (aSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
4836     aSurf = (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurf))->BasisSurface();
4837   }
4838   gp_Pnt2d pnt = aPC->Value((fp+lp)/2);
4839   Standard_Real u,v;
4840   pnt.Coord(u,v);
4841   //
4842   if (aSurf->IsUPeriodic()) {
4843     Standard_Real aPer = aSurf->UPeriod();
4844     Standard_Integer nshift = (Standard_Integer) ((u-umin)/aPer);
4845     if (u < umin+aPer*nshift) nshift--;
4846     umin += aPer*nshift;
4847     umax += aPer*nshift;
4848   }
4849   if (aSurf->IsVPeriodic()) {
4850     Standard_Real aPer = aSurf->VPeriod();
4851     Standard_Integer nshift = (Standard_Integer) ((v-vmin)/aPer);
4852     if (v < vmin+aPer*nshift) nshift--;
4853     vmin += aPer*nshift;
4854     vmax += aPer*nshift;
4855   }
4856   //
4857   //--------------------------------------------------------
4858   Standard_Boolean bRet;
4859   Standard_Integer j, aNbIntervals;
4860   Standard_Real aT, dT;
4861   gp_Pnt2d aP2D;
4862   //
4863   Geom2dAdaptor_Curve aGAC(aPC);
4864   aNbIntervals=aGAC.NbIntervals(GeomAbs_CN);
4865   //
4866   TColStd_Array1OfReal aTI(1, aNbIntervals+1);
4867   aGAC.Intervals(aTI,GeomAbs_CN);
4868   //
4869   bRet=Standard_False;
4870   //
4871   aT=aGAC.FirstParameter();
4872   for (j=1; j<=aNbIntervals; ++j) {
4873     dT=(aTI(j+1)-aTI(j))/NPoints;
4874     //
4875     for (i=1; i<NPoints; i++) {
4876       aT=aT+dT;
4877       aGAC.D0(aT, aP2D);
4878       aP2D.Coord(u,v);
4879     if (umin-u > tolU || u-umax > tolU ||
4880           vmin-v > tolV || v-vmax > tolV) {
4881         return bRet;
4882   }
4883 }
4884   }
4885   return !bRet;
4886 }
4887 //=======================================================================
4888 //function : CorrectPlaneBoundaries
4889 //purpose  :
4890 //=======================================================================
4891  void CorrectPlaneBoundaries(Standard_Real& aUmin,
4892                              Standard_Real& aUmax,
4893                              Standard_Real& aVmin,
4894                              Standard_Real& aVmax)
4895 {
4896   if (!(Precision::IsInfinite(aUmin) ||
4897         Precision::IsInfinite(aUmax))) {
4898     Standard_Real dU;
4899     //
4900     dU=0.1*(aUmax-aUmin);
4901     aUmin=aUmin-dU;
4902     aUmax=aUmax+dU;
4903   }
4904   if (!(Precision::IsInfinite(aVmin) ||
4905         Precision::IsInfinite(aVmax))) {
4906     Standard_Real dV;
4907     //
4908     dV=0.1*(aVmax-aVmin);
4909     aVmin=aVmin-dV;
4910     aVmax=aVmax+dV;
4911   }
4912 }