0027305: Using undefined variables, which causes devide by zero
[occt.git] / src / IntPatch / IntPatch_ImpPrmIntersection.cxx
1 // Created on: 1992-05-07
2 // Created by: Jacques GOUSSARD
3 // Copyright (c) 1992-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_HSurface.hxx>
20 #include <Adaptor3d_TopolTool.hxx>
21 #include <IntPatch_ArcFunction.hxx>
22 #include <IntPatch_ImpPrmIntersection.hxx>
23 #include <IntPatch_Line.hxx>
24 #include <IntPatch_Point.hxx>
25 #include <IntPatch_RLine.hxx>
26 #include <IntPatch_RstInt.hxx>
27 #include <IntPatch_SequenceOfLine.hxx>
28 #include <IntPatch_TheIWalking.hxx>
29 #include <IntPatch_TheIWLineOfTheIWalking.hxx>
30 #include <IntPatch_ThePathPointOfTheSOnBounds.hxx>
31 #include <IntPatch_TheSegmentOfTheSOnBounds.hxx>
32 #include <IntPatch_TheSurfFunction.hxx>
33 #include <IntPatch_WLine.hxx>
34 #include <IntSurf.hxx>
35 #include <IntSurf_InteriorPoint.hxx>
36 #include <IntSurf_LineOn2S.hxx>
37 #include <IntSurf_PathPoint.hxx>
38 #include <IntSurf_PntOn2S.hxx>
39 #include <IntSurf_SequenceOfPathPoint.hxx>
40 #include <Standard_ConstructionError.hxx>
41 #include <Standard_DomainError.hxx>
42 #include <Standard_NumericError.hxx>
43 #include <Standard_OutOfRange.hxx>
44 #include <Standard_TypeMismatch.hxx>
45 #include <StdFail_NotDone.hxx>
46 #include <TColStd_Array1OfInteger.hxx>
47
48 #ifndef OCCT_DEBUG
49 #define No_Standard_RangeError
50 #define No_Standard_OutOfRange
51 #endif
52
53 #include <math_Vector.hxx>
54 #include <math_Matrix.hxx>
55 #include <TopTrans_CurveTransition.hxx>
56 #include <TopAbs_State.hxx>
57 #include <TopAbs_Orientation.hxx>
58 #include <TColStd_Array1OfInteger.hxx>
59 #include <TColStd_Array1OfReal.hxx>
60
61 #include <IntSurf_SequenceOfInteriorPoint.hxx>
62 #include <IntSurf_QuadricTool.hxx>
63 #include <GeomAbs_SurfaceType.hxx>
64 #include <IntAna2d_AnaIntersection.hxx>
65 #include <gp_Lin2d.hxx>
66 #include <ElCLib.hxx>
67
68 #include <Bnd_Box2d.hxx>
69 #include <IntPatch_PointLine.hxx>
70
71 #include <Extrema_GenLocateExtPS.hxx>
72 #include <math_FunctionSetRoot.hxx>
73
74 static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLine,
75                                         const Standard_Boolean       IsReversed,
76                                         const IntSurf_Quadric&       theQuad,
77                                         const Handle(Adaptor3d_TopolTool)& thePDomain,
78                                         const Handle(Adaptor3d_HSurface)&  theQSurf,
79                                         const Handle(Adaptor3d_HSurface)&  theOtherSurf,
80                                         const Standard_Real                theArcTol,
81                                         const Standard_Real                theTolTang,
82                                         IntPatch_SequenceOfLine&           theLines);
83 static 
84   void ComputeTangency (const IntPatch_TheSOnBounds& solrst,
85   IntSurf_SequenceOfPathPoint& seqpdep,
86   const Handle(Adaptor3d_TopolTool)& Domain,
87   IntPatch_TheSurfFunction& Func,
88   const Handle(Adaptor3d_HSurface)& PSurf,
89   TColStd_Array1OfInteger& Destination);
90 static 
91   void Recadre(const Standard_Boolean ,
92   GeomAbs_SurfaceType typeS1,
93   GeomAbs_SurfaceType typeS2,
94   IntPatch_Point&  pt,
95   const Handle(IntPatch_TheIWLineOfTheIWalking)& iwline,
96   Standard_Integer Param,
97   Standard_Real U1,
98   Standard_Real V1,
99   Standard_Real U2,
100   Standard_Real V2);
101
102 static 
103   Standard_Boolean IsCoincide(IntPatch_TheSurfFunction& theFunc,
104                               const Handle(IntPatch_PointLine)& theLine,
105                               const Handle(Adaptor2d_HCurve2d)& theArc,
106                               const Standard_Boolean isTheSurface1Using,
107                               const Standard_Real theToler3D,
108                               const Standard_Real theToler2D,
109                               const Standard_Real thePeriod);
110
111 enum PrePoint_Type
112 {
113   PrePoint_NONE,
114   PrePoint_SEAMU,
115   PrePoint_SEAMV,
116   PrePoint_SEAMUV,
117   PrePoint_POLESEAMU,
118   PrePoint_POLE
119 };
120
121 static PrePoint_Type IsSeamOrPole(const Handle(Adaptor3d_HSurface)& theQSurf,
122                                   const Handle(IntSurf_LineOn2S)& theLine,
123                                   const Standard_Boolean IsReversed,
124                                   const Standard_Integer theRefIndex,
125                                   const Standard_Real theDeltaMax)
126 {
127   if((theRefIndex < 1) || (theRefIndex >= theLine->NbPoints()))
128     return PrePoint_NONE;
129
130   //Parameters on Quadric and on parametric for reference point
131   Standard_Real aUQRef, aVQRef, aUPRef, aVPRef;
132   Standard_Real aUQNext, aVQNext, aUPNext, aVPNext;
133
134   if(IsReversed)
135   {
136     theLine->Value(theRefIndex).Parameters  (aUPRef, aVPRef, aUQRef, aVQRef);
137     theLine->Value(theRefIndex+1).Parameters(aUPNext, aVPNext, aUQNext, aVQNext);
138   }
139   else
140   {
141     theLine->Value(theRefIndex).Parameters  (aUQRef, aVQRef, aUPRef, aVPRef);
142     theLine->Value(theRefIndex+1).Parameters(aUQNext, aVQNext, aUPNext, aVPNext);
143   }
144
145   const GeomAbs_SurfaceType aType = theQSurf->GetType();
146
147   const Standard_Real aDeltaU = Abs(aUQRef - aUQNext);
148
149   if((aType != GeomAbs_Torus) && (aDeltaU < theDeltaMax))
150     return PrePoint_NONE;
151
152   switch(aType)
153   {
154   case GeomAbs_Cylinder:
155     return PrePoint_SEAMU;
156
157   case GeomAbs_Torus:
158     {
159       const Standard_Real aDeltaV = Abs(aVQRef - aVQNext);
160
161       if((aDeltaU >= theDeltaMax) && (aDeltaV >= theDeltaMax))
162         return PrePoint_SEAMUV;
163
164       if(aDeltaU >= theDeltaMax)
165         return PrePoint_SEAMU;
166
167       if(aDeltaV >= theDeltaMax)
168         return PrePoint_SEAMV;
169     }
170
171     break;
172   case GeomAbs_Sphere:
173   case GeomAbs_Cone:
174     return PrePoint_POLESEAMU;
175   default:
176     break;
177   }
178
179   return PrePoint_NONE;
180 }
181
182 // The function for searching intersection point, which 
183 // lies in the seam-edge of the quadric definetely.
184 class FuncPreciseSeam: public math_FunctionSetWithDerivatives
185 {
186 public:
187   FuncPreciseSeam(const Handle(Adaptor3d_HSurface)& theQSurf, const Handle(Adaptor3d_HSurface)& thePSurf, const Standard_Boolean isTheUSeam): myQSurf(theQSurf), myPSurf(thePSurf), myIsUSeam(isTheUSeam) {};
188   
189   Standard_EXPORT virtual Standard_Integer NbVariables() const
190   {
191     return 3;
192   };
193
194   Standard_EXPORT virtual Standard_Integer NbEquations() const
195   {
196     return 3;
197   }
198
199   Standard_EXPORT virtual Standard_Boolean Value (const math_Vector& theX, math_Vector& theF)
200   {
201     try
202     {
203       const Standard_Integer anIndX = theX.Lower(), anIndF = theF.Lower();
204       const gp_Pnt aP1(myPSurf->Value(theX(anIndX), theX(anIndX+1)));
205       const gp_Pnt aP2(myIsUSeam? myQSurf->Value(0.0, theX(anIndX+2)) : myQSurf->Value(theX(anIndX+2), 0.0));
206
207       (aP1.XYZ()-aP2.XYZ()).Coord(theF(anIndF), theF(anIndF+1), theF(anIndF+2));
208     }
209     catch(Standard_Failure)
210     {
211       return Standard_False;
212     }
213
214     return Standard_True;
215   };
216
217   Standard_EXPORT virtual Standard_Boolean Derivatives (const math_Vector& theX, math_Matrix& theD)
218   {
219     try
220     {
221       const Standard_Integer anIndX = theX.Lower(), anIndRD = theD.LowerRow(), anIndCD = theD.LowerCol();
222       gp_Pnt aPt;
223       gp_Vec aD1u, aD1v, aD2u, aD2v;
224       myPSurf->D1(theX(anIndX), theX(anIndX+1), aPt, aD1u, aD1v);
225       if(myIsUSeam)
226         myQSurf->D1(0.0, theX(anIndX+2), aPt, aD2u, aD2v);
227       else
228         myQSurf->D1(theX(anIndX+2), 0.0, aPt, aD2u, aD2v);
229
230       // d/dX1
231       aD1u.Coord(theD(anIndRD, anIndCD), theD(anIndRD+1, anIndCD), theD(anIndRD+2, anIndCD));
232
233       // d/dX1
234       aD1v.Coord(theD(anIndRD, anIndCD+1), theD(anIndRD+1, anIndCD+1), theD(anIndRD+2, anIndCD+1));
235
236       // d/dX3
237       if(myIsUSeam)
238         aD2v.Reversed().Coord(theD(anIndRD, anIndCD+2), theD(anIndRD+1, anIndCD+2), theD(anIndRD+2, anIndCD+2));
239       else
240         aD2u.Reversed().Coord(theD(anIndRD, anIndCD+2), theD(anIndRD+1, anIndCD+2), theD(anIndRD+2, anIndCD+2));
241     }
242     catch(Standard_Failure)
243     {
244       return Standard_False;
245     }
246
247     return Standard_True;
248   };
249
250   Standard_EXPORT virtual Standard_Boolean Values (const math_Vector& theX, math_Vector& theF, math_Matrix& theD)
251   {
252     try
253     {
254       const Standard_Integer anIndX = theX.Lower(), anIndF = theF.Lower(), anIndRD = theD.LowerRow(), anIndCD = theD.LowerCol();
255       gp_Pnt aP1, aP2;
256       gp_Vec aD1u, aD1v, aD2u, aD2v;
257       myPSurf->D1(theX(anIndX), theX(anIndX+1), aP1, aD1u, aD1v);
258       if(myIsUSeam)
259         myQSurf->D1(0.0, theX(anIndX+2), aP2, aD2u, aD2v);
260       else
261         myQSurf->D1(theX(anIndX+2), 0.0, aP2, aD2u, aD2v);
262
263       //Value
264       (aP1.XYZ()-aP2.XYZ()).Coord(theF(anIndF), theF(anIndF+1), theF(anIndF+2));
265
266       // d/dX1
267       aD1u.Coord(theD(anIndRD, anIndCD), theD(anIndRD+1, anIndCD), theD(anIndRD+2, anIndCD));
268
269       // d/dX1
270       aD1v.Coord(theD(anIndRD, anIndCD+1), theD(anIndRD+1, anIndCD+1), theD(anIndRD+2, anIndCD+1));
271
272       // d/dX3
273       if(myIsUSeam)
274         aD2v.Reversed().Coord(theD(anIndRD, anIndCD+2), theD(anIndRD+1, anIndCD+2), theD(anIndRD+2, anIndCD+2));
275       else
276         aD2u.Reversed().Coord(theD(anIndRD, anIndCD+2), theD(anIndRD+1, anIndCD+2), theD(anIndRD+2, anIndCD+2));
277     }
278     catch(Standard_Failure)
279     {
280       return Standard_False;
281     }
282
283     return Standard_True;
284   }
285
286 protected:
287   FuncPreciseSeam operator=(FuncPreciseSeam&);
288
289 private:
290   const Handle(Adaptor3d_HSurface)& myQSurf;
291   const Handle(Adaptor3d_HSurface)& myPSurf;
292   const Standard_Boolean myIsUSeam;
293 };
294
295 //=======================================================================
296 //function : IntPatch_ImpPrmIntersection
297 //purpose  : 
298 //=======================================================================
299 IntPatch_ImpPrmIntersection::IntPatch_ImpPrmIntersection ()
300   : done(Standard_False),
301   empt(Standard_False),
302   myIsStartPnt(Standard_False),
303   myUStart(0.0),
304   myVStart(0.0)
305 { }
306
307
308 //=======================================================================
309 //function : IntPatch_ImpPrmIntersection
310 //purpose  : 
311 //=======================================================================
312
313 IntPatch_ImpPrmIntersection::IntPatch_ImpPrmIntersection
314   (const Handle(Adaptor3d_HSurface)&    Surf1,
315   const Handle(Adaptor3d_TopolTool)&   D1,
316   const Handle(Adaptor3d_HSurface)&    Surf2,
317   const Handle(Adaptor3d_TopolTool)&   D2,
318   const Standard_Real    TolArc,
319   const Standard_Real    TolTang,
320   const Standard_Real    Fleche,
321   const Standard_Real    Pas)
322   : done(Standard_False),
323   empt(Standard_False),
324   myIsStartPnt(Standard_False),
325   myUStart(0.0),
326   myVStart(0.0)
327 {
328   Perform(Surf1,D1,Surf2,D2,TolArc,TolTang,Fleche,Pas);
329 }
330
331
332 //=======================================================================
333 //function : SetStartPoint
334 //purpose  : 
335 //=======================================================================
336
337 void IntPatch_ImpPrmIntersection::SetStartPoint(const Standard_Real U,
338   const Standard_Real V)
339 {
340   myIsStartPnt = Standard_True;
341   myUStart = U; myVStart = V;
342 }
343
344 //=======================================================================
345 //function : ComputeTangency
346 //purpose  : 
347 //=======================================================================
348 void ComputeTangency (const IntPatch_TheSOnBounds& solrst,
349   IntSurf_SequenceOfPathPoint& seqpdep,
350   const Handle(Adaptor3d_TopolTool)& Domain,
351   IntPatch_TheSurfFunction& Func,
352   const Handle(Adaptor3d_HSurface)& PSurf,
353   TColStd_Array1OfInteger& Destination)
354 {
355   Standard_Integer i,k, NbPoints, seqlength;
356   Standard_Real theparam,test;
357   Standard_Boolean fairpt, ispassing;
358   TopAbs_Orientation arcorien,vtxorien;
359   Handle(Adaptor2d_HCurve2d) thearc;
360   Handle(Adaptor3d_HVertex) vtx,vtxbis;
361   //Standard_Boolean ispassing;
362   IntPatch_ThePathPointOfTheSOnBounds PStart;
363   IntSurf_PathPoint PPoint;
364   gp_Vec vectg;
365   gp_Dir2d dirtg;
366   gp_Pnt ptbid;
367   gp_Vec d1u,d1v,v1,v2;
368   gp_Pnt2d p2d;
369   gp_Vec2d d2d;
370   //
371   double aX[2], aF[1], aD[1][2];
372   math_Vector X(aX, 1, 2);
373   math_Vector F(aF, 1, 1);
374   math_Matrix D(aD, 1, 1, 1, 2); 
375   //
376   seqlength = 0;
377   NbPoints = solrst.NbPoints();
378   for (i=1; i<= NbPoints; i++) {
379     if (Destination(i) == 0) {
380       PStart = solrst.Point(i);
381       thearc   = PStart.Arc();
382       theparam = PStart.Parameter();
383       arcorien = Domain->Orientation(thearc);
384       ispassing = (arcorien == TopAbs_INTERNAL || 
385         arcorien == TopAbs_EXTERNAL);
386
387       thearc->D0(theparam,p2d);
388       X(1) = p2d.X(); 
389       X(2) = p2d.Y();
390       PPoint.SetValue(PStart.Value(),X(1),X(2));
391
392       Func.Values(X,F,D);
393       if (Func.IsTangent()) {
394         PPoint.SetTangency(Standard_True);
395         Destination(i) = seqlength+1;
396         if (!PStart.IsNew()) {
397           vtx = PStart.Vertex();
398           for (k=i+1; k<=NbPoints; k++) {
399             if (Destination(k) ==0) {
400               PStart = solrst.Point(k);
401               if (!PStart.IsNew()) {
402                 vtxbis = PStart.Vertex();
403                 if (Domain->Identical(vtx,vtxbis)) {
404                   thearc   = PStart.Arc();
405                   theparam = PStart.Parameter();
406                   arcorien = Domain->Orientation(thearc);
407                   ispassing = ispassing && (arcorien == TopAbs_INTERNAL ||
408                     arcorien == TopAbs_EXTERNAL);
409
410                   thearc->D0(theparam,p2d);
411                   PPoint.AddUV(p2d.X(),p2d.Y());
412                   Destination(k) = seqlength+1;
413                 }
414               }
415             }
416           }
417         }
418         PPoint.SetPassing(ispassing);
419         seqpdep.Append(PPoint);
420         seqlength++;
421       }
422       else { // on a un point de depart potentiel
423
424         vectg = Func.Direction3d();
425         dirtg = Func.Direction2d();
426
427         PSurf->D1(X(1),X(2),ptbid,d1u,d1v);
428         thearc->D1(theparam,p2d,d2d);
429         v2.SetLinearForm(d2d.X(),d1u,d2d.Y(),d1v);
430         v1 = d1u.Crossed(d1v);
431
432         test = vectg.Dot(v1.Crossed(v2));
433         if (PStart.IsNew()) {
434           if ((test < 0. && arcorien == TopAbs_FORWARD) ||
435             (test > 0. && arcorien == TopAbs_REVERSED)) {
436               vectg.Reverse();
437               dirtg.Reverse();
438           }
439           PPoint.SetDirections(vectg,dirtg);
440           PPoint.SetPassing(ispassing);
441           Destination(i) = seqlength+1;
442           seqpdep.Append(PPoint);
443           seqlength++;
444         }
445         else { // traiter la transition complexe
446           gp_Dir bidnorm(1.,1.,1.);
447           Standard_Real tole = 1.e-8;
448           TopAbs_Orientation LocTrans;
449           TopTrans_CurveTransition comptrans;
450           comptrans.Reset(vectg,bidnorm,0.);
451           if (arcorien == TopAbs_FORWARD || 
452             arcorien == TopAbs_REVERSED) {
453               // pour essai
454
455               vtx = PStart.Vertex();
456               vtxorien = Domain->Orientation(vtx);
457               if (Abs(test) <= tole) {
458                 LocTrans = TopAbs_EXTERNAL; // et pourquoi pas INTERNAL
459               }
460               else {
461                 if (((test > 0.)&& arcorien == TopAbs_FORWARD) ||
462                   ((test < 0.)&& arcorien == TopAbs_REVERSED)){
463                     LocTrans = TopAbs_FORWARD;
464                 }
465                 else {
466                   LocTrans = TopAbs_REVERSED;
467                 }
468                 if (arcorien == TopAbs_REVERSED) {v2.Reverse();}
469               }
470
471               comptrans.Compare(tole,v2,bidnorm,0.,LocTrans,vtxorien);
472           }
473           Destination(i) = seqlength+1;
474           for (k= i+1; k<=NbPoints; k++) {
475             if (Destination(k) == 0) {
476               PStart = solrst.Point(k);
477               if (!PStart.IsNew()) {
478                 vtxbis = PStart.Vertex();
479                 if (Domain->Identical(vtx,vtxbis)) {
480                   thearc   = PStart.Arc();
481                   theparam = PStart.Parameter();
482                   arcorien = Domain->Orientation(thearc);
483
484                   PPoint.AddUV(X(1),X(2));
485
486                   thearc->D1(theparam,p2d,d2d);
487                   PPoint.AddUV(p2d.X(),p2d.Y());
488
489                   if (arcorien == TopAbs_FORWARD || 
490                     arcorien == TopAbs_REVERSED) {
491                       ispassing = Standard_False;
492                       v2.SetLinearForm(d2d.X(),d1u,d2d.Y(),d1v);
493
494                       test = vectg.Dot(v1.Crossed(v2));
495                       vtxorien = Domain->Orientation(PStart.Vertex());
496                       if (Abs(test) <= tole) {
497                         LocTrans = TopAbs_EXTERNAL; // et pourquoi pas INTERNAL
498                       }
499                       else {
500                         if (((test > 0.)&& arcorien == TopAbs_FORWARD) ||
501                           ((test < 0.)&& arcorien == TopAbs_REVERSED)){
502                             LocTrans = TopAbs_FORWARD;
503                         }
504                         else {
505                           LocTrans = TopAbs_REVERSED;
506                         }
507                         if (arcorien == TopAbs_REVERSED) {v2.Reverse();}
508                       }
509
510                       comptrans.Compare(tole,v2,bidnorm,0.,LocTrans,vtxorien);
511                   }
512                   Destination(k) = seqlength+1;
513                 }
514               }
515             }
516           }
517           fairpt = Standard_True;
518           if (!ispassing) {
519             TopAbs_State Before = comptrans.StateBefore();
520             TopAbs_State After  = comptrans.StateAfter();
521             if ((Before == TopAbs_UNKNOWN)||(After == TopAbs_UNKNOWN)) {
522               fairpt = Standard_False;
523             }
524             else if (Before == TopAbs_IN) {
525               if (After == TopAbs_IN) {
526                 ispassing = Standard_True;
527               }
528               else {
529                 vectg.Reverse();
530                 dirtg.Reverse();
531               }
532             }
533             else {
534               if (After !=TopAbs_IN) {
535                 fairpt = Standard_False;
536               }
537             }
538           }
539           if (fairpt) {
540             PPoint.SetDirections(vectg,dirtg);
541             PPoint.SetPassing(ispassing);
542             seqpdep.Append(PPoint);
543             seqlength++;
544           }
545           else { // il faut remettre en "ordre" si on ne garde pas le point.
546             for (k=i; k <=NbPoints ; k++) {
547               if (Destination(k)==seqlength + 1) {
548                 Destination(k) = -Destination(k);
549               }
550             }
551           }
552         }
553       }
554     }
555   }
556 }
557 //=======================================================================
558 //function : Recadre
559 //purpose  : 
560 //=======================================================================
561 void Recadre(const Standard_Boolean ,
562   GeomAbs_SurfaceType typeS1,
563   GeomAbs_SurfaceType typeS2,
564   IntPatch_Point&  pt,
565   const Handle(IntPatch_TheIWLineOfTheIWalking)& iwline,
566   Standard_Integer Param,
567   Standard_Real U1,
568   Standard_Real V1,
569   Standard_Real U2,
570   Standard_Real V2)
571 {
572   Standard_Real U1p,V1p,U2p,V2p;
573   iwline->Line()->Value(Param).Parameters(U1p,V1p,U2p,V2p);
574   switch(typeS1)
575   {
576   case GeomAbs_Torus:
577     while(V1<(V1p-1.5*M_PI)) V1+=M_PI+M_PI;
578     while(V1>(V1p+1.5*M_PI)) V1-=M_PI+M_PI;
579   case GeomAbs_Cylinder:
580   case GeomAbs_Cone:
581   case GeomAbs_Sphere:
582     while(U1<(U1p-1.5*M_PI)) U1+=M_PI+M_PI;
583     while(U1>(U1p+1.5*M_PI)) U1-=M_PI+M_PI;
584   default:
585     break;
586   }
587   switch(typeS2)
588   { 
589   case GeomAbs_Torus:
590     while(V2<(V2p-1.5*M_PI)) V2+=M_PI+M_PI;
591     while(V2>(V2p+1.5*M_PI)) V2-=M_PI+M_PI;
592   case GeomAbs_Cylinder:
593   case GeomAbs_Cone:
594   case GeomAbs_Sphere:
595     while(U2<(U2p-1.5*M_PI)) U2+=M_PI+M_PI;
596     while(U2>(U2p+1.5*M_PI)) U2-=M_PI+M_PI;
597   default:
598     break;
599   }
600   pt.SetParameters(U1,V1,U2,V2);
601 }
602
603 //=======================================================================
604 //function : Perform
605 //purpose  : 
606 //=======================================================================
607 void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Surf1,
608   const Handle(Adaptor3d_TopolTool)& D1,
609   const Handle(Adaptor3d_HSurface)& Surf2,
610   const Handle(Adaptor3d_TopolTool)& D2,
611   const Standard_Real TolArc,
612   const Standard_Real TolTang,
613   const Standard_Real Fleche,
614   const Standard_Real Pas)
615 {
616   Standard_Boolean reversed, procf, procl, dofirst, dolast;
617   Standard_Integer indfirst = 0, indlast = 0, ind2, NbSegm;
618   Standard_Integer NbPointIns, NbPointRst, Nblines, Nbpts, NbPointDep;
619   Standard_Real U1,V1,U2,V2,paramf,paraml,currentparam;
620
621   IntPatch_TheSegmentOfTheSOnBounds thesegm;
622   IntSurf_PathPoint PPoint;
623
624   Handle(IntPatch_RLine) rline;
625   Handle(IntPatch_WLine) wline;
626   IntPatch_ThePathPointOfTheSOnBounds PStart,PStartf,PStartl;
627   IntPatch_Point ptdeb,ptfin,ptbis;
628
629   IntPatch_IType typ;
630   IntSurf_Transition TLine,TArc;
631   IntSurf_TypeTrans trans1,trans2;
632   gp_Pnt valpt,ptbid;
633   gp_Vec tgline,tgrst,norm1,norm2,d1u,d1v;
634   gp_Dir DirNormale;
635   gp_Vec VecNormale;
636
637   gp_Pnt2d p2d;
638   gp_Vec2d d2d;
639
640   Handle(Adaptor2d_HCurve2d) currentarc;
641   GeomAbs_SurfaceType typeS1, typeS2;
642   IntSurf_Quadric Quad;
643   IntPatch_TheSurfFunction Func;
644   IntPatch_ArcFunction AFunc;
645   //
646   typeS1 = Surf1->GetType();
647   typeS2 = Surf2->GetType();
648
649   paramf =0.;
650   paraml =0.;
651   trans1 = IntSurf_Undecided;
652   trans2 = IntSurf_Undecided;
653   //
654   done = Standard_False;
655   empt = Standard_True;
656   slin.Clear();
657   spnt.Clear();
658   //
659   reversed = Standard_False;
660   switch (typeS1)
661   {
662   case GeomAbs_Plane:
663     Quad.SetValue(Surf1->Plane());
664     break;
665
666   case GeomAbs_Cylinder:
667     Quad.SetValue(Surf1->Cylinder());
668     break;
669
670   case GeomAbs_Sphere:
671     Quad.SetValue(Surf1->Sphere());
672     break;
673
674   case GeomAbs_Cone:
675     Quad.SetValue(Surf1->Cone());
676     break;
677
678   default:
679     {
680       reversed = Standard_True;
681       switch (typeS2)
682       {
683       case GeomAbs_Plane:
684         Quad.SetValue(Surf2->Plane());
685         break;
686
687       case GeomAbs_Cylinder:
688         Quad.SetValue(Surf2->Cylinder());
689         break;
690
691       case GeomAbs_Sphere:
692         Quad.SetValue(Surf2->Sphere());
693         break;
694
695       case GeomAbs_Cone:
696         Quad.SetValue(Surf2->Cone());
697         break;
698       default:
699         {
700           Standard_ConstructionError::Raise();
701           break;
702         }
703       } 
704     }
705     break;
706   }
707   //
708   Func.SetImplicitSurface(Quad);
709   Func.Set(IntSurf_QuadricTool::Tolerance(Quad));
710   AFunc.SetQuadric(Quad);
711   //
712   if (!reversed) {
713     Func.Set(Surf2);
714     AFunc.Set(Surf2);
715   }
716   else {
717     Func.Set(Surf1);
718     AFunc.Set(Surf1);
719   }
720   //
721   if (!reversed) {
722     solrst.Perform(AFunc,D2,TolArc,TolTang);
723   }
724   else {
725     solrst.Perform(AFunc,D1,TolArc,TolTang);
726   }
727   if (!solrst.IsDone()) {
728     return;
729   }
730   //
731   IntSurf_SequenceOfPathPoint seqpdep;
732   IntSurf_SequenceOfInteriorPoint seqpins;
733   //
734   NbPointRst = solrst.NbPoints();
735   TColStd_Array1OfInteger Destination(1,NbPointRst+1); Destination.Init(0);
736   if (NbPointRst) {
737     if (!reversed) {
738       ComputeTangency(solrst,seqpdep,D2,Func,Surf2,Destination);
739     }
740     else {
741       ComputeTangency(solrst,seqpdep,D1,Func,Surf1,Destination);
742     }
743   }
744   //
745   Standard_Boolean SearchIns = Standard_True;
746   if(Quad.TypeQuadric() == GeomAbs_Plane && solrst.NbSegments() > 0)
747   {
748     //For such kind of cases it is possible that whole surface is on one side of plane,
749     //plane only touches surface and does not cross it,
750     //so no inner points exist.
751     SearchIns = Standard_False;
752     Handle(Adaptor3d_TopolTool) T;
753     if(reversed)
754     {
755       T = D1;
756     }
757     else
758     {
759       T = D2;
760     }
761     Standard_Integer aNbSamples = 0;
762     aNbSamples = T->NbSamples();
763     gp_Pnt2d s2d;
764     gp_Pnt s3d;
765     Standard_Real aValf[1], aUVap[2];
766     math_Vector Valf(aValf,1,1), UVap(aUVap,1,2);
767     T->SamplePoint(1,s2d, s3d);
768     UVap(1)=s2d.X(); 
769     UVap(2)=s2d.Y();
770     Func.Value(UVap,Valf);
771     Standard_Real rvalf = Sign(1.,Valf(1));
772     for(Standard_Integer i = 2; i <= aNbSamples; ++i)
773     {
774       T->SamplePoint(i,s2d, s3d);
775       UVap(1)=s2d.X(); 
776       UVap(2)=s2d.Y();
777       Func.Value(UVap,Valf);
778       if(rvalf * Valf(1) < 0.)
779       {
780         SearchIns = Standard_True;
781         break;
782       }   
783     }
784   }
785   // Recherche des points interieurs
786   NbPointIns = 0;
787   if(SearchIns) {
788     if (!reversed) {
789       if (myIsStartPnt)
790         solins.Perform(Func,Surf2,myUStart,myVStart);
791       else
792         solins.Perform(Func,Surf2,D2,TolTang);
793     }
794     else {
795       if (myIsStartPnt)
796         solins.Perform(Func,Surf1,myUStart,myVStart);
797       else
798         solins.Perform(Func,Surf1,D1,TolTang);
799     }
800     NbPointIns = solins.NbPoints();
801     for (Standard_Integer i=1; i <= NbPointIns; i++) {
802       seqpins.Append(solins.Value(i));
803     }
804   }
805   //
806   NbPointDep=seqpdep.Length();
807   //
808   if (NbPointDep || NbPointIns) {
809     IntPatch_TheIWalking iwalk(TolTang,Fleche,Pas);
810     if (!reversed) {
811       iwalk.Perform(seqpdep,seqpins,Func,Surf2);
812     }
813     else {
814       iwalk.Perform(seqpdep,seqpins,Func,Surf1,Standard_True);
815     }
816     if(!iwalk.IsDone()) {
817       return;
818     }
819
820     Standard_Real Vmin, Vmax, TolV = 1.e-14;
821     if (!reversed) { //Surf1 is quadric
822       Vmin = Surf1->FirstVParameter();
823       Vmax = Surf1->LastVParameter();
824     }
825     else { //Surf2 is quadric
826       Vmin = Surf2->FirstVParameter();
827       Vmax = Surf2->LastVParameter();
828     }
829     //
830     Nblines = iwalk.NbLines();
831     for (Standard_Integer j=1; j<=Nblines; j++) {
832       const Handle(IntPatch_TheIWLineOfTheIWalking)&  iwline  = iwalk.Value(j);
833       const Handle(IntSurf_LineOn2S)&                 thelin  = iwline->Line();
834
835       Nbpts = thelin->NbPoints();
836       if(Nbpts>=2) { 
837         Standard_Integer k = 0;
838         tgline = iwline->TangentVector(k);      
839         if(k>=1 && k<=Nbpts) { } else { k=Nbpts>>1; } 
840         valpt = thelin->Value(k).Value();       
841
842         if (!reversed) {
843           thelin->Value(k).ParametersOnS2(U2,V2);
844           norm1 = Quad.Normale(valpt);
845           Surf2->D1(U2,V2,ptbid,d1u,d1v);
846           norm2 = d1u.Crossed(d1v);
847         }
848         else {
849           thelin->Value(k).ParametersOnS1(U2,V2);
850           norm2 = Quad.Normale(valpt);
851           Surf1->D1(U2,V2,ptbid,d1u,d1v);
852           norm1 = d1u.Crossed(d1v);
853         }
854         if (tgline.DotCross(norm2,norm1) > 0.) {
855           trans1 = IntSurf_Out;
856           trans2 = IntSurf_In;
857         }
858         else {
859           trans1 = IntSurf_In;
860           trans2 = IntSurf_Out;
861         }
862
863         //
864         Standard_Real AnU1,AnU2,AnV2;
865
866         GeomAbs_SurfaceType typQuad = Quad.TypeQuadric();
867         Standard_Boolean arecadr=Standard_False;
868         valpt = thelin->Value(1).Value();
869         Quad.Parameters(valpt,AnU1,V1);
870
871         if((V1 < Vmin) && (Vmin-V1 < TolV)) V1 = Vmin;
872         if((V1 > Vmax) && (V1-Vmax < TolV)) V1 = Vmax;
873
874         if(reversed) { 
875           thelin->SetUV(1,Standard_False,AnU1,V1); //-- on va lire u2,v2
876           thelin->Value(1).ParametersOnS1(AnU2,AnV2);
877         }
878         else { 
879           thelin->SetUV(1,Standard_True,AnU1,V1);  //-- on va lire u1,v1 
880           thelin->Value(1).ParametersOnS2(AnU2,AnV2);
881         }
882
883         if(typQuad==GeomAbs_Cylinder || 
884           typQuad==GeomAbs_Cone || 
885           typQuad==GeomAbs_Sphere) { 
886             arecadr=Standard_True; 
887         } 
888         //
889         for (k=2; k<=Nbpts; ++k) {
890           valpt = thelin->Value(k).Value();
891           Quad.Parameters(valpt,U1,V1);
892           //
893           if((V1 < Vmin) && (Vmin-V1 < TolV)) {
894             V1 = Vmin;
895           }
896           if((V1 > Vmax) && (V1-Vmax < TolV)) {
897             V1 = Vmax;
898           }
899           //
900           if(arecadr) {
901             //modified by NIZNHY-PKV Fri Mar 28 15:06:01 2008f
902             Standard_Real aCf, aTwoPI;
903             //
904             aCf=0.;
905             aTwoPI=M_PI+M_PI;
906             if ((U1-AnU1) >  1.5*M_PI) { 
907               while ((U1-AnU1) > (1.5*M_PI+aCf*aTwoPI)) {
908                 aCf=aCf+1.;
909               }
910               U1=U1-aCf*aTwoPI;
911             } 
912             //
913             else {
914               while ((U1-AnU1) < (-1.5*M_PI-aCf*aTwoPI)) {
915                 aCf=aCf+1.;
916               }
917               U1=U1+aCf*aTwoPI;
918             }
919             // was:
920             //if ((U1-AnU1) >  1.5*M_PI) { 
921             //  U1-=M_PI+M_PI;
922             //}
923             //else if ((U1-AnU1) < -1.5*M_PI) { 
924             //  U1+=M_PI+M_PI; 
925             //}
926             //modified by NIZNHY-PKV Fri Mar 28 15:06:11 2008t
927           }
928           //
929           if(reversed) { 
930             thelin->SetUV(k,Standard_False,U1,V1);
931
932             thelin->Value(k).ParametersOnS1(U2,V2);
933             switch(typeS1) { 
934             case GeomAbs_Cylinder:
935             case GeomAbs_Cone:
936             case GeomAbs_Sphere:
937             case GeomAbs_Torus:
938               while(U2<(AnU2-1.5*M_PI)) U2+=M_PI+M_PI;
939               while(U2>(AnU2+1.5*M_PI)) U2-=M_PI+M_PI;
940               break;
941             default: 
942               break;
943             }
944             if(typeS2==GeomAbs_Torus) { 
945               while(V2<(AnV2-1.5*M_PI)) V2+=M_PI+M_PI;
946               while(V2>(AnV2+1.5*M_PI)) V2-=M_PI+M_PI;
947             }
948             thelin->SetUV(k,Standard_True,U2,V2);
949           }
950           else { 
951             thelin->SetUV(k,Standard_True,U1,V1);
952
953             thelin->Value(k).ParametersOnS2(U2,V2);
954             switch(typeS2) { 
955             case GeomAbs_Cylinder:
956             case GeomAbs_Cone:
957             case GeomAbs_Sphere:
958             case GeomAbs_Torus:
959               while(U2<(AnU2-1.5*M_PI)) U2+=M_PI+M_PI;
960               while(U2>(AnU2+1.5*M_PI)) U2-=M_PI+M_PI;
961               break;
962             default: 
963               break;
964             }
965             if(typeS2==GeomAbs_Torus) { 
966               while(V2<(AnV2-1.5*M_PI)) V2+=M_PI+M_PI;
967               while(V2>(AnV2+1.5*M_PI)) V2-=M_PI+M_PI;
968             }
969             thelin->SetUV(k,Standard_False,U2,V2);
970
971           }
972
973           AnU1=U1;
974           AnU2=U2;
975           AnV2=V2;
976         }
977         // <-A
978         wline = new IntPatch_WLine(thelin,Standard_False,trans1,trans2);
979
980 #ifdef INTPATCH_IMPPRMINTERSECTION_DEBUG
981         wline->Dump(0);
982 #endif
983
984         if (   iwline->HasFirstPoint() 
985           && iwline->IsTangentAtBegining() == Standard_False) 
986         {
987           indfirst = iwline->FirstPointIndex();
988           PPoint = seqpdep(indfirst);
989           tgline = PPoint.Direction3d();
990           Standard_Integer themult = PPoint.Multiplicity();
991           for (Standard_Integer i=NbPointRst; i>=1; i--) {
992             if (Destination(i) == indfirst) {
993               if (!reversed) { //-- typeS1 = Pln || Cyl || Sph || Cone
994                 Quad.Parameters(PPoint.Value(),U1,V1);
995
996                 if((V1 < Vmin) && (Vmin-V1 < TolV)) V1 = Vmin;
997                 if((V1 > Vmax) && (V1-Vmax < TolV)) V1 = Vmax;
998
999                 PPoint.Parameters(themult,U2,V2);
1000                 Surf2->D1(U2,V2,ptbid,d1u,d1v); //-- @@@@
1001               }
1002               else {  //-- typeS1 != Pln && Cyl && Sph && Cone
1003                 Quad.Parameters(PPoint.Value(),U2,V2);
1004
1005                 if((V2 < Vmin) && (Vmin-V2 < TolV)) V2 = Vmin;
1006                 if((V2 > Vmax) && (V2-Vmax < TolV)) V2 = Vmax;
1007
1008                 PPoint.Parameters(themult,U1,V1);
1009                 Surf1->D1(U1,V1,ptbid,d1u,d1v); //-- @@@@
1010               }
1011
1012               VecNormale = d1u.Crossed(d1v);                      
1013               //-- Modif du 27 Septembre 94 (Recadrage des pts U,V) 
1014               ptdeb.SetValue(PPoint.Value(),TolArc,Standard_False);
1015               ptdeb.SetParameters(U1,V1,U2,V2);
1016               ptdeb.SetParameter(1.);
1017
1018               Recadre(reversed,typeS1,typeS2,ptdeb,iwline,1,U1,V1,U2,V2);
1019
1020               currentarc = solrst.Point(i).Arc();
1021               currentparam = solrst.Point(i).Parameter();
1022               currentarc->D1(currentparam,p2d,d2d);
1023               tgrst.SetLinearForm(d2d.X(),d1u,d2d.Y(),d1v);
1024
1025               Standard_Real squaremagnitudeVecNormale = VecNormale.SquareMagnitude();
1026               if(squaremagnitudeVecNormale > 1e-13) { 
1027                 DirNormale=VecNormale;
1028                 IntSurf::MakeTransition(tgline,tgrst,DirNormale,TLine,TArc);
1029               }
1030               else { 
1031                 TLine.SetValue(Standard_True,IntSurf_Undecided);
1032                 TArc.SetValue(Standard_True,IntSurf_Undecided);
1033               }
1034
1035               ptdeb.SetArc(reversed,currentarc,currentparam,TLine,TArc);
1036               if (!solrst.Point(i).IsNew()) {
1037                 ptdeb.SetVertex(reversed,solrst.Point(i).Vertex());
1038               }
1039               wline->AddVertex(ptdeb);
1040               if (themult == 0) {
1041                 wline->SetFirstPoint(wline->NbVertex());
1042               }
1043
1044               themult--;
1045             }
1046           }
1047         }
1048         else if (iwline->IsTangentAtBegining()) 
1049         {
1050           gp_Pnt psol = thelin->Value(1).Value();
1051           thelin->Value(1).ParametersOnS1(U1,V1);
1052           thelin->Value(1).ParametersOnS2(U2,V2);
1053           ptdeb.SetValue(psol,TolArc,Standard_True);
1054           ptdeb.SetParameters(U1,V1,U2,V2);
1055           ptdeb.SetParameter(1.);
1056           wline->AddVertex(ptdeb);
1057           wline->SetFirstPoint(wline->NbVertex());
1058         }
1059         else 
1060         { 
1061           gp_Pnt psol = thelin->Value(1).Value();
1062           thelin->Value(1).ParametersOnS1(U1,V1);
1063           thelin->Value(1).ParametersOnS2(U2,V2);
1064           ptdeb.SetValue(psol,TolArc,Standard_False);
1065           ptdeb.SetParameters(U1,V1,U2,V2);
1066           ptdeb.SetParameter(1.);
1067           wline->AddVertex(ptdeb);
1068           wline->SetFirstPoint(wline->NbVertex());
1069         }
1070
1071
1072         if (   iwline->HasLastPoint() 
1073           && iwline->IsTangentAtEnd() == Standard_False) 
1074         {
1075           indlast = iwline->LastPointIndex();
1076           PPoint = seqpdep(indlast);
1077           tgline = PPoint.Direction3d().Reversed();
1078           Standard_Integer themult = PPoint.Multiplicity();
1079           for (Standard_Integer i=NbPointRst; i >=1; i--) {
1080             if (Destination(i) == indlast) {
1081               if (!reversed) {
1082                 Quad.Parameters(PPoint.Value(),U1,V1);
1083
1084                 if((V1 < Vmin) && (Vmin-V1 < TolV)) V1 = Vmin;
1085                 if((V1 > Vmax) && (V1-Vmax < TolV)) V1 = Vmax;
1086
1087                 PPoint.Parameters(themult,U2,V2);
1088                 Surf2->D1(U2,V2,ptbid,d1u,d1v); //-- @@@@
1089                 VecNormale = d1u.Crossed(d1v);                    //-- @@@@
1090               }
1091               else {
1092                 Quad.Parameters(PPoint.Value(),U2,V2);
1093
1094                 if((V2 < Vmin) && (Vmin-V2 < TolV)) V2 = Vmin;
1095                 if((V2 > Vmax) && (V2-Vmax < TolV)) V2 = Vmax;
1096
1097                 PPoint.Parameters(themult,U1,V1);
1098                 Surf1->D1(U1,V1,ptbid,d1u,d1v); //-- @@@@
1099                 VecNormale = d1u.Crossed(d1v);                    //-- @@@@
1100               }
1101
1102               ptfin.SetValue(PPoint.Value(),TolArc,Standard_False);
1103               ptfin.SetParameters(U1,V1,U2,V2);
1104               ptfin.SetParameter(Nbpts);
1105
1106               Recadre(reversed,typeS1,typeS2,ptfin,iwline,Nbpts-1,U1,V1,U2,V2);
1107
1108               currentarc = solrst.Point(i).Arc();
1109               currentparam = solrst.Point(i).Parameter();
1110               currentarc->D1(currentparam,p2d,d2d);
1111               tgrst.SetLinearForm(d2d.X(),d1u,d2d.Y(),d1v);
1112
1113
1114               Standard_Real squaremagnitudeVecNormale = VecNormale.SquareMagnitude();
1115               if(squaremagnitudeVecNormale > 1e-13) { 
1116                 DirNormale=VecNormale;
1117                 IntSurf::MakeTransition(tgline,tgrst,DirNormale,TLine,TArc);
1118               }
1119               else { 
1120                 TLine.SetValue(Standard_True,IntSurf_Undecided);
1121                 TArc.SetValue(Standard_True,IntSurf_Undecided);
1122               }
1123
1124
1125               ptfin.SetArc(reversed,currentarc,currentparam,TLine,TArc);
1126               if (!solrst.Point(i).IsNew()) {
1127                 ptfin.SetVertex(reversed,solrst.Point(i).Vertex());
1128               }
1129               wline->AddVertex(ptfin);
1130               if (themult == 0) {
1131                 wline->SetLastPoint(wline->NbVertex());
1132               }
1133
1134               themult--;
1135             }
1136           }
1137         }
1138         else if (iwline->IsTangentAtEnd()) 
1139         {
1140           gp_Pnt psol = thelin->Value(Nbpts).Value();
1141           thelin->Value(Nbpts).ParametersOnS1(U1,V1);
1142           thelin->Value(Nbpts).ParametersOnS2(U2,V2);
1143           ptfin.SetValue(psol,TolArc,Standard_True);
1144           ptfin.SetParameters(U1,V1,U2,V2);
1145           ptfin.SetParameter(Nbpts);
1146           wline->AddVertex(ptfin);
1147           wline->SetLastPoint(wline->NbVertex());
1148         }
1149         else 
1150         { 
1151           gp_Pnt psol = thelin->Value(Nbpts).Value();
1152           thelin->Value(Nbpts).ParametersOnS1(U1,V1);
1153           thelin->Value(Nbpts).ParametersOnS2(U2,V2);
1154           ptfin.SetValue(psol,TolArc,Standard_False);
1155           ptfin.SetParameters(U1,V1,U2,V2);
1156           ptfin.SetParameter(Nbpts);
1157           wline->AddVertex(ptfin);
1158           wline->SetLastPoint(wline->NbVertex());
1159         }
1160         //
1161         // Il faut traiter les points de passage.
1162         slin.Append(wline);
1163       }// if(Nbpts>=2) {
1164     }// for (j=1; j<=Nblines; j++) {
1165
1166     // ON GERE LES RACCORDS ENTRE LIGNES. ELLE NE PEUVENT SE RACCORDER
1167     // QUE SUR DES POINTS DE TANGENCE
1168
1169
1170     Nblines = slin.Length();
1171     for (Standard_Integer j=1; j<=Nblines-1; j++) {
1172       dofirst = dolast = Standard_False;
1173       const  Handle(IntPatch_Line)& slinj = slin(j);
1174       Handle(IntPatch_WLine) wlin1 (Handle(IntPatch_WLine)::DownCast (slinj));
1175       if (wlin1->HasFirstPoint()) {
1176         ptdeb = wlin1->FirstPoint(indfirst);
1177         if (ptdeb.IsTangencyPoint()) {
1178           dofirst = Standard_True;
1179         }
1180       }
1181       if (wlin1->HasLastPoint()) {
1182         ptfin = wlin1->LastPoint(indlast);
1183         if (ptfin.IsTangencyPoint()) {
1184           dolast = Standard_True;
1185         }
1186       }
1187
1188       if (dofirst || dolast) {
1189         for (Standard_Integer k=j+1; k<=Nblines;k++) {
1190           const  Handle(IntPatch_Line)& slink = slin(k);
1191           Handle(IntPatch_WLine) wlin2 (Handle(IntPatch_WLine)::DownCast (slink));
1192           if (wlin2->HasFirstPoint()) {
1193             ptbis = wlin2->FirstPoint(ind2);
1194             if (ptbis.IsTangencyPoint()) {
1195               if (dofirst ) {
1196                 if (ptdeb.Value().Distance(ptbis.Value()) <= TolArc) {
1197                   ptdeb.SetMultiple(Standard_True);
1198                   if (!ptbis.IsMultiple()) {
1199                     ptbis.SetMultiple(Standard_True);
1200                     wlin2->Replace(ind2,ptbis);
1201                   }
1202                 }
1203               }
1204               if (dolast ) {
1205                 if (ptfin.Value().Distance(ptbis.Value()) <= TolArc) {
1206                   ptfin.SetMultiple(Standard_True);
1207                   if (!ptbis.IsMultiple()) {
1208                     ptbis.SetMultiple(Standard_True);
1209                     wlin2->Replace(ind2,ptbis);
1210                   }
1211                 }
1212               }
1213             }
1214           }
1215           if (wlin2->HasLastPoint()) {
1216             ptbis = wlin2->LastPoint(ind2);
1217             if (ptbis.IsTangencyPoint()) {
1218               if (dofirst ) {
1219                 if (ptdeb.Value().Distance(ptbis.Value()) <= TolArc) {
1220                   ptdeb.SetMultiple(Standard_True);
1221                   if (!ptbis.IsMultiple()) {
1222                     ptbis.SetMultiple(Standard_True);
1223                     wlin2->Replace(ind2,ptbis);
1224                   }
1225                 }
1226               }
1227               if (dolast ) {
1228                 if (ptfin.Value().Distance(ptbis.Value()) <= TolArc) {
1229                   ptfin.SetMultiple(Standard_True);
1230                   if (!ptbis.IsMultiple()) {
1231                     ptbis.SetMultiple(Standard_True);
1232                     wlin2->Replace(ind2,ptbis);
1233                   }
1234                 }
1235               }
1236             }
1237           }
1238         }
1239         if(dofirst) 
1240           wlin1->Replace(indfirst,ptdeb);
1241         if(dolast) 
1242           wlin1->Replace(indlast,ptfin);
1243       }
1244     }
1245   }// if (seqpdep.Length() != 0 || seqpins.Length() != 0) {
1246   //
1247   // Treatment the segments
1248   NbSegm = solrst.NbSegments();
1249   if (NbSegm) {
1250     for(Standard_Integer i=1; i<=NbSegm; i++) {
1251       thesegm = solrst.Segment(i);  
1252       //Check if segment is degenerated
1253       if(thesegm.HasFirstPoint() && thesegm.HasLastPoint())
1254       {
1255         Standard_Real tol2 = Precision::Confusion();
1256         tol2 *= tol2;
1257         const gp_Pnt& aPf = thesegm.FirstPoint().Value();
1258         const gp_Pnt& aPl = thesegm.LastPoint().Value();
1259         if(aPf.SquareDistance(aPl) <= tol2)
1260         {
1261           //segment can be degenerated - check inner point
1262           paramf = thesegm.FirstPoint().Parameter();
1263           paraml = thesegm.LastPoint().Parameter();
1264           gp_Pnt2d _p2d = 
1265             thesegm.Curve()->Value(.57735 * paramf + 0.42265 * paraml);
1266           gp_Pnt aPm;
1267           if(reversed)
1268           {
1269             Surf1->D0(_p2d.X(), _p2d.Y(), aPm);
1270           }
1271           else
1272           {
1273             Surf2->D0(_p2d.X(), _p2d.Y(), aPm);
1274           }
1275           if(aPm.SquareDistance(aPf) <= tol2)
1276           {
1277             //Degenerated
1278             continue;
1279           }
1280         }
1281       }
1282
1283
1284       //----------------------------------------------------------------------      
1285       // on cree une ligne d intersection contenant uniquement le segment.
1286       // VOIR POUR LA TRANSITION DE LA LIGNE
1287       // On ajoute aussi un polygone pour le traitement des intersections
1288       // entre ligne et restrictions de la surface implicite (PutVertexOnLine)
1289       //----------------------------------------------------------------------
1290       //-- Calcul de la transition sur la rline (12 fev 97)
1291       //-- reversed a le sens de OnFirst
1292       //--
1293       dofirst = dolast  = Standard_False;
1294       procf = Standard_False;
1295       procl = Standard_False;
1296       IntSurf_Transition TLineUnk,TArcUnk;
1297
1298       IntPatch_Point _thepointAtBeg;
1299       IntPatch_Point _thepointAtEnd;
1300
1301       Standard_Boolean TransitionOK=Standard_False;
1302
1303       if(thesegm.HasFirstPoint()) { 
1304         Standard_Real _u1,_v1,_u2,_v2;
1305
1306         dofirst = Standard_True;
1307         PStartf = thesegm.FirstPoint();
1308         paramf = PStartf.Parameter();
1309
1310         gp_Pnt2d _p2d     = thesegm.Curve()->Value(paramf);
1311         Handle(Adaptor3d_HVertex) _vtx;
1312         if(PStartf.IsNew()==Standard_False) 
1313           _vtx= PStartf.Vertex();
1314         const gp_Pnt& _Pp = PStartf.Value();
1315         _thepointAtBeg.SetValue(_Pp,PStartf.Tolerance(),Standard_False);
1316         if (!reversed) { //-- typeS1 = Pln || Cyl || Sph || Cone
1317           Quad.Parameters(_Pp,_u1,_v1);
1318           _u2=_p2d.X(); _v2=_p2d.Y();
1319         }
1320         else {  //-- typeS1 != Pln && Cyl && Sph && Cone
1321           Quad.Parameters(_Pp,_u2,_v2);
1322           _u1=_p2d.X(); _v1=_p2d.Y();
1323         }
1324         _thepointAtBeg.SetParameters(_u1,_v1,_u2,_v2);
1325         _thepointAtBeg.SetParameter(paramf);
1326         if(PStartf.IsNew()==Standard_False) 
1327           _thepointAtBeg.SetVertex(reversed,_vtx);
1328         _thepointAtBeg.SetArc(reversed,thesegm.Curve(),paramf,TLineUnk,TArcUnk);
1329
1330
1331         gp_Vec d1u1,d1v1,d1u2,d1v2; gp_Vec2d _d2d;
1332         Surf1->D1(_u1,_v1,ptbid,d1u1,d1v1);
1333         norm1 = d1u1.Crossed(d1v1);
1334         Surf2->D1(_u2,_v2,ptbid,d1u2,d1v2);
1335         norm2 = d1u2.Crossed(d1v2);
1336
1337         thesegm.Curve()->D1(paramf,_p2d,_d2d);
1338         if(reversed) { 
1339           tgline.SetLinearForm(_d2d.X(),d1u1,_d2d.Y(),d1v1);
1340         }
1341         else { 
1342           tgline.SetLinearForm(_d2d.X(),d1u2,_d2d.Y(),d1v2);
1343         }
1344         _u1=tgline.DotCross(norm2,norm1);
1345         TransitionOK=Standard_True;
1346         if (_u1 > 0.00000001) {
1347           trans1 = IntSurf_Out;
1348           trans2 = IntSurf_In;
1349         }
1350         else if(_u1 < -0.00000001) { 
1351           trans1 = IntSurf_In;
1352           trans2 = IntSurf_Out;
1353         }
1354         else { 
1355           TransitionOK=Standard_False;
1356         }
1357       }
1358       if(thesegm.HasLastPoint()) {  
1359         Standard_Real _u1,_v1,_u2,_v2;
1360
1361         dolast = Standard_True;
1362         PStartl = thesegm.LastPoint();
1363         paraml = PStartl.Parameter();
1364
1365         gp_Pnt2d _p2d = thesegm.Curve()->Value(paraml);
1366         Handle(Adaptor3d_HVertex) _vtx;
1367         if(PStartl.IsNew()==Standard_False) 
1368           _vtx = PStartl.Vertex();
1369         const gp_Pnt& _Pp = PStartl.Value();
1370         IntPatch_Point _thepoint;
1371         _thepointAtEnd.SetValue(_Pp,PStartl.Tolerance(),Standard_False);
1372         if (!reversed) { //-- typeS1 = Pln || Cyl || Sph || Cone
1373           Quad.Parameters(_Pp,_u1,_v1);
1374           _u2=_p2d.X(); _v2=_p2d.Y();
1375         }
1376         else {  //-- typeS1 != Pln && Cyl && Sph && Cone
1377           Quad.Parameters(_Pp,_u2,_v2);
1378           _u1=_p2d.X(); _v1=_p2d.Y();
1379         }
1380         _thepointAtEnd.SetParameters(_u1,_v1,_u2,_v2);
1381         _thepointAtEnd.SetParameter(paraml);
1382         if(PStartl.IsNew()==Standard_False)
1383           _thepointAtEnd.SetVertex(reversed,_vtx);
1384         _thepointAtEnd.SetArc(reversed,thesegm.Curve(),paraml,TLineUnk,TArcUnk);
1385
1386
1387
1388         gp_Vec d1u1,d1v1,d1u2,d1v2; gp_Vec2d _d2d;
1389         Surf1->D1(_u1,_v1,ptbid,d1u1,d1v1);
1390         norm1 = d1u1.Crossed(d1v1);
1391         Surf2->D1(_u2,_v2,ptbid,d1u2,d1v2);
1392         norm2 = d1u2.Crossed(d1v2);
1393
1394         thesegm.Curve()->D1(paraml,_p2d,_d2d);
1395         if(reversed) { 
1396           tgline.SetLinearForm(_d2d.X(),d1u1,_d2d.Y(),d1v1);
1397         }
1398         else { 
1399           tgline.SetLinearForm(_d2d.X(),d1u2,_d2d.Y(),d1v2);
1400         }
1401         _u1=tgline.DotCross(norm2,norm1);
1402         TransitionOK=Standard_True;
1403         if (_u1 > 0.00000001) {
1404           trans1 = IntSurf_Out;
1405           trans2 = IntSurf_In;
1406         }
1407         else if(_u1 < -0.00000001) { 
1408           trans1 = IntSurf_In;
1409           trans2 = IntSurf_Out;
1410         }
1411         else { 
1412           TransitionOK=Standard_False;
1413         }       
1414       }
1415       if(TransitionOK==Standard_False) { 
1416         //-- rline = new IntPatch_RLine (thesegm.Curve(),reversed,Standard_False);
1417         rline =  new IntPatch_RLine (Standard_False);
1418         if(reversed) { 
1419           rline->SetArcOnS1(thesegm.Curve());
1420         }
1421         else { 
1422           rline->SetArcOnS2(thesegm.Curve());
1423         }
1424       }
1425       else { 
1426         //-- rline = new IntPatch_RLine (thesegm.Curve(),reversed,Standard_False,trans1,trans2);
1427         rline =  new IntPatch_RLine (Standard_False,trans1,trans2);
1428         if(reversed) { 
1429           rline->SetArcOnS1(thesegm.Curve());
1430         }
1431         else { 
1432           rline->SetArcOnS2(thesegm.Curve());
1433         }
1434       }
1435
1436       //------------------------------
1437       //-- Ajout des points 
1438       //--
1439       if (thesegm.HasFirstPoint()) {
1440         rline->AddVertex(_thepointAtBeg);
1441         rline->SetFirstPoint(rline->NbVertex());
1442       }
1443
1444       if (thesegm.HasLastPoint()) {
1445         rline->AddVertex(_thepointAtEnd);
1446         rline->SetLastPoint(rline->NbVertex());
1447       }
1448
1449       // Polygone sur restriction solution
1450       if (dofirst && dolast) {
1451         Standard_Real prm;
1452         gp_Pnt ptpoly;
1453         IntSurf_PntOn2S p2s;
1454         Handle(IntSurf_LineOn2S) Thelin = new IntSurf_LineOn2S ();
1455         Handle(Adaptor2d_HCurve2d) arcsegm = thesegm.Curve();
1456         Standard_Integer nbsample = 100;
1457
1458         if (!reversed) {
1459           for (Standard_Integer j=1; j<=nbsample; j++) {
1460             prm = paramf + (j-1)*(paraml-paramf)/(nbsample-1);
1461             arcsegm->D0(prm,p2d);
1462             Surf2->D0(p2d.X(),p2d.Y(),ptpoly);
1463
1464             Quad.Parameters(ptpoly,U1,V1);
1465             p2s.SetValue(ptpoly,U1,V1,p2d.X(),p2d.Y());
1466             Thelin->Add(p2s);
1467           }
1468         }
1469         else {
1470           for (Standard_Integer j=1; j<=nbsample; j++) {
1471             prm = paramf + (j-1)*(paraml-paramf)/(nbsample-1);
1472             arcsegm->D0(prm,p2d);
1473             Surf1->D0(p2d.X(),p2d.Y(),ptpoly);
1474
1475             Quad.Parameters(ptpoly,U2,V2);
1476             p2s.SetValue(ptpoly,p2d.X(),p2d.Y(),U2,V2);
1477             Thelin->Add(p2s);
1478           }
1479         }
1480         rline->Add(Thelin);
1481       }
1482
1483       if (dofirst || dolast) {
1484         Nblines = slin.Length();
1485         for (Standard_Integer j=1; j<=Nblines; j++) {
1486           const Handle(IntPatch_Line)& slinj = slin(j);
1487           typ = slinj->ArcType();
1488           if (typ == IntPatch_Walking) {
1489             Nbpts = Handle(IntPatch_WLine)::DownCast (slinj)->NbVertex();
1490           }
1491           else {
1492             Nbpts = Handle(IntPatch_RLine)::DownCast (slinj)->NbVertex();
1493           }
1494           for (Standard_Integer k=1; k<=Nbpts;k++) {
1495             if (typ == IntPatch_Walking) {
1496               ptdeb = Handle(IntPatch_WLine)::DownCast (slinj)->Vertex(k);
1497             }
1498             else {
1499               ptdeb = Handle(IntPatch_RLine)::DownCast (slinj)->Vertex(k);
1500             }
1501             if (dofirst) {
1502
1503               if (ptdeb.Value().Distance(PStartf.Value()) <=TolArc) {
1504                 ptdeb.SetMultiple(Standard_True);
1505                 if (typ == IntPatch_Walking) {
1506                   Handle(IntPatch_WLine)::DownCast (slinj)->Replace(k,ptdeb);
1507                 }
1508                 else {
1509                   Handle(IntPatch_RLine)::DownCast (slinj)->Replace(k,ptdeb);
1510                 }
1511                 ptdeb.SetParameter(paramf);
1512                 rline->AddVertex(ptdeb);
1513                 if (!procf){
1514                   procf=Standard_True;
1515                   rline->SetFirstPoint(rline->NbVertex());
1516                 }
1517               }
1518             }
1519             if (dolast) {
1520               if(dofirst) { //-- on recharge le ptdeb
1521                 if (typ == IntPatch_Walking) {
1522                   ptdeb = Handle(IntPatch_WLine)::DownCast (slinj)->Vertex(k);
1523                 }
1524                 else {
1525                   ptdeb = Handle(IntPatch_RLine)::DownCast (slinj)->Vertex(k);
1526                 }
1527               }
1528               if (ptdeb.Value().Distance(PStartl.Value()) <=TolArc) {
1529                 ptdeb.SetMultiple(Standard_True);
1530                 if (typ == IntPatch_Walking) {
1531                   Handle(IntPatch_WLine)::DownCast (slinj)->Replace(k,ptdeb);
1532                 }
1533                 else {
1534                   Handle(IntPatch_RLine)::DownCast (slinj)->Replace(k,ptdeb);
1535                 }
1536                 ptdeb.SetParameter(paraml);
1537                 rline->AddVertex(ptdeb);
1538                 if (!procl){
1539                   procl=Standard_True;
1540                   rline->SetLastPoint(rline->NbVertex());
1541                 }
1542               }
1543             }
1544           }
1545         }
1546       }
1547       slin.Append(rline);
1548     }
1549   }// if (NbSegm) 
1550   //
1551   // on traite les restrictions de la surface implicite
1552
1553   for (Standard_Integer i=1, aNbLin = slin.Length(); i<=aNbLin; i++)
1554   {
1555     Handle(IntPatch_Line)& aL = slin(i);
1556     
1557     if (!reversed)
1558       IntPatch_RstInt::PutVertexOnLine(aL,Surf1,D1,Surf2,Standard_True,TolTang);
1559     else
1560       IntPatch_RstInt::PutVertexOnLine(aL,Surf2,D2,Surf1,Standard_False,TolTang);
1561
1562     if(aL->ArcType() == IntPatch_Walking)
1563     {
1564       const Handle(IntPatch_WLine) aWL = Handle(IntPatch_WLine)::DownCast(aL);
1565       slin.Append(aWL);
1566       slin.Remove(i);
1567       i--;
1568       aNbLin--;
1569     }
1570   }
1571
1572   // Now slin is filled as follows: lower indices correspond to Restriction line,
1573   // after (higher indices) - only Walking-line.
1574
1575   const Standard_Real aTol3d = Max(Func.Tolerance(), TolTang); 
1576   const Handle(Adaptor3d_HSurface)& aQSurf = (reversed) ? Surf2 : Surf1;
1577   const Handle(Adaptor3d_HSurface)& anOtherSurf = (reversed) ? Surf1 : Surf2;
1578
1579   for (Standard_Integer i = 1; i <= slin.Length(); i++)
1580   {
1581     const Handle(IntPatch_PointLine)& aL1 = Handle(IntPatch_PointLine)::DownCast(slin(i));
1582     const Handle(IntPatch_RLine)& aRL1 = Handle(IntPatch_RLine)::DownCast(aL1);
1583
1584     if(aRL1.IsNull())
1585     {
1586       //Walking-Walking cases are not supported
1587       break;
1588     }
1589
1590     const Handle(Adaptor2d_HCurve2d)& anArc = aRL1->IsArcOnS1() ? 
1591                                               aRL1->ArcOnS1() :
1592                                               aRL1->ArcOnS2();
1593     if(anArc->Curve2d().GetType() != GeomAbs_Line)
1594     {
1595       //Restriction line must be isoline.
1596       //Other cases are not supported by
1597       //existing algorithms.
1598
1599       break;
1600     }
1601
1602     Standard_Boolean isFirstDeleted = Standard_False;
1603
1604     for(Standard_Integer j = i + 1; j <= slin.Length(); j++)
1605     {
1606       Handle(IntPatch_PointLine) aL2 = Handle(IntPatch_PointLine)::DownCast(slin(j));
1607       Handle(IntPatch_RLine) aRL2 = Handle(IntPatch_RLine)::DownCast(aL2);
1608
1609       //Here aL1 (i-th line) is Restriction-line and aL2 (j-th line) is
1610       //Restriction or Walking
1611
1612       if(!aRL2.IsNull())
1613       {
1614         const Handle(Adaptor2d_HCurve2d)& anArc2 = aRL2->IsArcOnS1() ?
1615                                                    aRL2->ArcOnS1() :
1616                                                    aRL2->ArcOnS2();
1617         if(anArc2->Curve2d().GetType() != GeomAbs_Line)
1618         {
1619           //Restriction line must be isoline.
1620           //Other cases are not supported by
1621           //existing algorithms.
1622
1623           continue;
1624         }
1625       }
1626
1627       //aDir can be equal to one of following four values only
1628       //(because Reastriction line is boundary of rectangular surface):
1629       //either {0, 1} or {0, -1} or {1, 0} or {-1, 0}.
1630       const gp_Dir2d aDir = anArc->Curve2d().Line().Direction();
1631
1632       Standard_Real aTol2d = anOtherSurf->UResolution(aTol3d),
1633                     aPeriod = anOtherSurf->IsVPeriodic() ? anOtherSurf->VPeriod() : 0.0;
1634
1635       if(Abs(aDir.X()) < 0.5)
1636       {//Restriction directs along V-direction
1637         aTol2d = anOtherSurf->VResolution(aTol3d);
1638         aPeriod = anOtherSurf->IsUPeriodic() ? anOtherSurf->UPeriod() : 0.0;
1639       }
1640
1641       const Standard_Boolean isCoincide = IsCoincide(Func, aL2, anArc, aRL1->IsArcOnS1(),
1642                                                       aTol3d, aTol2d, aPeriod);
1643
1644       if(isCoincide)
1645       {
1646         if(aRL2.IsNull())
1647         {//Delete Walking-line
1648           slin.Remove(j);
1649           j--;
1650         }
1651         else
1652         {//Restriction-Restriction
1653           const Handle(Adaptor2d_HCurve2d)& anArc2 = aRL2->IsArcOnS1() ?
1654                                                      aRL2->ArcOnS1() :
1655                                                      aRL2->ArcOnS2();
1656
1657           const Standard_Real aRange2 = anArc2->LastParameter() - 
1658                                         anArc2->FirstParameter();
1659           const Standard_Real aRange1 = anArc->LastParameter() -
1660                                         anArc->FirstParameter();
1661
1662           if(aRange2 > aRange1)
1663           {
1664             isFirstDeleted = Standard_True;
1665             break;
1666           }
1667           else
1668           {//Delete j-th line
1669             slin.Remove(j);
1670             j--;
1671           }
1672         }
1673       }
1674     } //for(Standard_Integer j = i + 1; j <= slin.Length(); j++)
1675
1676     if(isFirstDeleted)
1677     {//Delete i-th line
1678       slin.Remove(i--);
1679     }
1680   }//for (Standard_Integer i = 1; i <= slin.Length(); i++)
1681
1682   empt = (slin.Length() == 0 && spnt.Length() == 0);
1683   done = Standard_True;
1684
1685
1686   if(slin.Length() == 0)
1687     return;
1688
1689   Standard_Boolean isDecomposeRequired =  (Quad.TypeQuadric() == GeomAbs_Cone) || 
1690                                           (Quad.TypeQuadric() == GeomAbs_Sphere) ||
1691                                           (Quad.TypeQuadric() == GeomAbs_Cylinder) ||
1692                                           (Quad.TypeQuadric() == GeomAbs_Torus);
1693
1694   if(!isDecomposeRequired)
1695     return;
1696
1697   // post processing for cones and spheres
1698
1699   const Handle(Adaptor3d_TopolTool)& PDomain = (reversed) ? D1 : D2;
1700
1701   IntPatch_SequenceOfLine dslin;
1702   Standard_Boolean isDecompose = Standard_False;
1703   for(Standard_Integer i = 1; i <= slin.Length(); i++ )
1704   {
1705     if(DecomposeResult( Handle(IntPatch_PointLine)::DownCast(slin(i)),
1706                                         reversed, Quad, PDomain, aQSurf,
1707                                         anOtherSurf, TolArc, aTol3d, dslin))
1708     {
1709       isDecompose = Standard_True;
1710     }
1711   }
1712
1713   if(!isDecompose)
1714     return;
1715
1716   slin.Clear();
1717   for(Standard_Integer i = 1; i <= dslin.Length(); i++ )
1718     slin.Append(dslin(i));
1719 }
1720
1721 // correct U parameter of the start point of line on Quadric
1722 // (change 0->2PI or vs, if necessary)
1723 static Standard_Real AdjustUFirst(Standard_Real U1,Standard_Real U2)
1724 {
1725   Standard_Real u = U1;
1726
1727   // case: no adjustment
1728   if( U1 > 0. && U1 < (2.*M_PI) )
1729     return u;
1730
1731   // case: near '0'
1732   if( U1 == 0. || fabs(U1) <= 1.e-9 ) {
1733     if( U2 > 0. && U2 < (2.*M_PI) )
1734       u = ( U2 < ((2.*M_PI)-U2) ) ? 0. : (2.*M_PI);
1735     else {
1736       Standard_Real uu = U2;
1737       if( U2 > (2.*M_PI) )
1738         while( uu > (2.*M_PI) )
1739           uu -= (2.*M_PI);
1740       else 
1741         while( uu < 0.)
1742           uu += (2.*M_PI);
1743
1744       u = ( uu < ((2.*M_PI)-uu) ) ? 0. : (2.*M_PI);
1745     }
1746   }
1747   // case: near '2PI'
1748   else if( U1 == (2.*M_PI) || fabs((2.*M_PI)-fabs(U1)) <= 1.e-9 ) {
1749     if( U2 > 0. && U2 < (2.*M_PI) )
1750       u = ( U2 < ((2.*M_PI)-U2) ) ? 0. : (2.*M_PI);
1751     else {
1752       Standard_Real uu = U2;
1753       if( U2 > (2.*M_PI) )
1754         while( uu > (2.*M_PI) )
1755           uu -= (2.*M_PI);
1756       else 
1757         while( uu < 0.)
1758           uu += (2.*M_PI);
1759
1760       u = ( uu < ((2.*M_PI)-uu) ) ? 0. : (2.*M_PI);
1761     }
1762   }
1763   // case: '<0. || >2PI'
1764   else {
1765     if(U1 < 0.)
1766       while(u < 0.)
1767         u += 2.*M_PI;
1768     if(U1 > (2.*M_PI))
1769       while(u > (2.*M_PI))
1770         u -= (2.*M_PI);
1771   }
1772
1773   return u;
1774 }
1775
1776 // collect vertices, reject equals
1777 static Handle(IntSurf_LineOn2S) GetVertices(const Handle(IntPatch_PointLine)& thePLine,
1778                                             const Standard_Real           TOL3D,
1779                                             const Standard_Real           TOL2D)
1780 {
1781   //  Standard_Real TOL3D = 1.e-12, TOL2D = 1.e-8;
1782
1783   Handle(IntSurf_LineOn2S) vertices = new IntSurf_LineOn2S();
1784
1785   Standard_Real U1 = 0., U2 = 0., V1 = 0., V2 = 0.;
1786   Standard_Integer i = 0, k = 0;
1787   Standard_Integer NbVrt = thePLine->NbVertex();
1788
1789   TColStd_Array1OfInteger anVrts(1,NbVrt);
1790   anVrts.Init(0);
1791
1792   // check equal vertices
1793   for(i = 1; i <= NbVrt; i++) {
1794
1795     if( anVrts(i) == -1 ) continue;
1796
1797     const IntPatch_Point& Pi = thePLine->Vertex(i);
1798
1799     for(k = (i+1); k <= NbVrt; k++) {
1800
1801       if( anVrts(k) == -1 ) continue;
1802
1803       const IntPatch_Point& Pk = thePLine->Vertex(k);
1804
1805       if(Pi.Value().Distance(Pk.Value()) <= TOL3D) {
1806         // suggest the points are equal;
1807         // test 2d parameters on surface
1808         Standard_Boolean sameU1 = Standard_False;
1809         Standard_Boolean sameV1 = Standard_False;
1810         Standard_Boolean sameU2 = Standard_False;
1811         Standard_Boolean sameV2 = Standard_False;
1812
1813         Pi.ParametersOnS1(U1,V1);
1814         Pk.ParametersOnS1(U2,V2);
1815         if(fabs(U1-U2) <= TOL2D) sameU1 = Standard_True;
1816         if(fabs(V1-V2) <= TOL2D) sameV1 = Standard_True;
1817
1818         Pi.ParametersOnS2(U1,V1);
1819         Pk.ParametersOnS2(U2,V2);
1820         if(fabs(U1-U2) <= TOL2D) sameU2 = Standard_True;
1821         if(fabs(V1-V2) <= TOL2D) sameV2 = Standard_True;
1822
1823         if((sameU1 && sameV1) && (sameU2 && sameV2))
1824           anVrts(k) = -1;
1825       }
1826     }
1827   }
1828
1829   // copy further processed vertices
1830   for(i = 1; i <= NbVrt; i++) {
1831     if( anVrts(i) == -1 ) continue;
1832     vertices->Add(thePLine->Vertex(i).PntOn2S());
1833   }
1834   return vertices;
1835 }
1836
1837 static void SearchVertices(const Handle(IntSurf_LineOn2S)& Line,
1838   const Handle(IntSurf_LineOn2S)& Vertices,
1839   TColStd_Array1OfInteger&        PTypes)
1840 {
1841   Standard_Integer nbp = Line->NbPoints(), nbv = Vertices->NbPoints();
1842   Standard_Integer ip = 0, iv = 0;
1843   for(ip = 1; ip <= nbp; ip++) {
1844     const IntSurf_PntOn2S& aP = Line->Value(ip);
1845     Standard_Integer type = 0;
1846     for(iv = 1; iv <= nbv; iv++) {
1847       const IntSurf_PntOn2S& aV = Vertices->Value(iv);
1848       if(aP.IsSame(aV, Precision::Confusion(), Precision::PConfusion())) {
1849         type = iv; 
1850         break;
1851       }
1852     }
1853     PTypes(ip) = type;
1854   }
1855 }
1856
1857 static inline Standard_Boolean IsSeamParameter(const Standard_Real U,
1858   const Standard_Real TOL2D)
1859 {
1860   return (fabs(U) <= TOL2D || fabs(2.*M_PI - U) <= TOL2D);
1861 }
1862
1863 static inline Standard_Real AdjustU(const Standard_Real U)
1864 {
1865   Standard_Real u = U, DBLPI = 2.*M_PI;
1866   if(u < 0. || u > DBLPI) {
1867     if(u < 0.)
1868       while(u < 0.)
1869         u += DBLPI;
1870     else
1871       while(u > DBLPI)
1872         u -= DBLPI;
1873   }
1874   return u;
1875 }
1876
1877 static inline void Correct2DBounds(const Standard_Real UF,
1878   const Standard_Real UL,
1879   const Standard_Real VF,
1880   const Standard_Real VL,
1881   const Standard_Real TOL2D,
1882   Standard_Real&      U,
1883   Standard_Real&      V)
1884 {
1885   Standard_Real Eps = 1.e-16;
1886   Standard_Real dUF = fabs(U - UF);
1887   Standard_Real dUL = fabs(U - UL);
1888   Standard_Real dVF = fabs(V - VF);
1889   Standard_Real dVL = fabs(V - VL);
1890   if(dUF <= TOL2D && dUF > Eps) U = UF;
1891   if(dUL <= TOL2D && dUL > Eps) U = UL;
1892   if(dVF <= TOL2D && dVF > Eps) V = VF;
1893   if(dVL <= TOL2D && dVL > Eps) V = VL;
1894 }
1895
1896 static void AdjustLine(Handle(IntSurf_LineOn2S)& Line,
1897   const Standard_Boolean    IsReversed,
1898   const Handle(Adaptor3d_HSurface)&         QSurf,
1899   const Standard_Real       TOL2D)
1900 {
1901   Standard_Real VF = QSurf->FirstVParameter();
1902   Standard_Real VL = QSurf->LastVParameter();
1903   Standard_Real UF = QSurf->FirstUParameter();
1904   Standard_Real UL = QSurf->LastUParameter();
1905
1906   Standard_Integer nbp = Line->NbPoints(), ip = 0;
1907   Standard_Real U = 0., V = 0.;
1908   for(ip = 1; ip <= nbp; ip++) {
1909     if(IsReversed) {
1910       Line->Value(ip).ParametersOnS2(U,V);
1911       U = AdjustU(U);
1912       Correct2DBounds(UF,UL,VF,VL,TOL2D,U,V);
1913       Line->SetUV(ip,Standard_False,U,V);
1914     }
1915     else {
1916       Line->Value(ip).ParametersOnS1(U,V);
1917       U = AdjustU(U);
1918       Correct2DBounds(UF,UL,VF,VL,TOL2D,U,V);
1919       Line->SetUV(ip,Standard_True,U,V);
1920     }
1921   }
1922 }
1923
1924 static Standard_Boolean InsertSeamVertices(Handle(IntSurf_LineOn2S)&       Line,
1925   const Standard_Boolean          IsReversed,
1926   Handle(IntSurf_LineOn2S)&       Vertices,
1927   const TColStd_Array1OfInteger&  PTypes,
1928   const Standard_Real             TOL2D)
1929 {
1930   Standard_Boolean result = Standard_False;
1931   Standard_Integer ip = 0, nbp = Line->NbPoints();
1932   Standard_Real U = 0., V = 0.;
1933   for(ip = 1; ip <= nbp; ip++) {
1934     Standard_Integer ipt = PTypes(ip);
1935     if(ipt != 0) {
1936       const IntSurf_PntOn2S& aP = Line->Value(ip);
1937       if(IsReversed)
1938         aP.ParametersOnS2(U,V); // S2 - quadric
1939       else
1940         aP.ParametersOnS1(U,V); // S1 - quadric
1941       U = AdjustU(U);
1942       if(IsSeamParameter(U,TOL2D)) {
1943         if(ip == 1 || ip == nbp) {
1944           Standard_Real U1 = 0., V1 = 0.;
1945           Standard_Integer ipp = (ip == 1) ? (ip+1) : (ip-1);
1946           if(IsReversed)
1947             Line->Value(ipp).ParametersOnS2(U1,V1); // S2 - quadric
1948           else
1949             Line->Value(ipp).ParametersOnS1(U1,V1); // S1 - quadric
1950           Standard_Real u = AdjustUFirst(U,U1);
1951           if(fabs(u-U) >= 1.5*M_PI) {
1952             Standard_Real U2 = 0., V2 = 0.;
1953             if(IsReversed) {
1954               Line->Value(ip).ParametersOnS1(U2,V2); // prm
1955               Line->SetUV(ip,Standard_False,u,V);
1956               Line->SetUV(ip,Standard_True,U2,V2);
1957             }
1958             else {
1959               Line->Value(ip).ParametersOnS2(U2,V2); // prm
1960               Line->SetUV(ip,Standard_True,u,V);
1961               Line->SetUV(ip,Standard_False,U2,V2);
1962             }
1963           }
1964         }
1965         else {
1966           Standard_Integer ipp = ip - 1;
1967           Standard_Integer ipn = ip + 1;
1968           Standard_Real U1 = 0., V1 = 0., U2 = 0., V2 = 0.;
1969           if(IsReversed) {
1970             Line->Value(ipp).ParametersOnS2(U1,V1); // quad
1971             Line->Value(ipn).ParametersOnS2(U2,V2); // quad
1972           }
1973           else {
1974             Line->Value(ipp).ParametersOnS1(U1,V1); // quad
1975             Line->Value(ipn).ParametersOnS1(U2,V2); // quad
1976           }
1977           U1 = AdjustU(U1);
1978           U2 = AdjustU(U2);
1979           Standard_Boolean pnearZero = (fabs(U1) < fabs(2.*M_PI-U1)) ? Standard_True : Standard_False;
1980           Standard_Boolean cnearZero = (fabs(U) < fabs(2.*M_PI-U)) ? Standard_True : Standard_False;
1981           if(pnearZero == cnearZero) {
1982             if(!IsSeamParameter(U2,TOL2D) && !IsSeamParameter(U1,TOL2D)) {
1983               Standard_Real nU = (cnearZero) ? (2.*M_PI) : 0.;
1984               IntSurf_PntOn2S nP;
1985               nP.SetValue(aP.Value());
1986               Standard_Real U3 = 0., V3 = 0.;
1987               if(IsReversed) {
1988                 Line->Value(ip).ParametersOnS1(U3,V3); // prm
1989                 nP.SetValue(Standard_False,nU,V);
1990                 nP.SetValue(Standard_True,U3,V3);
1991               }
1992               else {
1993                 Line->Value(ip).ParametersOnS2(U3,V3); // prm
1994                 nP.SetValue(Standard_True,nU,V);
1995                 nP.SetValue(Standard_False,U3,V3);
1996               }
1997               Line->InsertBefore(ipn,nP);
1998               Vertices->Add(nP);
1999               result = Standard_True;
2000               break;
2001             }
2002           }
2003           else {
2004             if(!IsSeamParameter(U2,TOL2D) && !IsSeamParameter(U1,TOL2D)) {
2005               Standard_Real nU = (cnearZero) ? (2.*M_PI) : 0.;
2006               IntSurf_PntOn2S nP;
2007               nP.SetValue(aP.Value());
2008               Standard_Real U3 = 0., V3 = 0.;
2009               if(IsReversed) {
2010                 Line->Value(ip).ParametersOnS1(U3,V3); // prm
2011                 nP.SetValue(Standard_False,nU,V);
2012                 nP.SetValue(Standard_True,U3,V3);
2013               }
2014               else {
2015                 Line->Value(ip).ParametersOnS2(U3,V3); // prm
2016                 nP.SetValue(Standard_True,nU,V);
2017                 nP.SetValue(Standard_False,U3,V3);
2018               }
2019               Line->InsertBefore(ip,nP);
2020               Vertices->Add(nP);
2021               result = Standard_True;
2022               break;
2023             }
2024             else {
2025               // Line->InsertBefore(ip,Line->Value(ipn));
2026               // Line->RemovePoint(ip+2);
2027               // result = Standard_True;
2028               // cout << "swap vertex " << endl;
2029               // break;
2030             }
2031           }
2032         }
2033       }
2034     }
2035   }
2036   return result;
2037 }
2038
2039 static void ToSmooth( const Handle(IntSurf_LineOn2S)& Line,
2040                       const Standard_Boolean          IsReversed,
2041                       const IntSurf_Quadric&          Quad,
2042                       const Standard_Boolean          IsFirst,
2043                       Standard_Real&                  D3D)
2044 {
2045   if(Line->NbPoints() <= 10)
2046     return;
2047
2048   D3D = 0.;
2049   Standard_Integer NbTestPnts = Line->NbPoints() / 5;
2050   if(NbTestPnts < 5) NbTestPnts = 5;
2051
2052   Standard_Integer startp = (IsFirst) ? 2 : (Line->NbPoints() - NbTestPnts - 2);
2053   Standard_Integer ip = 0;
2054   Standard_Real Uc = 0., Vc = 0., Un = 0., Vn = 0., DDU = 0., DDV = 0.;
2055
2056   for(ip = startp; ip <= NbTestPnts; ip++) {
2057     if(IsReversed) {
2058       Line->Value(ip).ParametersOnS2(Uc,Vc); // S2 - quadric
2059       Line->Value(ip+1).ParametersOnS2(Un,Vn);
2060     }
2061     else {
2062       Line->Value(ip).ParametersOnS1(Uc,Vc); // S1 - quadric
2063       Line->Value(ip+1).ParametersOnS1(Un,Vn);
2064     }
2065     DDU += fabs(fabs(Uc)-fabs(Un));
2066     DDV += fabs(fabs(Vc)-fabs(Vn));
2067
2068     if(ip > startp) {
2069       Standard_Real DP = Line->Value(ip).Value().Distance(Line->Value(ip-1).Value());
2070       D3D += DP;
2071     }
2072   }
2073
2074   DDU /= (Standard_Real) NbTestPnts + 1;
2075   DDV /= (Standard_Real) NbTestPnts + 1;
2076
2077   D3D /= (Standard_Real) NbTestPnts + 1;
2078
2079
2080   Standard_Integer Index1 = (IsFirst) ? 1 : (Line->NbPoints());
2081   Standard_Integer Index2 = (IsFirst) ? 2 : (Line->NbPoints()-1);
2082   Standard_Integer Index3 = (IsFirst) ? 3 : (Line->NbPoints()-2);
2083
2084   Standard_Boolean doU = Standard_False;
2085
2086   Standard_Real U1 = 0., U2 = 0., V1 = 0., V2 = 0., U3 = 0., V3 = 0.;
2087
2088   if(IsReversed) {
2089     Line->Value(Index1).ParametersOnS2(U1,V1); // S2 - quadric
2090     Line->Value(Index2).ParametersOnS2(U2,V2);
2091     Line->Value(Index3).ParametersOnS2(U3,V3);
2092   }
2093   else {
2094     Line->Value(Index1).ParametersOnS1(U1,V1); // S1 - quadric
2095     Line->Value(Index2).ParametersOnS1(U2,V2);
2096     Line->Value(Index3).ParametersOnS1(U3,V3);
2097   }
2098
2099   if(!doU && Quad.TypeQuadric() == GeomAbs_Sphere) {
2100     if(fabs(fabs(U1)-fabs(U2)) > (M_PI/16.)) doU = Standard_True;
2101
2102     if(doU && (fabs(U1) <= 1.e-9 || fabs(U1-2.*M_PI) <= 1.e-9)) {
2103       if(fabs(V1-M_PI/2.) <= 1.e-9 || fabs(V1+M_PI/2.) <= 1.e-9) {}
2104       else {
2105         doU = Standard_False;
2106       }
2107     }
2108   }
2109
2110   if(Quad.TypeQuadric() == GeomAbs_Cone) {
2111     Standard_Real Uapx = 0., Vapx = 0.;
2112     Quad.Parameters(Quad.Cone().Apex(),Uapx,Vapx);
2113
2114     if(fabs(fabs(U1)-fabs(U2)) > M_PI/32.) doU = Standard_True;
2115
2116     if(doU && (fabs(U1) <= 1.e-9 || fabs(U1-2.*M_PI) <= 1.e-9)) {
2117       if(fabs(V1-Vapx) <= 1.e-9) {}
2118       else {
2119         doU = Standard_False;
2120       }
2121     }
2122   }
2123
2124   if(doU) {
2125     Standard_Real dU = Min((DDU/10.),5.e-8);
2126     Standard_Real U = (U2 > U3) ? (U2 + dU) : (U2 - dU);
2127     if(IsReversed)
2128       Line->SetUV(Index1,Standard_False,U,V1);
2129     else
2130       Line->SetUV(Index1,Standard_True,U,V1);
2131     U1 = U;
2132   }
2133
2134
2135 static Standard_Boolean TestMiddleOnPrm(const IntSurf_PntOn2S& aP,
2136                                         const IntSurf_PntOn2S& aV,
2137                                         const Standard_Boolean IsReversed,
2138                                         const Standard_Real    ArcTol,
2139                                         const Handle(Adaptor3d_TopolTool)&  PDomain)
2140
2141 {
2142   Standard_Boolean result = Standard_False;
2143   Standard_Real Up = 0., Vp = 0., Uv = 0., Vv = 0.;
2144   if(IsReversed) {
2145     aP.ParametersOnS1(Up,Vp); //S1 - parametric
2146     aV.ParametersOnS1(Uv,Vv);
2147   }
2148   else {
2149     aP.ParametersOnS2(Up,Vp); // S2 - parametric
2150     aV.ParametersOnS2(Uv,Vv);
2151   }
2152   Standard_Real Um = (Up + Uv)*0.5, Vm = (Vp + Vv)*0.5;
2153   gp_Pnt2d a2DPntM(Um,Vm);
2154   TopAbs_State PosM = PDomain->Classify(a2DPntM,ArcTol);
2155   if(PosM == TopAbs_ON || PosM == TopAbs_IN )
2156     result = Standard_True;
2157   return result;
2158 }
2159
2160 static void VerifyVertices( const Handle(IntSurf_LineOn2S)&    Line,
2161                             const Standard_Boolean             IsReversed,
2162                             const Handle(IntSurf_LineOn2S)&    Vertices,
2163                             const Standard_Real                TOL2D,
2164                             const Standard_Real       ArcTol,
2165                             const Handle(Adaptor3d_TopolTool)& PDomain,
2166                             IntSurf_PntOn2S&          VrtF,
2167                             Standard_Boolean&         AddFirst,
2168                             IntSurf_PntOn2S&          VrtL,
2169                             Standard_Boolean&         AddLast)
2170 {
2171   Standard_Integer nbp = Line->NbPoints(), nbv = Vertices->NbPoints();
2172   Standard_Integer FIndexSame = 0, FIndexNear = 0, LIndexSame = 0, LIndexNear = 0;
2173   const IntSurf_PntOn2S& aPF = Line->Value(1);
2174   const IntSurf_PntOn2S& aPL = Line->Value(nbp);
2175   Standard_Real UF = 0., VF = 0., UL = 0., VL = 0.;
2176   if(IsReversed) {
2177     aPF.ParametersOnS2(UF,VF);
2178     aPL.ParametersOnS2(UL,VL);
2179   }
2180   else {
2181     aPF.ParametersOnS1(UF,VF);
2182     aPL.ParametersOnS1(UL,VL);
2183   }
2184   gp_Pnt2d a2DPF(UF,VF);
2185   gp_Pnt2d a2DPL(UL,VL);
2186   Standard_Real DistMinF = 1.e+100, DistMinL = 1.e+100;
2187   Standard_Integer FConjugated = 0, LConjugated = 0;
2188
2189   Standard_Integer iv = 0;
2190
2191   for(iv = 1; iv <= nbv; iv++) {
2192     Standard_Real Uv = 0., Vv = 0.;
2193     if(IsReversed) {
2194       Vertices->Value(iv).ParametersOnS2(Uv,Vv);
2195       Uv = AdjustU(Uv);
2196       Vertices->SetUV(iv,Standard_False,Uv,Vv);
2197     }
2198     else {
2199       Vertices->Value(iv).ParametersOnS1(Uv,Vv);
2200       Uv = AdjustU(Uv);
2201       Vertices->SetUV(iv,Standard_True,Uv,Vv);
2202     }
2203   }
2204
2205   for(iv = 1; iv <= nbv; iv++) {
2206     const IntSurf_PntOn2S& aV = Vertices->Value(iv);
2207     if(aPF.IsSame(aV, Precision::Confusion(), Precision::PConfusion())) {
2208       FIndexSame = iv;
2209       break;
2210     }
2211     else {
2212       Standard_Real Uv = 0., Vv = 0.;
2213       if(IsReversed)
2214         aV.ParametersOnS2(Uv,Vv);
2215       else
2216         aV.ParametersOnS1(Uv,Vv);
2217       gp_Pnt2d a2DV(Uv,Vv);
2218       Standard_Real Dist = a2DV.Distance(a2DPF);
2219       if(Dist < DistMinF) {
2220         DistMinF = Dist;
2221         FIndexNear = iv;
2222         if(FConjugated != 0)
2223           FConjugated = 0;
2224       }
2225       if(IsSeamParameter(Uv,TOL2D)) {
2226         Standard_Real Ucv = (fabs(Uv) < fabs(2.*M_PI-Uv)) ? (2.*M_PI) : 0.;
2227         gp_Pnt2d a2DCV(Ucv,Vv);
2228         Standard_Real CDist = a2DCV.Distance(a2DPF);
2229         if(CDist < DistMinF) {
2230           DistMinF = CDist;
2231           FConjugated = iv;
2232           FIndexNear = iv;
2233         }
2234       }
2235     }
2236   }
2237
2238   for(iv = 1; iv <= nbv; iv++) {
2239     const IntSurf_PntOn2S& aV = Vertices->Value(iv);
2240     if(aPL.IsSame(aV, Precision::Confusion(), Precision::PConfusion())) {
2241       LIndexSame = iv;
2242       break;
2243     }
2244     else {
2245       Standard_Real Uv = 0., Vv = 0.;
2246       if(IsReversed)
2247         aV.ParametersOnS2(Uv,Vv);
2248       else
2249         aV.ParametersOnS1(Uv,Vv);
2250       gp_Pnt2d a2DV(Uv,Vv);
2251       Standard_Real Dist = a2DV.Distance(a2DPL);
2252       if(Dist < DistMinL) {
2253         DistMinL = Dist;
2254         LIndexNear = iv;
2255         if(LConjugated != 0)
2256           LConjugated = 0;
2257       }
2258       if(IsSeamParameter(Uv,TOL2D)) {
2259         Standard_Real Ucv = (fabs(Uv) < fabs(2.*M_PI-Uv)) ? (2.*M_PI) : 0.;
2260         gp_Pnt2d a2DCV(Ucv,Vv);
2261         Standard_Real CDist = a2DCV.Distance(a2DPL);
2262         if(CDist < DistMinL) {
2263           DistMinL = CDist;
2264           LConjugated = iv;
2265           LIndexNear = iv;
2266         }
2267       }
2268     }
2269   }
2270
2271   AddFirst = Standard_False;
2272   AddLast  = Standard_False;
2273
2274   if(FIndexSame == 0) {
2275     if(FIndexNear != 0) {
2276       const IntSurf_PntOn2S& aV = Vertices->Value(FIndexNear);
2277       Standard_Real Uv = 0., Vv = 0.;
2278       if(IsReversed)
2279         aV.ParametersOnS2(Uv,Vv);
2280       else
2281         aV.ParametersOnS1(Uv,Vv);
2282       if(IsSeamParameter(Uv,TOL2D)) {
2283         Standard_Real Ucv = (fabs(Uv) < fabs(2.*M_PI-Uv)) ? (2.*M_PI) : 0.;
2284         Standard_Boolean test = TestMiddleOnPrm(aPF,aV,IsReversed,ArcTol,PDomain);
2285         if(test) {
2286           VrtF.SetValue(aV.Value());
2287           if(IsReversed) {
2288             Standard_Real U2 = 0., V2 = 0.;
2289             aV.ParametersOnS1(U2,V2); // S1 - prm
2290             VrtF.SetValue(Standard_True,U2,V2);
2291             if(FConjugated == 0)
2292               VrtF.SetValue(Standard_False,Uv,Vv);
2293             else
2294               VrtF.SetValue(Standard_False,Ucv,Vv);
2295           }
2296           else {
2297             Standard_Real U2 = 0., V2 = 0.;
2298             aV.ParametersOnS2(U2,V2); // S2 - prm
2299             VrtF.SetValue(Standard_False,U2,V2);
2300             if(FConjugated == 0)
2301               VrtF.SetValue(Standard_True,Uv,Vv);
2302             else
2303               VrtF.SetValue(Standard_True,Ucv,Vv);
2304           }
2305           Standard_Real Dist3D = VrtF.Value().Distance(aPF.Value());
2306           if(Dist3D > 1.5e-7 && DistMinF > TOL2D) {
2307             AddFirst = Standard_True;
2308           }
2309         }
2310       }
2311       else {
2312         // to do: analyze internal vertex
2313       }
2314     }
2315   }
2316
2317   if(LIndexSame == 0) {
2318     if(LIndexNear != 0) {
2319       const IntSurf_PntOn2S& aV = Vertices->Value(LIndexNear);
2320       Standard_Real Uv = 0., Vv = 0.;
2321       if(IsReversed)
2322         aV.ParametersOnS2(Uv,Vv);
2323       else
2324         aV.ParametersOnS1(Uv,Vv);
2325       if(IsSeamParameter(Uv,TOL2D)) {
2326         Standard_Real Ucv = (fabs(Uv) < fabs(2.*M_PI-Uv)) ? (2.*M_PI) : 0.;
2327         Standard_Boolean test = TestMiddleOnPrm(aPL,aV,IsReversed,ArcTol,PDomain);
2328         if(test) {
2329           VrtL.SetValue(aV.Value());
2330           if(IsReversed) {
2331             Standard_Real U2 = 0., V2 = 0.;
2332             aV.ParametersOnS1(U2,V2); // S1 - prm
2333             VrtL.SetValue(Standard_True,U2,V2);
2334             if(LConjugated == 0)
2335               VrtL.SetValue(Standard_False,Uv,Vv);
2336             else
2337               VrtL.SetValue(Standard_False,Ucv,Vv);
2338           }
2339           else {
2340             Standard_Real U2 = 0., V2 = 0.;
2341             aV.ParametersOnS2(U2,V2); // S2 - prm
2342             VrtL.SetValue(Standard_False,U2,V2);
2343             if(LConjugated == 0)
2344               VrtL.SetValue(Standard_True,Uv,Vv);
2345             else
2346               VrtL.SetValue(Standard_True,Ucv,Vv);
2347           }
2348           Standard_Real Dist3D = VrtL.Value().Distance(aPL.Value());
2349           if(Dist3D > 1.5e-7 && DistMinL > TOL2D) {
2350             AddLast = Standard_True;
2351           }
2352         }
2353       }
2354       else {
2355         // to do: analyze internal vertex
2356       }
2357     }
2358   }
2359 }
2360
2361 static Standard_Boolean AddVertices(Handle(IntSurf_LineOn2S)& Line,
2362   const IntSurf_PntOn2S&    VrtF,
2363   const Standard_Boolean    AddFirst,
2364   const IntSurf_PntOn2S&    VrtL,
2365   const Standard_Boolean    AddLast,
2366   const Standard_Real       D3DF,
2367   const Standard_Real       D3DL)
2368 {
2369   Standard_Boolean result = Standard_False;
2370   if(AddFirst) {
2371     Standard_Real DF = Line->Value(1).Value().Distance(VrtF.Value());
2372     if((D3DF*2.) > DF && DF > 1.5e-7) {
2373       Line->InsertBefore(1,VrtF);
2374       result = Standard_True;
2375     }
2376   }
2377   if(AddLast) {
2378     Standard_Real DL = Line->Value(Line->NbPoints()).Value().Distance(VrtL.Value());
2379     if((D3DL*2.) > DL && DL > 1.5e-7) {
2380       Line->Add(VrtL);
2381       result = Standard_True;
2382     }
2383   }
2384   return result;
2385 }
2386
2387
2388 static void PutIntVertices(const Handle(IntPatch_PointLine)&    Line,
2389   Handle(IntSurf_LineOn2S)& Result,
2390   Standard_Boolean          theIsReversed,
2391   Handle(IntSurf_LineOn2S)& Vertices,
2392   const Standard_Real       ArcTol)
2393 {
2394   Standard_Integer nbp = Result->NbPoints(), nbv = Vertices->NbPoints();
2395
2396   if(nbp < 3)
2397     return;
2398
2399   const Handle(IntPatch_RLine) aRLine = Handle(IntPatch_RLine)::DownCast(Line);
2400
2401   Standard_Integer ip = 0, iv = 0;
2402   gp_Pnt aPnt;
2403   IntPatch_Point thePnt;
2404   Standard_Real U1 = 0., V1 = 0., U2 = 0., V2 = 0.;
2405
2406   for(ip = 2; ip <= (nbp-1); ip++) {
2407     const IntSurf_PntOn2S& aP = Result->Value(ip);
2408     for(iv = 1; iv <= nbv; iv++) {
2409       const IntSurf_PntOn2S& aV = Vertices->Value(iv);
2410       if(aP.IsSame(aV, Precision::Confusion(), Precision::PConfusion())) {
2411         aPnt = Result->Value(ip).Value();
2412         Result->Value(ip).ParametersOnS1(U1,V1);
2413         Result->Value(ip).ParametersOnS2(U2,V2);
2414         thePnt.SetValue(aPnt,ArcTol,Standard_False);
2415         thePnt.SetParameters(U1,V1,U2,V2);
2416         
2417         Standard_Real aParam = (Standard_Real)ip;
2418
2419         if(!aRLine.IsNull())
2420         {
2421           //In fact, aRLine is always on the parametric surface.
2422           //If (theIsReversed == TRUE) then (U1, V1) - point on
2423           //parametric surface, otherwise - point on quadric.
2424           const Handle(Adaptor2d_HCurve2d)& anArc = aRLine->IsArcOnS1() ?
2425                                                     aRLine->ArcOnS1() :
2426                                                     aRLine->ArcOnS2();
2427
2428           const gp_Lin2d aLin(anArc->Curve2d().Line());
2429           gp_Pnt2d aPSurf;
2430
2431           if(theIsReversed)
2432           {
2433             aPSurf.SetCoord(U1, V1);
2434           }
2435           else
2436           {
2437             aPSurf.SetCoord(U2, V2);
2438           }
2439
2440           aParam = ElCLib::Parameter(aLin, aPSurf);
2441         }
2442         
2443         thePnt.SetParameter(aParam);
2444         Line->AddVertex(thePnt);
2445       }
2446     }
2447   }
2448 }
2449
2450 static Standard_Boolean HasInternals(Handle(IntSurf_LineOn2S)& Line,
2451   Handle(IntSurf_LineOn2S)& Vertices)
2452 {
2453   Standard_Integer nbp = Line->NbPoints(), nbv = Vertices->NbPoints();
2454   Standard_Integer ip = 0, iv = 0;
2455   Standard_Boolean result = Standard_False;
2456
2457   if(nbp < 3)
2458     return result;
2459
2460   for(ip = 2; ip <= (nbp-1); ip++) {
2461     const IntSurf_PntOn2S& aP = Line->Value(ip);
2462     for(iv = 1; iv <= nbv; iv++) {
2463       const IntSurf_PntOn2S& aV = Vertices->Value(iv);
2464       if(aP.IsSame(aV, Precision::Confusion(), Precision::PConfusion())) {
2465         result = Standard_True;
2466         break;
2467       }
2468     }
2469     if(result)
2470       break;
2471   }
2472
2473   return result;
2474 }
2475 static Handle(IntPatch_WLine) MakeSplitWLine (Handle(IntPatch_WLine)&        WLine,
2476   Standard_Boolean         Tang,
2477   IntSurf_TypeTrans        Trans1,
2478   IntSurf_TypeTrans        Trans2,
2479   Standard_Real            ArcTol,
2480   Standard_Integer         ParFirst,
2481   Standard_Integer         ParLast)
2482 {
2483   Handle(IntSurf_LineOn2S) SLine = WLine->Curve();
2484   Handle(IntSurf_LineOn2S) sline = new IntSurf_LineOn2S();
2485
2486   Standard_Integer ip = 0;
2487   for(ip = ParFirst; ip <= ParLast; ip++)
2488     sline->Add(SLine->Value(ip));
2489
2490   Handle(IntPatch_WLine) wline = new IntPatch_WLine(sline,Tang,Trans1,Trans2);
2491
2492   gp_Pnt aSPnt;
2493   IntPatch_Point TPntF,TPntL;
2494   Standard_Real uu1 = 0., vv1 = 0., uu2 = 0., vv2 = 0.;
2495
2496   aSPnt = sline->Value(1).Value();
2497   sline->Value(1).ParametersOnS1(uu1,vv1);
2498   sline->Value(1).ParametersOnS2(uu2,vv2);
2499   TPntF.SetValue(aSPnt,ArcTol,Standard_False);
2500   TPntF.SetParameters(uu1,vv1,uu2,vv2);
2501   TPntF.SetParameter(1.);
2502   wline->AddVertex(TPntF);
2503   wline->SetFirstPoint(1);
2504
2505   aSPnt =  sline->Value(sline->NbPoints()).Value();
2506   sline->Value(sline->NbPoints()).ParametersOnS1(uu1,vv1);
2507   sline->Value(sline->NbPoints()).ParametersOnS2(uu2,vv2);
2508   TPntL.SetValue(aSPnt,ArcTol,Standard_False);
2509   TPntL.SetParameters(uu1,vv1,uu2,vv2);
2510   TPntL.SetParameter((Standard_Real)sline->NbPoints());
2511   wline->AddVertex(TPntL);
2512   wline->SetLastPoint(wline->NbVertex());
2513
2514   return wline;
2515 }
2516
2517 static Standard_Boolean SplitOnSegments(Handle(IntPatch_WLine)&        WLine,
2518   Standard_Boolean         Tang,
2519   IntSurf_TypeTrans        Trans1,
2520   IntSurf_TypeTrans        Trans2,
2521   Standard_Real            ArcTol,
2522   IntPatch_SequenceOfLine& Segments)
2523 {
2524   Standard_Boolean result = Standard_False;
2525   Segments.Clear();
2526
2527   Standard_Integer nbv = WLine->NbVertex();
2528   if(nbv > 3) {
2529     Standard_Integer iv = 0;
2530     for(iv = 1; iv < nbv; iv++) {
2531       Standard_Integer firstPar = 
2532                     (Standard_Integer) WLine->Vertex(iv).ParameterOnLine();
2533       Standard_Integer lastPar  = 
2534                     (Standard_Integer) WLine->Vertex(iv+1).ParameterOnLine();
2535       if((lastPar - firstPar) <= 1)
2536         continue;
2537       else {
2538         Handle(IntPatch_WLine) splitwline = MakeSplitWLine(WLine,Tang,Trans1,Trans2,
2539                                                                 ArcTol,firstPar,lastPar);
2540         Segments.Append(splitwline);
2541         if(!result)
2542           result = Standard_True;
2543       }
2544     }
2545   }
2546   return result;
2547 }
2548
2549 //=======================================================================
2550 //function : DecomposeResult
2551 //purpose  : Split <theLine> in the places where it passes through seam edge
2552 //            or singularity (apex of cone or pole of sphere).
2553 //            This passage is detected by jump of U-parameter
2554 //            from point to point.
2555 //=======================================================================
2556 static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLine,
2557                                         const Standard_Boolean       IsReversed,
2558                                         const IntSurf_Quadric&       theQuad,
2559                                         const Handle(Adaptor3d_TopolTool)& thePDomain,
2560                                         const Handle(Adaptor3d_HSurface)&  theQSurf, //quadric
2561                                         const Handle(Adaptor3d_HSurface)&  thePSurf, //parametric
2562                                         const Standard_Real                theArcTol,
2563                                         const Standard_Real                theTolTang,
2564                                         IntPatch_SequenceOfLine&           theLines)
2565 {
2566   if(theLine->ArcType() == IntPatch_Restriction)
2567   {
2568     const Handle(IntPatch_RLine)& aRL = Handle(IntPatch_RLine)::DownCast(theLine);
2569     if(!aRL.IsNull())
2570     {
2571       const Handle(Adaptor2d_HCurve2d)& anArc = aRL->IsArcOnS1() ?
2572                                         aRL->ArcOnS1() :
2573                                         aRL->ArcOnS2();
2574       if(anArc->Curve2d().GetType() != GeomAbs_Line)
2575       {
2576         //Restriction line must be isoline.
2577         //Other cases are not supported by
2578         //existing algorithms.
2579
2580         return Standard_False;
2581       }
2582     }
2583   }
2584   
2585   const Standard_Real aDeltaUmax = M_PI_2;
2586   const Standard_Real aTOL3D = 1.e-10, 
2587                       aTOL2D = Precision::PConfusion(),
2588                       aTOL2DS = Precision::PConfusion();
2589
2590   const Handle(IntSurf_LineOn2S)& aSLine = theLine->Curve();
2591
2592   if(aSLine->NbPoints() <= 2)
2593   {
2594     return Standard_False;
2595   }
2596   
2597   //Deletes repeated vertices
2598   Handle(IntSurf_LineOn2S) aVLine = GetVertices(theLine,aTOL3D,aTOL2D);
2599   
2600   Handle(IntSurf_LineOn2S) aSSLine(aSLine);
2601
2602   if(aSSLine->NbPoints() <= 1)
2603     return Standard_False;
2604
2605   AdjustLine(aSSLine,IsReversed,theQSurf,aTOL2D);
2606
2607   if(theLine->ArcType() == IntPatch_Walking)
2608   {
2609     Standard_Boolean isInserted = Standard_True;
2610     while(isInserted)
2611     {
2612       const Standard_Integer aNbPnts = aSSLine->NbPoints();
2613       TColStd_Array1OfInteger aPTypes(1,aNbPnts);
2614       SearchVertices(aSSLine,aVLine,aPTypes);
2615       isInserted = InsertSeamVertices(aSSLine,IsReversed,aVLine,aPTypes,aTOL2D);
2616     }
2617   }
2618
2619   const Standard_Integer aLindex = aSSLine->NbPoints();
2620   Standard_Integer aFindex = 1, aBindex = 0;
2621
2622   // build WLine parts (if any)
2623   Standard_Boolean flNextLine = Standard_True;
2624   Standard_Boolean hasBeenDecomposed = Standard_False;
2625   PrePoint_Type aPrePointExist = PrePoint_NONE;
2626
2627   IntSurf_PntOn2S PrePoint;
2628   while(flNextLine)
2629   {
2630     // reset variables
2631     flNextLine = Standard_False;
2632     Standard_Boolean isDecomposited = Standard_False;
2633     Standard_Real U1 = 0., U2 = 0., V1 = 0., V2 = 0.;
2634
2635     Handle(IntSurf_LineOn2S) sline = new IntSurf_LineOn2S();
2636
2637     //if((Lindex-Findex+1) <= 2 )
2638     if((aLindex <= aFindex) && !aPrePointExist)
2639     {
2640       //break of "while(flNextLine)" cycle
2641       break;
2642     }
2643
2644     if(aPrePointExist)
2645     {
2646       //The last point of the line is the pole of the quadric.
2647       //Therefore, Walking-line has been broken in this point.
2648       //However, new line must start from this point. Here we must
2649       //find its 2D-coordinates.
2650
2651       //For sphere and cone, some intersection point is satisfied to the system
2652       //  \cos(U_{q}) = S_{x}(U_{s},V_{s})/F(V_{q}) 
2653       //  \sin(U_{q}) = S_{y}(U_{s},V_{s})/F(V_{q}) 
2654
2655       //where 
2656       //  @S_{x}@, @S_{y}@ are X and Y-coordinates of thePSurf;
2657       //  @U_{s}@ and @V_{s}@ are UV-parameters on thePSurf;
2658       //  @U_{q}@ and @V_{q}@ are UV-parameters on theQSurf;
2659       //  @F(V_{q}) @ is some function, which value independs on @U_{q}@
2660       //              (form of this function depends on the type of the quadric).
2661
2662       //When we go through the pole, the function @F(V_{q}) @ changes sign.
2663       //Therefore, some cases are possible, when only @\cos(U_{q}) @ or
2664       //only @ \sin(U_{q}) @ change sign.
2665
2666       //Consequently, when the line goes throug the pole, @U_{q}@ can be
2667       //changed on @\pi /2 @ (but not less).
2668
2669       //Here, we forbid "jumping" between two neighbor Walking-point
2670       //with step greater than pi/4
2671       const Standard_Real aPeriod = M_PI_2, aHalfPeriod = M_PI_4;
2672       const IntSurf_PntOn2S& aRefPt = aSSLine->Value(aFindex);
2673
2674       const Standard_Real aURes = theQSurf->UResolution(theArcTol),
2675                           aVRes = theQSurf->UResolution(theArcTol);
2676
2677       const Standard_Real aTol2d = (aPrePointExist == PrePoint_POLE) ? 0.0 : 
2678               (aPrePointExist == PrePoint_SEAMV)? aVRes :
2679               (aPrePointExist == PrePoint_SEAMUV)? Max(aURes, aVRes) : aURes;
2680
2681       if(!PrePoint.IsSame(aRefPt, Precision::Confusion(), aTol2d))
2682       {
2683         Standard_Real aURef = 0.0, aVRef = 0.0;
2684         Standard_Real aUquad = 0.0, aVquad = 0.0;
2685
2686         //Take parameters on quadric
2687         if(IsReversed)
2688         {
2689           PrePoint.ParametersOnS2(aUquad, aVquad);
2690           aRefPt.ParametersOnS2(aURef, aVRef);
2691         }
2692         else
2693         {
2694           PrePoint.ParametersOnS1(aUquad, aVquad);
2695           aRefPt.ParametersOnS1(aURef, aVRef);
2696         }
2697
2698         if(theQSurf->IsUPeriodic())
2699         {
2700           Standard_Real aDeltaPar = aURef-aUquad;
2701           const Standard_Real anIncr = Sign(aPeriod, aDeltaPar);
2702           while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod))
2703           {
2704             aUquad += anIncr;
2705             aDeltaPar = aURef-aUquad;
2706           }
2707         }
2708
2709         if(theQSurf->IsVPeriodic())
2710         {
2711           Standard_Real aDeltaPar = aVRef-aVquad;
2712           const Standard_Real anIncr = Sign(aPeriod, aDeltaPar);
2713           while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod))
2714           {
2715             aVquad += anIncr;
2716             aDeltaPar = aVRef-aVquad;
2717           }
2718         }
2719
2720         PrePoint.SetValue(!IsReversed, aUquad, aVquad);
2721         sline->Add(PrePoint);
2722       }
2723       else
2724       {
2725         //break of "while(flNextLine)" cycle
2726         break;
2727       }
2728     }
2729
2730     aPrePointExist = PrePoint_NONE;
2731
2732     // analyze other points
2733     for(Standard_Integer k = aFindex; k <= aLindex; k++)
2734     {
2735       if( k == aFindex )
2736       {
2737         PrePoint = aSSLine->Value(k);
2738         sline->Add(PrePoint);
2739         continue;
2740       }
2741
2742       if(IsReversed)
2743       {
2744         aSSLine->Value(k).ParametersOnS2(U1,V1);   // S2 - quadric, set U,V by Pnt3D
2745       }
2746       else
2747       {
2748         aSSLine->Value(k).ParametersOnS1(U1,V1);    // S1 - quadric, set U,V by Pnt3D
2749       }
2750
2751       aPrePointExist = IsSeamOrPole(theQSurf, aSSLine, IsReversed, k-1, aDeltaUmax);
2752
2753       if(aPrePointExist != PrePoint_NONE)
2754       {
2755         aBindex = k;
2756         isDecomposited = Standard_True;
2757         ////
2758         const Standard_Real aPeriod = M_PI+M_PI, aHalfPeriod = M_PI;
2759         const IntSurf_PntOn2S& aRefPt = aSSLine->Value(aBindex-1);
2760
2761         //Not quadric point
2762         Standard_Real aU0 = 0.0, aV0 = 0.0;
2763         //Quadric point
2764         Standard_Real aUQuadRef = 0.0, aVQuadRef = 0.0;
2765
2766         if(IsReversed)
2767         {
2768           aRefPt.Parameters(aU0, aV0, aUQuadRef, aVQuadRef);
2769         }
2770         else
2771         {
2772           aRefPt.Parameters(aUQuadRef, aVQuadRef, aU0, aV0);
2773         }
2774
2775         if(aPrePointExist == PrePoint_SEAMUV)
2776         {
2777           aPrePointExist = PrePoint_NONE;
2778
2779           gp_Pnt aPQuad;
2780           Standard_Real aUquad = 0.0;
2781           Standard_Real aVquad = 0.0; 
2782
2783           theQSurf->D0(aUquad, aVquad, aPQuad);
2784
2785           Extrema_GenLocateExtPS anExtr(aPQuad, thePSurf->Surface(), aU0, aV0,
2786                                         Precision::PConfusion(),
2787                                         Precision::PConfusion());
2788
2789           if(!anExtr.IsDone())
2790           {
2791             break;
2792           }
2793
2794           if(anExtr.SquareDistance() < theTolTang*theTolTang)
2795           {
2796             anExtr.Point().Parameter(aU0, aV0);
2797             gp_Pnt aP0(anExtr.Point().Value());
2798
2799             IntSurf_PntOn2S aNewPoint;
2800             aNewPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), IsReversed, aU0, aV0);
2801
2802             if(!aNewPoint.IsSame(aRefPt, Precision::Confusion()))
2803             {
2804               //Adjust found U-paramter to previous point of the Walking-line
2805               Standard_Real aDeltaPar = aUQuadRef-aUquad;
2806               const Standard_Real anIncrU = Sign(aPeriod, aDeltaPar);
2807               while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod))
2808               {
2809                 aUquad += anIncrU;
2810                 aDeltaPar = aUQuadRef-aUquad;
2811               }
2812
2813               //Adjust found V-paramter to previous point of the Walking-line
2814               aDeltaPar = aVQuadRef-aVquad;
2815               const Standard_Real anIncrV = Sign(aPeriod, aDeltaPar);
2816               while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod))
2817               {
2818                 aVquad += anIncrV;
2819                 aDeltaPar = aVQuadRef-aVquad;
2820               }
2821
2822               aNewPoint.SetValue(!IsReversed, aUquad, aVquad);
2823               
2824               sline->Add(aNewPoint);
2825               aPrePointExist = PrePoint_SEAMUV;
2826               PrePoint = aNewPoint;
2827             }
2828           }
2829         }
2830         else if(aPrePointExist == PrePoint_SEAMV)
2831         {//WLine goes through seam
2832           aPrePointExist = PrePoint_NONE;
2833
2834           FuncPreciseSeam aF(theQSurf, thePSurf, Standard_False);
2835           math_Vector aTol(1, 3), aStartPoint(1,3),
2836                       anInfBound(1, 3), aSupBound(1, 3);
2837           
2838           //Parameters on parametric surface
2839           Standard_Real aUp = 0.0, aVp = 0.0;
2840           if(IsReversed)
2841           {
2842             aSSLine->Value(k).ParametersOnS1(aUp, aVp);
2843           }
2844           else
2845           {
2846             aSSLine->Value(k).ParametersOnS2(aUp, aVp);
2847           }
2848
2849           aTol(1) = thePSurf->UResolution(theArcTol);
2850           aTol(2) = thePSurf->VResolution(theArcTol);
2851           aTol(3) = theQSurf->UResolution(theArcTol);
2852           aStartPoint(1) = 0.5*(aU0 + aUp);
2853           aStartPoint(2) = 0.5*(aV0 + aVp);
2854           aStartPoint(3) = 0.5*(aUQuadRef + U1);
2855           anInfBound(1) = thePSurf->FirstUParameter();
2856           anInfBound(2) = thePSurf->FirstVParameter();
2857           anInfBound(3) = theQSurf->FirstUParameter();
2858           aSupBound(1) = thePSurf->LastUParameter();
2859           aSupBound(2) = thePSurf->LastVParameter();
2860           aSupBound(3) = theQSurf->LastUParameter();
2861
2862           math_FunctionSetRoot aSRF(aF, aTol);
2863           aSRF.Perform(aF, aStartPoint, anInfBound, aSupBound);
2864
2865           if(!aSRF.IsDone())
2866           {
2867             break;
2868           }
2869
2870           // Now aStartPoint is useless. Therefore, we use it for keeping
2871           // new point.
2872           aSRF.Root(aStartPoint);
2873           
2874           //On parametric
2875           aU0 = aStartPoint(1);
2876           aV0 = aStartPoint(2);
2877
2878           //On quadric
2879           Standard_Real aUquad = aStartPoint(3);
2880           Standard_Real aVquad = 0.0; 
2881           const gp_Pnt aPQuad(theQSurf->Value(aUquad, aVquad));
2882           const gp_Pnt aP0(thePSurf->Value(aU0, aV0));
2883
2884           {
2885             //Adjust found U-paramter to previous point of the Walking-line
2886             Standard_Real aDeltaPar = aVQuadRef-aVquad;
2887             const Standard_Real anIncr = Sign(aPeriod, aDeltaPar);
2888             while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod))
2889             {
2890               aVquad += anIncr;
2891               aDeltaPar = aVQuadRef-aVquad;
2892             }
2893           }
2894
2895           IntSurf_PntOn2S aNewPoint;
2896           if(IsReversed)
2897             aNewPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aU0, aV0, aUquad, aVquad);
2898           else
2899             aNewPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0);
2900
2901           if(!aNewPoint.IsSame(aRefPt, Precision::Confusion(), Precision::PConfusion()))
2902           {
2903             aNewPoint.SetValue(!IsReversed, aUquad, aVquad);
2904             sline->Add(aNewPoint);
2905             aPrePointExist = PrePoint_SEAMV;
2906             PrePoint = aNewPoint;
2907           }
2908           else
2909           {
2910             if(sline->NbPoints() == 1)
2911             {
2912               //FIRST point of the sline is the pole of the quadric.
2913               //Therefore, there is no point in decomposition.
2914
2915               PrePoint = aRefPt;
2916               aPrePointExist = PrePoint_SEAMV;
2917             }
2918           }
2919         }
2920         else if(aPrePointExist == PrePoint_POLESEAMU)
2921         {//Check if WLine goes through pole
2922           
2923           aPrePointExist = PrePoint_NONE;
2924
2925           //aPQuad is Pole
2926           gp_Pnt aPQuad;
2927           Standard_Real aUquad = 0.0;
2928           Standard_Real aVquad = 0.0; 
2929           
2930           if(theQuad.TypeQuadric() == GeomAbs_Sphere)
2931           {
2932             aVquad = Sign(M_PI_2, aVQuadRef);
2933           }
2934           else if(theQuad.TypeQuadric() == GeomAbs_Cone)
2935           {
2936             const Standard_Real aRadius = theQuad.Cone().RefRadius();
2937             const Standard_Real aSemiAngle = theQuad.Cone().SemiAngle();
2938             aVquad = -aRadius/sin(aSemiAngle);
2939           }
2940           else
2941           {
2942             Standard_TypeMismatch::Raise( "IntPatch_ImpPrmIntersection.cxx,"
2943                                           " DecomposeResult(...): "
2944                                           "Unsupported quadric with Pole");
2945           }
2946           
2947           theQSurf->D0(aUquad, aVquad, aPQuad);
2948
2949           Extrema_GenLocateExtPS anExtr(aPQuad, thePSurf->Surface(), aU0, aV0,
2950                                         Precision::PConfusion(),
2951                                         Precision::PConfusion());
2952
2953           if(!anExtr.IsDone())
2954           {
2955             break;
2956           }
2957
2958           if(anExtr.SquareDistance() < theTolTang*theTolTang)
2959           { //Pole is an intersection point
2960             //(lies in the quadric and the parametric surface)
2961
2962             anExtr.Point().Parameter(aU0, aV0);
2963             gp_Pnt aP0(anExtr.Point().Value());
2964
2965             IntSurf_PntOn2S aNewPoint;
2966             aNewPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), IsReversed, aU0, aV0);
2967
2968             if(!aNewPoint.IsSame(aRefPt, Precision::Confusion()))
2969             {
2970               //Found pole does not exist in the Walking-line
2971               //It must be added there (with correct 2D-parameters)
2972               
2973               //2D-parameters of theparametric surface have already been found (aU0, aV0).
2974               //Let find 2D-parameters on the quadric.
2975
2976               //The algorithm depends on the type of the quadric. Here we consider a Sphere only.
2977               //Analogical result can be made for another types (e.g. cone, but formulas will
2978               //be different) in case of need.
2979               
2980               //First of all, we need in adjusting thePSurf in the coordinate system of the Sphere
2981               //(in order to make the equation of the sphere maximal simple). However, as it will be
2982               //shown later, thePSurf is used in algorithm in order to get its derivatives. Therefore,
2983               //for improving performance, transformation of these vectors is enough (there is no point
2984               //in transformation of full surface).
2985               
2986               gp_Pnt aPtemp;
2987               gp_Vec aVecDu, aVecDv;
2988               thePSurf->D1(aU0, aV0, aPtemp, aVecDu, aVecDv);
2989
2990               //Transforms parametric surface in coordinate-system of the quadric
2991               gp_Trsf aTr;
2992               aTr.SetTransformation((theQuad.TypeQuadric() == GeomAbs_Sphere) ?
2993                                       theQuad.Sphere().Position() :
2994                                       theQuad.Cone().Position());
2995
2996               //Derivatives of transformed thePSurf
2997               aVecDu.Transform(aTr);
2998               aVecDv.Transform(aTr);
2999
3000               if(theQuad.TypeQuadric() == GeomAbs_Sphere)
3001               {
3002                 //The intersection point (including the pole)
3003                 //must be satisfied to the following system:
3004
3005                 //    \left\{\begin{matrix}
3006                 //    R*\cos (U_{q})*\cos (V_{q})=S_{x}(U_{s},V_{s})
3007                 //    R*\sin (U_{q})*\cos (V_{q})=S_{y}(U_{s},V_{s})
3008                 //    R*\sin (V_{q})=S_{z}(U_{s},V_{s})
3009                 //    \end{matrix}\right,
3010                 //where 
3011                 //  R is the radius of the sphere;
3012                 //  @S_{x}@, @S_{y}@ and @S_{z}@ are X, Y and Z-coordinates of thePSurf;
3013                 //  @U_{s}@ and @V_{s}@ are equal to aU0 and aV0 corespondingly;
3014                 //  @U_{q}@ and @V_{q}@ are equal to aUquad and aVquad corespondingly.
3015
3016                 //Consequently (from first two equations), 
3017                 //  \left\{\begin{matrix}
3018                 //  \cos (U_{q}) = \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})}
3019                 //  \sin (U_{q}) = \frac{S_{y}(U_{s},V_{s})}{R*\cos (V_{q})}
3020                 //  \end{matrix}\right.
3021
3022                 //For pole, 
3023                 //  V_{q}=\pm \pi /2 \Rightarrow \cos (V_{q}) = 0 (denominator is equal to 0).
3024
3025                 //Therefore, computation U_{q} directly is impossibly.
3026                 //
3027                 //Let @V_{q}@ tends to @\pm \pi /2@.
3028                 //Then (indeterminate form is evaluated in accordance of L'Hospital rule),
3029                 //  \cos (U_{q}) = \lim_{V_{q} \to (\pi /2-0)} 
3030                 //  \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})}= 
3031                 //  -\lim_{V_{q} \to (\pi /2-0)}
3032                 //  \frac{\frac{\partial S_{x}}
3033                 //  {\partial U_{s}}*\frac{\mathrm{d} U_{s}} 
3034                 //  {\mathrm{d} V_{q}}+\frac{\partial S_{x}} 
3035                 //  {\partial V_{s}}*\frac{\mathrm{d} V_{s}} 
3036                 //  {\mathrm{d} V_{q}}}{R*\sin (V_{q})} =  
3037                 //  -\frac{1}{R}*\frac{\mathrm{d} U_{s}}
3038                 //  {\mathrm{d} V_{q}}*(\frac{\partial S_{x}} 
3039                 //  {\partial U_{s}}+\frac{\partial S_{x}}
3040                 //  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}
3041                 //  {\mathrm{d} U_{s}}) =
3042                 //  -\frac{1}{R}*\frac{\mathrm{d} V_{s}}
3043                 //  {\mathrm{d} V_{q}}*(\frac{\partial S_{x}} 
3044                 //  {\partial U_{s}}*\frac{\mathrm{d} U_{s}}
3045                 //  {\mathrm{d} V_{s}}+\frac{\partial S_{x}}
3046                 //  {\partial V_{s}}).
3047
3048                 //Analogicaly for @\sin (U_{q})@ (@S_{x}@ is substituted to @S_{y}@).
3049
3050                 //Let mean, that
3051                 //  \cos (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \cos (U_{q}) \left | _{V_{q} \to (\pi /2-0)}
3052                 //  \sin (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \sin (U_{q}) \left | _{V_{q} \to (\pi /2-0)}
3053
3054                 //From the 3rd equation of the system, we obtain
3055                 //  \frac{\mathrm{d} (R*\sin (V_{q}))}{\mathrm{d} V_{q}} =
3056                 //  \frac{\mathrm{d} S_{z}(U_{s},V_{s})}{\mathrm{d} V_{q}}
3057                 //or
3058                 //  R*\cos (V_{q}) = \frac{\partial S_{z}}{\partial U_{s}}*
3059                 //  \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}}
3060                 //  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}}.
3061
3062                 //If @V_{q}=\pm \pi /2@, then
3063                 //  \frac{\partial S_{z}}{\partial U_{s}}*
3064                 //  \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}}
3065                 //  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}} = 0.
3066
3067                 //Consequently, if @\frac{\partial S_{z}}{\partial U_{s}} \neq 0 @ then
3068                 //  \frac{\mathrm{d} U_{s}}{\mathrm{d} V_{s}} =
3069                 //  -\frac{\frac{\partial S_{z}}{\partial V_{s}}}
3070                 //  {\frac{\partial S_{z}}{\partial U_{s}}}.
3071
3072                 //If @ \frac{\partial S_{z}}{\partial V_{s}} \neq 0 @ then
3073                 //  \frac{\mathrm{d} V_{s}}{\mathrm{d} U_{s}} =
3074                 //  -\frac{\frac{\partial S_{z}}{\partial U_{s}}}
3075                 //  {\frac{\partial S_{z}}{\partial V_{s}}}
3076
3077                 //Cases, when @ \frac{\partial S_{z}}{\partial U_{s}} = 
3078                 //\frac{\partial S_{z}}{\partial V_{s}} = 0 @ are not consider here.
3079                 //The reason is written below.
3080
3081                 //Vector with {@ \cos (U_{q}) @, @ \sin (U_{q}) @} coordinates.
3082                 //Ask to pay attention to the fact that this vector is always normalyzed.
3083                 gp_Vec2d aV1;
3084
3085                 if( (Abs(aVecDu.Z()) < Precision::PConfusion()) &&
3086                   (Abs(aVecDv.Z()) < Precision::PConfusion()))
3087                 {
3088                   //Example of this exception is intersection a plane with a sphere
3089                   //when the plane tangents the sphere in some pole (i.e. only one 
3090                   //intersection point, not line). In this case, U-coordinate of the
3091                   //sphere is undefined (can be realy anything).
3092                   //Another reason is that we have tangent zone around the pole
3093                   //(see bug #26576).
3094                   //Computation correct value of aUquad is impossible. Therefore,
3095                   //we should throw an exception in this case.
3096                   //Also, any Walking line cannot be created in this case.
3097                   //Hovewer, Restriction line is not created by intersection algorithm.
3098                   //It is already exists (above we check simply, if this line is
3099                   //intersection line).
3100                   //Therefore, we can try to find the aUquad-parameter on (existing)
3101                   //Restriction line. Here, we will do it with
3102                   //extrapolation algorithm.
3103                   //Use interpolation algorithm is wrong because aUquad parameter
3104                   //jumps while the line going though the pole.
3105
3106                   if((theLine->ArcType() == IntPatch_Walking) ||
3107                     (aBindex < 3))
3108                   {
3109                     //We must have at least two previous points
3110                     //in order to do linear extrapolation.
3111                     Standard_NumericError::
3112                       Raise("IntPatch_ImpPrmIntersection.cxx, DecomposeResult(...): "
3113                       "Cannot find UV-coordinate for quadric in the pole");
3114                   }
3115                   else
3116                   {
3117 #ifdef INTPATCH_IMPPRMINTERSECTION_DEBUG
3118                     cout << "Cannot find UV-coordinate for quadric in the pole."
3119                       " See considered comment above. IntPatch_ImpPrmIntersection.cxx,"
3120                       " DecomposeResult(...)" << endl;
3121 #endif
3122
3123                     //    *----------*------------x
3124                     // QuadPrev   QuadRef     Quad (must be found)
3125
3126                     const IntSurf_PntOn2S& aPt2S = aSSLine->Value(aBindex-2);
3127                     //Quadric point
3128                     Standard_Real aUQuadPrev = 0.0, aVQuadPrev = 0.0;
3129                     if(IsReversed)
3130                     {
3131                       aPt2S.ParametersOnS2(aUQuadPrev, aVQuadPrev);
3132                     }
3133                     else
3134                     {
3135                       aPt2S.ParametersOnS1(aUQuadPrev, aVQuadPrev);
3136                     }
3137
3138                     Standard_NumericError_Raise_if(
3139                       Abs(aVQuadPrev - aVQuadRef) < gp::Resolution(),
3140                       "Division by zero");
3141
3142                     aUquad = 
3143                       aUQuadPrev + (aUQuadRef - aUQuadPrev)*
3144                       (aVquad - aVQuadPrev)/(aVQuadRef - aVQuadPrev);
3145                   }
3146                 }
3147                 else
3148                 {
3149                   if(Abs(aVecDu.Z()) > Abs(aVecDv.Z()))
3150                   {
3151                     const Standard_Real aDusDvs = aVecDv.Z()/aVecDu.Z();
3152
3153                     aV1.SetCoord( aVecDu.X()*aDusDvs - aVecDv.X(),
3154                       aVecDu.Y()*aDusDvs - aVecDv.Y());
3155                   }
3156                   else
3157                   {
3158                     const Standard_Real aDvsDus = aVecDu.Z()/aVecDv.Z();
3159                     aV1.SetCoord( aVecDv.X()*aDvsDus - aVecDu.X(),
3160                       aVecDv.Y()*aDvsDus - aVecDu.Y());
3161                   }
3162
3163                   aV1.Normalize();
3164
3165                   if(Abs(aV1.X()) > Abs(aV1.Y()))
3166                     aUquad = Sign(asin(aV1.Y()), aVquad);
3167                   else
3168                     aUquad = Sign(acos(aV1.X()), aVquad);
3169                 }
3170
3171                 {
3172                   //Adjust found U-paramter to previous point of the Walking-line
3173                   Standard_Real aDeltaPar = aUQuadRef-aUquad;
3174                   const Standard_Real anIncr = Sign(aPeriod, aDeltaPar);
3175                   while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod))
3176                   {
3177                     aUquad += anIncr;
3178                     aDeltaPar = aUQuadRef-aUquad;
3179                   }
3180                 }
3181               }
3182
3183               aNewPoint.SetValue(!IsReversed, aUquad, aVquad);
3184               
3185               sline->Add(aNewPoint);
3186               aPrePointExist = PrePoint_POLE;
3187               PrePoint = aNewPoint;
3188             } // if(!aNewPoint.IsSame(aRefPt, Precision::Confusion()))
3189             else
3190             {
3191               aPrePointExist = PrePoint_NONE;
3192
3193               if(sline->NbPoints() == 1)
3194               {
3195                 //FIRST point of the sline is the pole of the quadric.
3196                 //Therefore, there is no point in decomposition.
3197
3198                 PrePoint = aRefPt;
3199                 aPrePointExist = PrePoint_POLE;
3200               }
3201             }
3202           } //if(anExtr.SquareDistance() < aTol*aTol)
3203           else
3204           {//Pole is not an intersection point
3205             aPrePointExist = PrePoint_SEAMU;
3206           }
3207         }
3208
3209         if(aPrePointExist == PrePoint_SEAMU)
3210         {//WLine goes through seam
3211
3212           aPrePointExist = PrePoint_NONE;
3213
3214           FuncPreciseSeam aF(theQSurf, thePSurf, Standard_True);
3215           math_Vector aTol(1, 3), aStartPoint(1,3),
3216                       anInfBound(1, 3), aSupBound(1, 3);
3217           
3218           //Parameters on parametric surface
3219           Standard_Real aUp = 0.0, aVp = 0.0;
3220           if(IsReversed)
3221           {
3222             aSSLine->Value(k).ParametersOnS1(aUp, aVp);
3223           }
3224           else
3225           {
3226             aSSLine->Value(k).ParametersOnS2(aUp, aVp);
3227           }
3228
3229           aTol(1) = thePSurf->UResolution(theArcTol);
3230           aTol(2) = thePSurf->VResolution(theArcTol);
3231           aTol(3) = theQSurf->VResolution(theArcTol);
3232           aStartPoint(1) = 0.5*(aU0 + aUp);
3233           aStartPoint(2) = 0.5*(aV0 + aVp);
3234           aStartPoint(3) = 0.5*(aVQuadRef + V1);
3235           anInfBound(1) = thePSurf->FirstUParameter();
3236           anInfBound(2) = thePSurf->FirstVParameter();
3237           anInfBound(3) = theQSurf->FirstVParameter();
3238           aSupBound(1) = thePSurf->LastUParameter();
3239           aSupBound(2) = thePSurf->LastVParameter();
3240           aSupBound(3) = theQSurf->LastVParameter();
3241
3242           math_FunctionSetRoot aSRF(aF, aTol);
3243           aSRF.Perform(aF, aStartPoint, anInfBound, aSupBound);
3244
3245           if(!aSRF.IsDone())
3246           {
3247             break;
3248           }
3249
3250           // Now aStartPoint is useless. Therefore, we use it for keeping
3251           // new point.
3252           aSRF.Root(aStartPoint);
3253           
3254           //On parametric
3255           aU0 = aStartPoint(1);
3256           aV0 = aStartPoint(2);
3257
3258           //On quadric
3259           Standard_Real aUquad = 0.0;
3260           Standard_Real aVquad = aStartPoint(3); 
3261           const gp_Pnt aPQuad(theQSurf->Value(aUquad, aVquad));
3262           const gp_Pnt aP0(thePSurf->Value(aU0, aV0));
3263
3264           {
3265             //Adjust found U-paramter to previous point of the Walking-line
3266             Standard_Real aDeltaPar = aUQuadRef-aUquad;
3267             const Standard_Real anIncr = Sign(aPeriod, aDeltaPar);
3268             while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod))
3269             {
3270               aUquad += anIncr;
3271               aDeltaPar = aUQuadRef-aUquad;
3272             }
3273           }
3274
3275           IntSurf_PntOn2S aNewPoint;
3276           if(IsReversed)
3277             aNewPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aU0, aV0, aUquad, aVquad);
3278           else
3279             aNewPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0);
3280
3281           if(!aNewPoint.IsSame(aRefPt, Precision::Confusion(), Precision::PConfusion()))