f9af783a78e93fbc09b97253940f1cb4b122a934
[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-2012 OPEN CASCADE SAS
5 //
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
10 //
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13 //
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
20
21
22 #include <GeomInt_IntSS.ixx>
23
24 #include <Precision.hxx>
25 #include <gp_Pnt.hxx>
26 #include <gp_Pnt2d.hxx>
27 #include <gp_Pln.hxx>
28
29 #include <TColStd_Array1OfReal.hxx>
30 #include <TColStd_Array1OfInteger.hxx>
31 #include <TColStd_Array1OfListOfInteger.hxx>
32 #include <TColStd_SequenceOfReal.hxx>
33 #include <TColStd_ListOfInteger.hxx>
34 #include <TColStd_ListIteratorOfListOfInteger.hxx>
35
36 #include <TColgp_Array1OfPnt.hxx>
37 #include <TColgp_Array1OfPnt2d.hxx>
38
39 #include <Adaptor3d_TopolTool.hxx>
40 #include <IntPatch_Line.hxx>
41 #include <IntPatch_WLine.hxx>
42 #include <IntPatch_GLine.hxx>
43 #include <IntPatch_ALineToWLine.hxx>
44 #include <IntPatch_IType.hxx>
45 #include <NCollection_IncAllocator.hxx>
46 #include <NCollection_List.hxx>
47 #include <NCollection_LocalArray.hxx>
48 #include <NCollection_StdAllocator.hxx>
49 #include <vector>
50
51 #include <AppParCurves_MultiBSpCurve.hxx>
52
53 #include <GeomAbs_SurfaceType.hxx>
54 #include <GeomAbs_CurveType.hxx>
55
56 #include <Geom_Curve.hxx>
57 #include <Geom_Line.hxx>
58 #include <Geom_Parabola.hxx>
59 #include <Geom_Hyperbola.hxx>
60 #include <Geom_TrimmedCurve.hxx>
61 #include <Geom_Circle.hxx>
62 #include <Geom_Ellipse.hxx>
63 #include <Geom_BSplineCurve.hxx>
64
65 #include <Geom2d_Curve.hxx>
66 #include <Geom2d_BSplineCurve.hxx>
67 #include <Geom2d_TrimmedCurve.hxx>
68
69 #include <GeomLib_CheckBSplineCurve.hxx>
70 #include <GeomLib_Check2dBSplineCurve.hxx>
71 #include <GeomProjLib.hxx>
72
73
74 #include <ElSLib.hxx>
75
76 #include <GeomInt_WLApprox.hxx>
77 #include <Extrema_ExtPS.hxx>
78
79 static
80   void BuildPCurves (Standard_Real f,
81                      Standard_Real l,
82                      Standard_Real& Tol,
83                      const Handle (Geom_Surface)& S,
84                      const Handle (Geom_Curve)&   C,
85                      Handle (Geom2d_Curve)& C2d);
86
87 static
88   void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
89                   const Handle(GeomAdaptor_HSurface)& HS2,
90                   const gp_Pnt& Ptref,
91                   Standard_Real& U1,
92                   Standard_Real& V1,
93                   Standard_Real& U2,
94                   Standard_Real& V2);
95
96 static 
97   Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
98                                    const Standard_Integer ideb,
99                                    const Standard_Integer ifin);
100
101 static
102   Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
103                                             const Standard_Integer                         ideb,
104                                             const Standard_Integer                         ifin,
105                                             const Standard_Boolean                         onFirst);
106
107 static 
108   Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine);
109
110
111 static
112   Standard_Boolean DecompositionOfWLine  (const Handle(IntPatch_WLine)& theWLine,
113                                           const Handle(GeomAdaptor_HSurface)&            theSurface1, 
114                                           const Handle(GeomAdaptor_HSurface)&            theSurface2,
115                                           const Standard_Real aTolSum, 
116                                           const GeomInt_LineConstructor&                 theLConstructor,
117                                           IntPatch_SequenceOfLine&                       theNewLines);
118
119 static
120   Standard_Real AdjustPeriodic(const Standard_Real theParameter,
121                                const Standard_Real parmin,
122                                const Standard_Real parmax,
123                                const Standard_Real thePeriod,
124                                Standard_Real&      theOffset) ;
125 static 
126   Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
127                                      const Standard_Real theFirstBoundary,
128                                      const Standard_Real theSecondBoundary,
129                                      const Standard_Real theResolution,
130                                      Standard_Boolean&   IsOnFirstBoundary);
131
132 static 
133   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
134                              const gp_Pnt2d&     theLastPoint,
135                              const Standard_Real theUmin, 
136                              const Standard_Real theUmax,
137                              const Standard_Real theVmin,
138                              const Standard_Real theVmax,
139                              gp_Pnt2d&           theNewPoint);
140
141
142 static
143   void AdjustUPeriodic (const Handle(Geom_Surface)& aS,
144                         Handle(Geom2d_Curve)& aC2D);
145 static
146   void GetQuadric(const Handle(GeomAdaptor_HSurface)& HS1,
147                   IntSurf_Quadric& quad1);
148
149
150 //
151
152
153 class ProjectPointOnSurf
154 {
155  public:
156   ProjectPointOnSurf() : myIsDone (Standard_False) {}
157   void Init(const Handle(Geom_Surface)& Surface,
158             const Standard_Real Umin,
159             const Standard_Real Usup,
160             const Standard_Real Vmin,
161             const Standard_Real Vsup);
162   void Init ();
163   void Perform(const gp_Pnt& P);
164   Standard_Boolean IsDone () const { return myIsDone; }
165   void LowerDistanceParameters(Standard_Real& U, Standard_Real& V ) const;
166   Standard_Real LowerDistance() const;
167  protected:
168   Standard_Boolean myIsDone;
169   Standard_Integer myIndex;
170   Extrema_ExtPS myExtPS;
171   GeomAdaptor_Surface myGeomAdaptor;
172 };
173
174 //=======================================================================
175 //function : Init
176 //purpose  : 
177 //=======================================================================
178   void ProjectPointOnSurf::Init (const Handle(Geom_Surface)& Surface,
179                                  const Standard_Real       Umin,
180                                  const Standard_Real       Usup,
181                                  const Standard_Real       Vmin,
182                                  const Standard_Real       Vsup )
183 {
184   const Standard_Real Tolerance = Precision::PConfusion();
185   //
186   myGeomAdaptor.Load(Surface, Umin,Usup,Vmin,Vsup);
187   myExtPS.Initialize(myGeomAdaptor, Umin, Usup, Vmin, Vsup, Tolerance, Tolerance);
188   myIsDone = Standard_False;
189 }
190 //=======================================================================
191 //function : Init
192 //purpose  : 
193 //=======================================================================
194   void ProjectPointOnSurf::Init ()
195 {
196   myIsDone = myExtPS.IsDone() && (myExtPS.NbExt() > 0);
197   if (myIsDone) {
198     // evaluate the lower distance and its index;
199     Standard_Real Dist2Min = myExtPS.SquareDistance(1);
200     myIndex = 1;
201     for (Standard_Integer i = 2; i <= myExtPS.NbExt(); i++)
202     {
203       const Standard_Real Dist2 = myExtPS.SquareDistance(i);
204       if (Dist2 < Dist2Min) {
205         Dist2Min = Dist2;
206         myIndex = i;
207       }
208     }
209   }
210 }
211 //=======================================================================
212 //function : Perform
213 //purpose  : 
214 //=======================================================================
215 void ProjectPointOnSurf::Perform(const gp_Pnt& P)
216 {
217   myExtPS.Perform(P);
218   Init ();
219 }
220 //=======================================================================
221 //function : LowerDistanceParameters
222 //purpose  : 
223 //=======================================================================
224 void ProjectPointOnSurf::LowerDistanceParameters (Standard_Real& U,
225                                                   Standard_Real& V ) const
226 {
227   StdFail_NotDone_Raise_if(!myIsDone, "GeomInt_IntSS::ProjectPointOnSurf::LowerDistanceParameters");
228   (myExtPS.Point(myIndex)).Parameter(U,V);
229 }
230 //=======================================================================
231 //function : LowerDistance
232 //purpose  : 
233 //=======================================================================
234 Standard_Real ProjectPointOnSurf::LowerDistance() const
235 {
236   StdFail_NotDone_Raise_if(!myIsDone, "GeomInt_IntSS::ProjectPointOnSurf::LowerDistance");
237   return sqrt(myExtPS.SquareDistance(myIndex));
238 }
239 //
240
241 //=======================================================================
242 //function : MakeCurve
243 //purpose  : 
244 //=======================================================================
245   void GeomInt_IntSS::MakeCurve(const Standard_Integer Index,
246                                 const Handle(Adaptor3d_TopolTool) & dom1,
247                                 const Handle(Adaptor3d_TopolTool) & dom2,
248                                 const Standard_Real Tol,
249                                 const Standard_Boolean Approx,
250                                 const Standard_Boolean ApproxS1,
251                                 const Standard_Boolean ApproxS2)
252
253 {
254   Standard_Boolean myApprox1, myApprox2, myApprox;
255   Standard_Real Tolpc, myTolApprox;
256   IntPatch_IType typl;
257   Handle(Geom2d_BSplineCurve) H1;
258   Handle(Geom_Surface) aS1, aS2;
259   //
260   Tolpc = Tol;
261   myApprox=Approx;
262   myApprox1=ApproxS1;
263   myApprox2=ApproxS2;
264   myTolApprox=0.0000001;
265   //
266   aS1=myHS1->ChangeSurface().Surface();
267   aS2=myHS2->ChangeSurface().Surface();
268   //
269   Handle(IntPatch_Line) L = myIntersector.Line(Index);
270   typl = L->ArcType();
271   //
272   if(typl==IntPatch_Walking) {
273     Handle(IntPatch_Line) anewL = 
274       ComputePurgedWLine(Handle(IntPatch_WLine)::DownCast(L));
275     //
276     if(anewL.IsNull()) {
277       return;
278     }
279     L = anewL;
280   }
281   //
282   // Line Constructor
283   myLConstruct.Perform(L);
284   if (!myLConstruct.IsDone() || myLConstruct.NbParts() <= 0) {
285     return;
286   }
287   // Do the Curve
288   Standard_Boolean ok;
289   Standard_Integer i, j,  aNbParts;
290   Standard_Real fprm, lprm;
291   Handle(Geom_Curve) newc;
292
293   switch (typl) {
294   //########################################  
295   // Line, Parabola, Hyperbola
296   //########################################  
297   case IntPatch_Lin:
298   case IntPatch_Parabola: 
299   case IntPatch_Hyperbola: {
300     if (typl == IntPatch_Lin) {
301       newc=new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
302     }
303     else if (typl == IntPatch_Parabola) {
304       newc=new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
305     }
306     else if (typl == IntPatch_Hyperbola) {
307       newc=new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
308     }
309     //
310     aNbParts=myLConstruct.NbParts();
311     for (i=1; i<=aNbParts; i++) {
312       myLConstruct.Part(i, fprm, lprm);
313       
314       if (!Precision::IsNegativeInfinite(fprm) && 
315           !Precision::IsPositiveInfinite(lprm)) {
316         // /cto/900/F1
317         //if (typl == IntPatch_Lin) {
318           //Standard_Real dPrm=1.3e-6;
319           //fprm=fprm-dPrm;
320           //lprm=lprm+dPrm;
321         //}
322         //
323         Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
324         sline.Append(aCT3D);
325         //
326         if(myApprox1) { 
327           Handle (Geom2d_Curve) C2d;
328           BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
329           if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
330             myTolReached2d=Tolpc;
331           }
332           slineS1.Append(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
333         }
334         else { 
335           slineS1.Append(H1);
336         }
337         //
338         if(myApprox2) { 
339     Handle (Geom2d_Curve) C2d;
340     BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
341     if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
342       myTolReached2d=Tolpc;
343     }
344     //
345     slineS2.Append(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
346   }
347   else { 
348     slineS2.Append(H1);
349   }
350       } // if (!Precision::IsNegativeInfinite(fprm) &&  !Precision::IsPositiveInfinite(lprm))
351
352       else {
353         GeomAbs_SurfaceType typS1 = myHS1->Surface().GetType();
354         GeomAbs_SurfaceType typS2 = myHS2->Surface().GetType();
355         if( typS1 == GeomAbs_SurfaceOfExtrusion ||
356           typS1 == GeomAbs_OffsetSurface ||
357           typS1 == GeomAbs_SurfaceOfRevolution ||
358           typS2 == GeomAbs_SurfaceOfExtrusion ||
359           typS2 == GeomAbs_OffsetSurface ||
360           typS2 == GeomAbs_SurfaceOfRevolution) {
361             sline.Append(newc);
362             slineS1.Append(H1);
363             slineS2.Append(H1);
364             continue;
365         }
366         Standard_Boolean bFNIt, bLPIt;
367         Standard_Real aTestPrm, dT=100.;
368         Standard_Real u1, v1, u2, v2, TolX;
369         //
370         bFNIt=Precision::IsNegativeInfinite(fprm);
371         bLPIt=Precision::IsPositiveInfinite(lprm);
372
373         aTestPrm=0.;
374
375         if (bFNIt && !bLPIt) {
376           aTestPrm=lprm-dT;
377         }
378         else if (!bFNIt && bLPIt) {
379           aTestPrm=fprm+dT;
380         }
381         //
382         gp_Pnt ptref(newc->Value(aTestPrm));
383         //
384         TolX = Precision::Confusion();
385         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
386         ok = (dom1->Classify(gp_Pnt2d(u1, v1), TolX) != TopAbs_OUT);
387         if(ok) { 
388           ok = (dom2->Classify(gp_Pnt2d(u2,v2),TolX) != TopAbs_OUT); 
389         }
390         if (ok) {
391           sline.Append(newc);
392           slineS1.Append(H1);
393           slineS2.Append(H1);
394         }
395       }
396     }// end of for (i=1; i<=myLConstruct.NbParts(); i++)
397    }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
398    break;
399
400   //########################################  
401   // Circle and Ellipse
402   //########################################  
403   case IntPatch_Circle: 
404   case IntPatch_Ellipse: {
405
406     if (typl == IntPatch_Circle) {
407       newc = new Geom_Circle
408         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
409     }
410     else { 
411       newc = new Geom_Ellipse
412         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
413     }
414     //
415     Standard_Real aPeriod, aNul, aRealEpsilon;
416     //
417     aRealEpsilon=RealEpsilon();
418     aNul=0.;
419     aPeriod=M_PI+M_PI;
420     //
421     aNbParts=myLConstruct.NbParts();
422     //
423     for (i=1; i<=aNbParts; i++) {
424       myLConstruct.Part(i, fprm, lprm);
425       //
426       if (Abs(fprm) > aRealEpsilon || Abs(lprm-aPeriod) > aRealEpsilon) {
427         //==============================================
428         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
429         //
430         sline.Append(aTC3D);
431         //
432         fprm=aTC3D->FirstParameter();
433         lprm=aTC3D->LastParameter ();
434         ////    
435         if(myApprox1) { 
436           Handle (Geom2d_Curve) C2d;
437           BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
438           if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
439             myTolReached2d=Tolpc;
440           }
441           slineS1.Append(C2d);
442         }
443         else { //// 
444           slineS1.Append(H1);
445         }
446         //
447         if(myApprox2) { 
448           Handle (Geom2d_Curve) C2d;
449           BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
450           if(Tolpc>myTolReached2d || myTolReached2d==0) { 
451             myTolReached2d=Tolpc;
452           }
453           slineS2.Append(C2d);
454         }
455         else { 
456           slineS2.Append(H1);  
457         }
458         //==============================================        
459       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
460       //
461       else {//  on regarde si on garde
462         //
463         if (aNbParts==1) {
464           if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
465             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
466             //
467             sline.Append(aTC3D);
468             fprm=aTC3D->FirstParameter();
469             lprm=aTC3D->LastParameter ();
470             
471             if(myApprox1) { 
472               Handle (Geom2d_Curve) C2d;
473               BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
474               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
475                 myTolReached2d=Tolpc;
476               }
477               slineS1.Append(C2d);
478             }
479             else { //// 
480               slineS1.Append(H1);
481             }
482
483             if(myApprox2) { 
484               Handle (Geom2d_Curve) C2d;
485               BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
486               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
487                 myTolReached2d=Tolpc;
488               }
489               slineS2.Append(C2d);
490             }
491             else { 
492               slineS2.Append(H1);
493             }
494             break;
495           }
496         }
497         //
498         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, TolX;
499         //
500         aTwoPIdiv17=2.*M_PI/17.;
501         //
502         for (j=0; j<=17; j++) {
503           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
504           TolX = Precision::Confusion();
505           
506           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
507           ok = (dom1->Classify(gp_Pnt2d(u1,v1),TolX) != TopAbs_OUT);
508           if(ok) { 
509             ok = (dom2->Classify(gp_Pnt2d(u2,v2),TolX) != TopAbs_OUT);
510           }
511           if (ok) {
512             sline.Append(newc);
513             //==============================================
514             if(myApprox1) { 
515               Handle (Geom2d_Curve) C2d;
516               BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
517               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
518                 myTolReached2d=Tolpc;
519               }
520               slineS1.Append(C2d);
521             }
522             else { 
523               slineS1.Append(H1);  
524             }
525             
526             if(myApprox2) { 
527               Handle (Geom2d_Curve) C2d;
528               BuildPCurves(fprm, lprm, Tolpc,myHS2->ChangeSurface().Surface(), newc, C2d);
529               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
530                 myTolReached2d=Tolpc;
531               } 
532               slineS2.Append(C2d);
533             }
534             else { 
535               slineS2.Append(H1);  
536             }
537             break;
538           }//  end of if (ok) {
539         }//  end of for (Standard_Integer j=0; j<=17; j++)
540       }//  end of else { on regarde si on garde
541     }// for (i=1; i<=myLConstruct.NbParts(); i++)
542   }// IntPatch_Circle: IntPatch_Ellipse:
543     break;
544     
545     //########################################  
546     // Analytic
547     //######################################## 
548   case IntPatch_Analytic: {
549     IntSurf_Quadric quad1,quad2;
550     //
551     GetQuadric(myHS1, quad1);
552     GetQuadric(myHS2, quad2);
553     //=========
554     IntPatch_ALineToWLine convert (quad1, quad2);
555       
556     if (!myApprox) {
557       Handle(Geom2d_BSplineCurve) aH1, aH2;
558       //
559       aNbParts=myLConstruct.NbParts();
560       for (i=1; i<=aNbParts; i++) {
561         myLConstruct.Part(i, fprm, lprm);
562         Handle(IntPatch_WLine) WL = 
563           convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
564         //
565         if(myApprox1) {
566           aH1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
567         }
568         
569         if(myApprox2) {
570           aH2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
571         }
572         sline.Append(MakeBSpline(WL,1,WL->NbPnts()));
573         slineS1.Append(aH1);
574         slineS2.Append(aH2);
575       }
576     } // if (!myApprox)
577
578     else { // myApprox=TRUE
579       GeomInt_WLApprox theapp3d;
580       Standard_Real tol2d = myTolApprox;
581       //        
582       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
583       
584       aNbParts=myLConstruct.NbParts();
585       for (i=1; i<=aNbParts; i++) {
586         myLConstruct.Part(i, fprm, lprm);
587         Handle(IntPatch_WLine) WL = 
588           convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
589
590         theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
591         if (!theapp3d.IsDone()) {
592           //
593           Handle(Geom2d_BSplineCurve) aH1, aH2;
594
595           if(myApprox1) {
596             aH1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
597           }
598           
599           if(myApprox2) {
600             aH2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
601           }
602           sline.Append(MakeBSpline(WL,1,WL->NbPnts()));
603           slineS1.Append(aH1);
604           slineS2.Append(aH1);
605         }
606         //
607         else {
608           if(myApprox1 || myApprox2) { 
609             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) { 
610               myTolReached2d = theapp3d.TolReached2d();
611             }
612           }
613           
614           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) { 
615             myTolReached3d = theapp3d.TolReached3d();
616           }
617           
618           Standard_Integer aNbMultiCurves, nbpoles;
619           aNbMultiCurves=theapp3d.NbMultiCurves();
620           for (j=1; j<=aNbMultiCurves; j++) {
621             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
622             
623             nbpoles = mbspc.NbPoles();
624             TColgp_Array1OfPnt tpoles(1, nbpoles);
625             mbspc.Curve(1, tpoles);
626             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
627                                                                mbspc.Knots(),
628                                                                mbspc.Multiplicities(),
629                                                                mbspc.Degree());
630             
631             GeomLib_CheckBSplineCurve Check(BS, myTolCheck, myTolAngCheck);
632             Check.FixTangent(Standard_True,Standard_True);
633             // 
634             sline.Append(BS);
635             //
636             if(myApprox1) { 
637               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
638               mbspc.Curve(2,tpoles2d);
639               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
640                                                                       mbspc.Knots(),
641                                                                       mbspc.Multiplicities(),
642                                                                       mbspc.Degree());
643
644               GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
645               newCheck.FixTangent(Standard_True,Standard_True);
646               slineS1.Append(BS2);              
647             }
648             else {
649               slineS1.Append(H1);
650             }
651             
652             if(myApprox2) { 
653               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
654               Standard_Integer TwoOrThree;
655               TwoOrThree=myApprox1 ? 3 : 2;
656               mbspc.Curve(TwoOrThree, tpoles2d);
657               Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
658                                                                        mbspc.Knots(),
659                                                                        mbspc.Multiplicities(),
660                                                                        mbspc.Degree());
661                 
662               GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
663               newCheck.FixTangent(Standard_True,Standard_True);
664               //        
665               slineS2.Append(BS2);
666             }
667             else { 
668               slineS2.Append(H1);
669             }
670             // 
671           }// for (j=1; j<=aNbMultiCurves; j++) {
672         }// else from if (!theapp3d.IsDone())
673       }// for (i=1; i<=aNbParts; i++) {
674     }// else { // myApprox=TRUE
675   }// case IntPatch_Analytic:
676     break;
677     
678     //########################################  
679     // Walking
680     //######################################## 
681   case IntPatch_Walking:{
682     Handle(IntPatch_WLine) WL = 
683       Handle(IntPatch_WLine)::DownCast(L);
684     //
685     Standard_Integer ifprm, ilprm;
686     //
687     if (!myApprox) {
688       aNbParts=myLConstruct.NbParts();
689       for (i=1; i<=aNbParts; i++) {
690         myLConstruct.Part(i, fprm, lprm);
691         ifprm=(Standard_Integer)fprm;
692         ilprm=(Standard_Integer)lprm;
693         //        
694         Handle(Geom2d_BSplineCurve) aH1, aH2;
695
696         if(myApprox1) {
697           aH1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
698         }
699         if(myApprox2) {
700           aH2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
701         }
702         //
703         Handle(Geom_Curve) aBSp=MakeBSpline(WL, ifprm, ilprm);
704         //      
705         sline.Append(aBSp);
706         slineS1.Append(aH1);
707         slineS2.Append(aH2);
708       }
709     }
710     //
711     else {
712       Standard_Boolean bIsDecomposited;
713       Standard_Integer nbiter, aNbSeqOfL;
714       GeomInt_WLApprox theapp3d;
715       IntPatch_SequenceOfLine aSeqOfL;
716       Standard_Real tol2d, aTolSS;
717       //        
718       tol2d = myTolApprox;
719       aTolSS=2.e-7;
720       if(myHS1 == myHS2) { 
721         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False);      
722       }
723       else { 
724         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
725       }
726       //
727       bIsDecomposited = 
728         DecompositionOfWLine(WL, myHS1, myHS2, aTolSS, myLConstruct, aSeqOfL);
729       //
730       aNbParts=myLConstruct.NbParts();
731       aNbSeqOfL=aSeqOfL.Length();
732       //
733       nbiter = (bIsDecomposited) ? aNbSeqOfL : aNbParts;
734       //
735       for(i = 1; i <= nbiter; i++) {
736         if(bIsDecomposited) {
737           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
738           ifprm = 1;
739           ilprm = WL->NbPnts();
740         }
741         else {
742           myLConstruct.Part(i, fprm, lprm);
743           ifprm = (Standard_Integer)fprm;
744           ilprm = (Standard_Integer)lprm;
745         }
746         //-- lbr : 
747         //-- Si une des surfaces est un plan , on approxime en 2d
748         //-- sur cette surface et on remonte les points 2d en 3d.
749         GeomAbs_SurfaceType typs1, typs2;
750         typs1 = myHS1->Surface().GetType();
751         typs2 = myHS2->Surface().GetType();
752         //
753         if(typs1 == GeomAbs_Plane) { 
754           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,
755                            Standard_True, myApprox2,
756                            ifprm,  ilprm);
757         }         
758         else if(typs2 == GeomAbs_Plane) { 
759           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,
760                            myApprox1,Standard_True,
761                            ifprm,  ilprm);
762         }
763         else { 
764           //
765           if (myHS1 != myHS2){
766             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
767                 (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
768              
769               theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
770               //Standard_Boolean bUseSurfaces;
771               //bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
772               //if (bUseSurfaces) {
773                 //theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False);
774               //}
775             }
776           }
777           //
778           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,
779                            myApprox1,myApprox2,
780                            ifprm,  ilprm);
781         }
782          
783         if (!theapp3d.IsDone()) {
784           //      
785           Handle(Geom2d_BSplineCurve) aH1, aH2;
786           //      
787           Handle(Geom_Curve) aBSp=MakeBSpline(WL, ifprm, ilprm);
788           if(myApprox1) {
789             aH1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
790           }
791           if(myApprox2) {
792             aH2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
793           }
794           //
795           sline.Append(aBSp);
796           slineS1.Append(aH1);
797           slineS2.Append(aH2);
798         }//if (!theapp3d.IsDone())
799         
800         else {
801           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) { 
802             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) { 
803               myTolReached2d = theapp3d.TolReached2d();
804             }
805           }
806           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) { 
807             myTolReached3d = myTolReached2d;
808           }
809           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) { 
810             myTolReached3d = theapp3d.TolReached3d();
811           }
812           
813           Standard_Integer aNbMultiCurves, nbpoles;
814           //
815           aNbMultiCurves=theapp3d.NbMultiCurves(); 
816           for (j=1; j<=aNbMultiCurves; j++) {
817             if(typs1 == GeomAbs_Plane) {
818               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
819               nbpoles = mbspc.NbPoles();
820               
821               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
822               TColgp_Array1OfPnt   tpoles(1,nbpoles);
823               
824               mbspc.Curve(1,tpoles2d);
825               const gp_Pln&  Pln = myHS1->Surface().Plane();
826               //
827               Standard_Integer ik; 
828               for(ik = 1; ik<= nbpoles; ik++) { 
829                 tpoles.SetValue(ik,
830                                 ElSLib::Value(tpoles2d.Value(ik).X(),
831                                               tpoles2d.Value(ik).Y(),
832                                               Pln));
833               }
834               //
835               Handle(Geom_BSplineCurve) BS = 
836                 new Geom_BSplineCurve(tpoles,
837                                       mbspc.Knots(),
838                                       mbspc.Multiplicities(),
839                                       mbspc.Degree());
840               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
841               Check.FixTangent(Standard_True, Standard_True);
842               //        
843               sline.Append(BS);
844               //
845               if(myApprox1) { 
846                 Handle(Geom2d_BSplineCurve) BS1 = 
847                   new Geom2d_BSplineCurve(tpoles2d,
848                                           mbspc.Knots(),
849                                           mbspc.Multiplicities(),
850                                           mbspc.Degree());
851                 GeomLib_Check2dBSplineCurve Check1(BS1,myTolCheck,myTolAngCheck);
852                 Check1.FixTangent(Standard_True,Standard_True);
853                 //      
854                 AdjustUPeriodic (aS1, BS1);  
855                 //
856                 slineS1.Append(BS1);
857               }
858               else {
859                 slineS1.Append(H1);
860               }
861
862               if(myApprox2) { 
863                 mbspc.Curve(2, tpoles2d);
864                 
865                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
866                                                                           mbspc.Knots(),
867                                                                           mbspc.Multiplicities(),
868                                                                           mbspc.Degree());
869                 GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
870                 newCheck.FixTangent(Standard_True,Standard_True);
871                 //
872                 AdjustUPeriodic (aS2, BS2);  
873                 //
874                 slineS2.Append(BS2); 
875               }
876               else { 
877                 slineS2.Append(H1); 
878               }
879             }//if(typs1 == GeomAbs_Plane) 
880             //
881             else if(typs2 == GeomAbs_Plane) { 
882               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
883               nbpoles = mbspc.NbPoles();
884               
885               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
886               TColgp_Array1OfPnt   tpoles(1,nbpoles);
887               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
888               const gp_Pln&  Pln = myHS2->Surface().Plane();
889               //
890               Standard_Integer ik; 
891               for(ik = 1; ik<= nbpoles; ik++) { 
892                 tpoles.SetValue(ik,
893                                 ElSLib::Value(tpoles2d.Value(ik).X(),
894                                               tpoles2d.Value(ik).Y(),
895                                               Pln));
896                 
897               }
898               //
899               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
900                                                                  mbspc.Knots(),
901                                                                  mbspc.Multiplicities(),
902                                                                  mbspc.Degree());
903               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
904               Check.FixTangent(Standard_True,Standard_True);
905               //        
906               sline.Append(BS);
907               //
908               if(myApprox2) {
909                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
910                                                                         mbspc.Knots(),
911                                                                         mbspc.Multiplicities(),
912                                                                         mbspc.Degree());
913                 GeomLib_Check2dBSplineCurve Check1(BS1,myTolCheck,myTolAngCheck);
914                 Check1.FixTangent(Standard_True,Standard_True);
915                 //
916                 //
917                 AdjustUPeriodic (aS2, BS1);  
918                 //
919                 slineS2.Append(BS1);
920               }
921               else {
922                 slineS2.Append(H1);
923               }
924               
925               if(myApprox1) { 
926                 mbspc.Curve(1,tpoles2d);
927                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
928                                                                         mbspc.Knots(),
929                                                                         mbspc.Multiplicities(),
930                                                                         mbspc.Degree());
931                 GeomLib_Check2dBSplineCurve Check2(BS2,myTolCheck,myTolAngCheck);
932                 Check2.FixTangent(Standard_True,Standard_True);
933                 // 
934                 //
935                 AdjustUPeriodic (aS1, BS2);  
936                 //      
937                 slineS1.Append(BS2);
938               }
939               else { 
940                 slineS1.Append(H1);
941               }
942             } // else if(typs2 == GeomAbs_Plane)
943             //
944             else { // typs1!=GeomAbs_Plane && typs2!=GeomAbs_Plane
945               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
946               nbpoles = mbspc.NbPoles();
947               TColgp_Array1OfPnt tpoles(1,nbpoles);
948               mbspc.Curve(1,tpoles);
949               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
950                                                                  mbspc.Knots(),
951                                                                  mbspc.Multiplicities(),
952                                                                  mbspc.Degree());
953               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
954               Check.FixTangent(Standard_True,Standard_True);
955               //        
956         //Check IsClosed()
957         Standard_Real aDist = Max(BS->StartPoint().XYZ().SquareModulus(),
958                                   BS->EndPoint().XYZ().SquareModulus());
959         Standard_Real eps = Epsilon(aDist);
960         if(BS->StartPoint().SquareDistance(BS->EndPoint()) < 2.*eps &&
961           !BS->IsClosed() && !BS->IsPeriodic())
962         {
963           //force Closed()
964           gp_Pnt aPm((BS->Pole(1).XYZ() + BS->Pole(BS->NbPoles()).XYZ()) / 2.);
965           BS->SetPole(1, aPm);
966           BS->SetPole(BS->NbPoles(), aPm);
967         }
968               sline.Append(BS);
969               
970               if(myApprox1) { 
971                 TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
972                 mbspc.Curve(2,tpoles2d);
973                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
974                                                                         mbspc.Knots(),
975                                                                         mbspc.Multiplicities(),
976                                                                         mbspc.Degree());
977                 GeomLib_Check2dBSplineCurve newCheck(BS1,myTolCheck,myTolAngCheck);
978                 newCheck.FixTangent(Standard_True,Standard_True);
979                 //
980                 AdjustUPeriodic (aS1, BS1);  
981                 //
982                 slineS1.Append(BS1);
983               }
984               else {
985                 slineS1.Append(H1);
986               }
987               if(myApprox2) { 
988                 TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
989                 mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
990                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
991                                                                         mbspc.Knots(),
992                                                                         mbspc.Multiplicities(),
993                                                                         mbspc.Degree());
994                 GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
995                 newCheck.FixTangent(Standard_True,Standard_True);
996                 //
997                 AdjustUPeriodic (aS2, BS2);  
998                 //
999                 slineS2.Append(BS2);
1000               }
1001               else { 
1002                 slineS2.Append(H1);
1003               }
1004             }// else { // typs1!=GeomAbs_Plane && typs2!=GeomAbs_Plane
1005           }// for (j=1; j<=aNbMultiCurves; j++
1006         }
1007       }
1008     }
1009   }
1010     break;
1011     
1012   case IntPatch_Restriction: 
1013     break;
1014
1015   }
1016 }
1017 //
1018 //=======================================================================
1019 //function : AdjustUPeriodic
1020 //purpose  : 
1021 //=======================================================================
1022  void AdjustUPeriodic (const Handle(Geom_Surface)& aS, Handle(Geom2d_Curve)& aC2D)
1023 {
1024   if (aC2D.IsNull() || !aS->IsUPeriodic())
1025     return;
1026   //
1027   const Standard_Real aEps=Precision::PConfusion();//1.e-9
1028   const Standard_Real aEpsilon=Epsilon(10.);//1.77e-15 
1029   //
1030   Standard_Real umin,umax,vmin,vmax;
1031   aS->Bounds(umin,umax,vmin,vmax);
1032   const Standard_Real aPeriod = aS->UPeriod();
1033   
1034   const Standard_Real aT1=aC2D->FirstParameter();
1035   const Standard_Real aT2=aC2D->LastParameter();
1036   const Standard_Real aTx=aT1+0.467*(aT2-aT1);
1037   const gp_Pnt2d aPx=aC2D->Value(aTx);
1038   //
1039   Standard_Real aUx=aPx.X();
1040   if (fabs(aUx)<aEpsilon)
1041     aUx=0.;
1042   if (fabs(aUx-aPeriod)<aEpsilon)
1043     aUx=aPeriod;
1044   //
1045   Standard_Real dU=0.;
1046   while(aUx <(umin-aEps)) {
1047     aUx+=aPeriod;
1048     dU+=aPeriod;
1049   }
1050   while(aUx>(umax+aEps)) {
1051     aUx-=aPeriod;
1052     dU-=aPeriod;
1053   }
1054   // 
1055   if (dU!=0.) {
1056     gp_Vec2d aV2D(dU, 0.);
1057     aC2D->Translate(aV2D);
1058   }
1059 }
1060 //
1061 //
1062 //=======================================================================
1063 //function : BuildPCurves
1064 //purpose  : 
1065 //=======================================================================
1066  void BuildPCurves (Standard_Real f, Standard_Real l,
1067                     Standard_Real& Tol,
1068                     const Handle(Geom_Surface)& S,
1069                     const Handle(Geom_Curve)& C,
1070                     Handle(Geom2d_Curve)& C2d)
1071 {
1072   if (C2d.IsNull())
1073   {
1074     C2d = GeomProjLib::Curve2d(C,f,l,S,Tol);
1075     if (S->IsUPeriodic())
1076     {
1077       Standard_Real umin,umax,vmin,vmax;
1078       S->Bounds(umin,umax,vmin,vmax);
1079       // Recadre dans le domaine UV de la face
1080       const Standard_Real period = S->UPeriod();
1081       gp_Pnt2d Pf=C2d->Value(f);
1082       //
1083       Standard_Real U0=Pf.X();
1084       //
1085       const Standard_Real aEpsilon=Epsilon(10.);//1.77e-15 
1086       if (fabs(U0)<aEpsilon)
1087         U0=0.;
1088       if (fabs(U0-period)<aEpsilon)
1089         U0=period;
1090       //
1091       const Standard_Real aEps=Precision::PConfusion();//1.e-9
1092       Standard_Real du=0.;
1093       while(U0 <(umin-aEps)) { 
1094         U0+=period;
1095         du+=period;
1096       }
1097       while(U0>(umax+aEps)) { 
1098         U0-=period;
1099         du-=period;
1100       }
1101       if (du!=0.) {
1102         gp_Vec2d T1(du, 0.);
1103         C2d->Translate(T1);
1104       }
1105     }
1106   }
1107 }
1108 //
1109 //=======================================================================
1110 //function : GetQuadric
1111 //purpose  : 
1112 //=======================================================================
1113  void GetQuadric(const Handle(GeomAdaptor_HSurface)& HS1, IntSurf_Quadric& quad1)
1114 {
1115   switch (HS1->Surface().GetType())
1116   {
1117     case GeomAbs_Plane:    quad1.SetValue(HS1->Surface().Plane()); break;
1118     case GeomAbs_Cylinder: quad1.SetValue(HS1->Surface().Cylinder()); break;
1119     case GeomAbs_Cone:     quad1.SetValue(HS1->Surface().Cone()); break;
1120     case GeomAbs_Sphere:   quad1.SetValue(HS1->Surface().Sphere()); break;
1121     default: Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
1122   }
1123 }
1124
1125 //=======================================================================
1126 //function : Parameters
1127 //purpose  : 
1128 //=======================================================================
1129  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
1130                  const Handle(GeomAdaptor_HSurface)& HS2,
1131                  const gp_Pnt& Ptref,
1132                  Standard_Real& U1,
1133                  Standard_Real& V1,
1134                  Standard_Real& U2,
1135                  Standard_Real& V2)
1136 {
1137   IntSurf_Quadric quad1,quad2;
1138   //
1139   GetQuadric(HS1, quad1);
1140   GetQuadric(HS2, quad2);
1141   //
1142   quad1.Parameters(Ptref,U1,V1);
1143   quad2.Parameters(Ptref,U2,V2);
1144 }
1145
1146 //=======================================================================
1147 //function : MakeBSpline
1148 //purpose  : 
1149 //=======================================================================
1150 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
1151                                  const Standard_Integer ideb,
1152                                  const Standard_Integer ifin)
1153 {
1154   const Standard_Integer nbpnt = ifin-ideb+1;
1155   TColgp_Array1OfPnt poles(1,nbpnt);
1156   TColStd_Array1OfReal knots(1,nbpnt);
1157   TColStd_Array1OfInteger mults(1,nbpnt);
1158   Standard_Integer i = 1, ipidebm1 = ideb;
1159   for(; i<=nbpnt; ipidebm1++, i++)
1160   {
1161     poles(i) = WL->Point(ipidebm1).Value();
1162     mults(i) = 1;
1163     knots(i) = i-1;
1164   }
1165   mults(1) = mults(nbpnt) = 2;
1166   return new Geom_BSplineCurve(poles,knots,mults,1);
1167 }
1168
1169 //=======================================================================
1170 //function : MakeBSpline2d
1171 //purpose  : 
1172 //=======================================================================
1173 Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
1174                                           const Standard_Integer ideb,
1175                                           const Standard_Integer ifin,
1176                                           const Standard_Boolean onFirst)
1177 {
1178   const Standard_Integer nbpnt = ifin-ideb+1;
1179   TColgp_Array1OfPnt2d poles(1,nbpnt);
1180   TColStd_Array1OfReal knots(1,nbpnt);
1181   TColStd_Array1OfInteger mults(1,nbpnt);
1182   Standard_Integer i = 1, ipidebm1 = ideb;
1183   for(; i <= nbpnt; ipidebm1++, i++)
1184   {
1185     Standard_Real U, V;
1186     if(onFirst)
1187           theWLine->Point(ipidebm1).ParametersOnS1(U, V);
1188     else
1189           theWLine->Point(ipidebm1).ParametersOnS2(U, V);
1190     poles(i).SetCoord(U, V);
1191     mults(i) = 1;
1192     knots(i) = i-1;
1193   }
1194   mults(1) = mults(nbpnt) = 2;
1195   return new Geom2d_BSplineCurve(poles,knots,mults,1);
1196 }
1197
1198 //=========================================================================
1199 // static function : ComputePurgedWLine
1200 // purpose : Removes equal points (leave one of equal points) from theWLine
1201 //           and recompute vertex parameters.
1202 //           Returns new WLine or null WLine if the number
1203 //           of the points is less than 2.
1204 //=========================================================================
1205 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine) 
1206 {
1207   Handle(IntPatch_WLine) aResult;
1208   Handle(IntPatch_WLine) aLocalWLine;
1209   Handle(IntPatch_WLine) aTmpWLine = theWLine;
1210
1211   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1212   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
1213   Standard_Integer i, k, v, nb, nbvtx;
1214   nbvtx = theWLine->NbVertex();
1215   nb = theWLine->NbPnts();
1216
1217   for(i = 1; i <= nb; i++) {
1218     aLineOn2S->Add(theWLine->Point(i));
1219   }
1220
1221   for(v = 1; v <= nbvtx; v++) {
1222     aLocalWLine->AddVertex(theWLine->Vertex(v));
1223   }
1224   
1225   for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
1226     Standard_Integer aStartIndex = i + 1;
1227     Standard_Integer anEndIndex = i + 5;
1228     nb = aLineOn2S->NbPoints();
1229     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
1230
1231     if((aStartIndex >= nb) || (anEndIndex <= 1)) {
1232       continue;
1233     }
1234     k = aStartIndex;
1235
1236     while(k <= anEndIndex) {
1237       
1238       if(i != k) {
1239         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
1240         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
1241         
1242         if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
1243           aTmpWLine = aLocalWLine;
1244           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
1245
1246           for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
1247             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
1248             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
1249
1250             if(avertexindex >= k) {
1251               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
1252             }
1253             aLocalWLine->AddVertex(aVertex);
1254           }
1255           aLineOn2S->RemovePoint(k);
1256           anEndIndex--;
1257           continue;
1258         }
1259       }
1260       k++;
1261     }
1262   }
1263
1264   if(aLineOn2S->NbPoints() > 1) {
1265     aResult = aLocalWLine;
1266   }
1267   return aResult;
1268 }
1269 //=======================================================================
1270 //function : DecompositionOfWLine
1271 //purpose  : 
1272 //=======================================================================
1273 Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
1274                                       const Handle(GeomAdaptor_HSurface)&            theSurface1, 
1275                                       const Handle(GeomAdaptor_HSurface)&            theSurface2,
1276                                       const Standard_Real aTolSum, 
1277                                       const GeomInt_LineConstructor&                 theLConstructor,
1278                                       IntPatch_SequenceOfLine&                       theNewLines)
1279 {
1280   typedef NCollection_List<Standard_Integer> ListOfInteger;
1281   //have to use std::vector, not NCollection_Vector in order to use copy constructor of
1282   //ListOfInteger which will be created with specific allocator instance
1283   typedef std::vector<ListOfInteger, NCollection_StdAllocator<
1284       ListOfInteger> > ArrayOfListOfInteger;
1285
1286   Standard_Boolean bIsPrevPointOnBoundary, bIsCurrentPointOnBoundary;
1287   Standard_Integer nblines, aNbPnts, aNbParts, pit, i, j, aNbListOfPointIndex;
1288   Standard_Real aTol, umin, umax, vmin, vmax;
1289
1290   //an inc allocator, it will contain wasted space (upon list's Clear()) but it should
1291   //still be faster than the standard allocator, and wasted memory should not be
1292   //significant and will be limited by time span of this function;
1293   //this is a separate allocator from the anIncAlloc below what provides better data
1294   //locality in the latter (by avoiding wastes which will only be in anIdxAlloc)
1295   Handle(NCollection_IncAllocator) anIdxAlloc = new NCollection_IncAllocator();
1296   ListOfInteger aListOfPointIndex (anIdxAlloc);
1297
1298   //GeomAPI_ProjectPointOnSurf aPrj1, aPrj2;
1299   ProjectPointOnSurf aPrj1, aPrj2;
1300   Handle(Geom_Surface) aSurf1,  aSurf2;
1301   //
1302   aNbParts=theLConstructor.NbParts();
1303   aNbPnts=theWLine->NbPnts();
1304   //
1305   if((!aNbPnts) || (!aNbParts)){
1306     return Standard_False;
1307   }
1308   //
1309   Handle(NCollection_IncAllocator) anIncAlloc = new NCollection_IncAllocator();
1310   NCollection_StdAllocator<ListOfInteger> anAlloc (anIncAlloc);
1311   const ListOfInteger aDummy (anIncAlloc); //empty list to be copy constructed from
1312   ArrayOfListOfInteger anArrayOfLines (aNbPnts + 1, aDummy, anAlloc);
1313
1314   NCollection_LocalArray<Standard_Integer> anArrayOfLineTypeArr (aNbPnts + 1);
1315   Standard_Integer* anArrayOfLineType = anArrayOfLineTypeArr;
1316   //
1317   nblines = 0;
1318   aTol = Precision::Confusion();
1319   //
1320   aSurf1 = theSurface1->ChangeSurface().Surface();
1321   aSurf1->Bounds(umin, umax, vmin, vmax);
1322   aPrj1.Init(aSurf1, umin, umax, vmin, vmax);
1323   //
1324   aSurf2 = theSurface2->ChangeSurface().Surface();
1325   aSurf2->Bounds(umin, umax, vmin, vmax);
1326   aPrj2.Init(aSurf2, umin, umax, vmin, vmax);
1327   //
1328   //
1329   bIsPrevPointOnBoundary=Standard_False;
1330   for(pit=1; pit<=aNbPnts; pit++) {
1331     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
1332     bIsCurrentPointOnBoundary=Standard_False;
1333     //
1334     // whether aPoint is on boundary or not
1335     //
1336     for(i=0; i<2; i++) {// exploration Surface 1,2 
1337       Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
1338       aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
1339       //
1340       for(j=0; j<2; j++) {// exploration of coordinate U,V
1341         Standard_Boolean isperiodic;
1342         //
1343         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1344         if(!isperiodic) {
1345           continue;
1346         }
1347         //
1348         Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
1349         Standard_Real aParameter, anoffset, anAdjustPar;
1350         Standard_Boolean bIsOnFirstBoundary, bIsPointOnBoundary;
1351         //
1352         aResolution = (!j) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
1353         aPeriod     = (!j) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1354         alowerboundary = (!j) ? umin : vmin;
1355         aupperboundary = (!j) ? umax : vmax;
1356         U=0.;V=0.;//?
1357         //
1358         if(!i){
1359           aPoint.ParametersOnS1(U, V);
1360         }
1361         else{
1362           aPoint.ParametersOnS2(U, V);
1363         }
1364         //
1365         aParameter = (!j) ? U : V;
1366         anoffset=0.;
1367         anAdjustPar=AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1368         //
1369         bIsOnFirstBoundary=Standard_True;
1370         //
1371         bIsPointOnBoundary=
1372           IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
1373         
1374         if(bIsPointOnBoundary) {
1375           bIsCurrentPointOnBoundary = Standard_True;
1376           break;
1377         }
1378       }// for(j=0; j<2; j++)
1379       
1380       if(bIsCurrentPointOnBoundary){
1381         break;
1382       }
1383     }// for(i=0; i<2; i++) 
1384     //
1385     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
1386
1387       if(!aListOfPointIndex.IsEmpty()) {
1388         nblines++;
1389         anArrayOfLines[nblines] = aListOfPointIndex;
1390         anArrayOfLineType[nblines] = bIsPrevPointOnBoundary;
1391         aListOfPointIndex.Clear();
1392       }
1393       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
1394     }
1395     aListOfPointIndex.Append(pit);
1396   } // for(pit=1; pit<=aNbPnts; pit++)
1397   //
1398   aNbListOfPointIndex=aListOfPointIndex.Extent();
1399   if(aNbListOfPointIndex) {
1400     nblines++;
1401     anArrayOfLines[nblines] = aListOfPointIndex;
1402     anArrayOfLineType[nblines] = bIsPrevPointOnBoundary;
1403     aListOfPointIndex.Clear();
1404   }
1405   //
1406   if(nblines <= 1){
1407     return Standard_False;
1408   }
1409   //
1410   // Correct wlines.begin
1411   Standard_Integer aLineType;
1412   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
1413   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S (new NCollection_IncAllocator());
1414   //
1415   for(i = 1; i <= nblines; i++) {
1416     aLineType=anArrayOfLineType[i];
1417     if(aLineType) {
1418       continue;
1419     }
1420     //
1421     const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1422     if(aListOfIndex.Extent() < 2) {
1423       continue;
1424     }
1425     //
1426     TColStd_ListOfInteger aListOfFLIndex;
1427     Standard_Integer aneighbourindex, aLineTypeNeib;
1428     //
1429     for(j = 0; j < 2; j++) {// neighbour line choice 
1430       aneighbourindex = (!j) ? (i-1) : (i+1);
1431       if((aneighbourindex < 1) || (aneighbourindex > nblines)){
1432         continue;
1433       }
1434       //
1435       aLineTypeNeib=anArrayOfLineType[aneighbourindex];
1436       if(!aLineTypeNeib){
1437         continue;
1438       }
1439       //
1440       const ListOfInteger& aNeighbour = anArrayOfLines[aneighbourindex];
1441       Standard_Integer anIndex = (!j) ? aNeighbour.Last() : aNeighbour.First();
1442       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
1443       // check if need use derivative.begin .end [absence]
1444       //
1445       IntSurf_PntOn2S aNewP = aPoint;
1446       Standard_Integer surfit, parit;
1447       //
1448       for(surfit = 0; surfit < 2; ++surfit) {
1449
1450         Handle(GeomAdaptor_HSurface) aGASurface = (!surfit) ? theSurface1 : theSurface2;
1451         
1452         aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
1453         Standard_Real U=0., V=0.;
1454
1455         if(!surfit) {
1456           aNewP.ParametersOnS1(U, V);
1457         }
1458         else {
1459           aNewP.ParametersOnS2(U, V);
1460         }
1461         //
1462         Standard_Integer nbboundaries = 0;
1463         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
1464         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
1465         //
1466         for(parit = 0; parit < 2; parit++) {
1467           Standard_Boolean isperiodic = (!parit) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1468
1469           Standard_Real aResolution = (!parit) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
1470           Standard_Real alowerboundary = (!parit) ? umin : vmin;
1471           Standard_Real aupperboundary = (!parit) ? umax : vmax;
1472
1473           Standard_Real aParameter = (!parit) ? U : V;
1474           Standard_Boolean bIsOnFirstBoundary = Standard_True;
1475   
1476           if(!isperiodic) {
1477             if(IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary)) {
1478               bIsUBoundary = (!parit);
1479               bIsFirstBoundary = bIsOnFirstBoundary;
1480               nbboundaries++;
1481             }
1482           }
1483           else {
1484             Standard_Real aPeriod     = (!parit) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1485             Standard_Real anoffset = 0.;
1486             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1487
1488             if(IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary)) {
1489               bIsUBoundary = (parit == 0);
1490               bIsFirstBoundary = bIsOnFirstBoundary;
1491               nbboundaries++;
1492             }
1493           }
1494         }
1495         //
1496         Standard_Boolean bComputeLineEnd = Standard_False;
1497         
1498         if(nbboundaries == 2) {
1499           bComputeLineEnd = Standard_True;
1500         }
1501         else if(nbboundaries == 1) {
1502           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1503
1504           if(isperiodic) {
1505             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
1506             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
1507             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1508             Standard_Real aParameter = (bIsUBoundary) ? U : V;
1509             Standard_Real anoffset = 0.;
1510             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1511
1512             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
1513             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
1514             anotherPar += anoffset;
1515             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
1516             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
1517             Standard_Real nU1, nV1;
1518
1519             if(surfit == 0)
1520               aNeighbourPoint.ParametersOnS1(nU1, nV1);
1521             else
1522               aNeighbourPoint.ParametersOnS2(nU1, nV1);
1523             
1524             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
1525             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
1526             bComputeLineEnd = Standard_True;
1527             Standard_Boolean bCheckAngle1 = Standard_False;
1528             Standard_Boolean bCheckAngle2 = Standard_False;
1529             gp_Vec2d aNewVec;
1530             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
1531             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
1532             //
1533             if(((adist1 - adist2) > Precision::PConfusion()) && 
1534                (adist2 < (aPeriod / 4.))) {
1535               bCheckAngle1 = Standard_True;
1536               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
1537
1538               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
1539                 aNewP.SetValue((surfit == 0), anewU, anewV);
1540                 bCheckAngle1 = Standard_False;
1541               }
1542             }
1543             else if(adist1 < (aPeriod / 4.)) {
1544               bCheckAngle2 = Standard_True;
1545               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
1546
1547               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
1548                 bCheckAngle2 = Standard_False;
1549               }
1550             }
1551             //
1552             if(bCheckAngle1 || bCheckAngle2) {
1553               // assume there are at least two points in line (see "if" above)
1554               Standard_Integer anindexother = aneighbourpointindex;
1555               
1556               while((anindexother <= aListOfIndex.Last()) && (anindexother >= aListOfIndex.First())) {
1557                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
1558                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
1559                 Standard_Real nU2, nV2;
1560                 
1561                 if(surfit == 0)
1562                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1563                 else
1564                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1565                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
1566
1567                 if(aVecOld.SquareMagnitude() <= (gp::Resolution() * gp::Resolution())) {
1568                   continue;
1569                 }
1570                 else {
1571                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
1572
1573                   if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
1574
1575                     if(bCheckAngle1) {
1576                       Standard_Real U1, U2, V1, V2;
1577                       IntSurf_PntOn2S atmppoint = aNewP;
1578                       atmppoint.SetValue((surfit == 0), anewU, anewV);
1579                       atmppoint.Parameters(U1, V1, U2, V2);
1580                       gp_Pnt P1 = theSurface1->Value(U1, V1);
1581                       gp_Pnt P2 = theSurface2->Value(U2, V2);
1582                       gp_Pnt P0 = aPoint.Value();
1583
1584                       if(P0.IsEqual(P1, aTol) &&
1585                          P0.IsEqual(P2, aTol) &&
1586                          P1.IsEqual(P2, aTol)) {
1587                         bComputeLineEnd = Standard_False;
1588                         aNewP.SetValue((surfit == 0), anewU, anewV);
1589                       }
1590                     }
1591                     
1592                     if(bCheckAngle2) {
1593                       bComputeLineEnd = Standard_False;
1594                     }
1595                   }
1596                   break;
1597                 }
1598               } // end while(anindexother...)
1599             }
1600           }
1601         }
1602         //
1603         if(bComputeLineEnd) {
1604           Standard_Integer aneighbourpointindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
1605           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
1606           Standard_Real nU1, nV1;
1607
1608           if(surfit == 0)
1609             aNeighbourPoint.ParametersOnS1(nU1, nV1);
1610           else
1611             aNeighbourPoint.ParametersOnS2(nU1, nV1);
1612           gp_Pnt2d ap1(nU1, nV1);
1613           gp_Pnt2d ap2(nU1, nV1);
1614           Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
1615
1616           while((aneighbourpointindex2 <= aListOfIndex.Last()) && (aneighbourpointindex2 >= aListOfIndex.First())) {
1617             aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
1618             const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
1619             Standard_Real nU2, nV2;
1620
1621             if(surfit == 0)
1622               aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1623             else
1624               aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1625             ap2.SetX(nU2);
1626             ap2.SetY(nV2);
1627
1628             if(ap1.SquareDistance(ap2) > (gp::Resolution() * gp::Resolution())) {
1629               break;
1630             }
1631           }       
1632           gp_Pnt2d anewpoint;
1633           Standard_Boolean found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
1634
1635           if(found) {
1636             // check point
1637             Standard_Real aCriteria =aTolSum;// BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
1638             //GeomAPI_ProjectPointOnSurf& aProjector = (surfit == 0) ? aPrj2 : aPrj1;
1639             ProjectPointOnSurf& aProjector = (surfit == 0) ? aPrj2 : aPrj1;
1640             Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
1641
1642             gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
1643             aProjector.Perform(aP3d);
1644
1645             if(aProjector.IsDone()) {
1646               if(aProjector.LowerDistance() < aCriteria) {
1647                 Standard_Real foundU = U, foundV = V;
1648                 aProjector.LowerDistanceParameters(foundU, foundV);
1649
1650                 if(surfit == 0)
1651                   aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
1652                 else
1653                   aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
1654               }
1655             }
1656           }
1657         }
1658       }
1659       aSeqOfPntOn2S->Add(aNewP);
1660       aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
1661     }
1662     anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
1663   }
1664   // Correct wlines.end
1665
1666   // Split wlines.begin
1667   for(j = 1; j <= theLConstructor.NbParts(); j++) {
1668     Standard_Real fprm = 0., lprm = 0.;
1669     theLConstructor.Part(j, fprm, lprm);
1670     Standard_Integer ifprm = (Standard_Integer)fprm;
1671     Standard_Integer ilprm = (Standard_Integer)lprm;
1672     //
1673     Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1674     //
1675     for(i = 1; i <= nblines; i++) {
1676       if(anArrayOfLineType[i] != 0) {
1677         continue;
1678       }
1679       const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1680
1681       if(aListOfIndex.Extent() < 2) {
1682         continue;
1683       }
1684       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
1685       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
1686       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
1687
1688       if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
1689         bhasfirstpoint = (i != 1);
1690       }
1691
1692       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
1693         bhaslastpoint = (i != nblines);
1694       }
1695       Standard_Boolean bIsFirstInside = ((ifprm >= aListOfIndex.First()) && (ifprm <= aListOfIndex.Last()));
1696       Standard_Boolean bIsLastInside =  ((ilprm >= aListOfIndex.First()) && (ilprm <= aListOfIndex.Last()));
1697
1698       if(!bIsFirstInside && !bIsLastInside) {
1699         if((ifprm < aListOfIndex.First()) && (ilprm > aListOfIndex.Last())) {
1700           // append whole line, and boundaries if neccesary
1701           if(bhasfirstpoint) {
1702             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
1703             aLineOn2S->Add(aP);
1704           }
1705           ListOfInteger::Iterator anIt(aListOfIndex);
1706
1707           for(; anIt.More(); anIt.Next()) {
1708             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1709             aLineOn2S->Add(aP);
1710           }
1711
1712           if(bhaslastpoint) {
1713             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
1714             aLineOn2S->Add(aP);
1715           }
1716
1717           // check end of split line (end is almost always)
1718           Standard_Integer aneighbour = i + 1;
1719           Standard_Boolean bIsEndOfLine = Standard_True;
1720
1721           if(aneighbour <= nblines) {
1722             const ListOfInteger& aListOfNeighbourIndex = anArrayOfLines[aneighbour];
1723
1724             if((anArrayOfLineType[aneighbour] != 0) &&
1725                (aListOfNeighbourIndex.IsEmpty())) {
1726               bIsEndOfLine = Standard_False;
1727             }
1728           }
1729
1730           if(bIsEndOfLine) {
1731             if(aLineOn2S->NbPoints() > 1) {
1732               Handle(IntPatch_WLine) aNewWLine = 
1733                 new IntPatch_WLine(aLineOn2S, Standard_False);
1734               theNewLines.Append(aNewWLine);
1735             }
1736             aLineOn2S = new IntSurf_LineOn2S();
1737           }
1738         }
1739         continue;
1740       }
1741       // end if(!bIsFirstInside && !bIsLastInside)
1742
1743       if(bIsFirstInside && bIsLastInside) {
1744         // append inside points between ifprm and ilprm
1745         ListOfInteger::Iterator anIt(aListOfIndex);
1746
1747         for(; anIt.More(); anIt.Next()) {
1748           if((anIt.Value() < ifprm) || (anIt.Value() > ilprm))
1749             continue;
1750           const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1751           aLineOn2S->Add(aP);
1752         }
1753       }
1754       else {
1755
1756         if(bIsFirstInside) {
1757           // append points from ifprm to last point + boundary point
1758           ListOfInteger::Iterator anIt(aListOfIndex);
1759
1760           for(; anIt.More(); anIt.Next()) {
1761             if(anIt.Value() < ifprm)
1762               continue;
1763             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1764             aLineOn2S->Add(aP);
1765           }
1766
1767           if(bhaslastpoint) {
1768             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
1769             aLineOn2S->Add(aP);
1770           }
1771           // check end of split line (end is almost always)
1772           Standard_Integer aneighbour = i + 1;
1773           Standard_Boolean bIsEndOfLine = Standard_True;
1774
1775           if(aneighbour <= nblines) {
1776             const ListOfInteger& aListOfNeighbourIndex = anArrayOfLines[aneighbour];
1777
1778             if((anArrayOfLineType[aneighbour] != 0) &&
1779                (aListOfNeighbourIndex.IsEmpty())) {
1780               bIsEndOfLine = Standard_False;
1781             }
1782           }
1783
1784           if(bIsEndOfLine) {
1785             if(aLineOn2S->NbPoints() > 1) {
1786               Handle(IntPatch_WLine) aNewWLine = 
1787                 new IntPatch_WLine(aLineOn2S, Standard_False);
1788               theNewLines.Append(aNewWLine);
1789             }
1790             aLineOn2S = new IntSurf_LineOn2S();
1791           }
1792         }
1793         // end if(bIsFirstInside)
1794
1795         if(bIsLastInside) {
1796           // append points from first boundary point to ilprm
1797           if(bhasfirstpoint) {
1798             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
1799             aLineOn2S->Add(aP);
1800           }
1801           ListOfInteger::Iterator anIt(aListOfIndex);
1802
1803           for(; anIt.More(); anIt.Next()) {
1804             if(anIt.Value() > ilprm)
1805               continue;
1806             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1807             aLineOn2S->Add(aP);
1808           }
1809         }
1810         //end if(bIsLastInside)
1811       }
1812     }
1813
1814     if(aLineOn2S->NbPoints() > 1) {
1815       Handle(IntPatch_WLine) aNewWLine = 
1816         new IntPatch_WLine(aLineOn2S, Standard_False);
1817       theNewLines.Append(aNewWLine);
1818     }
1819   }
1820   // Split wlines.end
1821   //
1822   // cda002/I3
1823   Standard_Real fprm, lprm;
1824   Standard_Integer ifprm, ilprm, aNbPoints, aIndex;
1825   //
1826   aNbParts=theLConstructor.NbParts();
1827   //
1828   for(j = 1; j <= aNbParts; j++) {
1829     theLConstructor.Part(j, fprm, lprm);
1830     ifprm=(Standard_Integer)fprm;
1831     ilprm=(Standard_Integer)lprm;
1832     //
1833     if ((ilprm-ifprm)==1) {
1834       for(i = 1; i <= nblines; i++) {
1835         aLineType=anArrayOfLineType[i];
1836         if(aLineType) {
1837           continue;
1838         }
1839         //
1840         const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1841         aNbPoints=aListOfIndex.Extent();
1842         if(aNbPoints==1) {
1843           aIndex=aListOfIndex.First();
1844           if (aIndex==ifprm || aIndex==ilprm) {
1845             Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1846             const IntSurf_PntOn2S& aP1 = theWLine->Point(ifprm);
1847             const IntSurf_PntOn2S& aP2 = theWLine->Point(ilprm);
1848             aLineOn2S->Add(aP1);
1849             aLineOn2S->Add(aP2);
1850             Handle(IntPatch_WLine) aNewWLine = 
1851               new IntPatch_WLine(aLineOn2S, Standard_False);
1852             theNewLines.Append(aNewWLine);
1853           }
1854         }
1855       }
1856     }
1857   }
1858   //
1859   return Standard_True;
1860 }
1861
1862 //=======================================================================
1863 //function : AdjustPeriodic
1864 //purpose  : 
1865 //=======================================================================
1866 Standard_Real AdjustPeriodic(const Standard_Real theParameter,
1867                              const Standard_Real parmin,
1868                              const Standard_Real parmax,
1869                              const Standard_Real thePeriod,
1870                              Standard_Real&      theOffset)
1871 {
1872   Standard_Real aresult = theParameter;
1873   theOffset = 0.;
1874   while(aresult < parmin) {
1875     aresult += thePeriod;
1876     theOffset += thePeriod;
1877   }
1878   while(aresult > parmax) {
1879     aresult -= thePeriod;
1880     theOffset -= thePeriod;
1881   }
1882   return aresult;
1883 }
1884
1885 //=======================================================================
1886 //function : IsPointOnBoundary
1887 //purpose  : 
1888 //=======================================================================
1889 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
1890                                    const Standard_Real theFirstBoundary,
1891                                    const Standard_Real theSecondBoundary,
1892                                    const Standard_Real theResolution,
1893                                    Standard_Boolean&   IsOnFirstBoundary)
1894 {
1895   IsOnFirstBoundary = Standard_True;
1896   if(fabs(theParameter - theFirstBoundary) < theResolution)
1897     return Standard_True;
1898   if(fabs(theParameter - theSecondBoundary) < theResolution)
1899   {
1900     IsOnFirstBoundary = Standard_False;
1901     return Standard_True;
1902   }
1903   return Standard_False;
1904 }
1905 //=======================================================================
1906 //function : FindPoint
1907 //purpose  : 
1908 //=======================================================================
1909 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
1910                            const gp_Pnt2d&     theLastPoint,
1911                            const Standard_Real theUmin,
1912                            const Standard_Real theUmax,
1913                            const Standard_Real theVmin,
1914                            const Standard_Real theVmax,
1915                            gp_Pnt2d&           theNewPoint)
1916 {
1917   gp_Vec2d aVec(theFirstPoint, theLastPoint);
1918   Standard_Integer i = 0, j = 0;
1919
1920   for(i = 0; i < 4; i++) {
1921     gp_Vec2d anOtherVec;
1922     gp_Vec2d anOtherVecNormal;
1923     gp_Pnt2d aprojpoint = theLastPoint;    
1924
1925     if((i % 2) == 0) {
1926       anOtherVec.SetX(0.);
1927       anOtherVec.SetY(1.);
1928       anOtherVecNormal.SetX(1.);
1929       anOtherVecNormal.SetY(0.);
1930
1931       if(i < 2)
1932         aprojpoint.SetX(theUmin);
1933       else
1934         aprojpoint.SetX(theUmax);
1935     }
1936     else {
1937       anOtherVec.SetX(1.);
1938       anOtherVec.SetY(0.);
1939       anOtherVecNormal.SetX(0.);
1940       anOtherVecNormal.SetY(1.);
1941
1942       if(i < 2)
1943         aprojpoint.SetY(theVmin);
1944       else
1945         aprojpoint.SetY(theVmax);
1946     }
1947     gp_Vec2d anormvec = aVec;
1948     anormvec.Normalize();
1949     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
1950
1951     if(fabs(adot1) < Precision::Angular())
1952       continue;
1953     Standard_Real adist = 0.;
1954
1955     if((i % 2) == 0) {
1956       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
1957     }
1958     else {
1959       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
1960     }
1961     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
1962
1963     for(j = 0; j < 2; j++) {
1964       anoffset = (j == 0) ? anoffset : -anoffset;
1965       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
1966       gp_Vec2d acurvec(theLastPoint, acurpoint);
1967
1968       //
1969       Standard_Real aDotX, anAngleX, aPC;
1970       //
1971       aDotX=aVec.Dot(acurvec);
1972       anAngleX=aVec.Angle(acurvec);
1973       aPC=Precision::PConfusion();
1974       //
1975       if(aDotX > 0. && fabs(anAngleX) < aPC) {
1976       //
1977         if((i % 2) == 0) {
1978           if((acurpoint.Y() >= theVmin) &&
1979              (acurpoint.Y() <= theVmax)) {
1980             theNewPoint = acurpoint;
1981             return Standard_True;
1982           }
1983         }
1984         else {
1985           if((acurpoint.X() >= theUmin) &&
1986              (acurpoint.X() <= theUmax)) {
1987             theNewPoint = acurpoint;
1988             return Standard_True;
1989           }
1990         }
1991       }
1992     }
1993   }
1994   return Standard_False;
1995 }
1996
1997 /*
1998 static
1999   void DumpWLine (const Handle(IntPatch_WLine)& theWLine);
2000 //=======================================================================
2001 //function : DumpWLine
2002 //purpose  : 
2003 //=======================================================================
2004   void DumpWLine (const Handle(IntPatch_WLine)& theWLine)
2005 {
2006   Standard_Integer i, nbp;
2007   Standard_Real u1,v1,u2,v2;
2008   //
2009   nbp=theWLine->NbPnts();
2010   printf(" nbp=%d\n", nbp);
2011   for(i=1;i<=nbp;i++) {
2012     const IntSurf_PntOn2S& aPoint = theWLine->Point(i);
2013     const gp_Pnt& aP=aPoint.Value();
2014     aPoint.Parameters(u1,v1,u2,v2);
2015     printf("point i%d %lf %lf %lf [ %lf %lf ] [ %lf %lf ]\n", 
2016            i, aP.X(), aP.Y(), aP.Z(), u1, v1, u2, v2);
2017   }
2018 }
2019 */