91832c5e044ec31adce89e02300ca59a37148ace
[occt.git] / src / GeomInt / GeomInt_IntSS_1.cxx
1 // Created on: 1995-01-27
2 // Created by: Jacques GOUSSARD
3 // Copyright (c) 1995-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17
18 #include <Adaptor2d_HCurve2d.hxx>
19 #include <Adaptor3d_TopolTool.hxx>
20 #include <AppParCurves_MultiBSpCurve.hxx>
21 #include <Approx_CurveOnSurface.hxx>
22 #include <Bnd_Box2d.hxx>
23 #include <ElSLib.hxx>
24 #include <Extrema_ExtPS.hxx>
25 #include <Geom2d_BSplineCurve.hxx>
26 #include <Geom2d_Curve.hxx>
27 #include <Geom2d_Line.hxx>
28 #include <Geom2d_TrimmedCurve.hxx>
29 #include <Geom2dAdaptor.hxx>
30 #include <Geom2dAdaptor_Curve.hxx>
31 #include <Geom2dInt_GInter.hxx>
32 #include <Geom_BSplineCurve.hxx>
33 #include <Geom_Circle.hxx>
34 #include <Geom_Curve.hxx>
35 #include <Geom_Ellipse.hxx>
36 #include <Geom_Hyperbola.hxx>
37 #include <Geom_Line.hxx>
38 #include <Geom_Parabola.hxx>
39 #include <Geom_Surface.hxx>
40 #include <Geom_TrimmedCurve.hxx>
41 #include <GeomAbs_CurveType.hxx>
42 #include <GeomAbs_SurfaceType.hxx>
43 #include <GeomAdaptor.hxx>
44 #include <GeomAdaptor_HSurface.hxx>
45 #include <GeomInt.hxx>
46 #include <GeomInt_IntSS.hxx>
47 #include <GeomInt_WLApprox.hxx>
48 #include <GeomLib_Check2dBSplineCurve.hxx>
49 #include <GeomLib_CheckBSplineCurve.hxx>
50 #include <GeomProjLib.hxx>
51 #include <gp_Lin2d.hxx>
52 #include <gp_Pln.hxx>
53 #include <gp_Pnt.hxx>
54 #include <gp_Pnt2d.hxx>
55 #include <IntPatch_ALine.hxx>
56 #include <IntPatch_ALineToWLine.hxx>
57 #include <IntPatch_GLine.hxx>
58 #include <IntPatch_IType.hxx>
59 #include <IntPatch_Line.hxx>
60 #include <IntPatch_RLine.hxx>
61 #include <IntPatch_WLine.hxx>
62 #include <IntRes2d_IntersectionPoint.hxx>
63 #include <IntRes2d_IntersectionSegment.hxx>
64 #include <NCollection_IncAllocator.hxx>
65 #include <NCollection_List.hxx>
66 #include <NCollection_LocalArray.hxx>
67 #include <NCollection_StdAllocator.hxx>
68 #include <Precision.hxx>
69 #include <Standard_OutOfRange.hxx>
70 #include <StdFail_NotDone.hxx>
71 #include <TColgp_Array1OfPnt.hxx>
72 #include <TColgp_Array1OfPnt2d.hxx>
73 #include <TColStd_Array1OfInteger.hxx>
74 #include <TColStd_Array1OfListOfInteger.hxx>
75 #include <TColStd_Array1OfReal.hxx>
76 #include <TColStd_ListIteratorOfListOfInteger.hxx>
77 #include <TColStd_ListOfInteger.hxx>
78 #include <TColStd_SequenceOfReal.hxx>
79
80 #include <algorithm>
81 #include <vector>
82 static
83   void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
84                   const Handle(GeomAdaptor_HSurface)& HS2,
85                   const gp_Pnt& Ptref,
86                   Standard_Real& U1,
87                   Standard_Real& V1,
88                   Standard_Real& U2,
89                   Standard_Real& V2);
90
91 static 
92   Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
93                                    const Standard_Integer ideb,
94                                    const Standard_Integer ifin);
95
96 static
97   Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
98                                             const Standard_Integer                         ideb,
99                                             const Standard_Integer                         ifin,
100                                             const Standard_Boolean                         onFirst);
101
102 static 
103   Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine);
104
105
106 static
107   Standard_Boolean DecompositionOfWLine  (const Handle(IntPatch_WLine)& theWLine,
108                                           const Handle(GeomAdaptor_HSurface)&            theSurface1, 
109                                           const Handle(GeomAdaptor_HSurface)&            theSurface2,
110                                           const Standard_Real aTolSum, 
111                                           const GeomInt_LineConstructor&                 theLConstructor,
112                                           IntPatch_SequenceOfLine&                       theNewLines);
113
114 static
115   Standard_Real AdjustPeriodic(const Standard_Real theParameter,
116                                const Standard_Real parmin,
117                                const Standard_Real parmax,
118                                const Standard_Real thePeriod,
119                                Standard_Real&      theOffset) ;
120 static 
121   Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
122                                      const Standard_Real theFirstBoundary,
123                                      const Standard_Real theSecondBoundary,
124                                      const Standard_Real theResolution,
125                                      Standard_Boolean&   IsOnFirstBoundary);
126
127 static 
128   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
129                              const gp_Pnt2d&     theLastPoint,
130                              const Standard_Real theUmin, 
131                              const Standard_Real theUmax,
132                              const Standard_Real theVmin,
133                              const Standard_Real theVmax,
134                              gp_Pnt2d&           theNewPoint);
135
136
137 static
138   void AdjustUPeriodic (const Handle(Geom_Surface)& aS,
139                         const Handle(Geom2d_Curve)& aC2D);
140 static
141   void GetQuadric(const Handle(GeomAdaptor_HSurface)& HS1,
142                   IntSurf_Quadric& quad1);
143
144 class ProjectPointOnSurf
145 {
146  public:
147   ProjectPointOnSurf() : myIsDone (Standard_False),myIndex(0) {}
148   void Init(const Handle(Geom_Surface)& Surface,
149             const Standard_Real Umin,
150             const Standard_Real Usup,
151             const Standard_Real Vmin,
152             const Standard_Real Vsup);
153   void Init ();
154   void Perform(const gp_Pnt& P);
155   Standard_Boolean IsDone () const { return myIsDone; }
156   void LowerDistanceParameters(Standard_Real& U, Standard_Real& V ) const;
157   Standard_Real LowerDistance() const;
158  protected:
159   Standard_Boolean myIsDone;
160   Standard_Integer myIndex;
161   Extrema_ExtPS myExtPS;
162   GeomAdaptor_Surface myGeomAdaptor;
163 };
164
165 //=======================================================================
166 //function : Init
167 //purpose  : 
168 //=======================================================================
169   void ProjectPointOnSurf::Init (const Handle(Geom_Surface)& Surface,
170                                  const Standard_Real       Umin,
171                                  const Standard_Real       Usup,
172                                  const Standard_Real       Vmin,
173                                  const Standard_Real       Vsup )
174 {
175   const Standard_Real Tolerance = Precision::PConfusion();
176   //
177   myGeomAdaptor.Load(Surface, Umin,Usup,Vmin,Vsup);
178   myExtPS.Initialize(myGeomAdaptor, Umin, Usup, Vmin, Vsup, Tolerance, Tolerance);
179   myIsDone = Standard_False;
180 }
181 //=======================================================================
182 //function : Init
183 //purpose  : 
184 //=======================================================================
185   void ProjectPointOnSurf::Init ()
186 {
187   myIsDone = myExtPS.IsDone() && (myExtPS.NbExt() > 0);
188   if (myIsDone) {
189     // evaluate the lower distance and its index;
190     Standard_Real Dist2Min = myExtPS.SquareDistance(1);
191     myIndex = 1;
192     for (Standard_Integer i = 2; i <= myExtPS.NbExt(); i++)
193     {
194       const Standard_Real Dist2 = myExtPS.SquareDistance(i);
195       if (Dist2 < Dist2Min) {
196         Dist2Min = Dist2;
197         myIndex = i;
198       }
199     }
200   }
201 }
202 //=======================================================================
203 //function : Perform
204 //purpose  : 
205 //=======================================================================
206 void ProjectPointOnSurf::Perform(const gp_Pnt& P)
207 {
208   myExtPS.Perform(P);
209   Init ();
210 }
211 //=======================================================================
212 //function : LowerDistanceParameters
213 //purpose  : 
214 //=======================================================================
215 void ProjectPointOnSurf::LowerDistanceParameters (Standard_Real& U,
216                                                   Standard_Real& V ) const
217 {
218   StdFail_NotDone_Raise_if(!myIsDone, "GeomInt_IntSS::ProjectPointOnSurf::LowerDistanceParameters");
219   (myExtPS.Point(myIndex)).Parameter(U,V);
220 }
221 //=======================================================================
222 //function : LowerDistance
223 //purpose  : 
224 //=======================================================================
225 Standard_Real ProjectPointOnSurf::LowerDistance() const
226 {
227   StdFail_NotDone_Raise_if(!myIsDone, "GeomInt_IntSS::ProjectPointOnSurf::LowerDistance");
228   return sqrt(myExtPS.SquareDistance(myIndex));
229 }
230 //
231
232 //=======================================================================
233 //function : MakeCurve
234 //purpose  : 
235 //=======================================================================
236 void GeomInt_IntSS::MakeCurve(const Standard_Integer Index,
237                               const Handle(Adaptor3d_TopolTool) & dom1,
238                               const Handle(Adaptor3d_TopolTool) & dom2,
239                               const Standard_Real Tol,
240                               const Standard_Boolean Approx,
241                               const Standard_Boolean ApproxS1,
242                               const Standard_Boolean ApproxS2)
243
244 {
245   Standard_Boolean myApprox1, myApprox2, myApprox;
246   Standard_Real Tolpc, myTolApprox;
247   IntPatch_IType typl;
248   Handle(Geom2d_BSplineCurve) H1;
249   Handle(Geom_Surface) aS1, aS2;
250   //
251   Tolpc = Tol;
252   myApprox=Approx;
253   myApprox1=ApproxS1;
254   myApprox2=ApproxS2;
255   myTolApprox=0.0000001;
256   //
257   aS1=myHS1->ChangeSurface().Surface();
258   aS2=myHS2->ChangeSurface().Surface();
259   //
260   Handle(IntPatch_Line) L = myIntersector.Line(Index);
261   typl = L->ArcType();
262   //
263   if(typl==IntPatch_Walking) {
264     Handle(IntPatch_Line) anewL = 
265       ComputePurgedWLine(Handle(IntPatch_WLine)::DownCast(L));
266     //
267     if(anewL.IsNull()) {
268       return;
269     }
270     L = anewL;
271   }
272   //
273   // Line Constructor
274   myLConstruct.Perform(L);
275   if (!myLConstruct.IsDone() || myLConstruct.NbParts() <= 0) {
276     return;
277   }
278   // Do the Curve
279   Standard_Boolean ok;
280   Standard_Integer i, j,  aNbParts;
281   Standard_Real fprm, lprm;
282   Handle(Geom_Curve) newc;
283
284   switch (typl) {
285     //########################################  
286     // Line, Parabola, Hyperbola
287     //########################################  
288   case IntPatch_Lin:
289   case IntPatch_Parabola: 
290   case IntPatch_Hyperbola: {
291     if (typl == IntPatch_Lin) {
292       newc=new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
293     }
294     else if (typl == IntPatch_Parabola) {
295       newc=new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
296     }
297     else if (typl == IntPatch_Hyperbola) {
298       newc=new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
299     }
300     //
301     aNbParts=myLConstruct.NbParts();
302     for (i=1; i<=aNbParts; i++) {
303       myLConstruct.Part(i, fprm, lprm);
304
305       if (!Precision::IsNegativeInfinite(fprm) && 
306         !Precision::IsPositiveInfinite(lprm)) {
307           Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
308           sline.Append(aCT3D);
309           //
310           if(myApprox1) { 
311             Handle (Geom2d_Curve) C2d;
312             BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
313             if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
314               myTolReached2d=Tolpc;
315             }
316             slineS1.Append(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
317           }
318           else { 
319             slineS1.Append(H1);
320           }
321           //
322           if(myApprox2) { 
323             Handle (Geom2d_Curve) C2d;
324             BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
325             if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
326               myTolReached2d=Tolpc;
327             }
328             //
329             slineS2.Append(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
330           }
331           else { 
332             slineS2.Append(H1);
333           }
334       } // if (!Precision::IsNegativeInfinite(fprm) &&  !Precision::IsPositiveInfinite(lprm))
335
336       else {
337         GeomAbs_SurfaceType typS1 = myHS1->Surface().GetType();
338         GeomAbs_SurfaceType typS2 = myHS2->Surface().GetType();
339         if( typS1 == GeomAbs_SurfaceOfExtrusion ||
340           typS1 == GeomAbs_OffsetSurface ||
341           typS1 == GeomAbs_SurfaceOfRevolution ||
342           typS2 == GeomAbs_SurfaceOfExtrusion ||
343           typS2 == GeomAbs_OffsetSurface ||
344           typS2 == GeomAbs_SurfaceOfRevolution) {
345             sline.Append(newc);
346             slineS1.Append(H1);
347             slineS2.Append(H1);
348             continue;
349         }
350         Standard_Boolean bFNIt, bLPIt;
351         Standard_Real aTestPrm, dT=100.;
352         Standard_Real u1, v1, u2, v2, TolX;
353         //
354         bFNIt=Precision::IsNegativeInfinite(fprm);
355         bLPIt=Precision::IsPositiveInfinite(lprm);
356
357         aTestPrm=0.;
358
359         if (bFNIt && !bLPIt) {
360           aTestPrm=lprm-dT;
361         }
362         else if (!bFNIt && bLPIt) {
363           aTestPrm=fprm+dT;
364         }
365         //
366         gp_Pnt ptref(newc->Value(aTestPrm));
367         //
368         TolX = Precision::Confusion();
369         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
370         ok = (dom1->Classify(gp_Pnt2d(u1, v1), TolX) != TopAbs_OUT);
371         if(ok) { 
372           ok = (dom2->Classify(gp_Pnt2d(u2,v2),TolX) != TopAbs_OUT); 
373         }
374         if (ok) {
375           sline.Append(newc);
376           slineS1.Append(H1);
377           slineS2.Append(H1);
378         }
379       }
380     }// end of for (i=1; i<=myLConstruct.NbParts(); i++)
381   }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
382   break;
383
384                            //########################################  
385                            // Circle and Ellipse
386                            //########################################  
387   case IntPatch_Circle: 
388   case IntPatch_Ellipse: {
389
390     if (typl == IntPatch_Circle) {
391       newc = new Geom_Circle
392         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
393     }
394     else { 
395       newc = new Geom_Ellipse
396         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
397     }
398     //
399     Standard_Real aPeriod, aRealEpsilon;
400     //
401     aRealEpsilon=RealEpsilon();
402     aPeriod=M_PI+M_PI;
403     //
404     aNbParts=myLConstruct.NbParts();
405     //
406     for (i=1; i<=aNbParts; i++) {
407       myLConstruct.Part(i, fprm, lprm);
408       //
409       if (Abs(fprm) > aRealEpsilon || Abs(lprm-aPeriod) > aRealEpsilon) {
410         //==============================================
411         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
412         //
413         sline.Append(aTC3D);
414         //
415         fprm=aTC3D->FirstParameter();
416         lprm=aTC3D->LastParameter ();
417         ////    
418         if(myApprox1) { 
419           Handle (Geom2d_Curve) C2d;
420           BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
421           if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
422             myTolReached2d=Tolpc;
423           }
424           slineS1.Append(C2d);
425         }
426         else { //// 
427           slineS1.Append(H1);
428         }
429         //
430         if(myApprox2) { 
431           Handle (Geom2d_Curve) C2d;
432           BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
433           if(Tolpc>myTolReached2d || myTolReached2d==0) { 
434             myTolReached2d=Tolpc;
435           }
436           slineS2.Append(C2d);
437         }
438         else { 
439           slineS2.Append(H1);  
440         }
441         //==============================================        
442       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
443       //
444       else {//  on regarde si on garde
445         //
446         if (aNbParts==1) {
447           if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
448             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
449             //
450             sline.Append(aTC3D);
451             fprm=aTC3D->FirstParameter();
452             lprm=aTC3D->LastParameter ();
453
454             if(myApprox1) { 
455               Handle (Geom2d_Curve) C2d;
456               BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
457               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
458                 myTolReached2d=Tolpc;
459               }
460               slineS1.Append(C2d);
461             }
462             else { //// 
463               slineS1.Append(H1);
464             }
465
466             if(myApprox2) { 
467               Handle (Geom2d_Curve) C2d;
468               BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
469               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
470                 myTolReached2d=Tolpc;
471               }
472               slineS2.Append(C2d);
473             }
474             else { 
475               slineS2.Append(H1);
476             }
477             break;
478           }
479         }
480         //
481         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, TolX;
482         //
483         aTwoPIdiv17=2.*M_PI/17.;
484         //
485         for (j=0; j<=17; j++) {
486           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
487           TolX = Precision::Confusion();
488
489           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
490           ok = (dom1->Classify(gp_Pnt2d(u1,v1),TolX) != TopAbs_OUT);
491           if(ok) { 
492             ok = (dom2->Classify(gp_Pnt2d(u2,v2),TolX) != TopAbs_OUT);
493           }
494           if (ok) {
495             sline.Append(newc);
496             //==============================================
497             if(myApprox1) { 
498               Handle (Geom2d_Curve) C2d;
499               BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
500               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
501                 myTolReached2d=Tolpc;
502               }
503               slineS1.Append(C2d);
504             }
505             else { 
506               slineS1.Append(H1);  
507             }
508
509             if(myApprox2) { 
510               Handle (Geom2d_Curve) C2d;
511               BuildPCurves(fprm, lprm, Tolpc,myHS2->ChangeSurface().Surface(), newc, C2d);
512               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
513                 myTolReached2d=Tolpc;
514               } 
515               slineS2.Append(C2d);
516             }
517             else { 
518               slineS2.Append(H1);  
519             }
520             break;
521           }//  end of if (ok) {
522         }//  end of for (Standard_Integer j=0; j<=17; j++)
523       }//  end of else { on regarde si on garde
524     }// for (i=1; i<=myLConstruct.NbParts(); i++)
525   }// IntPatch_Circle: IntPatch_Ellipse
526   break;
527
528                          //########################################  
529                          // Analytic
530                          //######################################## 
531   case IntPatch_Analytic: {
532     IntSurf_Quadric quad1,quad2;
533     //
534     GetQuadric(myHS1, quad1);
535     GetQuadric(myHS2, quad2);
536     //=========
537     IntPatch_ALineToWLine convert (quad1, quad2);
538
539     if (!myApprox) {
540       Handle(Geom2d_BSplineCurve) aH1, aH2;
541       //
542       aNbParts=myLConstruct.NbParts();
543       for (i=1; i<=aNbParts; i++) {
544         myLConstruct.Part(i, fprm, lprm);
545         Handle(IntPatch_WLine) WL = 
546           convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
547         //
548         if(myApprox1) {
549           aH1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
550         }
551
552         if(myApprox2) {
553           aH2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
554         }
555         sline.Append(MakeBSpline(WL,1,WL->NbPnts()));
556         slineS1.Append(aH1);
557         slineS2.Append(aH2);
558       }
559     } // if (!myApprox)
560
561     else { // myApprox=TRUE
562       GeomInt_WLApprox theapp3d;
563       Standard_Real tol2d = myTolApprox;
564       //        
565       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
566
567       aNbParts=myLConstruct.NbParts();
568       for (i=1; i<=aNbParts; i++) {
569         myLConstruct.Part(i, fprm, lprm);
570         Handle(IntPatch_WLine) WL = 
571           convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
572
573         theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
574         if (!theapp3d.IsDone()) {
575           //
576           Handle(Geom2d_BSplineCurve) aH1, aH2;
577
578           if(myApprox1) {
579             aH1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
580           }
581
582           if(myApprox2) {
583             aH2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
584           }
585           sline.Append(MakeBSpline(WL,1,WL->NbPnts()));
586           slineS1.Append(aH1);
587           slineS2.Append(aH1);
588         }
589         //
590         else {
591           if(myApprox1 || myApprox2) { 
592             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) { 
593               myTolReached2d = theapp3d.TolReached2d();
594             }
595           }
596
597           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) { 
598             myTolReached3d = theapp3d.TolReached3d();
599           }
600
601           Standard_Integer aNbMultiCurves, nbpoles;
602           aNbMultiCurves=theapp3d.NbMultiCurves();
603           for (j=1; j<=aNbMultiCurves; j++) {
604             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
605
606             nbpoles = mbspc.NbPoles();
607             TColgp_Array1OfPnt tpoles(1, nbpoles);
608             mbspc.Curve(1, tpoles);
609             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
610               mbspc.Knots(),
611               mbspc.Multiplicities(),
612               mbspc.Degree());
613
614             GeomLib_CheckBSplineCurve Check(BS, myTolCheck, myTolAngCheck);
615             Check.FixTangent(Standard_True,Standard_True);
616             // 
617             sline.Append(BS);
618             //
619             if(myApprox1) { 
620               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
621               mbspc.Curve(2,tpoles2d);
622               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
623                 mbspc.Knots(),
624                 mbspc.Multiplicities(),
625                 mbspc.Degree());
626
627               GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
628               newCheck.FixTangent(Standard_True,Standard_True);
629               slineS1.Append(BS2);              
630             }
631             else {
632               slineS1.Append(H1);
633             }
634
635             if(myApprox2) { 
636               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
637               Standard_Integer TwoOrThree;
638               TwoOrThree=myApprox1 ? 3 : 2;
639               mbspc.Curve(TwoOrThree, tpoles2d);
640               Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
641                 mbspc.Knots(),
642                 mbspc.Multiplicities(),
643                 mbspc.Degree());
644
645               GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
646               newCheck.FixTangent(Standard_True,Standard_True);
647               //        
648               slineS2.Append(BS2);
649             }
650             else { 
651               slineS2.Append(H1);
652             }
653             // 
654           }// for (j=1; j<=aNbMultiCurves; j++) {
655         }// else from if (!theapp3d.IsDone())
656       }// for (i=1; i<=aNbParts; i++) {
657     }// else { // myApprox=TRUE
658   }// case IntPatch_Analytic:
659   break;
660
661                           //########################################  
662                           // Walking
663                           //######################################## 
664   case IntPatch_Walking:{
665     Handle(IntPatch_WLine) WL = 
666       Handle(IntPatch_WLine)::DownCast(L);
667     //
668     Standard_Integer ifprm, ilprm;
669     //
670     if (!myApprox) {
671       aNbParts=myLConstruct.NbParts();
672       for (i=1; i<=aNbParts; i++) {
673         myLConstruct.Part(i, fprm, lprm);
674         ifprm=(Standard_Integer)fprm;
675         ilprm=(Standard_Integer)lprm;
676         //        
677         Handle(Geom2d_BSplineCurve) aH1, aH2;
678
679         if(myApprox1) {
680           aH1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
681         }
682         if(myApprox2) {
683           aH2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
684         }
685         //
686         Handle(Geom_Curve) aBSp=MakeBSpline(WL, ifprm, ilprm);
687         //      
688         sline.Append(aBSp);
689         slineS1.Append(aH1);
690         slineS2.Append(aH2);
691       }
692     }
693     //
694     else {
695       Standard_Boolean bIsDecomposited;
696       Standard_Integer nbiter, aNbSeqOfL;
697       GeomInt_WLApprox theapp3d;
698       IntPatch_SequenceOfLine aSeqOfL;
699       Standard_Real tol2d, aTolSS;
700       //        
701       tol2d = myTolApprox;
702       aTolSS=2.e-7;
703       if(myHS1 == myHS2) { 
704         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False);      
705       }
706       else { 
707         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
708       }
709       //
710       bIsDecomposited = 
711         DecompositionOfWLine(WL, myHS1, myHS2, aTolSS, myLConstruct, aSeqOfL);
712       //
713       aNbParts=myLConstruct.NbParts();
714       aNbSeqOfL=aSeqOfL.Length();
715       //
716       nbiter = (bIsDecomposited) ? aNbSeqOfL : aNbParts;
717       //
718       for(i = 1; i <= nbiter; i++) {
719         if(bIsDecomposited) {
720           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
721           ifprm = 1;
722           ilprm = WL->NbPnts();
723         }
724         else {
725           myLConstruct.Part(i, fprm, lprm);
726           ifprm = (Standard_Integer)fprm;
727           ilprm = (Standard_Integer)lprm;
728         }
729         //-- lbr : 
730         //-- Si une des surfaces est un plan , on approxime en 2d
731         //-- sur cette surface et on remonte les points 2d en 3d.
732         GeomAbs_SurfaceType typs1, typs2;
733         typs1 = myHS1->Surface().GetType();
734         typs2 = myHS2->Surface().GetType();
735         //
736         if(typs1 == GeomAbs_Plane) { 
737           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,
738             Standard_True, myApprox2,
739             ifprm,  ilprm);
740         }         
741         else if(typs2 == GeomAbs_Plane) { 
742           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,
743             myApprox1,Standard_True,
744             ifprm,  ilprm);
745         }
746         else { 
747           //
748           if (myHS1 != myHS2){
749             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
750               (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
751
752                 theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
753                 //Standard_Boolean bUseSurfaces;
754                 //bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
755                 //if (bUseSurfaces) {
756                 //theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False);
757                 //}
758             }
759           }
760           //
761           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,
762             myApprox1,myApprox2,
763             ifprm,  ilprm);
764         }
765
766         if (!theapp3d.IsDone()) {
767           //      
768           Handle(Geom2d_BSplineCurve) aH1, aH2;
769           //      
770           Handle(Geom_Curve) aBSp=MakeBSpline(WL, ifprm, ilprm);
771           if(myApprox1) {
772             aH1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
773           }
774           if(myApprox2) {
775             aH2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
776           }
777           //
778           sline.Append(aBSp);
779           slineS1.Append(aH1);
780           slineS2.Append(aH2);
781         }//if (!theapp3d.IsDone())
782
783         else {
784           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) { 
785             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) { 
786               myTolReached2d = theapp3d.TolReached2d();
787             }
788           }
789           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) { 
790             myTolReached3d = myTolReached2d;
791           }
792           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) { 
793             myTolReached3d = theapp3d.TolReached3d();
794           }
795
796           Standard_Integer aNbMultiCurves, nbpoles;
797           //
798           aNbMultiCurves=theapp3d.NbMultiCurves(); 
799           for (j=1; j<=aNbMultiCurves; j++) {
800             if(typs1 == GeomAbs_Plane) {
801               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
802               nbpoles = mbspc.NbPoles();
803
804               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
805               TColgp_Array1OfPnt   tpoles(1,nbpoles);
806
807               mbspc.Curve(1,tpoles2d);
808               const gp_Pln&  Pln = myHS1->Surface().Plane();
809               //
810               Standard_Integer ik; 
811               for(ik = 1; ik<= nbpoles; ik++) { 
812                 tpoles.SetValue(ik,
813                   ElSLib::Value(tpoles2d.Value(ik).X(),
814                   tpoles2d.Value(ik).Y(),
815                   Pln));
816               }
817               //
818               Handle(Geom_BSplineCurve) BS = 
819                 new Geom_BSplineCurve(tpoles,
820                 mbspc.Knots(),
821                 mbspc.Multiplicities(),
822                 mbspc.Degree());
823               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
824               Check.FixTangent(Standard_True, Standard_True);
825               //        
826               sline.Append(BS);
827               //
828               if(myApprox1) { 
829                 Handle(Geom2d_BSplineCurve) BS1 = 
830                   new Geom2d_BSplineCurve(tpoles2d,
831                   mbspc.Knots(),
832                   mbspc.Multiplicities(),
833                   mbspc.Degree());
834                 GeomLib_Check2dBSplineCurve Check1(BS1,myTolCheck,myTolAngCheck);
835                 Check1.FixTangent(Standard_True,Standard_True);
836                 //      
837                 AdjustUPeriodic (aS1, BS1);  
838                 //
839                 slineS1.Append(BS1);
840               }
841               else {
842                 slineS1.Append(H1);
843               }
844
845               if(myApprox2) { 
846                 mbspc.Curve(2, tpoles2d);
847
848                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
849                   mbspc.Knots(),
850                   mbspc.Multiplicities(),
851                   mbspc.Degree());
852                 GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
853                 newCheck.FixTangent(Standard_True,Standard_True);
854                 //
855                 AdjustUPeriodic (aS2, BS2);  
856                 //
857                 slineS2.Append(BS2); 
858               }
859               else { 
860                 slineS2.Append(H1); 
861               }
862             }//if(typs1 == GeomAbs_Plane) 
863             //
864             else if(typs2 == GeomAbs_Plane) { 
865               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
866               nbpoles = mbspc.NbPoles();
867
868               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
869               TColgp_Array1OfPnt   tpoles(1,nbpoles);
870               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
871               const gp_Pln&  Pln = myHS2->Surface().Plane();
872               //
873               Standard_Integer ik; 
874               for(ik = 1; ik<= nbpoles; ik++) { 
875                 tpoles.SetValue(ik,
876                   ElSLib::Value(tpoles2d.Value(ik).X(),
877                   tpoles2d.Value(ik).Y(),
878                   Pln));
879
880               }
881               //
882               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
883                 mbspc.Knots(),
884                 mbspc.Multiplicities(),
885                 mbspc.Degree());
886               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
887               Check.FixTangent(Standard_True,Standard_True);
888               //        
889               sline.Append(BS);
890               //
891               if(myApprox2) {
892                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
893                   mbspc.Knots(),
894                   mbspc.Multiplicities(),
895                   mbspc.Degree());
896                 GeomLib_Check2dBSplineCurve Check1(BS1,myTolCheck,myTolAngCheck);
897                 Check1.FixTangent(Standard_True,Standard_True);
898                 //
899                 //
900                 AdjustUPeriodic (aS2, BS1);  
901                 //
902                 slineS2.Append(BS1);
903               }
904               else {
905                 slineS2.Append(H1);
906               }
907
908               if(myApprox1) { 
909                 mbspc.Curve(1,tpoles2d);
910                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
911                   mbspc.Knots(),
912                   mbspc.Multiplicities(),
913                   mbspc.Degree());
914                 GeomLib_Check2dBSplineCurve Check2(BS2,myTolCheck,myTolAngCheck);
915                 Check2.FixTangent(Standard_True,Standard_True);
916                 // 
917                 //
918                 AdjustUPeriodic (aS1, BS2);  
919                 //      
920                 slineS1.Append(BS2);
921               }
922               else { 
923                 slineS1.Append(H1);
924               }
925             } // else if(typs2 == GeomAbs_Plane)
926             //
927             else { // typs1!=GeomAbs_Plane && typs2!=GeomAbs_Plane
928               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
929               nbpoles = mbspc.NbPoles();
930               TColgp_Array1OfPnt tpoles(1,nbpoles);
931               mbspc.Curve(1,tpoles);
932               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
933                 mbspc.Knots(),
934                 mbspc.Multiplicities(),
935                 mbspc.Degree());
936               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
937               Check.FixTangent(Standard_True,Standard_True);
938               //        
939               //Check IsClosed()
940               Standard_Real aDist = Max(BS->StartPoint().XYZ().SquareModulus(),
941                 BS->EndPoint().XYZ().SquareModulus());
942               Standard_Real eps = Epsilon(aDist);
943               if(BS->StartPoint().SquareDistance(BS->EndPoint()) < 2.*eps)
944               {
945                 // Avoid creating B-splines containing two coincident poles only
946                 if (mbspc.Degree() == 1 && nbpoles == 2)
947                   continue;
948
949                 if (!BS->IsClosed() && !BS->IsPeriodic())
950                 {
951                   //force Closed()
952                   gp_Pnt aPm((BS->Pole(1).XYZ() + BS->Pole(BS->NbPoles()).XYZ()) / 2.);
953                   BS->SetPole(1, aPm);
954                   BS->SetPole(BS->NbPoles(), aPm);
955                 }
956               }
957               sline.Append(BS);
958
959               if(myApprox1) { 
960                 TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
961                 mbspc.Curve(2,tpoles2d);
962                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
963                   mbspc.Knots(),
964                   mbspc.Multiplicities(),
965                   mbspc.Degree());
966                 GeomLib_Check2dBSplineCurve newCheck(BS1,myTolCheck,myTolAngCheck);
967                 newCheck.FixTangent(Standard_True,Standard_True);
968                 //
969                 AdjustUPeriodic (aS1, BS1);  
970                 //
971                 slineS1.Append(BS1);
972               }
973               else {
974                 slineS1.Append(H1);
975               }
976               if(myApprox2) { 
977                 TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
978                 mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
979                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
980                   mbspc.Knots(),
981                   mbspc.Multiplicities(),
982                   mbspc.Degree());
983                 GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
984                 newCheck.FixTangent(Standard_True,Standard_True);
985                 //
986                 AdjustUPeriodic (aS2, BS2);  
987                 //
988                 slineS2.Append(BS2);
989               }
990               else { 
991                 slineS2.Append(H1);
992               }
993             }// else { // typs1!=GeomAbs_Plane && typs2!=GeomAbs_Plane
994           }// for (j=1; j<=aNbMultiCurves; j++
995         }
996       }
997     }
998   }
999   break;
1000
1001   case IntPatch_Restriction: 
1002     {
1003       Handle(IntPatch_RLine) RL = 
1004         Handle(IntPatch_RLine)::DownCast(L);
1005       Handle(Geom_Curve) aC3d;
1006       Handle(Geom2d_Curve) aC2d1, aC2d2;
1007       Standard_Real aTolReached;
1008       TreatRLine(RL, myHS1, myHS2, aC3d,
1009                   aC2d1, aC2d2, aTolReached);
1010
1011       if(aC3d.IsNull())
1012         break;
1013
1014       Bnd_Box2d aBox1, aBox2;
1015
1016       const Standard_Real aU1f = myHS1->FirstUParameter(),
1017                           aV1f = myHS1->FirstVParameter(),
1018                           aU1l = myHS1->LastUParameter(),
1019                           aV1l = myHS1->LastVParameter();
1020       const Standard_Real aU2f = myHS2->FirstUParameter(),
1021                           aV2f = myHS2->FirstVParameter(),
1022                           aU2l = myHS2->LastUParameter(),
1023                           aV2l = myHS2->LastVParameter();
1024
1025       aBox1.Add(gp_Pnt2d(aU1f, aV1f));
1026       aBox1.Add(gp_Pnt2d(aU1l, aV1l));
1027       aBox2.Add(gp_Pnt2d(aU2f, aV2f));
1028       aBox2.Add(gp_Pnt2d(aU2l, aV2l));
1029
1030       GeomInt_VectorOfReal anArrayOfParameters;
1031
1032       //We consider here that the intersection line is same-parameter-line
1033       anArrayOfParameters.Append(aC3d->FirstParameter());
1034       anArrayOfParameters.Append(aC3d->LastParameter());
1035
1036       TrimILineOnSurfBoundaries(aC2d1, aC2d2, aBox1, aBox2, anArrayOfParameters);
1037
1038       const Standard_Integer aNbIntersSolutionsm1 = anArrayOfParameters.Length() - 1;
1039
1040       //Trim RLine found.
1041       for(Standard_Integer anInd = 0; anInd < aNbIntersSolutionsm1; anInd++)
1042       {
1043         const Standard_Real aParF = anArrayOfParameters(anInd),
1044                             aParL = anArrayOfParameters(anInd+1);
1045
1046         if((aParL - aParF) <= Precision::PConfusion())
1047           continue;
1048
1049         const Standard_Real aPar = 0.5*(aParF + aParL);
1050         gp_Pnt2d aPt;
1051
1052         Handle(Geom2d_Curve) aCurv2d1, aCurv2d2;
1053         if(!aC2d1.IsNull())
1054         {
1055           aC2d1->D0(aPar, aPt);
1056
1057           if(aBox1.IsOut(aPt))
1058             continue;
1059
1060           if(myApprox1)
1061             aCurv2d1 = new Geom2d_TrimmedCurve(aC2d1, aParF, aParL);
1062         }
1063
1064         if(!aC2d2.IsNull())
1065         {
1066           aC2d2->D0(aPar, aPt);
1067
1068           if(aBox2.IsOut(aPt))
1069             continue;
1070
1071           if(myApprox2)
1072             aCurv2d2 = new Geom2d_TrimmedCurve(aC2d2, aParF, aParL);
1073         }
1074
1075         Handle(Geom_Curve) aCurv3d = new Geom_TrimmedCurve(aC3d, aParF, aParL);
1076
1077         sline.Append(aCurv3d);
1078         slineS1.Append(aCurv2d1);
1079         slineS2.Append(aCurv2d2);
1080       }
1081     }
1082     break;
1083   }
1084 }
1085 //
1086 //=======================================================================
1087 //function : AdjustUPeriodic
1088 //purpose  : 
1089 //=======================================================================
1090  static void AdjustUPeriodic (const Handle(Geom_Surface)& aS, const Handle(Geom2d_Curve)& aC2D)
1091 {
1092   if (aC2D.IsNull() || !aS->IsUPeriodic())
1093     return;
1094   //
1095   const Standard_Real aEps=Precision::PConfusion();//1.e-9
1096   const Standard_Real aEpsilon=Epsilon(10.);//1.77e-15 
1097   //
1098   Standard_Real umin,umax,vmin,vmax;
1099   aS->Bounds(umin,umax,vmin,vmax);
1100   const Standard_Real aPeriod = aS->UPeriod();
1101   
1102   const Standard_Real aT1=aC2D->FirstParameter();
1103   const Standard_Real aT2=aC2D->LastParameter();
1104   const Standard_Real aTx=aT1+0.467*(aT2-aT1);
1105   const gp_Pnt2d aPx=aC2D->Value(aTx);
1106   //
1107   Standard_Real aUx=aPx.X();
1108   if (fabs(aUx)<aEpsilon)
1109     aUx=0.;
1110   if (fabs(aUx-aPeriod)<aEpsilon)
1111     aUx=aPeriod;
1112   //
1113   Standard_Real dU=0.;
1114   while(aUx <(umin-aEps)) {
1115     aUx+=aPeriod;
1116     dU+=aPeriod;
1117   }
1118   while(aUx>(umax+aEps)) {
1119     aUx-=aPeriod;
1120     dU-=aPeriod;
1121   }
1122   // 
1123   if (dU!=0.) {
1124     gp_Vec2d aV2D(dU, 0.);
1125     aC2D->Translate(aV2D);
1126   }
1127 }
1128 //
1129 //
1130 //=======================================================================
1131 //function : GetQuadric
1132 //purpose  : 
1133 //=======================================================================
1134  static void GetQuadric(const Handle(GeomAdaptor_HSurface)& HS1, IntSurf_Quadric& quad1)
1135 {
1136   switch (HS1->Surface().GetType())
1137   {
1138     case GeomAbs_Plane:    quad1.SetValue(HS1->Surface().Plane()); break;
1139     case GeomAbs_Cylinder: quad1.SetValue(HS1->Surface().Cylinder()); break;
1140     case GeomAbs_Cone:     quad1.SetValue(HS1->Surface().Cone()); break;
1141     case GeomAbs_Sphere:   quad1.SetValue(HS1->Surface().Sphere()); break;
1142     case GeomAbs_Torus:    quad1.SetValue(HS1->Surface().Torus()); break;
1143     default: Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
1144   }
1145 }
1146
1147 //=======================================================================
1148 //function : Parameters
1149 //purpose  : 
1150 //=======================================================================
1151  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
1152                  const Handle(GeomAdaptor_HSurface)& HS2,
1153                  const gp_Pnt& Ptref,
1154                  Standard_Real& U1,
1155                  Standard_Real& V1,
1156                  Standard_Real& U2,
1157                  Standard_Real& V2)
1158 {
1159   IntSurf_Quadric quad1,quad2;
1160   //
1161   GetQuadric(HS1, quad1);
1162   GetQuadric(HS2, quad2);
1163   //
1164   quad1.Parameters(Ptref,U1,V1);
1165   quad2.Parameters(Ptref,U2,V2);
1166 }
1167
1168 //=======================================================================
1169 //function : MakeBSpline
1170 //purpose  : 
1171 //=======================================================================
1172 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
1173                                  const Standard_Integer ideb,
1174                                  const Standard_Integer ifin)
1175 {
1176   const Standard_Integer nbpnt = ifin-ideb+1;
1177   TColgp_Array1OfPnt poles(1,nbpnt);
1178   TColStd_Array1OfReal knots(1,nbpnt);
1179   TColStd_Array1OfInteger mults(1,nbpnt);
1180   Standard_Integer i = 1, ipidebm1 = ideb;
1181   for(; i<=nbpnt; ipidebm1++, i++)
1182   {
1183     poles(i) = WL->Point(ipidebm1).Value();
1184     mults(i) = 1;
1185     knots(i) = i-1;
1186   }
1187   mults(1) = mults(nbpnt) = 2;
1188   return new Geom_BSplineCurve(poles,knots,mults,1);
1189 }
1190
1191 //=======================================================================
1192 //function : MakeBSpline2d
1193 //purpose  : 
1194 //=======================================================================
1195 Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
1196                                           const Standard_Integer ideb,
1197                                           const Standard_Integer ifin,
1198                                           const Standard_Boolean onFirst)
1199 {
1200   const Standard_Integer nbpnt = ifin-ideb+1;
1201   TColgp_Array1OfPnt2d poles(1,nbpnt);
1202   TColStd_Array1OfReal knots(1,nbpnt);
1203   TColStd_Array1OfInteger mults(1,nbpnt);
1204   Standard_Integer i = 1, ipidebm1 = ideb;
1205   for(; i <= nbpnt; ipidebm1++, i++)
1206   {
1207     Standard_Real U, V;
1208     if(onFirst)
1209           theWLine->Point(ipidebm1).ParametersOnS1(U, V);
1210     else
1211           theWLine->Point(ipidebm1).ParametersOnS2(U, V);
1212     poles(i).SetCoord(U, V);
1213     mults(i) = 1;
1214     knots(i) = i-1;
1215   }
1216   mults(1) = mults(nbpnt) = 2;
1217   return new Geom2d_BSplineCurve(poles,knots,mults,1);
1218 }
1219
1220 //=========================================================================
1221 // static function : ComputePurgedWLine
1222 // purpose : Removes equal points (leave one of equal points) from theWLine
1223 //           and recompute vertex parameters.
1224 //           Returns new WLine or null WLine if the number
1225 //           of the points is less than 2.
1226 //=========================================================================
1227 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine) 
1228 {
1229   Handle(IntPatch_WLine) aResult;
1230   Handle(IntPatch_WLine) aLocalWLine;
1231   Handle(IntPatch_WLine) aTmpWLine = theWLine;
1232
1233   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1234   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
1235   Standard_Integer i, k, v, nb, nbvtx;
1236   nbvtx = theWLine->NbVertex();
1237   nb = theWLine->NbPnts();
1238
1239   for(i = 1; i <= nb; i++) {
1240     aLineOn2S->Add(theWLine->Point(i));
1241   }
1242
1243   for(v = 1; v <= nbvtx; v++) {
1244     aLocalWLine->AddVertex(theWLine->Vertex(v));
1245   }
1246   
1247   for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
1248     Standard_Integer aStartIndex = i + 1;
1249     Standard_Integer anEndIndex = i + 5;
1250     nb = aLineOn2S->NbPoints();
1251     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
1252
1253     if((aStartIndex >= nb) || (anEndIndex <= 1)) {
1254       continue;
1255     }
1256     k = aStartIndex;
1257
1258     while(k <= anEndIndex) {
1259       
1260       if(i != k) {
1261         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
1262         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
1263         
1264         if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
1265           aTmpWLine = aLocalWLine;
1266           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
1267
1268           for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
1269             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
1270             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
1271
1272             if(avertexindex >= k) {
1273               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
1274             }
1275             aLocalWLine->AddVertex(aVertex);
1276           }
1277           aLineOn2S->RemovePoint(k);
1278           anEndIndex--;
1279           continue;
1280         }
1281       }
1282       k++;
1283     }
1284   }
1285
1286   if(aLineOn2S->NbPoints() > 1) {
1287     aResult = aLocalWLine;
1288   }
1289   return aResult;
1290 }
1291 //=======================================================================
1292 //function : DecompositionOfWLine
1293 //purpose  : 
1294 //=======================================================================
1295 Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
1296                                       const Handle(GeomAdaptor_HSurface)&            theSurface1, 
1297                                       const Handle(GeomAdaptor_HSurface)&            theSurface2,
1298                                       const Standard_Real aTolSum, 
1299                                       const GeomInt_LineConstructor&                 theLConstructor,
1300                                       IntPatch_SequenceOfLine&                       theNewLines)
1301 {
1302   typedef NCollection_List<Standard_Integer> ListOfInteger;
1303   //have to use std::vector, not NCollection_Vector in order to use copy constructor of
1304   //ListOfInteger which will be created with specific allocator instance
1305   typedef std::vector<ListOfInteger, NCollection_StdAllocator<
1306       ListOfInteger> > ArrayOfListOfInteger;
1307
1308   Standard_Boolean bIsPrevPointOnBoundary, bIsCurrentPointOnBoundary;
1309   Standard_Integer nblines, aNbPnts, aNbParts, pit, i, j, aNbListOfPointIndex;
1310   Standard_Real aTol, umin, umax, vmin, vmax;
1311
1312   //an inc allocator, it will contain wasted space (upon list's Clear()) but it should
1313   //still be faster than the standard allocator, and wasted memory should not be
1314   //significant and will be limited by time span of this function;
1315   //this is a separate allocator from the anIncAlloc below what provides better data
1316   //locality in the latter (by avoiding wastes which will only be in anIdxAlloc)
1317   Handle(NCollection_IncAllocator) anIdxAlloc = new NCollection_IncAllocator();
1318   ListOfInteger aListOfPointIndex (anIdxAlloc);
1319
1320   //GeomAPI_ProjectPointOnSurf aPrj1, aPrj2;
1321   ProjectPointOnSurf aPrj1, aPrj2;
1322   Handle(Geom_Surface) aSurf1,  aSurf2;
1323   //
1324   aNbParts=theLConstructor.NbParts();
1325   aNbPnts=theWLine->NbPnts();
1326   //
1327   if((!aNbPnts) || (!aNbParts)){
1328     return Standard_False;
1329   }
1330   //
1331   Handle(NCollection_IncAllocator) anIncAlloc = new NCollection_IncAllocator();
1332   NCollection_StdAllocator<ListOfInteger> anAlloc (anIncAlloc);
1333   const ListOfInteger aDummy (anIncAlloc); //empty list to be copy constructed from
1334   ArrayOfListOfInteger anArrayOfLines (aNbPnts + 1, aDummy, anAlloc);
1335
1336   NCollection_LocalArray<Standard_Integer> anArrayOfLineTypeArr (aNbPnts + 1);
1337   Standard_Integer* anArrayOfLineType = anArrayOfLineTypeArr;
1338   //
1339   nblines = 0;
1340   aTol = Precision::Confusion();
1341   //
1342   aSurf1 = theSurface1->ChangeSurface().Surface();
1343   aSurf1->Bounds(umin, umax, vmin, vmax);
1344   aPrj1.Init(aSurf1, umin, umax, vmin, vmax);
1345   //
1346   aSurf2 = theSurface2->ChangeSurface().Surface();
1347   aSurf2->Bounds(umin, umax, vmin, vmax);
1348   aPrj2.Init(aSurf2, umin, umax, vmin, vmax);
1349   //
1350   //
1351   bIsPrevPointOnBoundary=Standard_False;
1352   for(pit=1; pit<=aNbPnts; pit++) {
1353     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
1354     bIsCurrentPointOnBoundary=Standard_False;
1355     //
1356     // whether aPoint is on boundary or not
1357     //
1358     for(i=0; i<2; i++) {// exploration Surface 1,2 
1359       Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
1360       aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
1361       //
1362       for(j=0; j<2; j++) {// exploration of coordinate U,V
1363         Standard_Boolean isperiodic;
1364         //
1365         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1366         if(!isperiodic) {
1367           continue;
1368         }
1369         //
1370         Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
1371         Standard_Real aParameter, anoffset, anAdjustPar;
1372         Standard_Boolean bIsOnFirstBoundary, bIsPointOnBoundary;
1373         //
1374         aResolution = (!j) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
1375         aPeriod     = (!j) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1376         alowerboundary = (!j) ? umin : vmin;
1377         aupperboundary = (!j) ? umax : vmax;
1378         U=0.;V=0.;//?
1379         //
1380         if(!i){
1381           aPoint.ParametersOnS1(U, V);
1382         }
1383         else{
1384           aPoint.ParametersOnS2(U, V);
1385         }
1386         //
1387         aParameter = (!j) ? U : V;
1388         anoffset=0.;
1389         anAdjustPar=AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1390         //
1391         bIsOnFirstBoundary=Standard_True;
1392         //
1393         bIsPointOnBoundary=
1394           IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
1395         
1396         if(bIsPointOnBoundary) {
1397           bIsCurrentPointOnBoundary = Standard_True;
1398           break;
1399         }
1400       }// for(j=0; j<2; j++)
1401       
1402       if(bIsCurrentPointOnBoundary){
1403         break;
1404       }
1405     }// for(i=0; i<2; i++) 
1406     //
1407     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
1408
1409       if(!aListOfPointIndex.IsEmpty()) {
1410         nblines++;
1411         anArrayOfLines[nblines] = aListOfPointIndex;
1412         anArrayOfLineType[nblines] = bIsPrevPointOnBoundary;
1413         aListOfPointIndex.Clear();
1414       }
1415       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
1416     }
1417     aListOfPointIndex.Append(pit);
1418   } // for(pit=1; pit<=aNbPnts; pit++)
1419   //
1420   aNbListOfPointIndex=aListOfPointIndex.Extent();
1421   if(aNbListOfPointIndex) {
1422     nblines++;
1423     anArrayOfLines[nblines].Assign (aListOfPointIndex);
1424     anArrayOfLineType[nblines] = bIsPrevPointOnBoundary;
1425     aListOfPointIndex.Clear();
1426   }
1427   //
1428   if(nblines <= 1){
1429     return Standard_False;
1430   }
1431   //
1432   // Correct wlines.begin
1433   Standard_Integer aLineType;
1434   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
1435   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S (new NCollection_IncAllocator());
1436   //
1437   for(i = 1; i <= nblines; i++) {
1438     aLineType=anArrayOfLineType[i];
1439     if(aLineType) {
1440       continue;
1441     }
1442     //
1443     const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1444     if(aListOfIndex.Extent() < 2) {
1445       continue;
1446     }
1447     //
1448     TColStd_ListOfInteger aListOfFLIndex;
1449     Standard_Integer aneighbourindex, aLineTypeNeib;
1450     //
1451     for(j = 0; j < 2; j++) {// neighbour line choice 
1452       aneighbourindex = (!j) ? (i-1) : (i+1);
1453       if((aneighbourindex < 1) || (aneighbourindex > nblines)){
1454         continue;
1455       }
1456       //
1457       aLineTypeNeib=anArrayOfLineType[aneighbourindex];
1458       if(!aLineTypeNeib){
1459         continue;
1460       }
1461       //
1462       const ListOfInteger& aNeighbour = anArrayOfLines[aneighbourindex];
1463       Standard_Integer anIndex = (!j) ? aNeighbour.Last() : aNeighbour.First();
1464       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
1465       // check if need use derivative.begin .end [absence]
1466       //
1467       IntSurf_PntOn2S aNewP = aPoint;
1468       Standard_Integer surfit, parit;
1469       //
1470       for(surfit = 0; surfit < 2; ++surfit) {
1471
1472         Handle(GeomAdaptor_HSurface) aGASurface = (!surfit) ? theSurface1 : theSurface2;
1473         
1474         umin = aGASurface->FirstUParameter();
1475         umax = aGASurface->LastUParameter();
1476         vmin = aGASurface->FirstVParameter();
1477         vmax = aGASurface->LastVParameter();
1478         Standard_Real U=0., V=0.;
1479
1480         if(!surfit) {
1481           aNewP.ParametersOnS1(U, V);
1482         }
1483         else {
1484           aNewP.ParametersOnS2(U, V);
1485         }
1486         //
1487         Standard_Integer nbboundaries = 0;
1488         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
1489         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
1490         //
1491         for(parit = 0; parit < 2; parit++) {
1492           Standard_Boolean isperiodic = (!parit) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1493
1494           Standard_Real aResolution = (!parit) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
1495           Standard_Real alowerboundary = (!parit) ? umin : vmin;
1496           Standard_Real aupperboundary = (!parit) ? umax : vmax;
1497
1498           Standard_Real aParameter = (!parit) ? U : V;
1499           Standard_Boolean bIsOnFirstBoundary = Standard_True;
1500   
1501           if(!isperiodic) {
1502             if(IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary)) {
1503               bIsUBoundary = (!parit);
1504               bIsFirstBoundary = bIsOnFirstBoundary;
1505               nbboundaries++;
1506             }
1507           }
1508           else {
1509             Standard_Real aPeriod     = (!parit) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1510             Standard_Real anoffset = 0.;
1511             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1512
1513             if(IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary)) {
1514               bIsUBoundary = (parit == 0);
1515               bIsFirstBoundary = bIsOnFirstBoundary;
1516               nbboundaries++;
1517             }
1518           }
1519         }
1520         //
1521         Standard_Boolean bComputeLineEnd = Standard_False;
1522         
1523         if(nbboundaries == 2) {
1524           bComputeLineEnd = Standard_True;
1525         }
1526         else if(nbboundaries == 1) {
1527           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1528
1529           if(isperiodic) {
1530             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
1531             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
1532             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1533             Standard_Real aParameter = (bIsUBoundary) ? U : V;
1534             Standard_Real anoffset = 0.;
1535             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1536
1537             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
1538             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
1539             anotherPar += anoffset;
1540             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
1541             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
1542             Standard_Real nU1, nV1;
1543
1544             if(surfit == 0)
1545               aNeighbourPoint.ParametersOnS1(nU1, nV1);
1546             else
1547               aNeighbourPoint.ParametersOnS2(nU1, nV1);
1548             
1549             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
1550             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
1551             bComputeLineEnd = Standard_True;
1552             Standard_Boolean bCheckAngle1 = Standard_False;
1553             Standard_Boolean bCheckAngle2 = Standard_False;
1554             gp_Vec2d aNewVec;
1555             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
1556             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
1557             //
1558             if(((adist1 - adist2) > Precision::PConfusion()) && 
1559                (adist2 < (aPeriod / 4.))) {
1560               bCheckAngle1 = Standard_True;
1561               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
1562
1563               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
1564                 aNewP.SetValue((surfit == 0), anewU, anewV);
1565                 bCheckAngle1 = Standard_False;
1566               }
1567             }
1568             else if(adist1 < (aPeriod / 4.)) {
1569               bCheckAngle2 = Standard_True;
1570               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
1571
1572               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
1573                 bCheckAngle2 = Standard_False;
1574               }
1575             }
1576             //
1577             if(bCheckAngle1 || bCheckAngle2) {
1578               // assume there are at least two points in line (see "if" above)
1579               Standard_Integer anindexother = aneighbourpointindex;
1580               
1581               while((anindexother <= aListOfIndex.Last()) && (anindexother >= aListOfIndex.First())) {
1582                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
1583                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
1584                 Standard_Real nU2, nV2;
1585                 
1586                 if(surfit == 0)
1587                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1588                 else
1589                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1590                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
1591
1592                 if(aVecOld.SquareMagnitude() <= (gp::Resolution() * gp::Resolution())) {
1593                   continue;
1594                 }
1595                 else {
1596                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
1597
1598                   if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
1599
1600                     if(bCheckAngle1) {
1601                       Standard_Real U1, U2, V1, V2;
1602                       IntSurf_PntOn2S atmppoint = aNewP;
1603                       atmppoint.SetValue((surfit == 0), anewU, anewV);
1604                       atmppoint.Parameters(U1, V1, U2, V2);
1605                       gp_Pnt P1 = theSurface1->Value(U1, V1);
1606                       gp_Pnt P2 = theSurface2->Value(U2, V2);
1607                       gp_Pnt P0 = aPoint.Value();
1608
1609                       if(P0.IsEqual(P1, aTol) &&
1610                          P0.IsEqual(P2, aTol) &&
1611                          P1.IsEqual(P2, aTol)) {
1612                         bComputeLineEnd = Standard_False;
1613                         aNewP.SetValue((surfit == 0), anewU, anewV);
1614                       }
1615                     }
1616                     
1617                     if(bCheckAngle2) {
1618                       bComputeLineEnd = Standard_False;
1619                     }
1620                   }
1621                   break;
1622                 }
1623               } // end while(anindexother...)
1624             }
1625           }
1626         }
1627         //
1628         if(bComputeLineEnd) {
1629           Standard_Integer aneighbourpointindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
1630           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
1631           Standard_Real nU1, nV1;
1632
1633           if(surfit == 0)
1634             aNeighbourPoint.ParametersOnS1(nU1, nV1);
1635           else
1636             aNeighbourPoint.ParametersOnS2(nU1, nV1);
1637           gp_Pnt2d ap1(nU1, nV1);
1638           gp_Pnt2d ap2(nU1, nV1);
1639           Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
1640
1641           while((aneighbourpointindex2 <= aListOfIndex.Last()) && (aneighbourpointindex2 >= aListOfIndex.First())) {
1642             aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
1643             const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
1644             Standard_Real nU2, nV2;
1645
1646             if(surfit == 0)
1647               aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1648             else
1649               aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1650             ap2.SetX(nU2);
1651             ap2.SetY(nV2);
1652
1653             if(ap1.SquareDistance(ap2) > (gp::Resolution() * gp::Resolution())) {
1654               break;
1655             }
1656           }       
1657           gp_Pnt2d anewpoint;
1658           Standard_Boolean found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
1659
1660           if(found) {
1661             // check point
1662             Standard_Real aCriteria =aTolSum;// BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
1663             //GeomAPI_ProjectPointOnSurf& aProjector = (surfit == 0) ? aPrj2 : aPrj1;
1664             ProjectPointOnSurf& aProjector = (surfit == 0) ? aPrj2 : aPrj1;
1665             Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
1666
1667             gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
1668             aProjector.Perform(aP3d);
1669
1670             if(aProjector.IsDone()) {
1671               if(aProjector.LowerDistance() < aCriteria) {
1672                 Standard_Real foundU = U, foundV = V;
1673                 aProjector.LowerDistanceParameters(foundU, foundV);
1674
1675                 if(surfit == 0)
1676                   aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
1677                 else
1678                   aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
1679               }
1680             }
1681           }
1682         }
1683       }
1684       aSeqOfPntOn2S->Add(aNewP);
1685       aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
1686     }
1687     anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
1688   }
1689   // Correct wlines.end
1690
1691   // Split wlines.begin
1692   for(j = 1; j <= theLConstructor.NbParts(); j++) {
1693     Standard_Real fprm = 0., lprm = 0.;
1694     theLConstructor.Part(j, fprm, lprm);
1695     Standard_Integer ifprm = (Standard_Integer)fprm;
1696     Standard_Integer ilprm = (Standard_Integer)lprm;
1697     //
1698     Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1699     //
1700     for(i = 1; i <= nblines; i++) {
1701       if(anArrayOfLineType[i] != 0) {
1702         continue;
1703       }
1704       const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1705
1706       if(aListOfIndex.Extent() < 2) {
1707         continue;
1708       }
1709       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
1710       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
1711       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
1712
1713       if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
1714         bhasfirstpoint = (i != 1);
1715       }
1716
1717       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
1718         bhaslastpoint = (i != nblines);
1719       }
1720       Standard_Boolean bIsFirstInside = ((ifprm >= aListOfIndex.First()) && (ifprm <= aListOfIndex.Last()));
1721       Standard_Boolean bIsLastInside =  ((ilprm >= aListOfIndex.First()) && (ilprm <= aListOfIndex.Last()));
1722
1723       if(!bIsFirstInside && !bIsLastInside) {
1724         if((ifprm < aListOfIndex.First()) && (ilprm > aListOfIndex.Last())) {
1725           // append whole line, and boundaries if neccesary
1726           if(bhasfirstpoint) {
1727             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
1728             aLineOn2S->Add(aP);
1729           }
1730           ListOfInteger::Iterator anIt(aListOfIndex);
1731
1732           for(; anIt.More(); anIt.Next()) {
1733             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1734             aLineOn2S->Add(aP);
1735           }
1736
1737           if(bhaslastpoint) {
1738             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
1739             aLineOn2S->Add(aP);
1740           }
1741
1742           // check end of split line (end is almost always)
1743           Standard_Integer aneighbour = i + 1;
1744           Standard_Boolean bIsEndOfLine = Standard_True;
1745
1746           if(aneighbour <= nblines) {
1747             const ListOfInteger& aListOfNeighbourIndex = anArrayOfLines[aneighbour];
1748
1749             if((anArrayOfLineType[aneighbour] != 0) &&
1750                (aListOfNeighbourIndex.IsEmpty())) {
1751               bIsEndOfLine = Standard_False;
1752             }
1753           }
1754
1755           if(bIsEndOfLine) {
1756             if(aLineOn2S->NbPoints() > 1) {
1757               Handle(IntPatch_WLine) aNewWLine = 
1758                 new IntPatch_WLine(aLineOn2S, Standard_False);
1759               theNewLines.Append(aNewWLine);
1760             }
1761             aLineOn2S = new IntSurf_LineOn2S();
1762           }
1763         }
1764         continue;
1765       }
1766       // end if(!bIsFirstInside && !bIsLastInside)
1767
1768       if(bIsFirstInside && bIsLastInside) {
1769         // append inside points between ifprm and ilprm
1770         ListOfInteger::Iterator anIt(aListOfIndex);
1771
1772         for(; anIt.More(); anIt.Next()) {
1773           if((anIt.Value() < ifprm) || (anIt.Value() > ilprm))
1774             continue;
1775           const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1776           aLineOn2S->Add(aP);
1777         }
1778       }
1779       else {
1780
1781         if(bIsFirstInside) {
1782           // append points from ifprm to last point + boundary point
1783           ListOfInteger::Iterator anIt(aListOfIndex);
1784
1785           for(; anIt.More(); anIt.Next()) {
1786             if(anIt.Value() < ifprm)
1787               continue;
1788             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1789             aLineOn2S->Add(aP);
1790           }
1791
1792           if(bhaslastpoint) {
1793             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
1794             aLineOn2S->Add(aP);
1795           }
1796           // check end of split line (end is almost always)
1797           Standard_Integer aneighbour = i + 1;
1798           Standard_Boolean bIsEndOfLine = Standard_True;
1799
1800           if(aneighbour <= nblines) {
1801             const ListOfInteger& aListOfNeighbourIndex = anArrayOfLines[aneighbour];
1802
1803             if((anArrayOfLineType[aneighbour] != 0) &&
1804                (aListOfNeighbourIndex.IsEmpty())) {
1805               bIsEndOfLine = Standard_False;
1806             }
1807           }
1808
1809           if(bIsEndOfLine) {
1810             if(aLineOn2S->NbPoints() > 1) {
1811               Handle(IntPatch_WLine) aNewWLine = 
1812                 new IntPatch_WLine(aLineOn2S, Standard_False);
1813               theNewLines.Append(aNewWLine);
1814             }
1815             aLineOn2S = new IntSurf_LineOn2S();
1816           }
1817         }
1818         // end if(bIsFirstInside)
1819
1820         if(bIsLastInside) {
1821           // append points from first boundary point to ilprm
1822           if(bhasfirstpoint) {
1823             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
1824             aLineOn2S->Add(aP);
1825           }
1826           ListOfInteger::Iterator anIt(aListOfIndex);
1827
1828           for(; anIt.More(); anIt.Next()) {
1829             if(anIt.Value() > ilprm)
1830               continue;
1831             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1832             aLineOn2S->Add(aP);
1833           }
1834         }
1835         //end if(bIsLastInside)
1836       }
1837     }
1838
1839     if(aLineOn2S->NbPoints() > 1) {
1840       Handle(IntPatch_WLine) aNewWLine = 
1841         new IntPatch_WLine(aLineOn2S, Standard_False);
1842       theNewLines.Append(aNewWLine);
1843     }
1844   }
1845   // Split wlines.end
1846   //
1847   // cda002/I3
1848   Standard_Real fprm, lprm;
1849   Standard_Integer ifprm, ilprm, aNbPoints, aIndex;
1850   //
1851   aNbParts=theLConstructor.NbParts();
1852   //
1853   for(j = 1; j <= aNbParts; j++) {
1854     theLConstructor.Part(j, fprm, lprm);
1855     ifprm=(Standard_Integer)fprm;
1856     ilprm=(Standard_Integer)lprm;
1857     //
1858     if ((ilprm-ifprm)==1) {
1859       for(i = 1; i <= nblines; i++) {
1860         aLineType=anArrayOfLineType[i];
1861         if(aLineType) {
1862           continue;
1863         }
1864         //
1865         const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1866         aNbPoints=aListOfIndex.Extent();
1867         if(aNbPoints==1) {
1868           aIndex=aListOfIndex.First();
1869           if (aIndex==ifprm || aIndex==ilprm) {
1870             Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1871             const IntSurf_PntOn2S& aP1 = theWLine->Point(ifprm);
1872             const IntSurf_PntOn2S& aP2 = theWLine->Point(ilprm);
1873             aLineOn2S->Add(aP1);
1874             aLineOn2S->Add(aP2);
1875             Handle(IntPatch_WLine) aNewWLine = 
1876               new IntPatch_WLine(aLineOn2S, Standard_False);
1877             theNewLines.Append(aNewWLine);
1878           }
1879         }
1880       }
1881     }
1882   }
1883   //
1884   return Standard_True;
1885 }
1886
1887 //=======================================================================
1888 //function : AdjustPeriodic
1889 //purpose  : 
1890 //=======================================================================
1891 Standard_Real AdjustPeriodic(const Standard_Real theParameter,
1892                              const Standard_Real parmin,
1893                              const Standard_Real parmax,
1894                              const Standard_Real thePeriod,
1895                              Standard_Real&      theOffset)
1896 {
1897   Standard_Real aresult = theParameter;
1898   theOffset = 0.;
1899   while(aresult < parmin) {
1900     aresult += thePeriod;
1901     theOffset += thePeriod;
1902   }
1903   while(aresult > parmax) {
1904     aresult -= thePeriod;
1905     theOffset -= thePeriod;
1906   }
1907   return aresult;
1908 }
1909
1910 //=======================================================================
1911 //function : IsPointOnBoundary
1912 //purpose  : 
1913 //=======================================================================
1914 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
1915                                    const Standard_Real theFirstBoundary,
1916                                    const Standard_Real theSecondBoundary,
1917                                    const Standard_Real theResolution,
1918                                    Standard_Boolean&   IsOnFirstBoundary)
1919 {
1920   IsOnFirstBoundary = Standard_True;
1921   if(fabs(theParameter - theFirstBoundary) < theResolution)
1922     return Standard_True;
1923   if(fabs(theParameter - theSecondBoundary) < theResolution)
1924   {
1925     IsOnFirstBoundary = Standard_False;
1926     return Standard_True;
1927   }
1928   return Standard_False;
1929 }
1930 //=======================================================================
1931 //function : FindPoint
1932 //purpose  : 
1933 //=======================================================================
1934 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
1935                            const gp_Pnt2d&     theLastPoint,
1936                            const Standard_Real theUmin,
1937                            const Standard_Real theUmax,
1938                            const Standard_Real theVmin,
1939                            const Standard_Real theVmax,
1940                            gp_Pnt2d&           theNewPoint)
1941 {
1942   gp_Vec2d aVec(theFirstPoint, theLastPoint);
1943   Standard_Integer i = 0, j = 0;
1944
1945   for(i = 0; i < 4; i++) {
1946     gp_Vec2d anOtherVec;
1947     gp_Vec2d anOtherVecNormal;
1948     gp_Pnt2d aprojpoint = theLastPoint;    
1949
1950     if((i % 2) == 0) {
1951       anOtherVec.SetX(0.);
1952       anOtherVec.SetY(1.);
1953       anOtherVecNormal.SetX(1.);
1954       anOtherVecNormal.SetY(0.);
1955
1956       if(i < 2)
1957         aprojpoint.SetX(theUmin);
1958       else
1959         aprojpoint.SetX(theUmax);
1960     }
1961     else {
1962       anOtherVec.SetX(1.);
1963       anOtherVec.SetY(0.);
1964       anOtherVecNormal.SetX(0.);
1965       anOtherVecNormal.SetY(1.);
1966
1967       if(i < 2)
1968         aprojpoint.SetY(theVmin);
1969       else
1970         aprojpoint.SetY(theVmax);
1971     }
1972     gp_Vec2d anormvec = aVec;
1973     anormvec.Normalize();
1974     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
1975
1976     if(fabs(adot1) < Precision::Angular())
1977       continue;
1978     Standard_Real adist = 0.;
1979
1980     if((i % 2) == 0) {
1981       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
1982     }
1983     else {
1984       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
1985     }
1986     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
1987
1988     for(j = 0; j < 2; j++) {
1989       anoffset = (j == 0) ? anoffset : -anoffset;
1990       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
1991       gp_Vec2d acurvec(theLastPoint, acurpoint);
1992
1993       //
1994       Standard_Real aDotX, anAngleX, aPC;
1995       //
1996       aDotX=aVec.Dot(acurvec);
1997       anAngleX=aVec.Angle(acurvec);
1998       aPC=Precision::PConfusion();
1999       //
2000       if(aDotX > 0. && fabs(anAngleX) < aPC) {
2001       //
2002         if((i % 2) == 0) {
2003           if((acurpoint.Y() >= theVmin) &&
2004              (acurpoint.Y() <= theVmax)) {
2005             theNewPoint = acurpoint;
2006             return Standard_True;
2007           }
2008         }
2009         else {
2010           if((acurpoint.X() >= theUmin) &&
2011              (acurpoint.X() <= theUmax)) {
2012             theNewPoint = acurpoint;
2013             return Standard_True;
2014           }
2015         }
2016       }
2017     }
2018   }
2019   return Standard_False;
2020 }
2021
2022 //=======================================================================
2023 //function : ParametersOfNearestPointOnSurface
2024 //purpose  : 
2025 //=======================================================================
2026 static Standard_Boolean ParametersOfNearestPointOnSurface(const Extrema_ExtPS theExtr,
2027                                                           Standard_Real& theU,
2028                                                           Standard_Real& theV)
2029 {
2030   if(!theExtr.IsDone() || !theExtr.NbExt())
2031     return Standard_False;
2032
2033   Standard_Integer anIndex = 1;
2034   Standard_Real aMinSQDist = theExtr.SquareDistance(anIndex);
2035   for(Standard_Integer i = 2; i <= theExtr.NbExt(); i++)
2036   {
2037     Standard_Real aSQD = theExtr.SquareDistance(i);
2038     if (aSQD < aMinSQDist)
2039     {
2040       aMinSQDist = aSQD;
2041       anIndex = i;
2042     }
2043   }
2044
2045   theExtr.Point(anIndex).Parameter(theU, theV);
2046
2047   return Standard_True;
2048 }
2049
2050 //=======================================================================
2051 //function : GetSegmentBoundary
2052 //purpose  : 
2053 //=======================================================================
2054 static void GetSegmentBoundary( const IntRes2d_IntersectionSegment& theSegm,
2055                                 const Handle(Geom2d_Curve)& theCurve,
2056                                 GeomInt_VectorOfReal& theArrayOfParameters)
2057 {
2058   Standard_Real aU1 = theCurve->FirstParameter(), aU2 = theCurve->LastParameter();
2059
2060   if(theSegm.HasFirstPoint())
2061   {
2062     const IntRes2d_IntersectionPoint& anIPF = theSegm.FirstPoint();
2063     aU1 = anIPF.ParamOnFirst();
2064   }
2065
2066   if(theSegm.HasLastPoint())
2067   {
2068     const IntRes2d_IntersectionPoint& anIPL = theSegm.LastPoint();
2069     aU2 = anIPL.ParamOnFirst();
2070   }
2071
2072   theArrayOfParameters.Append(aU1);
2073   theArrayOfParameters.Append(aU2);
2074 }
2075
2076 //=======================================================================
2077 //function : IntersectCurveAndBoundary
2078 //purpose  : 
2079 //=======================================================================
2080 static void IntersectCurveAndBoundary(const Handle(Geom2d_Curve)& theC2d,
2081                                       const Handle(Geom2d_Curve)* const theArrBounds,
2082                                       const Standard_Integer theNumberOfCurves,
2083                                       const Standard_Real theTol,
2084                                       GeomInt_VectorOfReal& theArrayOfParameters)
2085 {
2086   if(theC2d.IsNull())
2087     return;
2088
2089   Geom2dAdaptor_Curve anAC1(theC2d);
2090   for(Standard_Integer aCurID = 0; aCurID < theNumberOfCurves; aCurID++)
2091   {
2092     if(theArrBounds[aCurID].IsNull())
2093       continue;
2094
2095     Geom2dAdaptor_Curve anAC2(theArrBounds[aCurID]);
2096     Geom2dInt_GInter anIntCC2d(anAC1, anAC2, theTol, theTol);
2097
2098     if(!anIntCC2d.IsDone() || anIntCC2d.IsEmpty())
2099       continue;
2100
2101     for (Standard_Integer aPntID = 1; aPntID <= anIntCC2d.NbPoints(); aPntID++)
2102     {
2103       const Standard_Real aParam = anIntCC2d.Point(aPntID).ParamOnFirst();
2104       theArrayOfParameters.Append(aParam);
2105     }
2106
2107     for (Standard_Integer aSegmID = 1; aSegmID <= anIntCC2d.NbSegments(); aSegmID++)
2108     {
2109       GetSegmentBoundary(anIntCC2d.Segment(aSegmID), theC2d, theArrayOfParameters);
2110     }
2111   }
2112 }
2113
2114 //=======================================================================
2115 //function : TreatRLine
2116 //purpose  : Approx of Restriction line
2117 //=======================================================================
2118 void GeomInt_IntSS::TreatRLine(const Handle(IntPatch_RLine)& theRL, 
2119                 const Handle(GeomAdaptor_HSurface)& theHS1, 
2120                 const Handle(GeomAdaptor_HSurface)& theHS2, 
2121                 Handle(Geom_Curve)& theC3d,
2122                 Handle(Geom2d_Curve)& theC2d1, 
2123                 Handle(Geom2d_Curve)& theC2d2, 
2124                 Standard_Real& theTolReached)
2125 {
2126   Handle(GeomAdaptor_HSurface) aGAHS;
2127   Handle(Adaptor2d_HCurve2d) anAHC2d;
2128   Standard_Real tf, tl;
2129   gp_Lin2d aL;
2130   // It is assumed that 2d curve is 2d line (rectangular surface domain)
2131   if(theRL->IsArcOnS1())
2132   {
2133     aGAHS = theHS1;
2134     anAHC2d = theRL->ArcOnS1();
2135     theRL->ParamOnS1(tf, tl);
2136     theC2d1 = Geom2dAdaptor::MakeCurve(anAHC2d->Curve2d());
2137   }
2138   else if (theRL->IsArcOnS2())
2139   {
2140     aGAHS = theHS2;
2141     anAHC2d = theRL->ArcOnS2();
2142     theRL->ParamOnS2(tf, tl);
2143     theC2d2 = Geom2dAdaptor::MakeCurve(anAHC2d->Curve2d());
2144   }
2145   else
2146   {
2147     return;
2148   }
2149   //
2150   //To provide sameparameter it is necessary to get 3d curve as
2151   //approximation of curve on surface.
2152   Standard_Integer aMaxDeg = 8;
2153   Standard_Integer aMaxSeg = 1000;
2154   Approx_CurveOnSurface anApp(anAHC2d, aGAHS, tf, tl, Precision::Confusion(),
2155                               GeomAbs_C1, aMaxDeg, aMaxSeg, 
2156                               Standard_True, Standard_False);
2157   if(!anApp.HasResult())
2158     return;
2159
2160   theC3d = anApp.Curve3d();
2161   theTolReached = anApp.MaxError3d();
2162   Standard_Real aTol = Precision::Confusion();
2163   if(theRL->IsArcOnS1())
2164   {
2165     Handle(Geom_Surface) aS = GeomAdaptor::MakeSurface(theHS2->Surface());
2166     BuildPCurves (tf, tl, aTol, 
2167                   aS, theC3d, theC2d2);
2168   }
2169   if(theRL->IsArcOnS2())
2170   {
2171     Handle(Geom_Surface) aS = GeomAdaptor::MakeSurface(theHS1->Surface());
2172     BuildPCurves (tf, tl, aTol, 
2173                   aS, theC3d, theC2d1);
2174   }
2175   theTolReached = Max(theTolReached, aTol);
2176 }
2177
2178 //=======================================================================
2179 //function : BuildPCurves
2180 //purpose  : 
2181 //=======================================================================
2182 void GeomInt_IntSS::BuildPCurves (Standard_Real f,
2183                                   Standard_Real l,
2184                                   Standard_Real& Tol,
2185                                   const Handle (Geom_Surface)& S,
2186                                   const Handle (Geom_Curve)&   C,
2187                                   Handle (Geom2d_Curve)& C2d)
2188 {
2189   if (!C2d.IsNull()) {
2190     return;
2191   }
2192   //
2193   Standard_Real umin,umax,vmin,vmax;
2194   // 
2195   S->Bounds(umin, umax, vmin, vmax);
2196   // in class ProjLib_Function the range of parameters is shrank by 1.e-09
2197   if((l - f) > 2.e-09) {
2198     C2d = GeomProjLib::Curve2d(C,f,l,S,umin,umax,vmin,vmax,Tol);
2199     //
2200     if (C2d.IsNull()) {
2201       // proj. a circle that goes through the pole on a sphere to the sphere     
2202       Tol += Precision::Confusion();
2203       C2d = GeomProjLib::Curve2d(C,f,l,S,Tol);
2204     }
2205   }
2206   else {
2207     if((l - f) > Epsilon(Abs(f)))
2208     {
2209       //The domain of C2d is [Epsilon(Abs(f)), 2.e-09]
2210       //On this small range C2d can be considered as segment 
2211       //of line.
2212
2213       Standard_Real aU=0., aV=0.;
2214       GeomAdaptor_Surface anAS;
2215       anAS.Load(S);
2216       Extrema_ExtPS anExtr;
2217       const gp_Pnt aP3d1 = C->Value(f);
2218       const gp_Pnt aP3d2 = C->Value(l);
2219
2220       anExtr.SetAlgo(Extrema_ExtAlgo_Grad);
2221       anExtr.Initialize(anAS, umin, umax, vmin, vmax,
2222                                 Precision::Confusion(), Precision::Confusion());
2223       anExtr.Perform(aP3d1);
2224
2225       if(ParametersOfNearestPointOnSurface(anExtr, aU, aV))
2226       {
2227         const gp_Pnt2d aP2d1(aU, aV);
2228
2229         anExtr.Perform(aP3d2);
2230
2231         if(ParametersOfNearestPointOnSurface(anExtr, aU, aV))
2232         {
2233           const gp_Pnt2d aP2d2(aU, aV);
2234
2235           if(aP2d1.Distance(aP2d2) > gp::Resolution())
2236           {
2237             TColgp_Array1OfPnt2d poles(1,2);
2238             TColStd_Array1OfReal knots(1,2);
2239             TColStd_Array1OfInteger mults(1,2);
2240             poles(1) = aP2d1;
2241             poles(2) = aP2d2;
2242             knots(1) = f;
2243             knots(2) = l;
2244             mults(1) = mults(2) = 2;
2245
2246             C2d = new Geom2d_BSplineCurve(poles,knots,mults,1);
2247
2248             //Check same parameter in middle point .begin
2249             const gp_Pnt PMid(C->Value(0.5*(f+l)));
2250             const gp_Pnt2d pmidcurve2d(0.5*(aP2d1.XY() + aP2d2.XY()));
2251             const gp_Pnt aPC(anAS.Value(pmidcurve2d.X(), pmidcurve2d.Y()));
2252             const Standard_Real aDist = PMid.Distance(aPC);
2253             Tol = Max(aDist, Tol);
2254             //Check same parameter in middle point .end
2255           }
2256         }
2257       }
2258     }
2259   }
2260   //
2261   if (S->IsUPeriodic() && !C2d.IsNull()) {
2262     // Recadre dans le domaine UV de la face
2263     Standard_Real aTm, U0, aEps, period, du, U0x;
2264     Standard_Boolean bAdjust;
2265     //
2266     aEps = Precision::PConfusion();
2267     period = S->UPeriod();
2268     //
2269     aTm = .5*(f + l);
2270     gp_Pnt2d pm = C2d->Value(aTm);
2271     U0 = pm.X();
2272     //
2273     bAdjust = 
2274       GeomInt::AdjustPeriodic(U0, umin, umax, period, U0x, du, aEps);
2275     if (bAdjust) {
2276       gp_Vec2d T1(du, 0.);
2277       C2d->Translate(T1);
2278     }
2279   }
2280 }
2281
2282 //=======================================================================
2283 //function : TrimILineOnSurfBoundaries
2284 //purpose  : This function finds intersection points of given curves with
2285 //            surface boundaries and fills theArrayOfParameters by parameters
2286 //            along the given curves corresponding of these points.
2287 //=======================================================================
2288 void GeomInt_IntSS::TrimILineOnSurfBoundaries(const Handle(Geom2d_Curve)& theC2d1,
2289                                               const Handle(Geom2d_Curve)& theC2d2,
2290                                               const Bnd_Box2d& theBound1,
2291                                               const Bnd_Box2d& theBound2,
2292                                               GeomInt_VectorOfReal& theArrayOfParameters)
2293 {
2294   //Rectangular boundaries of two surfaces: [0]:U=Ufirst, [1]:U=Ulast,
2295   //                                        [2]:V=Vfirst, [3]:V=Vlast 
2296   const Standard_Integer aNumberOfCurves = 4;
2297   Handle(Geom2d_Curve) aCurS1Bounds[aNumberOfCurves];
2298   Handle(Geom2d_Curve) aCurS2Bounds[aNumberOfCurves];
2299
2300   Standard_Real aU1f=0.0, aV1f=0.0, aU1l=0.0, aV1l=0.0;
2301   Standard_Real aU2f=0.0, aV2f=0.0, aU2l=0.0, aV2l=0.0;
2302
2303   theBound1.Get(aU1f, aV1f, aU1l, aV1l);
2304   theBound2.Get(aU2f, aV2f, aU2l, aV2l);
2305
2306   Standard_Real aDelta = aV1l-aV1f;
2307   if(Abs(aDelta) > RealSmall())
2308   {
2309     if(!Precision::IsInfinite(aU1f))
2310     {
2311       aCurS1Bounds[0] = new Geom2d_Line(gp_Pnt2d(aU1f, aV1f), gp_Dir2d(0.0, 1.0));
2312
2313       if(!Precision::IsInfinite(aDelta))
2314         aCurS1Bounds[0] = new Geom2d_TrimmedCurve(aCurS1Bounds[0], 0, aDelta);
2315     }
2316
2317     if(!Precision::IsInfinite(aU1l))
2318     {
2319       aCurS1Bounds[1] = new Geom2d_Line(gp_Pnt2d(aU1l, aV1f), gp_Dir2d(0.0, 1.0));
2320       if(!Precision::IsInfinite(aDelta))
2321         aCurS1Bounds[1] = new Geom2d_TrimmedCurve(aCurS1Bounds[1], 0, aDelta);
2322     }
2323   }
2324
2325   aDelta = aU1l-aU1f;
2326   if(Abs(aDelta) > RealSmall())
2327   {
2328     if(!Precision::IsInfinite(aV1f))
2329     {
2330       aCurS1Bounds[2] = new Geom2d_Line(gp_Pnt2d(aU1f, aV1f), gp_Dir2d(1.0, 0.0));
2331       if(!Precision::IsInfinite(aDelta))
2332         aCurS1Bounds[2] = new Geom2d_TrimmedCurve(aCurS1Bounds[2], 0, aDelta);
2333     }
2334
2335     if(!Precision::IsInfinite(aV1l))
2336     {
2337       aCurS1Bounds[3] = new Geom2d_Line(gp_Pnt2d(aU1l, aV1l), gp_Dir2d(1.0, 0.0));
2338       if(!Precision::IsInfinite(aDelta))
2339         aCurS1Bounds[3] = new Geom2d_TrimmedCurve(aCurS1Bounds[3], 0, aDelta);
2340     }
2341   }
2342
2343   aDelta = aV2l-aV2f;
2344   if(Abs(aDelta) > RealSmall())
2345   {
2346     if(!Precision::IsInfinite(aU2f))
2347     {
2348       aCurS2Bounds[0] = new Geom2d_Line(gp_Pnt2d(aU2f, aV2f), gp_Dir2d(0.0, 1.0));
2349       if(!Precision::IsInfinite(aDelta))
2350         aCurS2Bounds[0] = new Geom2d_TrimmedCurve(aCurS2Bounds[0], 0, aDelta);
2351     }
2352
2353     if(!Precision::IsInfinite(aU2l))
2354     {
2355       aCurS2Bounds[1] = new Geom2d_Line(gp_Pnt2d(aU2l, aV2f), gp_Dir2d(0.0, 1.0));
2356       if(!Precision::IsInfinite(aDelta))
2357         aCurS2Bounds[1] = new Geom2d_TrimmedCurve(aCurS2Bounds[1], 0, aDelta);
2358     }
2359   }
2360
2361   aDelta = aU2l-aU2f;
2362   if(Abs(aDelta) > RealSmall())
2363   {
2364     if(!Precision::IsInfinite(aV2f))
2365     {
2366       aCurS2Bounds[2] = new Geom2d_Line(gp_Pnt2d(aU2f, aV2f), gp_Dir2d(1.0, 0.0));
2367       if(!Precision::IsInfinite(aDelta))
2368         aCurS2Bounds[2] = new Geom2d_TrimmedCurve(aCurS2Bounds[2], 0, aDelta);
2369     }
2370
2371     if(!Precision::IsInfinite(aV2l))
2372     {
2373       aCurS2Bounds[3] = new Geom2d_Line(gp_Pnt2d(aU2l, aV2l), gp_Dir2d(1.0, 0.0));
2374       if(!Precision::IsInfinite(aDelta))
2375         aCurS2Bounds[3] = new Geom2d_TrimmedCurve(aCurS2Bounds[3], 0, aDelta);
2376     }
2377   }
2378
2379   const Standard_Real anIntTol = 10.0*Precision::Confusion();
2380
2381   IntersectCurveAndBoundary(theC2d1, aCurS1Bounds,
2382                         aNumberOfCurves, anIntTol, theArrayOfParameters);
2383
2384   IntersectCurveAndBoundary(theC2d2, aCurS2Bounds,
2385                         aNumberOfCurves, anIntTol, theArrayOfParameters);
2386
2387   std::sort(theArrayOfParameters.begin(), theArrayOfParameters.end());
2388 }