0026621: Boolean Cut does not work on two solids
[occt.git] / src / IntTools / IntTools_FaceFace.cxx
1 // Created on: 2000-11-23
2 // Created by: Michael KLOKOV
3 // Copyright (c) 2000-2014 OPEN CASCADE SAS
4 //
5 // This file is part of Open CASCADE Technology software library.
6 //
7 // This library is free software; you can redistribute it and/or modify it under
8 // the terms of the GNU Lesser General Public License version 2.1 as published
9 // by the Free Software Foundation, with special exception defined in the file
10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11 // distribution for complete text of the license and disclaimer of any warranty.
12 //
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
15
16 #include <IntTools_FaceFace.hxx>
17
18 #include <Precision.hxx>
19
20 #include <TColStd_HArray1OfReal.hxx>
21 #include <TColStd_Array1OfReal.hxx>
22 #include <TColStd_Array1OfInteger.hxx>
23 #include <TColStd_SequenceOfReal.hxx>
24 #include <TColStd_ListOfInteger.hxx>
25 #include <TColStd_ListIteratorOfListOfInteger.hxx>
26 #include <TColStd_Array1OfListOfInteger.hxx>
27
28 #include <gp_Lin2d.hxx>
29 #include <gp_Ax22d.hxx>
30 #include <gp_Circ2d.hxx>
31 #include <gp_Torus.hxx>
32 #include <gp_Cylinder.hxx>
33
34 #include <Bnd_Box.hxx>
35
36 #include <TColgp_HArray1OfPnt2d.hxx>
37 #include <TColgp_SequenceOfPnt2d.hxx>
38 #include <TColgp_Array1OfPnt.hxx>
39 #include <TColgp_Array1OfPnt2d.hxx>
40
41 #include <IntAna_QuadQuadGeo.hxx>
42
43 #include <IntSurf_PntOn2S.hxx>
44 #include <IntSurf_LineOn2S.hxx>
45 #include <IntSurf_PntOn2S.hxx>
46 #include <IntSurf_ListOfPntOn2S.hxx>
47 #include <IntRes2d_Domain.hxx>
48 #include <ProjLib_Plane.hxx>
49
50 #include <IntPatch_GLine.hxx>
51 #include <IntPatch_RLine.hxx>
52 #include <IntPatch_WLine.hxx>
53 #include <IntPatch_ALine.hxx>
54 #include <IntPatch_ALineToWLine.hxx>
55
56 #include <ElSLib.hxx>
57 #include <ElCLib.hxx>
58
59 #include <Extrema_ExtCC.hxx>
60 #include <Extrema_POnCurv.hxx>
61 #include <BndLib_AddSurface.hxx>
62
63 #include <Adaptor3d_SurfacePtr.hxx>
64 #include <Adaptor2d_HLine2d.hxx>
65
66 #include <GeomAbs_SurfaceType.hxx>
67 #include <GeomAbs_CurveType.hxx>
68
69 #include <Geom_Surface.hxx>
70 #include <Geom_Line.hxx>
71 #include <Geom_Circle.hxx>
72 #include <Geom_Ellipse.hxx>
73 #include <Geom_Parabola.hxx>
74 #include <Geom_Hyperbola.hxx>
75 #include <Geom_TrimmedCurve.hxx>
76 #include <Geom_BSplineCurve.hxx>
77 #include <Geom_RectangularTrimmedSurface.hxx>
78 #include <Geom_OffsetSurface.hxx>
79 #include <Geom_Curve.hxx>
80 #include <Geom_Conic.hxx>
81
82 #include <Geom2d_TrimmedCurve.hxx>
83 #include <Geom2d_BSplineCurve.hxx>
84 #include <Geom2d_Line.hxx>
85 #include <Geom2d_Curve.hxx>
86 #include <Geom2d_Circle.hxx>
87
88 #include <Geom2dAPI_InterCurveCurve.hxx>
89 #include <Geom2dInt_GInter.hxx>
90 #include <Geom2dAdaptor.hxx>
91 #include <GeomAdaptor_Curve.hxx>
92 #include <GeomAdaptor_HSurface.hxx>
93 #include <GeomAdaptor_Surface.hxx>
94 #include <GeomLib_CheckBSplineCurve.hxx>
95 #include <GeomLib_Check2dBSplineCurve.hxx>
96
97 #include <GeomInt_WLApprox.hxx>
98 #include <GeomProjLib.hxx>
99 #include <GeomAPI_ProjectPointOnSurf.hxx>
100 #include <Geom2dAdaptor_Curve.hxx>
101 #include <TopoDS.hxx>
102 #include <TopoDS_Edge.hxx>
103 #include <TopExp_Explorer.hxx>
104
105 #include <BRep_Tool.hxx>
106 #include <BRepTools.hxx>
107 #include <BRepAdaptor_Surface.hxx>
108
109 #include <IntTools_Curve.hxx>
110 #include <IntTools_Tools.hxx>
111 #include <IntTools_Tools.hxx>
112 #include <IntTools_TopolTool.hxx>
113 #include <IntTools_PntOnFace.hxx>
114 #include <IntTools_PntOn2Faces.hxx>
115 #include <IntTools_Context.hxx>
116 #include <IntSurf_ListIteratorOfListOfPntOn2S.hxx>
117 #include <GeomInt.hxx>
118
119 #include <Approx_CurveOnSurface.hxx>
120 #include <GeomAdaptor.hxx>
121 #include <GeomInt_IntSS.hxx>
122
123 static
124   void RefineVector(gp_Vec2d& aV2D);
125 #ifdef OCCT_DEBUG_DUMPWLINE
126 static
127   void DumpWLine(const Handle(IntPatch_WLine)& aWLine);
128 #endif
129 //
130 static
131   void TolR3d(const TopoDS_Face& ,
132               const TopoDS_Face& ,
133               Standard_Real& );
134 static 
135   Handle(Geom_Curve) MakeBSpline (const Handle(IntPatch_WLine)&,
136                                   const Standard_Integer,
137                                   const Standard_Integer);
138
139 static 
140   void Parameters(const Handle(GeomAdaptor_HSurface)&,
141                   const Handle(GeomAdaptor_HSurface)&,
142                   const gp_Pnt&,
143                   Standard_Real&,
144                   Standard_Real&,
145                   Standard_Real&,
146                   Standard_Real&);
147
148 static 
149   void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
150                                 const Standard_Real theTolerance,
151                                 Standard_Real&      theumin,
152                                 Standard_Real&      theumax, 
153                                 Standard_Real&      thevmin, 
154                                 Standard_Real&      thevmax);
155
156 static
157   Standard_Boolean NotUseSurfacesForApprox
158           (const TopoDS_Face& aF1,
159            const TopoDS_Face& aF2,
160            const Handle(IntPatch_WLine)& WL,
161            const Standard_Integer ifprm,
162            const Standard_Integer ilprm);
163
164 static 
165   Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)       &theWLine,
166                                             const Handle(GeomAdaptor_HSurface) &theS1,
167                                             const Handle(GeomAdaptor_HSurface) &theS2);
168
169 static 
170   Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
171                                             const Standard_Integer                         ideb,
172                                             const Standard_Integer                         ifin,
173                                             const Standard_Boolean                         onFirst);
174
175 static 
176   Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
177                                         const Handle(GeomAdaptor_HSurface)&            theSurface1, 
178                                         const Handle(GeomAdaptor_HSurface)&            theSurface2,
179                                         const TopoDS_Face&                             theFace1,
180                                         const TopoDS_Face&                             theFace2,
181                                         const GeomInt_LineConstructor&                 theLConstructor,
182                                         const Standard_Boolean                         theAvoidLConstructor,
183                                         IntPatch_SequenceOfLine&                       theNewLines,
184                                         Standard_Real&                                 theReachedTol3d,
185                                         const Handle(IntTools_Context)& );
186
187 static 
188   Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter, 
189                                           const Handle(Geom_Curve)& theCurve, 
190                                           const TopoDS_Face&        theFace1, 
191                                           const TopoDS_Face&        theFace2,
192                                           const Standard_Real       theOtherParameter,
193                                           const Standard_Boolean    bIncreasePar,
194                                           Standard_Real&            theNewParameter,
195                                           const Handle(IntTools_Context)& );
196
197 static 
198   Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve);
199
200 static 
201   Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
202                                      const Standard_Real theFirstBoundary,
203                                      const Standard_Real theSecondBoundary,
204                                      const Standard_Real theResolution,
205                                      Standard_Boolean&   IsOnFirstBoundary);
206 static
207   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
208                              const gp_Pnt2d&     theLastPoint,
209                              const Standard_Real theUmin, 
210                              const Standard_Real theUmax,
211                              const Standard_Real theVmin,
212                              const Standard_Real theVmax,
213                              gp_Pnt2d&           theNewPoint);
214
215
216 static 
217   Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
218                                        const Handle(GeomAdaptor_HSurface)&  theSurface2,
219                                        const TopoDS_Face&                   theFace1,
220                                        const TopoDS_Face&                   theFace2,
221                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
222                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
223                                        Handle(TColStd_HArray1OfReal)&       theResultRadius,
224                                        const Handle(IntTools_Context)& );
225
226 static
227   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
228                              const gp_Pnt2d&     theLastPoint,
229                              const Standard_Real theUmin, 
230                              const Standard_Real theUmax,
231                              const Standard_Real theVmin,
232                              const Standard_Real theVmax,
233                              const gp_Pnt2d&     theTanZoneCenter,
234                              const Standard_Real theZoneRadius,
235                              Handle(GeomAdaptor_HSurface) theGASurface,
236                              gp_Pnt2d&           theNewPoint);
237
238 static
239   Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
240                                    const gp_Pnt2d&     theTanZoneCenter,
241                                    const Standard_Real theZoneRadius,
242                                    Handle(GeomAdaptor_HSurface) theGASurface);
243
244 static
245   gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
246                              const gp_Pnt2d&     theOriginalPoint,
247                              Handle(GeomAdaptor_HSurface) theGASurface);
248 static
249   Standard_Boolean  ApproxWithPCurves(const gp_Cylinder& theCyl, 
250                                       const gp_Sphere& theSph);
251
252 static void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1, 
253                            const Handle(GeomAdaptor_HSurface)& theS2, 
254                            const Standard_Real TolAng, 
255                            const Standard_Real TolTang, 
256                            const Standard_Boolean theApprox1,
257                            const Standard_Boolean theApprox2,
258                            IntTools_SequenceOfCurves& theSeqOfCurve, 
259                            Standard_Boolean& theTangentFaces);
260
261 static Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS, 
262                                       const gp_Lin2d& theLin2d, 
263                                       const Standard_Real theTol,
264                                       Standard_Real& theP1, 
265                                       Standard_Real& theP2);
266 //
267 static
268   void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
269                         const Handle(GeomAdaptor_HSurface)& aHS2,
270                         Standard_Integer& iDegMin,
271                         Standard_Integer& iNbIter,
272                         Standard_Integer& iDegMax);
273
274 static
275   void Tolerances(const Handle(GeomAdaptor_HSurface)& aHS1,
276                   const Handle(GeomAdaptor_HSurface)& aHS2,
277                   Standard_Real& aTolTang);
278
279 static
280   Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
281                              const GeomAbs_SurfaceType aType2);
282 static
283   Standard_Integer IndexType(const GeomAbs_SurfaceType aType);
284
285 //
286 static
287   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC, 
288                                const TopoDS_Face& aFace);
289
290 static
291   Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
292                             const Standard_Real aT,
293                             GeomAPI_ProjectPointOnSurf& theProjPS);
294
295 static
296   Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
297                                 const Standard_Real theFirst,
298                                 const Standard_Real theLast,
299                                 GeomAPI_ProjectPointOnSurf& theProjPS,
300                                 const Standard_Real theEps);
301
302 static
303   Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
304                                 const Standard_Real theFirst,
305                                 const Standard_Real theLast,
306                                 const TopoDS_Face& theFace,
307                                 const Handle(IntTools_Context)& theContext);
308
309 static
310   void CorrectPlaneBoundaries(Standard_Real& aUmin,
311                               Standard_Real& aUmax, 
312                               Standard_Real& aVmin, 
313                               Standard_Real& aVmax);
314
315 //=======================================================================
316 //function : 
317 //purpose  : 
318 //=======================================================================
319 IntTools_FaceFace::IntTools_FaceFace()
320 {
321   myIsDone=Standard_False;
322   myTangentFaces=Standard_False;
323   //
324   myHS1 = new GeomAdaptor_HSurface ();
325   myHS2 = new GeomAdaptor_HSurface ();
326   myTolReached2d=0.; 
327   myTolReached3d=0.;
328   SetParameters(Standard_True, Standard_True, Standard_True, 1.e-07);
329   
330 }
331 //=======================================================================
332 //function : SetContext
333 //purpose  : 
334 //======================================================================= 
335 void IntTools_FaceFace::SetContext(const Handle(IntTools_Context)& aContext)
336 {
337   myContext=aContext;
338 }
339 //=======================================================================
340 //function : Context
341 //purpose  : 
342 //======================================================================= 
343 const Handle(IntTools_Context)& IntTools_FaceFace::Context()const
344 {
345   return myContext;
346 }
347 //=======================================================================
348 //function : Face1
349 //purpose  : 
350 //======================================================================= 
351 const TopoDS_Face& IntTools_FaceFace::Face1() const
352 {
353   return myFace1;
354 }
355 //=======================================================================
356 //function : Face2
357 //purpose  : 
358 //======================================================================= 
359 const TopoDS_Face& IntTools_FaceFace::Face2() const
360 {
361   return myFace2;
362 }
363 //=======================================================================
364 //function : TangentFaces
365 //purpose  : 
366 //======================================================================= 
367 Standard_Boolean IntTools_FaceFace::TangentFaces() const
368 {
369   return myTangentFaces;
370 }
371 //=======================================================================
372 //function : Points
373 //purpose  : 
374 //======================================================================= 
375 const IntTools_SequenceOfPntOn2Faces& IntTools_FaceFace::Points() const
376 {
377   return myPnts;
378 }
379 //=======================================================================
380 //function : IsDone
381 //purpose  : 
382 //======================================================================= 
383 Standard_Boolean IntTools_FaceFace::IsDone() const
384 {
385   return myIsDone;
386 }
387 //=======================================================================
388 //function : TolReached3d
389 //purpose  : 
390 //=======================================================================
391 Standard_Real IntTools_FaceFace::TolReached3d() const
392 {
393   return myTolReached3d;
394 }
395 //=======================================================================
396 //function : Lines
397 //purpose  : return lines of intersection
398 //=======================================================================
399 const IntTools_SequenceOfCurves& IntTools_FaceFace::Lines() const
400 {
401   StdFail_NotDone_Raise_if
402     (!myIsDone,
403      "IntTools_FaceFace::Lines() => myIntersector NOT DONE");
404   return mySeqOfCurve;
405 }
406 //=======================================================================
407 //function : TolReached2d
408 //purpose  : 
409 //=======================================================================
410 Standard_Real IntTools_FaceFace::TolReached2d() const
411 {
412   return myTolReached2d;
413 }
414 // =======================================================================
415 // function: SetParameters
416 //
417 // =======================================================================
418 void IntTools_FaceFace::SetParameters(const Standard_Boolean ToApproxC3d,
419                                       const Standard_Boolean ToApproxC2dOnS1,
420                                       const Standard_Boolean ToApproxC2dOnS2,
421                                       const Standard_Real ApproximationTolerance) 
422 {
423   myApprox = ToApproxC3d;
424   myApprox1 = ToApproxC2dOnS1;
425   myApprox2 = ToApproxC2dOnS2;
426   myTolApprox = ApproximationTolerance;
427 }
428 //=======================================================================
429 //function : SetList
430 //purpose  : 
431 //=======================================================================
432 void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
433 {
434   myListOfPnts = aListOfPnts;  
435 }
436
437
438 static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
439                                         const TopoDS_Face& theF2)
440 {
441   const Standard_Real Tolang = 1.e-8;
442   const Standard_Real aTolF1=BRep_Tool::Tolerance(theF1);
443   const Standard_Real aTolF2=BRep_Tool::Tolerance(theF2);
444   const Standard_Real aTolSum = aTolF1 + aTolF2;
445   Standard_Real aHigh = 0.0;
446
447   const BRepAdaptor_Surface aBAS1(theF1), aBAS2(theF2);
448   const GeomAbs_SurfaceType aType1=aBAS1.GetType();
449   const GeomAbs_SurfaceType aType2=aBAS2.GetType();
450   
451   gp_Pln aS1;
452   gp_Cylinder aS2;
453   if(aType1 == GeomAbs_Plane)
454   {
455     aS1=aBAS1.Plane();
456   }
457   else if(aType2 == GeomAbs_Plane)
458   {
459     aS1=aBAS2.Plane();
460   }
461   else
462   {
463     return Standard_True;
464   }
465
466   if(aType1 == GeomAbs_Cylinder)
467   {
468     aS2=aBAS1.Cylinder();
469     const Standard_Real VMin = aBAS1.FirstVParameter();
470     const Standard_Real VMax = aBAS1.LastVParameter();
471
472     if( Precision::IsNegativeInfinite(VMin) ||
473         Precision::IsPositiveInfinite(VMax))
474           return Standard_True;
475     else
476       aHigh = VMax - VMin;
477   }
478   else if(aType2 == GeomAbs_Cylinder)
479   {
480     aS2=aBAS2.Cylinder();
481
482     const Standard_Real VMin = aBAS2.FirstVParameter();
483     const Standard_Real VMax = aBAS2.LastVParameter();
484
485     if( Precision::IsNegativeInfinite(VMin) ||
486         Precision::IsPositiveInfinite(VMax))
487           return Standard_True;
488     else
489       aHigh = VMax - VMin;
490   }
491   else
492   {
493     return Standard_True;
494   }
495
496   IntAna_QuadQuadGeo inter;
497   inter.Perform(aS1,aS2,Tolang,aTolSum, aHigh);
498   if(inter.TypeInter() == IntAna_Ellipse)
499   {
500     const gp_Elips anEl = inter.Ellipse(1);
501     const Standard_Real aMajorR = anEl.MajorRadius();
502     const Standard_Real aMinorR = anEl.MinorRadius();
503     
504     return (aMajorR < 100000.0 * aMinorR);
505   }
506   else
507   {
508     return inter.IsDone();
509   }
510 }
511 //=======================================================================
512 //function : Perform
513 //purpose  : intersect surfaces of the faces
514 //=======================================================================
515 void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
516                                 const TopoDS_Face& aF2)
517 {
518   Standard_Boolean RestrictLine = Standard_False, hasCone = Standard_False;
519   
520   if (myContext.IsNull()) {
521     myContext=new IntTools_Context;
522   }
523
524   mySeqOfCurve.Clear();
525   myTolReached2d=0.;
526   myTolReached3d=0.;
527   myIsDone = Standard_False;
528   myNbrestr=0;//?
529
530   myFace1=aF1;
531   myFace2=aF2;
532
533   const BRepAdaptor_Surface aBAS1(myFace1, Standard_False);
534   const BRepAdaptor_Surface aBAS2(myFace2, Standard_False);
535   GeomAbs_SurfaceType aType1=aBAS1.GetType();
536   GeomAbs_SurfaceType aType2=aBAS2.GetType();
537
538   const Standard_Boolean bReverse=SortTypes(aType1, aType2);
539   if (bReverse)
540   {
541     myFace1=aF2;
542     myFace2=aF1;
543     aType1=aBAS2.GetType();
544     aType2=aBAS1.GetType();
545
546     if (myListOfPnts.Extent())
547     {
548       Standard_Real aU1,aV1,aU2,aV2;
549       IntSurf_ListIteratorOfListOfPntOn2S aItP2S;
550       //
551       aItP2S.Initialize(myListOfPnts);
552       for (; aItP2S.More(); aItP2S.Next())
553       {
554         IntSurf_PntOn2S& aP2S=aItP2S.Value();
555         aP2S.Parameters(aU1,aV1,aU2,aV2);
556         aP2S.SetValue(aU2,aV2,aU1,aV1);
557       }
558     }
559     //
560     Standard_Boolean anAproxTmp = myApprox1;
561     myApprox1 = myApprox2;
562     myApprox2 = anAproxTmp;
563   }
564
565
566   const Handle(Geom_Surface) S1=BRep_Tool::Surface(myFace1);
567   const Handle(Geom_Surface) S2=BRep_Tool::Surface(myFace2);
568
569   const Standard_Real aTolF1=BRep_Tool::Tolerance(myFace1);
570   const Standard_Real aTolF2=BRep_Tool::Tolerance(myFace2);
571
572   Standard_Real TolArc = aTolF1 + aTolF2;
573   Standard_Real TolTang = TolArc;
574
575   const Standard_Boolean isFace1Quad = (aType1 == GeomAbs_Cylinder ||
576                                         aType1 == GeomAbs_Cone ||
577                                         aType1 == GeomAbs_Torus);
578
579   const Standard_Boolean isFace2Quad = (aType2 == GeomAbs_Cylinder ||
580                                         aType2 == GeomAbs_Cone ||
581                                         aType2 == GeomAbs_Torus);
582
583   if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane)  {
584     Standard_Real umin, umax, vmin, vmax;
585     //
586     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
587     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
588     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
589     //
590     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
591     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
592     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
593     //
594     Standard_Real TolAng = 1.e-8;
595     //
596     PerformPlanes(myHS1, myHS2, TolAng, TolTang, myApprox1, myApprox2, 
597                   mySeqOfCurve, myTangentFaces);
598     //
599     myIsDone = Standard_True;
600     
601     if(!myTangentFaces) {
602       const Standard_Integer NbLinPP = mySeqOfCurve.Length();
603       if(NbLinPP) {
604         Standard_Real aTolFMax;
605         myTolReached3d = 1.e-7;
606         aTolFMax=Max(aTolF1, aTolF2);
607         if (aTolFMax>myTolReached3d) {
608           myTolReached3d=aTolFMax;
609         }
610         //
611         myTolReached2d = myTolReached3d;
612
613         if (bReverse) {
614           Handle(Geom2d_Curve) aC2D1, aC2D2;
615           const Standard_Integer aNbLin = mySeqOfCurve.Length();
616           for (Standard_Integer i = 1; i <= aNbLin; ++i) {
617             IntTools_Curve& aIC=mySeqOfCurve(i);
618             aC2D1=aIC.FirstCurve2d();
619             aC2D2=aIC.SecondCurve2d();
620             aIC.SetFirstCurve2d(aC2D2);
621             aIC.SetSecondCurve2d(aC2D1);
622           }
623         }
624       }
625     }
626     return;
627   }//if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane){
628
629   if ((aType1==GeomAbs_Plane) && isFace2Quad)
630   {
631     Standard_Real umin, umax, vmin, vmax;
632     // F1
633     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax); 
634     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
635     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
636     // F2
637     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
638     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
639     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
640     //
641     if( aType2==GeomAbs_Cone ) { 
642       TolArc = 0.0001; 
643       hasCone = Standard_True; 
644     }
645   }
646   else if ((aType2==GeomAbs_Plane) && isFace1Quad)
647   {
648     Standard_Real umin, umax, vmin, vmax;
649     //F1
650     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
651     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
652     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
653     // F2
654     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
655     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
656     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
657     //
658     if( aType1==GeomAbs_Cone ) {
659       TolArc = 0.0001; 
660       hasCone = Standard_True; 
661     }
662   }
663   else
664   {
665     Standard_Real umin, umax, vmin, vmax;
666     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
667     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
668     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
669     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
670     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
671     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
672   }
673
674   const Handle(IntTools_TopolTool) dom1 = new IntTools_TopolTool(myHS1);
675   const Handle(IntTools_TopolTool) dom2 = new IntTools_TopolTool(myHS2);
676
677   myLConstruct.Load(dom1, dom2, myHS1, myHS2);
678   
679
680   Tolerances(myHS1, myHS2, TolTang);
681
682   {
683     const Standard_Real UVMaxStep = 0.001;
684     const Standard_Real Deflection = (hasCone) ? 0.085 : 0.1;
685     myIntersector.SetTolerances(TolArc, TolTang, UVMaxStep, Deflection); 
686   }
687   
688   if((myHS1->IsUClosed() && !myHS1->IsUPeriodic()) || 
689      (myHS1->IsVClosed() && !myHS1->IsVPeriodic()) ||
690      (myHS2->IsUClosed() && !myHS2->IsUPeriodic()) || 
691      (myHS2->IsVClosed() && !myHS2->IsVPeriodic()))
692   {
693     RestrictLine = Standard_True;
694   }
695   //
696   if((aType1 != GeomAbs_BSplineSurface) &&
697       (aType1 != GeomAbs_BezierSurface)  &&
698      (aType1 != GeomAbs_OtherSurface)  &&
699      (aType2 != GeomAbs_BSplineSurface) &&
700       (aType2 != GeomAbs_BezierSurface)  &&
701      (aType2 != GeomAbs_OtherSurface))
702   {
703     RestrictLine = Standard_True;
704
705     if ((aType1 == GeomAbs_Torus) ||
706         (aType2 == GeomAbs_Torus))
707     {
708       myListOfPnts.Clear();
709     }
710   }
711
712   //
713   if(!RestrictLine)
714   {
715     TopExp_Explorer aExp;
716     for(Standard_Integer i = 0; (!RestrictLine) && (i < 2); i++)
717     {
718       const TopoDS_Face& aF=(!i) ? myFace1 : myFace2;
719       aExp.Init(aF, TopAbs_EDGE);
720       for(; aExp.More(); aExp.Next())
721       {
722         const TopoDS_Edge& aE=TopoDS::Edge(aExp.Current());
723
724         if(BRep_Tool::Degenerated(aE))
725         {
726           RestrictLine = Standard_True;
727           break;
728         }
729       }
730     }
731   }
732
733   const Standard_Boolean isGeomInt = isTreatAnalityc(aF1, aF2);
734   myIntersector.Perform(myHS1, dom1, myHS2, dom2, TolArc, TolTang, 
735                                   myListOfPnts, RestrictLine, isGeomInt);
736
737   myIsDone = myIntersector.IsDone();
738
739   if (myIsDone)
740   {
741     myTangentFaces=myIntersector.TangentFaces();
742     if (myTangentFaces) {
743       return;
744     }
745     //
746     if(RestrictLine) {
747       myListOfPnts.Clear(); // to use LineConstructor
748     }
749     //
750     const Standard_Integer aNbLin = myIntersector.NbLines();
751     for (Standard_Integer i=1; i <= aNbLin; ++i) {
752       MakeCurve(i, dom1, dom2);
753     }
754     //
755     ComputeTolReached3d();
756     //
757     if (bReverse) {
758       Handle(Geom2d_Curve) aC2D1, aC2D2;
759       //
760       const Standard_Integer aNbLin=mySeqOfCurve.Length();
761       for (Standard_Integer i=1; i<=aNbLin; ++i)
762       {
763         IntTools_Curve& aIC=mySeqOfCurve(i);
764         aC2D1=aIC.FirstCurve2d();
765         aC2D2=aIC.SecondCurve2d();
766         aIC.SetFirstCurve2d(aC2D2);
767         aIC.SetSecondCurve2d(aC2D1);
768       }
769     }
770
771     // Points
772     Standard_Boolean bValid2D1, bValid2D2;
773     Standard_Real U1,V1,U2,V2;
774     IntTools_PntOnFace aPntOnF1, aPntOnF2;
775     IntTools_PntOn2Faces aPntOn2Faces;
776     //
777     const Standard_Integer aNbPnts = myIntersector.NbPnts();
778     for (Standard_Integer i=1; i <= aNbPnts; ++i)
779     {
780       const IntSurf_PntOn2S& aISPnt=myIntersector.Point(i).PntOn2S();
781       const gp_Pnt& aPnt=aISPnt.Value();
782       aISPnt.Parameters(U1,V1,U2,V2);
783       //
784       // check the validity of the intersection point for the faces
785       bValid2D1 = myContext->IsPointInOnFace(myFace1, gp_Pnt2d(U1, V1));
786       if (!bValid2D1) {
787         continue;
788       }
789       //
790       bValid2D2 = myContext->IsPointInOnFace(myFace2, gp_Pnt2d(U2, V2));
791       if (!bValid2D2) {
792         continue;
793       }
794       //
795       // add the intersection point
796       aPntOnF1.Init(myFace1, aPnt, U1, V1);
797       aPntOnF2.Init(myFace2, aPnt, U2, V2);
798       //
799       if (!bReverse)
800       {
801         aPntOn2Faces.SetP1(aPntOnF1);
802         aPntOn2Faces.SetP2(aPntOnF2);
803       }
804       else
805       {
806         aPntOn2Faces.SetP2(aPntOnF1);
807         aPntOn2Faces.SetP1(aPntOnF2);
808       }
809
810       myPnts.Append(aPntOn2Faces);
811     }
812   }
813 }
814
815 //=======================================================================
816 //function : ComputeTolerance
817 //purpose  : 
818 //=======================================================================
819 Standard_Real IntTools_FaceFace::ComputeTolerance()
820 {
821   Standard_Integer i, j, aNbLin;
822   Standard_Real aFirst, aLast, aD, aDMax, aT;
823   Handle(Geom_Surface) aS1, aS2;
824   //
825   aDMax = 0;
826   aNbLin = mySeqOfCurve.Length();
827   //
828   aS1 = myHS1->ChangeSurface().Surface();
829   aS2 = myHS2->ChangeSurface().Surface();
830   //
831   for (i = 1; i <= aNbLin; ++i)
832   {
833     const IntTools_Curve& aIC = mySeqOfCurve(i);
834     const Handle(Geom_Curve)& aC3D = aIC.Curve();
835     if (aC3D.IsNull())
836     {
837       continue;
838     }
839     //
840     aFirst = aC3D->FirstParameter();
841     aLast  = aC3D->LastParameter();
842     //
843     const Handle(Geom2d_Curve)& aC2D1 = aIC.FirstCurve2d();
844     const Handle(Geom2d_Curve)& aC2D2 = aIC.SecondCurve2d();
845     //
846     for (j = 0; j < 2; ++j)
847     {
848       const Handle(Geom2d_Curve)& aC2D = !j ? aC2D1 : aC2D2;
849       const Handle(Geom_Surface)& aS = !j ? aS1 : aS2;
850       //
851       if (!aC2D.IsNull())
852       {
853         if (IntTools_Tools::ComputeTolerance
854             (aC3D, aC2D, aS, aFirst, aLast, aD, aT))
855         {
856           if (aD > aDMax)
857           {
858             aDMax = aD;
859           }
860         }
861       }
862       else
863       {
864         const TopoDS_Face& aF = !j ? myFace1 : myFace2;
865         aD = FindMaxDistance(aC3D, aFirst, aLast, aF, myContext);
866         if (aD > aDMax)
867         {
868           aDMax = aD;
869         }
870       }
871     }
872   }
873   //
874   return aDMax;
875 }
876
877 //=======================================================================
878 //function :ComputeTolReached3d 
879 //purpose  : 
880 //=======================================================================
881   void IntTools_FaceFace::ComputeTolReached3d()
882 {
883   Standard_Integer aNbLin;
884   GeomAbs_SurfaceType aType1, aType2;
885   //
886   aNbLin=myIntersector.NbLines();
887   if (!aNbLin) {
888     return;
889   }
890   //
891   aType1=myHS1->Surface().GetType();
892   aType2=myHS2->Surface().GetType();
893   //
894   if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder)
895   {
896     if (aNbLin==2)
897     { 
898       Handle(IntPatch_Line) aIL1, aIL2;
899       IntPatch_IType aTL1, aTL2;
900       //
901       aIL1=myIntersector.Line(1);
902       aIL2=myIntersector.Line(2);
903       aTL1=aIL1->ArcType();
904       aTL2=aIL2->ArcType();
905       if (aTL1==IntPatch_Lin && aTL2==IntPatch_Lin) {
906         Standard_Real aD, aDTresh, dTol;
907         gp_Lin aL1, aL2;
908         //
909         dTol=1.e-8;
910         aDTresh=1.5e-6;
911         //
912         aL1=Handle(IntPatch_GLine)::DownCast(aIL1)->Line();
913         aL2=Handle(IntPatch_GLine)::DownCast(aIL2)->Line();
914         aD=aL1.Distance(aL2);
915         aD=0.5*aD;
916         if (aD<aDTresh)
917         {//In order to avoid creation too thin face
918           myTolReached3d=aD+dTol;
919         }
920       }
921     }
922   }// if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
923   //
924
925   Standard_Real aDMax = ComputeTolerance();
926   if (aDMax > myTolReached3d)
927   {
928       myTolReached3d = aDMax;
929     }
930   }
931
932 //=======================================================================
933 //function : MakeCurve
934 //purpose  : 
935 //=======================================================================
936   void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
937                                     const Handle(Adaptor3d_TopolTool)& dom1,
938                                     const Handle(Adaptor3d_TopolTool)& dom2) 
939 {
940   Standard_Boolean bDone, rejectSurface, reApprox, bAvoidLineConstructor;
941   Standard_Boolean ok, bPCurvesOk;
942   Standard_Integer i, j, aNbParts;
943   Standard_Real fprm, lprm;
944   Standard_Real Tolpc;
945   Handle(IntPatch_Line) L;
946   IntPatch_IType typl;
947   Handle(Geom_Curve) newc;
948   //
949   const Standard_Real TOLCHECK   =0.0000001;
950   const Standard_Real TOLANGCHECK=0.1;
951   //
952   rejectSurface = Standard_False;
953   reApprox = Standard_False;
954   //
955   bPCurvesOk = Standard_True;
956  
957  reapprox:;
958   
959   Tolpc = myTolApprox;
960   bAvoidLineConstructor = Standard_False;
961   L = myIntersector.Line(Index);
962   typl = L->ArcType();
963   //
964   if(typl==IntPatch_Walking) {
965     Handle(IntPatch_Line) anewL;
966     //
967     Handle(IntPatch_WLine) aWLine (Handle(IntPatch_WLine)::DownCast(L));
968     anewL = ComputePurgedWLine(aWLine, myHS1, myHS2);
969     if(anewL.IsNull()) {
970       return;
971     }
972     L = anewL;
973
974     //Handle(IntPatch_WLine) aWLineX (Handle(IntPatch_WLine)::DownCast(L));
975     //DumpWLine(aWLineX);
976
977     //
978     if(!myListOfPnts.IsEmpty()) {
979       bAvoidLineConstructor = Standard_True;
980     }
981
982     Standard_Integer nbp = aWLine->NbPnts();
983     const IntSurf_PntOn2S& p1 = aWLine->Point(1);
984     const IntSurf_PntOn2S& p2 = aWLine->Point(nbp);
985
986     const gp_Pnt& P1 = p1.Value();
987     const gp_Pnt& P2 = p2.Value();
988
989     if(P1.SquareDistance(P2) < 1.e-14) {
990       bAvoidLineConstructor = Standard_False;
991     }
992   }
993
994   typl=L->ArcType();
995
996   //
997   // Line Constructor
998   if(!bAvoidLineConstructor) {
999     myLConstruct.Perform(L);
1000     //
1001     bDone=myLConstruct.IsDone();
1002     if(!bDone)
1003     {
1004       return;
1005     }
1006
1007     if(typl != IntPatch_Restriction)
1008     {
1009       aNbParts=myLConstruct.NbParts();
1010       if (aNbParts <= 0)
1011       {
1012         return;
1013       }
1014     }
1015   }
1016   // Do the Curve
1017   
1018   
1019   switch (typl) {
1020   //########################################  
1021   // Line, Parabola, Hyperbola
1022   //########################################  
1023   case IntPatch_Lin:
1024   case IntPatch_Parabola: 
1025   case IntPatch_Hyperbola: {
1026     if (typl == IntPatch_Lin) {
1027       newc = 
1028         new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
1029     }
1030
1031     else if (typl == IntPatch_Parabola) {
1032       newc = 
1033         new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
1034     }
1035     
1036     else if (typl == IntPatch_Hyperbola) {
1037       newc = 
1038         new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
1039     }
1040     //
1041     // myTolReached3d
1042     if (typl == IntPatch_Lin) {
1043       TolR3d (myFace1, myFace2, myTolReached3d);
1044     }
1045     //
1046     aNbParts=myLConstruct.NbParts();
1047     for (i=1; i<=aNbParts; i++) {
1048       Standard_Boolean bFNIt, bLPIt;
1049       //
1050       myLConstruct.Part(i, fprm, lprm);
1051         //
1052       bFNIt=Precision::IsNegativeInfinite(fprm);
1053       bLPIt=Precision::IsPositiveInfinite(lprm);
1054       //
1055       if (!bFNIt && !bLPIt) {
1056         //
1057         IntTools_Curve aCurve;
1058         //
1059         Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
1060         aCurve.SetCurve(aCT3D);
1061         if (typl == IntPatch_Parabola) {
1062           Standard_Real aTolF1, aTolF2, aTolBase;
1063           
1064           aTolF1 = BRep_Tool::Tolerance(myFace1);
1065           aTolF2 = BRep_Tool::Tolerance(myFace2);
1066           aTolBase=aTolF1+aTolF2;
1067           myTolReached3d=IntTools_Tools::CurveTolerance(aCT3D, aTolBase);
1068         }
1069         //
1070         aCurve.SetCurve(new Geom_TrimmedCurve(newc, fprm, lprm));
1071         if(myApprox1) { 
1072           Handle (Geom2d_Curve) C2d;
1073           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1074                 myHS1->ChangeSurface().Surface(), newc, C2d);
1075           if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
1076             myTolReached2d=Tolpc;
1077           }
1078             //            
1079             aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
1080           }
1081           else { 
1082             Handle(Geom2d_BSplineCurve) H1;
1083             //
1084             aCurve.SetFirstCurve2d(H1);
1085           }
1086         //
1087         if(myApprox2) { 
1088           Handle (Geom2d_Curve) C2d;
1089           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1090                     myHS2->ChangeSurface().Surface(), newc, C2d);
1091           if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
1092             myTolReached2d=Tolpc;
1093           }
1094           //
1095           aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
1096           }
1097         else { 
1098           Handle(Geom2d_BSplineCurve) H1;
1099           //
1100           aCurve.SetSecondCurve2d(H1);
1101         }
1102         mySeqOfCurve.Append(aCurve);
1103       } //if (!bFNIt && !bLPIt) {
1104       else {
1105         //  on regarde si on garde
1106         //
1107         Standard_Real aTestPrm, dT=100.;
1108         //
1109         aTestPrm=0.;
1110         if (bFNIt && !bLPIt) {
1111           aTestPrm=lprm-dT;
1112         }
1113         else if (!bFNIt && bLPIt) {
1114           aTestPrm=fprm+dT;
1115         }
1116         else {
1117           // i.e, if (bFNIt && bLPIt)
1118           aTestPrm=IntTools_Tools::IntermediatePoint(-dT, dT);
1119         }
1120         //
1121         gp_Pnt ptref(newc->Value(aTestPrm));
1122         //
1123         GeomAbs_SurfaceType typS1 = myHS1->GetType();
1124         GeomAbs_SurfaceType typS2 = myHS2->GetType();
1125         if( typS1 == GeomAbs_SurfaceOfExtrusion ||
1126             typS1 == GeomAbs_OffsetSurface ||
1127             typS1 == GeomAbs_SurfaceOfRevolution ||
1128             typS2 == GeomAbs_SurfaceOfExtrusion ||
1129             typS2 == GeomAbs_OffsetSurface ||
1130             typS2 == GeomAbs_SurfaceOfRevolution) {
1131           Handle(Geom2d_BSplineCurve) H1;
1132           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1133           continue;
1134         }
1135
1136         Standard_Real u1, v1, u2, v2, Tol;
1137         
1138         Tol = Precision::Confusion();
1139         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
1140         ok = (dom1->Classify(gp_Pnt2d(u1, v1), Tol) != TopAbs_OUT);
1141         if(ok) { 
1142           ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT); 
1143         }
1144         if (ok) {
1145           Handle(Geom2d_BSplineCurve) H1;
1146           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1147         }
1148       }
1149     }// for (i=1; i<=aNbParts; i++) {
1150   }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
1151     break;
1152
1153   //########################################  
1154   // Circle and Ellipse
1155   //########################################  
1156   case IntPatch_Circle: 
1157   case IntPatch_Ellipse: {
1158
1159     if (typl == IntPatch_Circle) {
1160       newc = new Geom_Circle
1161         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
1162     }
1163     else { //IntPatch_Ellipse
1164       newc = new Geom_Ellipse
1165         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
1166     }
1167     //
1168     // myTolReached3d
1169     TolR3d (myFace1, myFace2, myTolReached3d);
1170     //
1171     aNbParts=myLConstruct.NbParts();
1172     //
1173     Standard_Real aPeriod, aNul;
1174     TColStd_SequenceOfReal aSeqFprm,  aSeqLprm;
1175     
1176     aNul=0.;
1177     aPeriod=M_PI+M_PI;
1178
1179     for (i=1; i<=aNbParts; i++) {
1180       myLConstruct.Part(i, fprm, lprm);
1181
1182       if (fprm < aNul && lprm > aNul) {
1183         // interval that goes through 0. is divided on two intervals;
1184         while (fprm<aNul || fprm>aPeriod) fprm=fprm+aPeriod;
1185         while (lprm<aNul || lprm>aPeriod) lprm=lprm+aPeriod;
1186         //
1187         if((aPeriod - fprm) > Tolpc) {
1188           aSeqFprm.Append(fprm);
1189           aSeqLprm.Append(aPeriod);
1190         }
1191         else {
1192           gp_Pnt P1 = newc->Value(fprm);
1193           gp_Pnt P2 = newc->Value(aPeriod);
1194           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1195           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1196
1197           if(P1.Distance(P2) > aTolDist) {
1198             Standard_Real anewpar = fprm;
1199
1200             if(ParameterOutOfBoundary(fprm, newc, myFace1, myFace2, 
1201                                       lprm, Standard_False, anewpar, myContext)) {
1202               fprm = anewpar;
1203             }
1204             aSeqFprm.Append(fprm);
1205             aSeqLprm.Append(aPeriod);
1206           }
1207         }
1208
1209         //
1210         if((lprm - aNul) > Tolpc) {
1211           aSeqFprm.Append(aNul);
1212           aSeqLprm.Append(lprm);
1213         }
1214         else {
1215           gp_Pnt P1 = newc->Value(aNul);
1216           gp_Pnt P2 = newc->Value(lprm);
1217           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1218           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1219
1220           if(P1.Distance(P2) > aTolDist) {
1221             Standard_Real anewpar = lprm;
1222
1223             if(ParameterOutOfBoundary(lprm, newc, myFace1, myFace2, 
1224                                       fprm, Standard_True, anewpar, myContext)) {
1225               lprm = anewpar;
1226             }
1227             aSeqFprm.Append(aNul);
1228             aSeqLprm.Append(lprm);
1229           }
1230         }
1231       }
1232       else {
1233         // usual interval 
1234         aSeqFprm.Append(fprm);
1235         aSeqLprm.Append(lprm);
1236       }
1237     }
1238     
1239     //
1240     aNbParts=aSeqFprm.Length();
1241     for (i=1; i<=aNbParts; i++) {
1242       fprm=aSeqFprm(i);
1243       lprm=aSeqLprm(i);
1244       //
1245       Standard_Real aRealEpsilon=RealEpsilon();
1246       if (Abs(fprm) > aRealEpsilon || Abs(lprm-2.*M_PI) > aRealEpsilon) {
1247         //==============================================
1248         ////
1249         IntTools_Curve aCurve;
1250         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1251         aCurve.SetCurve(aTC3D);
1252         fprm=aTC3D->FirstParameter();
1253         lprm=aTC3D->LastParameter ();
1254         ////         
1255         if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {//// 
1256           if(myApprox1) { 
1257             Handle (Geom2d_Curve) C2d;
1258             GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc, 
1259                         myHS1->ChangeSurface().Surface(), newc, C2d);
1260             if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1261               myTolReached2d=Tolpc;
1262             }
1263             //
1264             aCurve.SetFirstCurve2d(C2d);
1265           }
1266           else { //// 
1267             Handle(Geom2d_BSplineCurve) H1;
1268             aCurve.SetFirstCurve2d(H1);
1269           }
1270
1271
1272           if(myApprox2) { 
1273             Handle (Geom2d_Curve) C2d;
1274             GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1275                     myHS2->ChangeSurface().Surface(),newc,C2d);
1276             if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1277               myTolReached2d=Tolpc;
1278             }
1279             //
1280             aCurve.SetSecondCurve2d(C2d);
1281           }
1282           else { 
1283             Handle(Geom2d_BSplineCurve) H1;
1284             aCurve.SetSecondCurve2d(H1);
1285           }
1286         }
1287         
1288         else { 
1289           Handle(Geom2d_BSplineCurve) H1;
1290           aCurve.SetFirstCurve2d(H1);
1291           aCurve.SetSecondCurve2d(H1);
1292         }
1293         mySeqOfCurve.Append(aCurve);
1294           //==============================================        
1295       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
1296
1297       else {
1298         //  on regarde si on garde
1299         //
1300         if (aNbParts==1) {
1301 //           if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
1302           if (Abs(fprm) <= aRealEpsilon && Abs(lprm-2.*M_PI) <= aRealEpsilon) {
1303             IntTools_Curve aCurve;
1304             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1305             aCurve.SetCurve(aTC3D);
1306             fprm=aTC3D->FirstParameter();
1307             lprm=aTC3D->LastParameter ();
1308             
1309             if(myApprox1) { 
1310               Handle (Geom2d_Curve) C2d;
1311               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1312                     myHS1->ChangeSurface().Surface(),newc,C2d);
1313               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1314                 myTolReached2d=Tolpc;
1315               }
1316               //
1317               aCurve.SetFirstCurve2d(C2d);
1318             }
1319             else { //// 
1320               Handle(Geom2d_BSplineCurve) H1;
1321               aCurve.SetFirstCurve2d(H1);
1322             }
1323
1324             if(myApprox2) { 
1325               Handle (Geom2d_Curve) C2d;
1326               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1327                     myHS2->ChangeSurface().Surface(),newc,C2d);
1328               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1329                 myTolReached2d=Tolpc;
1330               }
1331               //
1332               aCurve.SetSecondCurve2d(C2d);
1333             }
1334             else { 
1335               Handle(Geom2d_BSplineCurve) H1;
1336               aCurve.SetSecondCurve2d(H1);
1337             }
1338             mySeqOfCurve.Append(aCurve);
1339             break;
1340           }
1341         }
1342         //
1343         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, Tol;
1344
1345         aTwoPIdiv17=2.*M_PI/17.;
1346
1347         for (j=0; j<=17; j++) {
1348           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
1349           Tol = Precision::Confusion();
1350
1351           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
1352           ok = (dom1->Classify(gp_Pnt2d(u1,v1),Tol) != TopAbs_OUT);
1353           if(ok) { 
1354             ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1355           }
1356           if (ok) {
1357             IntTools_Curve aCurve;
1358             aCurve.SetCurve(newc);
1359             //==============================================
1360             if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {
1361               
1362               if(myApprox1) { 
1363                 Handle (Geom2d_Curve) C2d;
1364                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc, 
1365                         myHS1->ChangeSurface().Surface(), newc, C2d);
1366                 if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1367                   myTolReached2d=Tolpc;
1368                 }
1369                 // 
1370                 aCurve.SetFirstCurve2d(C2d);
1371               }
1372               else { 
1373                 Handle(Geom2d_BSplineCurve) H1;
1374                 aCurve.SetFirstCurve2d(H1);
1375               }
1376                 
1377               if(myApprox2) { 
1378                 Handle (Geom2d_Curve) C2d;
1379                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1380                         myHS2->ChangeSurface().Surface(), newc, C2d);
1381                 if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1382                   myTolReached2d=Tolpc;
1383                 }
1384                 //                 
1385                 aCurve.SetSecondCurve2d(C2d);
1386               }
1387                 
1388               else { 
1389                 Handle(Geom2d_BSplineCurve) H1;
1390                 aCurve.SetSecondCurve2d(H1);
1391               }
1392             }//  end of if (typl == IntPatch_Circle || typl == IntPatch_Ellipse)
1393              
1394             else { 
1395               Handle(Geom2d_BSplineCurve) H1;
1396               //         
1397               aCurve.SetFirstCurve2d(H1);
1398               aCurve.SetSecondCurve2d(H1);
1399             }
1400             //==============================================        
1401             //
1402             mySeqOfCurve.Append(aCurve);
1403             break;
1404
1405             }//  end of if (ok) {
1406           }//  end of for (Standard_Integer j=0; j<=17; j++)
1407         }//  end of else { on regarde si on garde
1408       }// for (i=1; i<=myLConstruct.NbParts(); i++)
1409     }// IntPatch_Circle: IntPatch_Ellipse:
1410     break;
1411     
1412   case IntPatch_Analytic: {
1413     IntSurf_Quadric quad1,quad2;
1414     GeomAbs_SurfaceType typs = myHS1->Surface().GetType();
1415     
1416     switch (typs) {
1417       case GeomAbs_Plane:
1418         quad1.SetValue(myHS1->Surface().Plane());
1419         break;
1420       case GeomAbs_Cylinder:
1421         quad1.SetValue(myHS1->Surface().Cylinder());
1422         break;
1423       case GeomAbs_Cone:
1424         quad1.SetValue(myHS1->Surface().Cone());
1425         break;
1426       case GeomAbs_Sphere:
1427         quad1.SetValue(myHS1->Surface().Sphere());
1428         break;
1429     case GeomAbs_Torus:
1430       quad1.SetValue(myHS1->Surface().Torus());
1431       break;
1432       default:
1433         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 1");
1434       }
1435       
1436     typs = myHS2->Surface().GetType();
1437     
1438     switch (typs) {
1439       case GeomAbs_Plane:
1440         quad2.SetValue(myHS2->Surface().Plane());
1441         break;
1442       case GeomAbs_Cylinder:
1443         quad2.SetValue(myHS2->Surface().Cylinder());
1444         break;
1445       case GeomAbs_Cone:
1446         quad2.SetValue(myHS2->Surface().Cone());
1447         break;
1448       case GeomAbs_Sphere:
1449         quad2.SetValue(myHS2->Surface().Sphere());
1450         break;
1451     case GeomAbs_Torus:
1452       quad2.SetValue(myHS2->Surface().Torus());
1453       break;
1454       default:
1455         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 2");
1456       }
1457     //
1458     //=========
1459     IntPatch_ALineToWLine convert (quad1, quad2);
1460       
1461     if (!myApprox) {
1462       aNbParts=myLConstruct.NbParts();
1463       for (i=1; i<=aNbParts; i++) {
1464         myLConstruct.Part(i, fprm, lprm);
1465         Handle(IntPatch_WLine) WL = 
1466           convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
1467         //
1468         Handle(Geom2d_BSplineCurve) H1;
1469         Handle(Geom2d_BSplineCurve) H2;
1470
1471         if(myApprox1) {
1472           H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1473         }
1474         
1475         if(myApprox2) {
1476           H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1477         }
1478         //          
1479         mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1480       }
1481     } // if (!myApprox)
1482
1483     else { // myApprox=TRUE
1484       GeomInt_WLApprox theapp3d;
1485       // 
1486       Standard_Real tol2d = myTolApprox;
1487       //         
1488       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
1489       
1490       aNbParts=myLConstruct.NbParts();
1491       for (i=1; i<=aNbParts; i++) {
1492         myLConstruct.Part(i, fprm, lprm);
1493         Handle(IntPatch_WLine) WL = 
1494           convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
1495
1496         theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
1497
1498         if (!theapp3d.IsDone()) {
1499           //
1500           Handle(Geom2d_BSplineCurve) H1;
1501           Handle(Geom2d_BSplineCurve) H2;
1502
1503           if(myApprox1) {
1504             H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1505           }
1506           
1507           if(myApprox2) {
1508             H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1509           }
1510           //          
1511           mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1512         }
1513
1514         else {
1515           if(myApprox1 || myApprox2) { 
1516             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) { 
1517               myTolReached2d = theapp3d.TolReached2d();
1518             }
1519           }
1520           
1521           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) { 
1522             myTolReached3d = theapp3d.TolReached3d();
1523           }
1524
1525           Standard_Integer aNbMultiCurves, nbpoles;
1526           aNbMultiCurves=theapp3d.NbMultiCurves();
1527           for (j=1; j<=aNbMultiCurves; j++) {
1528             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1529             nbpoles = mbspc.NbPoles();
1530             
1531             TColgp_Array1OfPnt tpoles(1, nbpoles);
1532             mbspc.Curve(1, tpoles);
1533             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1534                                                                mbspc.Knots(),
1535                                                                mbspc.Multiplicities(),
1536                                                                mbspc.Degree());
1537             
1538             GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1539             Check.FixTangent(Standard_True,Standard_True);
1540             // 
1541             IntTools_Curve aCurve;
1542             aCurve.SetCurve(BS);
1543             
1544             if(myApprox1) { 
1545               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1546               mbspc.Curve(2,tpoles2d);
1547               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1548                                                                       mbspc.Knots(),
1549                                                                       mbspc.Multiplicities(),
1550                                                                       mbspc.Degree());
1551
1552               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1553               newCheck.FixTangent(Standard_True,Standard_True);
1554               //                 
1555               aCurve.SetFirstCurve2d(BS2);
1556             }
1557             else {
1558               Handle(Geom2d_BSplineCurve) H1;
1559               aCurve.SetFirstCurve2d(H1);
1560             }
1561             
1562             if(myApprox2) { 
1563               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
1564               Standard_Integer TwoOrThree;
1565               TwoOrThree=myApprox1 ? 3 : 2;
1566               mbspc.Curve(TwoOrThree, tpoles2d);
1567               Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
1568                                                                        mbspc.Knots(),
1569                                                                        mbspc.Multiplicities(),
1570                                                                        mbspc.Degree());
1571                 
1572               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1573               newCheck.FixTangent(Standard_True,Standard_True);
1574               //         
1575               aCurve.SetSecondCurve2d(BS2);
1576             }
1577             else { 
1578               Handle(Geom2d_BSplineCurve) H2;
1579               aCurve.SetSecondCurve2d(H2);
1580             }
1581             // 
1582             mySeqOfCurve.Append(aCurve);
1583
1584           }// for (j=1; j<=aNbMultiCurves; j++) {
1585         }// else from if (!theapp3d.IsDone())
1586       }// for (i=1; i<=aNbParts; i++) {
1587     }// else { // myApprox=TRUE
1588   }// case IntPatch_Analytic:
1589     break;
1590
1591   case IntPatch_Walking:{
1592     Handle(IntPatch_WLine) WL = 
1593       Handle(IntPatch_WLine)::DownCast(L);
1594
1595 #ifdef OCCT_DEBUG
1596     //WL->Dump(0);
1597 #endif
1598
1599     //
1600     Standard_Integer ifprm, ilprm;
1601     //
1602     if (!myApprox) {
1603       aNbParts = 1;
1604       if(!bAvoidLineConstructor){
1605         aNbParts=myLConstruct.NbParts();
1606       }
1607       for (i=1; i<=aNbParts; ++i) {
1608         Handle(Geom2d_BSplineCurve) H1, H2;
1609         Handle(Geom_Curve) aBSp;
1610         //
1611         if(bAvoidLineConstructor) {
1612           ifprm = 1;
1613           ilprm = WL->NbPnts();
1614         }
1615         else {
1616           myLConstruct.Part(i, fprm, lprm);
1617           ifprm=(Standard_Integer)fprm;
1618           ilprm=(Standard_Integer)lprm;
1619         }
1620         //
1621         if(myApprox1) {
1622           H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1623         }
1624         //
1625         if(myApprox2) {
1626           H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1627         }
1628         //           
1629         aBSp=MakeBSpline(WL, ifprm, ilprm);
1630         IntTools_Curve aIC(aBSp, H1, H2);
1631         mySeqOfCurve.Append(aIC);
1632       }// for (i=1; i<=aNbParts; ++i) {
1633     }// if (!myApprox) {
1634     //
1635     else { // X
1636       Standard_Boolean bIsDecomposited;
1637       Standard_Integer nbiter, aNbSeqOfL;
1638       Standard_Real tol2d, aTolApproxImp;
1639       IntPatch_SequenceOfLine aSeqOfL;
1640       GeomInt_WLApprox theapp3d;
1641       Approx_ParametrizationType aParType = Approx_ChordLength;
1642       //
1643       Standard_Boolean anApprox1 = myApprox1;
1644       Standard_Boolean anApprox2 = myApprox2;
1645       //
1646       aTolApproxImp=1.e-5;
1647       tol2d = myTolApprox;
1648
1649       GeomAbs_SurfaceType typs1, typs2;
1650       typs1 = myHS1->Surface().GetType();
1651       typs2 = myHS2->Surface().GetType();
1652       Standard_Boolean anWithPC = Standard_True;
1653
1654       if(typs1 == GeomAbs_Cylinder && typs2 == GeomAbs_Sphere) {
1655         anWithPC = 
1656           ApproxWithPCurves(myHS1->Surface().Cylinder(), myHS2->Surface().Sphere());
1657       }
1658       else if (typs1 == GeomAbs_Sphere && typs2 == GeomAbs_Cylinder) {
1659         anWithPC = 
1660           ApproxWithPCurves(myHS2->Surface().Cylinder(), myHS1->Surface().Sphere());
1661       }
1662       //
1663       if(!anWithPC) {
1664         myTolApprox = aTolApproxImp;//1.e-5; 
1665         anApprox1 = Standard_False;
1666         anApprox2 = Standard_False;
1667         //         
1668         tol2d = myTolApprox;
1669       }
1670         
1671       if(myHS1 == myHS2) { 
1672         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1673         rejectSurface = Standard_True;
1674       }
1675       else { 
1676         if(reApprox && !rejectSurface)
1677           theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1678         else {
1679           Standard_Integer iDegMax, iDegMin, iNbIter;
1680           //
1681           ApproxParameters(myHS1, myHS2, iDegMin, iDegMax, iNbIter);
1682           theapp3d.SetParameters(myTolApprox, tol2d, iDegMin, iDegMax, iNbIter, Standard_True, aParType);
1683         }
1684       }
1685       //
1686       Standard_Real aReachedTol = Precision::Confusion();
1687       bIsDecomposited=DecompositionOfWLine(WL,
1688                                            myHS1, 
1689                                            myHS2, 
1690                                            myFace1, 
1691                                            myFace2, 
1692                                            myLConstruct, 
1693                                            bAvoidLineConstructor, 
1694                                            aSeqOfL, 
1695                                            aReachedTol,
1696                                            myContext);
1697       if ( bIsDecomposited && ( myTolReached3d < aReachedTol ) ) {
1698         myTolReached3d = aReachedTol;
1699       }
1700       //
1701       aNbSeqOfL=aSeqOfL.Length();
1702       //
1703       if (bIsDecomposited) {
1704         nbiter=aNbSeqOfL;
1705       }
1706       else {
1707         nbiter=1;
1708         aNbParts=1;
1709         if (!bAvoidLineConstructor) {
1710           aNbParts=myLConstruct.NbParts();
1711           nbiter=aNbParts;
1712         }
1713       }
1714       //
1715       for(i = 1; i <= nbiter; ++i) {
1716         if(bIsDecomposited) {
1717           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
1718           ifprm = 1;
1719           ilprm = WL->NbPnts();
1720         }
1721         else {
1722           if(bAvoidLineConstructor) {
1723             ifprm = 1;
1724             ilprm = WL->NbPnts();
1725           }
1726           else {
1727             myLConstruct.Part(i, fprm, lprm);
1728             ifprm = (Standard_Integer)fprm;
1729             ilprm = (Standard_Integer)lprm;
1730           }
1731         }
1732         //-- lbr : 
1733         //-- Si une des surfaces est un plan , on approxime en 2d
1734         //-- sur cette surface et on remonte les points 2d en 3d.
1735         if(typs1 == GeomAbs_Plane) { 
1736           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,Standard_True, myApprox2,ifprm,ilprm);
1737         }          
1738         else if(typs2 == GeomAbs_Plane) { 
1739           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,myApprox1,Standard_True,ifprm,ilprm);
1740         }
1741         else { 
1742           //
1743           if (myHS1 != myHS2){
1744             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
1745                 (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
1746              
1747               theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True, aParType);
1748               
1749               Standard_Boolean bUseSurfaces;
1750               bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
1751               if (bUseSurfaces) {
1752                 // ######
1753                 rejectSurface = Standard_True;
1754                 // ######
1755                 theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1756               }
1757             }
1758           }
1759           //
1760           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,anApprox1,anApprox2,ifprm,ilprm);
1761         }
1762           //           
1763         if (!theapp3d.IsDone()) {
1764           Handle(Geom2d_BSplineCurve) H1;
1765           Handle(Geom2d_BSplineCurve) H2;
1766           //           
1767           Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
1768           //
1769           if(myApprox1) {
1770             H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1771           }
1772           //
1773           if(myApprox2) {
1774             H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1775           }
1776           //           
1777           IntTools_Curve aIC(aBSp, H1, H2);
1778           mySeqOfCurve.Append(aIC);
1779         }
1780         
1781         else {
1782           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) { 
1783             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) { 
1784               myTolReached2d = theapp3d.TolReached2d();
1785             }
1786           }
1787           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) { 
1788             myTolReached3d = myTolReached2d;
1789             //
1790             if (typs1==GeomAbs_Torus || typs2==GeomAbs_Torus) {
1791               if (myTolReached3d<1.e-6) {
1792                 myTolReached3d = theapp3d.TolReached3d();
1793                 myTolReached3d=1.e-6;
1794               }
1795             }
1796           }
1797           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) { 
1798             myTolReached3d = theapp3d.TolReached3d();
1799           }
1800           
1801           Standard_Integer aNbMultiCurves, nbpoles;
1802           aNbMultiCurves=theapp3d.NbMultiCurves(); 
1803           for (j=1; j<=aNbMultiCurves; j++) {
1804             if(typs1 == GeomAbs_Plane) {
1805               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1806               nbpoles = mbspc.NbPoles();
1807               
1808               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1809               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1810               
1811               mbspc.Curve(1,tpoles2d);
1812               const gp_Pln&  Pln = myHS1->Surface().Plane();
1813               //
1814               Standard_Integer ik; 
1815               for(ik = 1; ik<= nbpoles; ik++) { 
1816                 tpoles.SetValue(ik,
1817                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1818                                               tpoles2d.Value(ik).Y(),
1819                                               Pln));
1820               }
1821               //
1822               Handle(Geom_BSplineCurve) BS = 
1823                 new Geom_BSplineCurve(tpoles,
1824                                       mbspc.Knots(),
1825                                       mbspc.Multiplicities(),
1826                                       mbspc.Degree());
1827               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1828               Check.FixTangent(Standard_True, Standard_True);
1829               //         
1830               IntTools_Curve aCurve;
1831               aCurve.SetCurve(BS);
1832
1833               if(myApprox1) { 
1834                 Handle(Geom2d_BSplineCurve) BS1 = 
1835                   new Geom2d_BSplineCurve(tpoles2d,
1836                                           mbspc.Knots(),
1837                                           mbspc.Multiplicities(),
1838                                           mbspc.Degree());
1839                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1840                 Check1.FixTangent(Standard_True,Standard_True);
1841                 //
1842                 // ############################################
1843                 if(!rejectSurface && !reApprox) {
1844                   Standard_Boolean isValid = IsCurveValid(BS1);
1845                   if(!isValid) {
1846                     reApprox = Standard_True;
1847                     goto reapprox;
1848                   }
1849                 }
1850                 // ############################################
1851                 aCurve.SetFirstCurve2d(BS1);
1852               }
1853               else {
1854                 Handle(Geom2d_BSplineCurve) H1;
1855                 aCurve.SetFirstCurve2d(H1);
1856               }
1857
1858               if(myApprox2) { 
1859                 mbspc.Curve(2, tpoles2d);
1860                 
1861                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
1862                                                                           mbspc.Knots(),
1863                                                                           mbspc.Multiplicities(),
1864                                                                           mbspc.Degree());
1865                 GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1866                 newCheck.FixTangent(Standard_True,Standard_True);
1867                 
1868                 // ###########################################
1869                 if(!rejectSurface && !reApprox) {
1870                   Standard_Boolean isValid = IsCurveValid(BS2);
1871                   if(!isValid) {
1872                     reApprox = Standard_True;
1873                     goto reapprox;
1874                   }
1875                 }
1876                 // ###########################################
1877                 // 
1878                 aCurve.SetSecondCurve2d(BS2);
1879               }
1880               else { 
1881                 Handle(Geom2d_BSplineCurve) H2;
1882                 //                 
1883                   aCurve.SetSecondCurve2d(H2);
1884               }
1885               //
1886               mySeqOfCurve.Append(aCurve);
1887
1888             }//if(typs1 == GeomAbs_Plane) {
1889             
1890             else if(typs2 == GeomAbs_Plane)
1891             {
1892               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1893               nbpoles = mbspc.NbPoles();
1894               
1895               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1896               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1897               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
1898               const gp_Pln&  Pln = myHS2->Surface().Plane();
1899               //
1900               Standard_Integer ik; 
1901               for(ik = 1; ik<= nbpoles; ik++) { 
1902                 tpoles.SetValue(ik,
1903                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1904                                               tpoles2d.Value(ik).Y(),
1905                                               Pln));
1906                 
1907               }
1908               //
1909               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1910                                                                  mbspc.Knots(),
1911                                                                  mbspc.Multiplicities(),
1912                                                                  mbspc.Degree());
1913               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1914               Check.FixTangent(Standard_True,Standard_True);
1915               //         
1916               IntTools_Curve aCurve;
1917               aCurve.SetCurve(BS);
1918
1919               if(myApprox2) {
1920                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
1921                                                                         mbspc.Knots(),
1922                                                                         mbspc.Multiplicities(),
1923                                                                         mbspc.Degree());
1924                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1925                 Check1.FixTangent(Standard_True,Standard_True);
1926                 //         
1927                 // ###########################################
1928                 if(!rejectSurface && !reApprox) {
1929                   Standard_Boolean isValid = IsCurveValid(BS1);
1930                   if(!isValid) {
1931                     reApprox = Standard_True;
1932                     goto reapprox;
1933                   }
1934                 }
1935                 // ###########################################
1936                 bPCurvesOk = CheckPCurve(BS1, myFace2);
1937                 aCurve.SetSecondCurve2d(BS1);
1938               }
1939               else {
1940                 Handle(Geom2d_BSplineCurve) H2;
1941                 aCurve.SetSecondCurve2d(H2);
1942               }
1943               
1944               if(myApprox1) { 
1945                 mbspc.Curve(1,tpoles2d);
1946                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1947                                                                         mbspc.Knots(),
1948                                                                         mbspc.Multiplicities(),
1949                                                                         mbspc.Degree());
1950                 GeomLib_Check2dBSplineCurve Check2(BS2,TOLCHECK,TOLANGCHECK);
1951                 Check2.FixTangent(Standard_True,Standard_True);
1952                 //
1953                 // ###########################################
1954                 if(!rejectSurface && !reApprox) {
1955                   Standard_Boolean isValid = IsCurveValid(BS2);
1956                   if(!isValid) {
1957                     reApprox = Standard_True;
1958                     goto reapprox;
1959                   }
1960                 }
1961                 // ###########################################
1962                 bPCurvesOk = bPCurvesOk && CheckPCurve(BS2, myFace1);
1963                 aCurve.SetFirstCurve2d(BS2);
1964               }
1965               else { 
1966                 Handle(Geom2d_BSplineCurve) H1;
1967                 //                 
1968                 aCurve.SetFirstCurve2d(H1);
1969               }
1970               //
1971               //if points of the pcurves are out of the faces bounds
1972               //create 3d and 2d curves without approximation
1973               if (!bPCurvesOk) {
1974                 Handle(Geom2d_BSplineCurve) H1, H2;
1975                 bPCurvesOk = Standard_True;
1976                 //           
1977                 Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
1978                 
1979                 if(myApprox1) {
1980                   H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1981                   bPCurvesOk = CheckPCurve(H1, myFace1);
1982                 }
1983                 
1984                 if(myApprox2) {
1985                   H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1986                   bPCurvesOk = bPCurvesOk && CheckPCurve(H2, myFace2);
1987                 }
1988                 //
1989                 //if pcurves created without approximation are out of the 
1990                 //faces bounds, use approximated 3d and 2d curves
1991                 if (bPCurvesOk) {
1992                   IntTools_Curve aIC(aBSp, H1, H2);
1993                   mySeqOfCurve.Append(aIC);
1994                 } else {
1995                   mySeqOfCurve.Append(aCurve);
1996                 }
1997               } else {
1998                 mySeqOfCurve.Append(aCurve);
1999               }
2000
2001             }// else if(typs2 == GeomAbs_Plane)
2002             //
2003             else { //typs2 != GeomAbs_Plane && typs1 != GeomAbs_Plane
2004               Standard_Boolean bIsValid1, bIsValid2;
2005               Handle(Geom_BSplineCurve) BS;
2006               Handle(Geom2d_BSplineCurve) aH2D;        
2007               IntTools_Curve aCurve; 
2008               //
2009               bIsValid1=Standard_True;
2010               bIsValid2=Standard_True;
2011               //
2012               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
2013               nbpoles = mbspc.NbPoles();
2014               TColgp_Array1OfPnt tpoles(1,nbpoles);
2015               mbspc.Curve(1,tpoles);
2016               BS=new Geom_BSplineCurve(tpoles,
2017                                                                  mbspc.Knots(),
2018                                                                  mbspc.Multiplicities(),
2019                                                                  mbspc.Degree());
2020               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
2021               Check.FixTangent(Standard_True,Standard_True);
2022               //                 
2023               aCurve.SetCurve(BS);
2024               aCurve.SetFirstCurve2d(aH2D);
2025               aCurve.SetSecondCurve2d(aH2D);
2026               //
2027               if(myApprox1) { 
2028                 if(anApprox1) {
2029                   Handle(Geom2d_BSplineCurve) BS1;
2030                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2031                   mbspc.Curve(2,tpoles2d);
2032                   //
2033                   BS1=new Geom2d_BSplineCurve(tpoles2d,
2034                                                                         mbspc.Knots(),
2035                                                                         mbspc.Multiplicities(),
2036                                                                         mbspc.Degree());
2037                   GeomLib_Check2dBSplineCurve newCheck(BS1,TOLCHECK,TOLANGCHECK);
2038                   newCheck.FixTangent(Standard_True,Standard_True);
2039                   //         
2040                   if (!reApprox) {
2041                     bIsValid1=CheckPCurve(BS1, myFace1);
2042                   }
2043                   //
2044                   aCurve.SetFirstCurve2d(BS1);
2045                 }
2046                 else {
2047                   Handle(Geom2d_BSplineCurve) BS1;
2048                   fprm = BS->FirstParameter();
2049                   lprm = BS->LastParameter();
2050
2051                   Handle(Geom2d_Curve) C2d;
2052                   Standard_Real aTol = myTolApprox;
2053                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
2054                             myHS1->ChangeSurface().Surface(), BS, C2d);
2055                   BS1 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2056                   aCurve.SetFirstCurve2d(BS1);
2057                 }
2058               } // if(myApprox1) { 
2059                 //                 
2060               if(myApprox2) { 
2061                 if(anApprox2) {
2062                   Handle(Geom2d_BSplineCurve) BS2;
2063                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2064                   mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
2065                   BS2=new Geom2d_BSplineCurve(tpoles2d,
2066                                                                         mbspc.Knots(),
2067                                                                         mbspc.Multiplicities(),
2068                                                                         mbspc.Degree());
2069                   GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
2070                   newCheck.FixTangent(Standard_True,Standard_True);
2071                 //                 
2072                   if (!reApprox) {
2073                     bIsValid2=CheckPCurve(BS2, myFace2);        
2074                   }
2075                   aCurve.SetSecondCurve2d(BS2);
2076                 }
2077                 else {
2078                   Handle(Geom2d_BSplineCurve) BS2;
2079                   fprm = BS->FirstParameter();
2080                   lprm = BS->LastParameter();
2081
2082                   Handle(Geom2d_Curve) C2d;
2083                   Standard_Real aTol = myTolApprox;
2084                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
2085                             myHS2->ChangeSurface().Surface(), BS, C2d);
2086                   BS2 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2087                   aCurve.SetSecondCurve2d(BS2);
2088                 }
2089               } //if(myApprox2) { 
2090               if (!bIsValid1 || !bIsValid2) {
2091                 myTolApprox=aTolApproxImp;//1.e-5;
2092                 tol2d = myTolApprox;
2093                 reApprox = Standard_True;
2094                 goto reapprox;
2095               }
2096                 //                 
2097               mySeqOfCurve.Append(aCurve);
2098             }
2099           }
2100         }
2101       }
2102     }// else { // X
2103   }// case IntPatch_Walking:{
2104     break;
2105     
2106   case IntPatch_Restriction: 
2107     {
2108       Handle(IntPatch_RLine) RL = 
2109         Handle(IntPatch_RLine)::DownCast(L);
2110       Handle(Geom_Curve) aC3d;
2111       Handle(Geom2d_Curve) aC2d1, aC2d2;
2112       Standard_Real aTolReached;
2113       GeomInt_IntSS::TreatRLine(RL, myHS1, myHS2, aC3d,
2114                                   aC2d1, aC2d2, aTolReached);
2115
2116       if(aC3d.IsNull())
2117         break;
2118
2119       Bnd_Box2d aBox1, aBox2;
2120
2121       const Standard_Real aU1f = myHS1->FirstUParameter(),
2122                           aV1f = myHS1->FirstVParameter(),
2123                           aU1l = myHS1->LastUParameter(),
2124                           aV1l = myHS1->LastVParameter();
2125       const Standard_Real aU2f = myHS2->FirstUParameter(),
2126                           aV2f = myHS2->FirstVParameter(),
2127                           aU2l = myHS2->LastUParameter(),
2128                           aV2l = myHS2->LastVParameter();
2129
2130       aBox1.Add(gp_Pnt2d(aU1f, aV1f));
2131       aBox1.Add(gp_Pnt2d(aU1l, aV1l));
2132       aBox2.Add(gp_Pnt2d(aU2f, aV2f));
2133       aBox2.Add(gp_Pnt2d(aU2l, aV2l));
2134
2135       GeomInt_VectorOfReal anArrayOfParameters;
2136         
2137       //We consider here that the intersection line is same-parameter-line
2138       anArrayOfParameters.Append(aC3d->FirstParameter());
2139       anArrayOfParameters.Append(aC3d->LastParameter());
2140
2141       GeomInt_IntSS::
2142         TrimILineOnSurfBoundaries(aC2d1, aC2d2, aBox1, aBox2, anArrayOfParameters);
2143
2144       const Standard_Integer aNbIntersSolutionsm1 = anArrayOfParameters.Length() - 1;
2145
2146       //Trim RLine found.
2147       for(Standard_Integer anInd = 0; anInd < aNbIntersSolutionsm1; anInd++)
2148       {
2149         const Standard_Real aParF = anArrayOfParameters(anInd),
2150                             aParL = anArrayOfParameters(anInd+1);
2151
2152         if((aParL - aParF) <= Precision::PConfusion())
2153           continue;
2154
2155         const Standard_Real aPar = 0.5*(aParF + aParL);
2156         gp_Pnt2d aPt;
2157
2158         Handle(Geom2d_Curve) aCurv2d1, aCurv2d2;
2159         if(!aC2d1.IsNull())
2160         {
2161           aC2d1->D0(aPar, aPt);
2162
2163           if(aBox1.IsOut(aPt))
2164             continue;
2165
2166           if(myApprox1)
2167             aCurv2d1 = new Geom2d_TrimmedCurve(aC2d1, aParF, aParL);
2168         }
2169
2170         if(!aC2d2.IsNull())
2171         {
2172           aC2d2->D0(aPar, aPt);
2173
2174           if(aBox2.IsOut(aPt))
2175             continue;
2176
2177           if(myApprox2)
2178             aCurv2d2 = new Geom2d_TrimmedCurve(aC2d2, aParF, aParL);
2179         }
2180
2181         Handle(Geom_Curve) aCurv3d = new Geom_TrimmedCurve(aC3d, aParF, aParL);
2182
2183         IntTools_Curve aIC(aCurv3d, aCurv2d1, aCurv2d2);
2184         mySeqOfCurve.Append(aIC);
2185       }
2186     }
2187     break;
2188   default:
2189     break;
2190
2191   }
2192 }
2193
2194 //=======================================================================
2195 //function : Parameters
2196 //purpose  : 
2197 //=======================================================================
2198  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
2199                  const Handle(GeomAdaptor_HSurface)& HS2,
2200                  const gp_Pnt& Ptref,
2201                  Standard_Real& U1,
2202                  Standard_Real& V1,
2203                  Standard_Real& U2,
2204                  Standard_Real& V2)
2205 {
2206
2207   IntSurf_Quadric quad1,quad2;
2208   GeomAbs_SurfaceType typs = HS1->Surface().GetType();
2209
2210   switch (typs) {
2211   case GeomAbs_Plane:
2212     quad1.SetValue(HS1->Surface().Plane());
2213     break;
2214   case GeomAbs_Cylinder:
2215     quad1.SetValue(HS1->Surface().Cylinder());
2216     break;
2217   case GeomAbs_Cone:
2218     quad1.SetValue(HS1->Surface().Cone());
2219     break;
2220   case GeomAbs_Sphere:
2221     quad1.SetValue(HS1->Surface().Sphere());
2222     break;
2223   case GeomAbs_Torus:
2224     quad1.SetValue(HS1->Surface().Torus());
2225     break;
2226   default:
2227     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2228   }
2229   
2230   typs = HS2->Surface().GetType();
2231   switch (typs) {
2232   case GeomAbs_Plane:
2233     quad2.SetValue(HS2->Surface().Plane());
2234     break;
2235   case GeomAbs_Cylinder:
2236     quad2.SetValue(HS2->Surface().Cylinder());
2237     break;
2238   case GeomAbs_Cone:
2239     quad2.SetValue(HS2->Surface().Cone());
2240     break;
2241   case GeomAbs_Sphere:
2242     quad2.SetValue(HS2->Surface().Sphere());
2243     break;
2244   case GeomAbs_Torus:
2245     quad2.SetValue(HS2->Surface().Torus());
2246     break;
2247   default:
2248     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2249   }
2250
2251   quad1.Parameters(Ptref,U1,V1);
2252   quad2.Parameters(Ptref,U2,V2);
2253 }
2254
2255 //=======================================================================
2256 //function : MakeBSpline
2257 //purpose  : 
2258 //=======================================================================
2259 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
2260                                  const Standard_Integer ideb,
2261                                  const Standard_Integer ifin)
2262 {
2263   Standard_Integer i,nbpnt = ifin-ideb+1;
2264   TColgp_Array1OfPnt poles(1,nbpnt);
2265   TColStd_Array1OfReal knots(1,nbpnt);
2266   TColStd_Array1OfInteger mults(1,nbpnt);
2267   Standard_Integer ipidebm1;
2268   for(i=1,ipidebm1=i+ideb-1; i<=nbpnt;ipidebm1++, i++) {
2269     poles(i) = WL->Point(ipidebm1).Value();
2270     mults(i) = 1;
2271     knots(i) = i-1;
2272   }
2273   mults(1) = mults(nbpnt) = 2;
2274   return
2275     new Geom_BSplineCurve(poles,knots,mults,1);
2276 }
2277 //
2278
2279 //=======================================================================
2280 //function : MakeBSpline2d
2281 //purpose  : 
2282 //=======================================================================
2283 Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
2284                                           const Standard_Integer ideb,
2285                                           const Standard_Integer ifin,
2286                                           const Standard_Boolean onFirst)
2287 {
2288   Standard_Integer i, nbpnt = ifin-ideb+1;
2289   TColgp_Array1OfPnt2d poles(1,nbpnt);
2290   TColStd_Array1OfReal knots(1,nbpnt);
2291   TColStd_Array1OfInteger mults(1,nbpnt);
2292   Standard_Integer ipidebm1;
2293
2294   for(i = 1, ipidebm1 = i+ideb-1; i <= nbpnt; ipidebm1++, i++) {
2295       Standard_Real U, V;
2296       if(onFirst)
2297         theWLine->Point(ipidebm1).ParametersOnS1(U, V);
2298       else
2299         theWLine->Point(ipidebm1).ParametersOnS2(U, V);
2300       poles(i).SetCoord(U, V);
2301       mults(i) = 1;
2302       knots(i) = i-1;
2303     }
2304     mults(1) = mults(nbpnt) = 2;
2305
2306   return new Geom2d_BSplineCurve(poles,knots,mults,1);
2307 }
2308 //=======================================================================
2309 //function : PrepareLines3D
2310 //purpose  : 
2311 //=======================================================================
2312   void IntTools_FaceFace::PrepareLines3D(const Standard_Boolean bToSplit)
2313 {
2314   Standard_Integer i, aNbCurves;
2315   GeomAbs_SurfaceType aType1, aType2;
2316   IntTools_SequenceOfCurves aNewCvs;
2317   //
2318   // 1. Treatment closed  curves
2319   aNbCurves=mySeqOfCurve.Length();
2320   for (i=1; i<=aNbCurves; ++i) {
2321     const IntTools_Curve& aIC=mySeqOfCurve(i);
2322     //
2323     if (bToSplit) {
2324       Standard_Integer j, aNbC;
2325       IntTools_SequenceOfCurves aSeqCvs;
2326       //
2327       aNbC=IntTools_Tools::SplitCurve(aIC, aSeqCvs);
2328       if (aNbC) {
2329         for (j=1; j<=aNbC; ++j) {
2330           const IntTools_Curve& aICNew=aSeqCvs(j);
2331           aNewCvs.Append(aICNew);
2332         }
2333       }
2334       else {
2335         aNewCvs.Append(aIC);
2336       }
2337     }
2338     else {
2339       aNewCvs.Append(aIC);
2340     }
2341   }
2342   //
2343   // 2. Plane\Cone intersection when we had 4 curves
2344   aType1=myHS1->GetType();
2345   aType2=myHS2->GetType();
2346   aNbCurves=aNewCvs.Length();
2347   //
2348   if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone) ||
2349       (aType2==GeomAbs_Plane && aType1==GeomAbs_Cone)) {
2350     if (aNbCurves==4) {
2351       GeomAbs_CurveType aCType1;
2352       //
2353       aCType1=aNewCvs(1).Type();
2354       if (aCType1==GeomAbs_Line) {
2355         IntTools_SequenceOfCurves aSeqIn, aSeqOut;
2356         //
2357         for (i=1; i<=aNbCurves; ++i) {
2358           const IntTools_Curve& aIC=aNewCvs(i);
2359           aSeqIn.Append(aIC);
2360         }
2361         //
2362         IntTools_Tools::RejectLines(aSeqIn, aSeqOut);
2363         //
2364         aNewCvs.Clear();
2365         aNbCurves=aSeqOut.Length(); 
2366         for (i=1; i<=aNbCurves; ++i) {
2367           const IntTools_Curve& aIC=aSeqOut(i);
2368           aNewCvs.Append(aIC);
2369         }
2370       }
2371     }
2372   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone)...
2373   //
2374   // 3. Fill  mySeqOfCurve
2375   mySeqOfCurve.Clear();
2376   aNbCurves=aNewCvs.Length();
2377   for (i=1; i<=aNbCurves; ++i) {
2378     const IntTools_Curve& aIC=aNewCvs(i);
2379     mySeqOfCurve.Append(aIC);
2380   }
2381 }
2382 //=======================================================================
2383 //function : CorrectSurfaceBoundaries
2384 //purpose  : 
2385 //=======================================================================
2386  void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
2387                               const Standard_Real theTolerance,
2388                               Standard_Real&      theumin,
2389                               Standard_Real&      theumax, 
2390                               Standard_Real&      thevmin, 
2391                               Standard_Real&      thevmax) 
2392 {
2393   Standard_Boolean enlarge, isuperiodic, isvperiodic;
2394   Standard_Real uinf, usup, vinf, vsup, delta;
2395   GeomAbs_SurfaceType aType;
2396   Handle(Geom_Surface) aSurface;
2397   //
2398   aSurface = BRep_Tool::Surface(theFace);
2399   aSurface->Bounds(uinf, usup, vinf, vsup);
2400   delta = theTolerance;
2401   enlarge = Standard_False;
2402   //
2403   GeomAdaptor_Surface anAdaptorSurface(aSurface);
2404   //
2405   if(aSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
2406     Handle(Geom_Surface) aBasisSurface = 
2407       (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface))->BasisSurface();
2408     
2409     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2410        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2411       return;
2412     }
2413   }
2414   //
2415   if(aSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2416     Handle(Geom_Surface) aBasisSurface = 
2417       (Handle(Geom_OffsetSurface)::DownCast(aSurface))->BasisSurface();
2418     
2419     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2420        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2421       return;
2422     }
2423   }
2424   //
2425   isuperiodic = anAdaptorSurface.IsUPeriodic();
2426   isvperiodic = anAdaptorSurface.IsVPeriodic();
2427   //
2428   aType=anAdaptorSurface.GetType();
2429   if((aType==GeomAbs_BezierSurface) ||
2430      (aType==GeomAbs_BSplineSurface) ||
2431      (aType==GeomAbs_SurfaceOfExtrusion) ||
2432      (aType==GeomAbs_SurfaceOfRevolution) ||
2433      (aType==GeomAbs_Cylinder)) {
2434     enlarge=Standard_True;
2435   }
2436   //
2437   if(!isuperiodic && enlarge) {
2438
2439     if(!Precision::IsInfinite(theumin) &&
2440         ((theumin - uinf) > delta))
2441       theumin -= delta;
2442     else {
2443       theumin = uinf;
2444     }
2445
2446     if(!Precision::IsInfinite(theumax) &&
2447         ((usup - theumax) > delta))
2448       theumax += delta;
2449     else
2450       theumax = usup;
2451   }
2452   //
2453   if(!isvperiodic && enlarge) {
2454     if(!Precision::IsInfinite(thevmin) &&
2455         ((thevmin - vinf) > delta)) {
2456       thevmin -= delta;
2457     }
2458     else { 
2459       thevmin = vinf;
2460     }
2461     if(!Precision::IsInfinite(thevmax) &&
2462         ((vsup - thevmax) > delta)) {
2463       thevmax += delta;
2464     }
2465     else {
2466       thevmax = vsup;
2467     }
2468   }
2469   //
2470   {
2471     Standard_Integer aNbP;
2472     Standard_Real aXP, dXfact, aXmid, aX1, aX2, aTolPA;
2473     //
2474     aTolPA=Precision::Angular();
2475     // U
2476     if (isuperiodic) {
2477       aXP=anAdaptorSurface.UPeriod();
2478       dXfact=theumax-theumin;
2479       if (dXfact-aTolPA>aXP) {
2480         aXmid=0.5*(theumax+theumin);
2481         aNbP=RealToInt(aXmid/aXP);
2482         if (aXmid<0.) {
2483           aNbP=aNbP-1;
2484         }
2485         aX1=aNbP*aXP;
2486         if (theumin>aTolPA) {
2487           aX1=theumin+aNbP*aXP;
2488         }
2489         aX2=aX1+aXP;
2490         if (theumin<aX1) {
2491           theumin=aX1;
2492         }
2493         if (theumax>aX2) {
2494           theumax=aX2;
2495         }
2496       }
2497     }
2498     // V
2499     if (isvperiodic) {
2500       aXP=anAdaptorSurface.VPeriod();
2501       dXfact=thevmax-thevmin;
2502       if (dXfact-aTolPA>aXP) {
2503         aXmid=0.5*(thevmax+thevmin);
2504         aNbP=RealToInt(aXmid/aXP);
2505         if (aXmid<0.) {
2506           aNbP=aNbP-1;
2507         }
2508         aX1=aNbP*aXP;
2509         if (thevmin>aTolPA) {
2510           aX1=thevmin+aNbP*aXP;
2511         }
2512         aX2=aX1+aXP;
2513         if (thevmin<aX1) {
2514           thevmin=aX1;
2515         }
2516         if (thevmax>aX2) {
2517           thevmax=aX2;
2518         }
2519       }
2520     }
2521   }
2522   //
2523   if(isuperiodic || isvperiodic) {
2524     Standard_Boolean correct = Standard_False;
2525     Standard_Boolean correctU = Standard_False;
2526     Standard_Boolean correctV = Standard_False;
2527     Bnd_Box2d aBox;
2528     TopExp_Explorer anExp;
2529
2530     for(anExp.Init(theFace, TopAbs_EDGE); anExp.More(); anExp.Next()) {
2531       if(BRep_Tool::IsClosed(TopoDS::Edge(anExp.Current()), theFace)) {
2532         correct = Standard_True;
2533         Standard_Real f, l;
2534         TopoDS_Edge anEdge = TopoDS::Edge(anExp.Current());
2535         
2536         for(Standard_Integer i = 0; i < 2; i++) {
2537           if(i==0) {
2538             anEdge.Orientation(TopAbs_FORWARD);
2539           }
2540           else {
2541             anEdge.Orientation(TopAbs_REVERSED);
2542           }
2543           Handle(Geom2d_Curve) aCurve = BRep_Tool::CurveOnSurface(anEdge, theFace, f, l);
2544           
2545           if(aCurve.IsNull()) {
2546             correct = Standard_False;
2547             break;
2548           }
2549           Handle(Geom2d_Line) aLine = Handle(Geom2d_Line)::DownCast(aCurve);
2550
2551           if(aLine.IsNull()) {
2552             correct = Standard_False;
2553             break;
2554           }
2555           gp_Dir2d anUDir(1., 0.);
2556           gp_Dir2d aVDir(0., 1.);
2557           Standard_Real anAngularTolerance = Precision::Angular();
2558
2559           correctU = correctU || aLine->Position().Direction().IsParallel(aVDir, anAngularTolerance);
2560           correctV = correctV || aLine->Position().Direction().IsParallel(anUDir, anAngularTolerance);
2561           
2562           gp_Pnt2d pp1 = aCurve->Value(f);
2563           aBox.Add(pp1);
2564           gp_Pnt2d pp2 = aCurve->Value(l);
2565           aBox.Add(pp2);
2566         }
2567         if(!correct)
2568           break;
2569       }
2570     }
2571
2572     if(correct) {
2573       Standard_Real umin, vmin, umax, vmax;
2574       aBox.Get(umin, vmin, umax, vmax);
2575
2576       if(isuperiodic && correctU) {
2577         
2578         if(theumin < umin)
2579           theumin = umin;
2580         
2581         if(theumax > umax) {
2582           theumax = umax;
2583         }
2584       }
2585       if(isvperiodic && correctV) {
2586         
2587         if(thevmin < vmin)
2588           thevmin = vmin;
2589         if(thevmax > vmax)
2590           thevmax = vmax;
2591       }
2592     }
2593   }
2594 }
2595 //
2596 //
2597 // The block is dedicated to determine whether WLine [ifprm, ilprm]
2598 // crosses the degenerated zone on each given surface or not.
2599 // If Yes -> We will not use info about surfaces during approximation
2600 // because inside degenerated zone of the surface the approx. algo.
2601 // uses wrong values of normal, etc., and resulting curve will have
2602 // oscillations that we would not like to have. 
2603  
2604
2605
2606 static
2607   Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2608                                      const Handle(Geom_Surface)& aS,
2609                                      const Standard_Integer iDir);
2610 static
2611   Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2612                                             const TopoDS_Face& aF1,
2613                                             const TopoDS_Face& aF2);
2614 //=======================================================================
2615 //function :  NotUseSurfacesForApprox
2616 //purpose  : 
2617 //=======================================================================
2618 Standard_Boolean NotUseSurfacesForApprox(const TopoDS_Face& aF1,
2619                                          const TopoDS_Face& aF2,
2620                                          const Handle(IntPatch_WLine)& WL,
2621                                          const Standard_Integer ifprm,
2622                                          const Standard_Integer ilprm)
2623 {
2624   Standard_Boolean bPInDZ;
2625
2626   Handle(IntSurf_LineOn2S) aLineOn2S=WL->Curve();
2627   
2628   const IntSurf_PntOn2S& aP2Sfprm=aLineOn2S->Value(ifprm);
2629   bPInDZ=IsPointInDegeneratedZone(aP2Sfprm, aF1, aF2);
2630   if (bPInDZ) {
2631     return bPInDZ;
2632   }
2633
2634   const IntSurf_PntOn2S& aP2Slprm=aLineOn2S->Value(ilprm);
2635   bPInDZ=IsPointInDegeneratedZone(aP2Slprm, aF1, aF2);
2636   
2637   return bPInDZ;
2638 }
2639 //=======================================================================
2640 //function : IsPointInDegeneratedZone
2641 //purpose  : 
2642 //=======================================================================
2643 Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2644                                           const TopoDS_Face& aF1,
2645                                           const TopoDS_Face& aF2)
2646                                           
2647 {
2648   Standard_Boolean bFlag=Standard_True;
2649   Standard_Real US11, US12, VS11, VS12, US21, US22, VS21, VS22;
2650   Standard_Real U1, V1, U2, V2, aDelta, aD;
2651   gp_Pnt2d aP2d;
2652
2653   Handle(Geom_Surface)aS1 = BRep_Tool::Surface(aF1);
2654   aS1->Bounds(US11, US12, VS11, VS12);
2655   GeomAdaptor_Surface aGAS1(aS1);
2656
2657   Handle(Geom_Surface)aS2 = BRep_Tool::Surface(aF2);
2658   aS1->Bounds(US21, US22, VS21, VS22);
2659   GeomAdaptor_Surface aGAS2(aS2);
2660   //
2661   //const gp_Pnt& aP=aP2S.Value();
2662   aP2S.Parameters(U1, V1, U2, V2);
2663   //
2664   aDelta=1.e-7;
2665   // Check on Surf 1
2666   aD=aGAS1.UResolution(aDelta);
2667   aP2d.SetCoord(U1, V1);
2668   if (fabs(U1-US11) < aD) {
2669     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2670     if (bFlag) {
2671       return bFlag;
2672     }
2673   }
2674   if (fabs(U1-US12) < aD) {
2675     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2676     if (bFlag) {
2677       return bFlag;
2678     }
2679   }
2680   aD=aGAS1.VResolution(aDelta);
2681   if (fabs(V1-VS11) < aDelta) {
2682     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2683     if (bFlag) {
2684       return bFlag;
2685     }
2686   }
2687   if (fabs(V1-VS12) < aDelta) {
2688     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2689     if (bFlag) {
2690       return bFlag;
2691     }
2692   }
2693   // Check on Surf 2
2694   aD=aGAS2.UResolution(aDelta);
2695   aP2d.SetCoord(U2, V2);
2696   if (fabs(U2-US21) < aDelta) {
2697     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2698     if (bFlag) {
2699       return bFlag;
2700     }
2701   }
2702   if (fabs(U2-US22) < aDelta) {
2703     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2704     if (bFlag) {
2705       return bFlag;
2706     }
2707   }
2708   aD=aGAS2.VResolution(aDelta);
2709   if (fabs(V2-VS21) < aDelta) {
2710     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2711     if (bFlag) {  
2712       return bFlag;
2713     }
2714   }
2715   if (fabs(V2-VS22) < aDelta) {
2716     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2717     if (bFlag) {
2718       return bFlag;
2719     }
2720   }
2721   return !bFlag;
2722 }
2723
2724 //=======================================================================
2725 //function : IsDegeneratedZone
2726 //purpose  : 
2727 //=======================================================================
2728 Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2729                                    const Handle(Geom_Surface)& aS,
2730                                    const Standard_Integer iDir)
2731 {
2732   Standard_Boolean bFlag=Standard_True;
2733   Standard_Real US1, US2, VS1, VS2, dY, dX, d1, d2, dD;
2734   Standard_Real aXm, aYm, aXb, aYb, aXe, aYe;
2735   aS->Bounds(US1, US2, VS1, VS2); 
2736
2737   gp_Pnt aPm, aPb, aPe;
2738   
2739   aXm=aP2d.X();
2740   aYm=aP2d.Y();
2741   
2742   aS->D0(aXm, aYm, aPm); 
2743   
2744   dX=1.e-5;
2745   dY=1.e-5;
2746   dD=1.e-12;
2747
2748   if (iDir==1) {
2749     aXb=aXm;
2750     aXe=aXm;
2751     aYb=aYm-dY;
2752     if (aYb < VS1) {
2753       aYb=VS1;
2754     }
2755     aYe=aYm+dY;
2756     if (aYe > VS2) {
2757       aYe=VS2;
2758     }
2759     aS->D0(aXb, aYb, aPb);
2760     aS->D0(aXe, aYe, aPe);
2761     
2762     d1=aPm.Distance(aPb);
2763     d2=aPm.Distance(aPe);
2764     if (d1 < dD && d2 < dD) {
2765       return bFlag;
2766     }
2767     return !bFlag;
2768   }
2769   //
2770   else if (iDir==2) {
2771     aYb=aYm;
2772     aYe=aYm;
2773     aXb=aXm-dX;
2774     if (aXb < US1) {
2775       aXb=US1;
2776     }
2777     aXe=aXm+dX;
2778     if (aXe > US2) {
2779       aXe=US2;
2780     }
2781     aS->D0(aXb, aYb, aPb);
2782     aS->D0(aXe, aYe, aPe);
2783     
2784     d1=aPm.Distance(aPb);
2785     d2=aPm.Distance(aPe);
2786     if (d1 < dD && d2 < dD) {
2787       return bFlag;
2788     }
2789     return !bFlag;
2790   }
2791   return !bFlag;
2792 }
2793
2794 // Check if aNextPnt lies inside of tube build on aBasePnt and aBaseVec.
2795 // In 2d space.
2796 static Standard_Boolean IsInsideIn2d(const gp_Pnt2d& aBasePnt,
2797                                      const gp_Vec2d& aBaseVec,
2798                                      const gp_Pnt2d& aNextPnt,
2799                                      const Standard_Real aSquareMaxDist)
2800 {
2801   gp_Vec2d aVec2d(aBasePnt, aNextPnt);
2802
2803   //d*d = (basevec^(nextpnt-basepnt))**2 / basevec**2
2804   Standard_Real aCross = aVec2d.Crossed(aBaseVec);
2805   Standard_Real aSquareDist = aCross * aCross
2806                             / aBaseVec.SquareMagnitude();
2807
2808   return (aSquareDist <= aSquareMaxDist);
2809 }
2810
2811 // Check if aNextPnt lies inside of tube build on aBasePnt and aBaseVec.
2812 // In 3d space.
2813 static Standard_Boolean IsInsideIn3d(const gp_Pnt& aBasePnt,
2814                                      const gp_Vec& aBaseVec,
2815                                      const gp_Pnt& aNextPnt,
2816                                      const Standard_Real aSquareMaxDist)
2817 {
2818   gp_Vec aVec(aBasePnt, aNextPnt);
2819
2820   //d*d = (basevec^(nextpnt-basepnt))**2 / basevec**2
2821   Standard_Real aSquareDist = aVec.CrossSquareMagnitude(aBaseVec)
2822                             / aBaseVec.SquareMagnitude();
2823
2824   return (aSquareDist <= aSquareMaxDist);
2825 }
2826
2827 //=========================================================================
2828 // static function : ComputePurgedWLine
2829 // purpose : Removes equal points (leave one of equal points) from theWLine
2830 //           and recompute vertex parameters.
2831 //           Removes exceed points using tube criteria:
2832 //           delete 7D point if it lies near to expected lines in 2d and 3d.
2833 //           Each task (2d, 2d, 3d) have its own tolerance and checked separately.
2834 //           Returns new WLine or null WLine if the number
2835 //           of the points is less than 2.
2836 //=========================================================================
2837 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)       &theWLine,
2838                                           const Handle(GeomAdaptor_HSurface) &theS1,
2839                                           const Handle(GeomAdaptor_HSurface) &theS2)
2840 {
2841   Standard_Integer i, k, v, nb, nbvtx;
2842   Handle(IntPatch_WLine) aResult;
2843   nbvtx = theWLine->NbVertex();
2844   nb = theWLine->NbPnts();
2845   if (nb==2) {
2846     const IntSurf_PntOn2S& p1 = theWLine->Point(1);
2847     const IntSurf_PntOn2S& p2 = theWLine->Point(2);
2848     if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
2849       return aResult;
2850     }
2851   }
2852   //
2853   Handle(IntPatch_WLine) aLocalWLine;
2854   Handle(IntPatch_WLine) aTmpWLine = theWLine;
2855   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
2856   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
2857   for(i = 1; i <= nb; i++) {
2858     aLineOn2S->Add(theWLine->Point(i));
2859   }
2860
2861   for(v = 1; v <= nbvtx; v++) {
2862     aLocalWLine->AddVertex(theWLine->Vertex(v));
2863   }
2864   for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
2865     Standard_Integer aStartIndex = i + 1;
2866     Standard_Integer anEndIndex = i + 5;
2867     nb = aLineOn2S->NbPoints();
2868     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
2869
2870     if((aStartIndex > nb) || (anEndIndex <= 1)) {
2871       continue;
2872     }
2873     k = aStartIndex;
2874
2875     while(k <= anEndIndex) {
2876       
2877       if(i != k) {
2878         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
2879         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
2880         
2881         Standard_Real UV[8];
2882         p1.Parameters(UV[0], UV[1], UV[2], UV[3]);
2883         p2.Parameters(UV[4], UV[5], UV[6], UV[7]);
2884
2885         Standard_Real aMax = Abs(UV[0]);
2886         for(Standard_Integer anIdx = 1; anIdx < 8; anIdx++)
2887         {
2888           if (aMax < Abs(UV[anIdx]))
2889             aMax = Abs(UV[anIdx]);
2890         }
2891
2892         if(p1.Value().IsEqual(p2.Value(), gp::Resolution()) ||
2893            Abs(UV[0] - UV[4]) + Abs(UV[1] - UV[5]) < 1.0e-16 * aMax ||
2894            Abs(UV[2] - UV[6]) + Abs(UV[3] - UV[7]) < 1.0e-16 * aMax )
2895         {
2896           aTmpWLine = aLocalWLine;
2897           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
2898           
2899           for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
2900             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
2901             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
2902
2903             if(avertexindex >= k) {
2904               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
2905             }
2906             aLocalWLine->AddVertex(aVertex);
2907           }
2908           aLineOn2S->RemovePoint(k);
2909           anEndIndex--;
2910           continue;
2911         }
2912       }
2913       k++;
2914     }
2915   }
2916
2917   if (aLineOn2S->NbPoints() <= 2)
2918   {
2919     if (aLineOn2S->NbPoints() == 2)
2920       return aLocalWLine;
2921     else
2922       return aResult;
2923   }
2924
2925   const Standard_Integer aMinNbBadDistr = 15;
2926   const Standard_Integer aNbSingleBezier = 30;
2927   // Avoid purge in case of C0 continuity:
2928   // Intersection approximator may produce invalid curve after purge, example:
2929   // bugs modalg_5 bug24731,
2930   // Do not run purger when base number of points is too small.
2931   if (theS1->UContinuity() == GeomAbs_C0 ||
2932       theS1->VContinuity() == GeomAbs_C0 ||
2933       theS2->UContinuity() == GeomAbs_C0 ||
2934       theS2->VContinuity() == GeomAbs_C0 ||
2935       nb < aNbSingleBezier)
2936   {
2937     return aLocalWLine;
2938   }
2939
2940   // 1 - Delete point.
2941   // 0 - Store point.
2942   // -1 - Vertex point (not delete).
2943   NCollection_Array1<Standard_Integer> aNewPointsHash(1, aLineOn2S->NbPoints());
2944   for(i = 1; i <= aLineOn2S->NbPoints(); i++)
2945     aNewPointsHash.SetValue(i, 0);
2946
2947   for(v = 1; v <= aLocalWLine->NbVertex(); v++) 
2948   {
2949     IntPatch_Point aVertex = aLocalWLine->Vertex(v);
2950     Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
2951     aNewPointsHash.SetValue(avertexindex, -1);
2952   }
2953
2954   // Workaround to handle case of small amount points after purge.
2955   // Test "boolean boptuc_complex B5" and similar.
2956   Standard_Integer aNbPnt = 0;
2957
2958   // Inital computations.
2959   Standard_Real UonS1[3], VonS1[3], UonS2[3], VonS2[3];
2960   aLineOn2S->Value(1).ParametersOnS1(UonS1[0], VonS1[0]);
2961   aLineOn2S->Value(2).ParametersOnS1(UonS1[1], VonS1[1]);
2962   aLineOn2S->Value(1).ParametersOnS2(UonS2[0], VonS2[0]);
2963   aLineOn2S->Value(2).ParametersOnS2(UonS2[1], VonS2[1]);
2964
2965   gp_Pnt2d aBase2dPnt1(UonS1[0], VonS1[0]);
2966   gp_Pnt2d aBase2dPnt2(UonS2[0], VonS2[0]);
2967   gp_Vec2d aBase2dVec1(UonS1[1] - UonS1[0], VonS1[1] - VonS1[0]);
2968   gp_Vec2d aBase2dVec2(UonS2[1] - UonS2[0], VonS2[1] - VonS2[0]);
2969   gp_Pnt   aBase3dPnt = aLineOn2S->Value(1).Value();
2970   gp_Vec   aBase3dVec(aLineOn2S->Value(1).Value(), aLineOn2S->Value(2).Value());
2971
2972   // Choose base tolerance and scale it to pipe algorithm.
2973   const Standard_Real aBaseTolerance = Precision::Approximation();
2974   Standard_Real aResS1Tol = Min(theS1->UResolution(aBaseTolerance),
2975                                 theS1->VResolution(aBaseTolerance));
2976   Standard_Real aResS2Tol = Min(theS2->UResolution(aBaseTolerance),
2977                                 theS2->VResolution(aBaseTolerance));
2978   Standard_Real aTol1 = aResS1Tol * aResS1Tol;
2979   Standard_Real aTol2 = aResS2Tol * aResS2Tol;
2980   Standard_Real aTol3d = aBaseTolerance * aBaseTolerance;
2981
2982   const Standard_Real aLimitCoeff = 0.99 * 0.99;
2983   for(i = 3; i <= aLineOn2S->NbPoints(); i++)
2984   {
2985     Standard_Boolean isDeleteState = Standard_False;
2986
2987     aLineOn2S->Value(i).ParametersOnS1(UonS1[2], VonS1[2]);
2988     aLineOn2S->Value(i).ParametersOnS2(UonS2[2], VonS2[2]);
2989     gp_Pnt2d aPnt2dOnS1(UonS1[2], VonS1[2]);
2990     gp_Pnt2d aPnt2dOnS2(UonS2[2], VonS2[2]);
2991     const gp_Pnt& aPnt3d = aLineOn2S->Value(i).Value();
2992
2993     if (aNewPointsHash(i - 1) != - 1 &&
2994         IsInsideIn2d(aBase2dPnt1, aBase2dVec1, aPnt2dOnS1, aTol1) &&
2995         IsInsideIn2d(aBase2dPnt2, aBase2dVec2, aPnt2dOnS2, aTol2) &&
2996         IsInsideIn3d(aBase3dPnt, aBase3dVec, aPnt3d, aTol3d) )
2997     {
2998       // Handle possible uneven parametrization on one of 2d subspaces.
2999       // Delete point only when expected lengths are close to each other (aLimitCoeff).
3000       // Example:
3001       // c2d1 - line
3002       // c3d - line
3003       // c2d2 - geometrically line, but have uneven parametrization -> c2d2 is bspline.
3004       gp_XY aPntOnS1[2]= { gp_XY(UonS1[1] - UonS1[0], VonS1[1] - VonS1[0])
3005                          , gp_XY(UonS1[2] - UonS1[1], VonS1[2] - VonS1[1])};
3006       gp_XY aPntOnS2[2]= { gp_XY(UonS2[1] - UonS2[0], VonS2[1] - VonS2[0])
3007                          , gp_XY(UonS2[2] - UonS2[1], VonS2[2] - VonS2[1])};
3008
3009       Standard_Real aStepOnS1 = aPntOnS1[0].SquareModulus() / aPntOnS1[1].SquareModulus();
3010       Standard_Real aStepOnS2 = aPntOnS2[0].SquareModulus() / aPntOnS2[1].SquareModulus();
3011
3012       Standard_Real aStepCoeff = Min(aStepOnS1, aStepOnS2) / Max(aStepOnS1, aStepOnS2);
3013
3014       if (aStepCoeff > aLimitCoeff)
3015       {
3016         // Set hash flag to "Delete" state.
3017         isDeleteState = Standard_True;
3018         aNewPointsHash.SetValue(i - 1, 1);
3019
3020         // Change middle point.
3021         UonS1[1] = UonS1[2];
3022         UonS2[1] = UonS2[2];
3023         VonS1[1] = VonS1[2];
3024         VonS2[1] = VonS2[2];
3025       }
3026     }
3027
3028     if (!isDeleteState)
3029     {
3030       // Compute new pipe parameters.
3031       UonS1[0] = UonS1[1];
3032       VonS1[0] = VonS1[1];
3033       UonS2[0] = UonS2[1];
3034       VonS2[0] = VonS2[1];
3035
3036       UonS1[1] = UonS1[2];
3037       VonS1[1] = VonS1[2];
3038       UonS2[1] = UonS2[2];
3039       VonS2[1] = VonS2[2];
3040
3041       aBase2dPnt1.SetCoord(UonS1[0], VonS1[0]);
3042       aBase2dPnt2.SetCoord(UonS2[0], VonS2[0]);
3043       aBase2dVec1.SetCoord(UonS1[1] - UonS1[0], VonS1[1] - VonS1[0]);
3044       aBase2dVec2.SetCoord(UonS2[1] - UonS2[0], VonS2[1] - VonS2[0]);
3045       aBase3dPnt = aLineOn2S->Value(i - 1).Value();
3046       aBase3dVec = gp_Vec(aLineOn2S->Value(i - 1).Value(), aLineOn2S->Value(i).Value());
3047
3048       aNbPnt++;
3049     }
3050   }
3051
3052   // Workaround to handle case of small amount of points after purge.
3053   // Test "boolean boptuc_complex B5" and similar.
3054   // This is possible since there are at least two points.
3055   if (aNewPointsHash(1) == -1 &&
3056       aNewPointsHash(2) == -1 &&
3057       aNbPnt <= 3)
3058   {
3059     // Delete first.
3060     aNewPointsHash(1) = 1;
3061   }
3062   if (aNewPointsHash(aLineOn2S->NbPoints() - 1) == -1 &&
3063       aNewPointsHash(aLineOn2S->NbPoints()    ) == -1 &&
3064       aNbPnt <= 3)
3065   {
3066     // Delete last.
3067     aNewPointsHash(aLineOn2S->NbPoints() ) = 1;
3068   }
3069
3070   // Purgre when too small amount of points left.
3071   if (aNbPnt <= 2)
3072   {
3073     for(i = aNewPointsHash.Lower(); i <= aNewPointsHash.Upper(); i++)
3074     {
3075       if (aNewPointsHash(i) != -1)
3076       {
3077         aNewPointsHash(i) = 1;
3078       }
3079     }
3080   }
3081
3082   // Handle possible bad distribution of points, 
3083   // which are will converted into one single bezier curve (less than 30 points).
3084   // Make distribution more even:
3085   // max step will be nearly to 0.1 of param distance.
3086   if (aNbPnt + 2 > aMinNbBadDistr &&
3087       aNbPnt + 2 < aNbSingleBezier )
3088   {
3089     for(Standard_Integer anIdx = 1; anIdx <= 8; anIdx++)
3090     {
3091       Standard_Integer aHashIdx = 
3092         Standard_Integer(anIdx * aLineOn2S->NbPoints() / 9);
3093
3094       //Store this point.
3095       aNewPointsHash(aHashIdx) = 0;
3096     }
3097   }
3098
3099   aTmpWLine = aLocalWLine;
3100   Handle(IntSurf_LineOn2S) aPurgedLineOn2S = new IntSurf_LineOn2S();
3101   aLocalWLine = new IntPatch_WLine(aPurgedLineOn2S, Standard_False);
3102   Standard_Integer anOldLineIdx = 1, aVertexIdx = 1;
3103   for(i = 1; i <= aNewPointsHash.Upper(); i++)
3104   {
3105     if (aNewPointsHash(i) == 0)
3106     {
3107       // Store this point.
3108       aPurgedLineOn2S->Add(aLineOn2S->Value(i));
3109       anOldLineIdx++;
3110     }
3111     else if (aNewPointsHash(i) == -1)
3112     {
3113       // Add vertex.
3114       IntPatch_Point aVertex = aTmpWLine->Vertex(aVertexIdx++);
3115       aVertex.SetParameter(anOldLineIdx++);
3116       aLocalWLine->AddVertex(aVertex);
3117       aPurgedLineOn2S->Add(aLineOn2S->Value(i));
3118     }
3119   }
3120
3121   if(aPurgedLineOn2S->NbPoints() > 1) {
3122     aResult = aLocalWLine;
3123   }
3124   return aResult;
3125 }
3126
3127 //=======================================================================
3128 //function : TolR3d
3129 //purpose  : 
3130 //=======================================================================
3131 void TolR3d(const TopoDS_Face& aF1,
3132             const TopoDS_Face& aF2,
3133             Standard_Real& myTolReached3d)
3134 {
3135   Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
3136       
3137   aTolTresh=2.999999e-3;
3138   aTolF1 = BRep_Tool::Tolerance(aF1);
3139   aTolF2 = BRep_Tool::Tolerance(aF2);
3140   aTolFMax=Max(aTolF1, aTolF2);
3141   
3142   if (aTolFMax>aTolTresh) {
3143     myTolReached3d=aTolFMax;
3144   }
3145 }
3146 //=======================================================================
3147 //function : IsPointOnBoundary
3148 //purpose  : 
3149 //=======================================================================
3150 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
3151                                    const Standard_Real theFirstBoundary,
3152                                    const Standard_Real theSecondBoundary,
3153                                    const Standard_Real theResolution,
3154                                    Standard_Boolean&   IsOnFirstBoundary) 
3155 {
3156   Standard_Boolean bRet;
3157   Standard_Integer i;
3158   Standard_Real adist;
3159   //
3160   bRet=Standard_False;
3161   for(i = 0; i < 2; ++i) {
3162     IsOnFirstBoundary = (i == 0);
3163     if (IsOnFirstBoundary) {
3164       adist = fabs(theParameter - theFirstBoundary);
3165     }
3166     else {
3167       adist = fabs(theParameter - theSecondBoundary);
3168     }
3169     if(adist < theResolution) {
3170       return !bRet;
3171     }
3172   }
3173   return bRet;
3174 }
3175 // ------------------------------------------------------------------------------------------------
3176 // static function: FindPoint
3177 // purpose:
3178 // ------------------------------------------------------------------------------------------------
3179 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
3180                            const gp_Pnt2d&     theLastPoint,
3181                            const Standard_Real theUmin, 
3182                            const Standard_Real theUmax,
3183                            const Standard_Real theVmin,
3184                            const Standard_Real theVmax,
3185                            gp_Pnt2d&           theNewPoint) {
3186   
3187   gp_Vec2d aVec(theFirstPoint, theLastPoint);
3188   Standard_Integer i = 0, j = 0;
3189
3190   for(i = 0; i < 4; i++) {
3191     gp_Vec2d anOtherVec;
3192     gp_Vec2d anOtherVecNormal;
3193     gp_Pnt2d aprojpoint = theLastPoint;    
3194
3195     if((i % 2) == 0) {
3196       anOtherVec.SetX(0.);
3197       anOtherVec.SetY(1.);
3198       anOtherVecNormal.SetX(1.);
3199       anOtherVecNormal.SetY(0.);
3200
3201       if(i < 2)
3202         aprojpoint.SetX(theUmin);
3203       else
3204         aprojpoint.SetX(theUmax);
3205     }
3206     else {
3207       anOtherVec.SetX(1.);
3208       anOtherVec.SetY(0.);
3209       anOtherVecNormal.SetX(0.);
3210       anOtherVecNormal.SetY(1.);
3211
3212       if(i < 2)
3213         aprojpoint.SetY(theVmin);
3214       else
3215         aprojpoint.SetY(theVmax);
3216     }
3217     gp_Vec2d anormvec = aVec;
3218     anormvec.Normalize();
3219     RefineVector(anormvec);
3220     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
3221
3222     if(fabs(adot1) < Precision::Angular())
3223       continue;
3224     Standard_Real adist = 0.;
3225     Standard_Boolean bIsOut = Standard_False;
3226
3227     if((i % 2) == 0) {
3228       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
3229       bIsOut = (i < 2) ? (theLastPoint.X() < theUmin) : (theLastPoint.X() > theUmax);
3230     }
3231     else {
3232       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
3233       bIsOut = (i < 2) ? (theLastPoint.Y() < theVmin) : (theLastPoint.Y() > theVmax);
3234     }
3235     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
3236
3237     for(j = 0; j < 2; j++) {
3238       anoffset = (j == 0) ? anoffset : -anoffset;
3239       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
3240       gp_Vec2d acurvec(theLastPoint, acurpoint);
3241       if ( bIsOut )
3242         acurvec.Reverse();
3243
3244       Standard_Real aDotX, anAngleX;
3245       //
3246       aDotX = aVec.Dot(acurvec);
3247       anAngleX = aVec.Angle(acurvec);
3248       //
3249       if(aDotX > 0. && fabs(anAngleX) < Precision::PConfusion()) {
3250         if((i % 2) == 0) {
3251           if((acurpoint.Y() >= theVmin) &&
3252              (acurpoint.Y() <= theVmax)) {
3253             theNewPoint = acurpoint;
3254             return Standard_True;
3255           }
3256         }
3257         else {
3258           if((acurpoint.X() >= theUmin) &&
3259              (acurpoint.X() <= theUmax)) {
3260             theNewPoint = acurpoint;
3261             return Standard_True;
3262           }
3263         }
3264       }
3265     }
3266   }
3267   return Standard_False;
3268 }
3269
3270
3271 // ------------------------------------------------------------------------------------------------
3272 // static function: FindPoint
3273 // purpose: Find point on the boundary of radial tangent zone
3274 // ------------------------------------------------------------------------------------------------
3275 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
3276                            const gp_Pnt2d&     theLastPoint,
3277                            const Standard_Real theUmin, 
3278                            const Standard_Real theUmax,
3279                            const Standard_Real theVmin,
3280                            const Standard_Real theVmax,
3281                            const gp_Pnt2d&     theTanZoneCenter,
3282                            const Standard_Real theZoneRadius,
3283                            Handle(GeomAdaptor_HSurface) theGASurface,
3284                            gp_Pnt2d&           theNewPoint) {
3285   theNewPoint = theLastPoint;
3286
3287   if ( !IsInsideTanZone( theLastPoint, theTanZoneCenter, theZoneRadius, theGASurface) )
3288     return Standard_False;
3289
3290   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3291   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3292
3293   Standard_Real aRadius = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3294   gp_Ax22d anAxis( theTanZoneCenter, gp_Dir2d(1, 0), gp_Dir2d(0, 1) );
3295   gp_Circ2d aCircle( anAxis, aRadius );
3296   
3297   //
3298   gp_Vec2d aDir( theLastPoint.XY() - theFirstPoint.XY() );
3299   Standard_Real aLength = aDir.Magnitude();
3300   if ( aLength <= gp::Resolution() )
3301     return Standard_False;
3302   gp_Lin2d aLine( theFirstPoint, aDir );
3303
3304   //
3305   Handle(Geom2d_Line) aCLine = new Geom2d_Line( aLine );
3306   Handle(Geom2d_TrimmedCurve) aC1 = new Geom2d_TrimmedCurve( aCLine, 0, aLength );
3307   Handle(Geom2d_Circle) aC2 = new Geom2d_Circle( aCircle );
3308
3309   Standard_Real aTol = aRadius * 0.001;
3310   aTol = ( aTol < Precision::PConfusion() ) ? Precision::PConfusion() : aTol;
3311
3312   Geom2dAPI_InterCurveCurve anIntersector;
3313   anIntersector.Init( aC1, aC2, aTol );
3314
3315   if ( anIntersector.NbPoints() == 0 )
3316     return Standard_False;
3317
3318   Standard_Boolean aFound = Standard_False;
3319   Standard_Real aMinDist = aLength * aLength;
3320   Standard_Integer i = 0;
3321   for ( i = 1; i <= anIntersector.NbPoints(); i++ ) {
3322     gp_Pnt2d aPInt = anIntersector.Point( i );
3323     if ( aPInt.SquareDistance( theFirstPoint ) < aMinDist ) {
3324       if ( ( aPInt.X() >= theUmin ) && ( aPInt.X() <= theUmax ) &&
3325            ( aPInt.Y() >= theVmin ) && ( aPInt.Y() <= theVmax ) ) {
3326         theNewPoint = aPInt;
3327         aFound = Standard_True;
3328       }
3329     }
3330   }
3331
3332   return aFound;
3333 }
3334
3335 // ------------------------------------------------------------------------------------------------
3336 // static function: IsInsideTanZone
3337 // purpose: Check if point is inside a radial tangent zone
3338 // ------------------------------------------------------------------------------------------------
3339 Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
3340                                  const gp_Pnt2d&     theTanZoneCenter,
3341                                  const Standard_Real theZoneRadius,
3342                                  Handle(GeomAdaptor_HSurface) theGASurface) {
3343
3344   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3345   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3346   Standard_Real aRadiusSQR = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3347   aRadiusSQR *= aRadiusSQR;
3348   if ( thePoint.SquareDistance( theTanZoneCenter ) <= aRadiusSQR )
3349     return Standard_True;
3350   return Standard_False;
3351 }
3352
3353 // ------------------------------------------------------------------------------------------------
3354 // static function: CheckTangentZonesExist
3355 // purpose: Check if tangent zone exists
3356 // ------------------------------------------------------------------------------------------------
3357 Standard_Boolean CheckTangentZonesExist( const Handle(GeomAdaptor_HSurface)& theSurface1,
3358                                         const Handle(GeomAdaptor_HSurface)&  theSurface2 ) 
3359 {
3360   if ( ( theSurface1->GetType() != GeomAbs_Torus ) ||
3361       ( theSurface2->GetType() != GeomAbs_Torus ) )
3362     return Standard_False;
3363
3364   gp_Torus aTor1 = theSurface1->Torus();
3365   gp_Torus aTor2 = theSurface2->Torus();
3366
3367   if ( aTor1.Location().Distance( aTor2.Location() ) > Precision::Confusion() )
3368     return Standard_False;
3369
3370   if ( ( fabs( aTor1.MajorRadius() - aTor2.MajorRadius() ) > Precision::Confusion() ) ||
3371        ( fabs( aTor1.MinorRadius() - aTor2.MinorRadius() ) > Precision::Confusion() ) )
3372     return Standard_False;
3373
3374   if ( ( aTor1.MajorRadius() < aTor1.MinorRadius() ) ||
3375        ( aTor2.MajorRadius() < aTor2.MinorRadius() ) )
3376     return Standard_False;
3377   return Standard_True;
3378 }
3379
3380 // ------------------------------------------------------------------------------------------------
3381 // static function: ComputeTangentZones
3382 // purpose: 
3383 // ------------------------------------------------------------------------------------------------
3384 Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
3385                                      const Handle(GeomAdaptor_HSurface)&  theSurface2,
3386                                      const TopoDS_Face&                   theFace1,
3387    &nbs