0028956: Coding Rules - eliminate several trivial compiler warnings
[occt.git] / src / IntCurveSurface / IntCurveSurface_Inter.gxx
1 // Created on: 1993-04-09
2 // Created by: Laurent BUCHARD
3 // Copyright (c) 1993-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 //#ifndef OCCT_DEBUG
18 //#define No_Standard_RangeError
19 //#define No_Standard_OutOfRange
20 //#endif
21
22
23 #define  TOLTANGENCY         0.00000001
24 #define  TOLERANCE_ANGULAIRE 1.e-12//0.00000001
25 #define  TOLERANCE           0.00000001
26
27 #define NBSAMPLESONCIRCLE  32
28 #define NBSAMPLESONELLIPSE 32
29 #define NBSAMPLESONPARAB   16
30 #define NBSAMPLESONHYPR    32
31
32 #include <ElSLib.hxx>
33 #include <Intf_SectionPoint.hxx>
34 #include <Intf_TangentZone.hxx>
35 #include <Intf_Tool.hxx>
36 #include <math_FunctionSetRoot.hxx>
37 #include <IntCurveSurface_IntersectionPoint.hxx>
38 #include <IntCurveSurface_IntersectionSegment.hxx>
39 #include <IntAna_IntConicQuad.hxx>
40 #include <IntAna_Quadric.hxx>
41 #include <gp_Vec.hxx>
42 #include <gp_Dir.hxx>
43 #include <Precision.hxx>
44 #include <IntAna_IntLinTorus.hxx>
45
46 #include <GeomAbs_Shape.hxx>
47 #include <GeomAbs_CurveType.hxx>
48 #include <TColStd_Array1OfReal.hxx>
49 #include <TColgp_Array1OfPnt.hxx>
50 #include <gp_Pnt.hxx>
51 #include <gp_Pln.hxx>
52 #include <gp_Lin.hxx>
53 #include <GProp_PEquation.hxx>
54 #include <GProp_PGProps.hxx>
55 #include <GProp_PrincipalProps.hxx>
56 #include <Extrema_ExtElC.hxx>
57 #include <Extrema_POnCurv.hxx>
58
59
60
61 #include <ProjLib_Plane.hxx>
62 #include <IntAna2d_AnaIntersection.hxx>
63 #include <gp_Lin2d.hxx>
64 #include <IntAna2d_IntPoint.hxx>
65 #include <gp_Parab2d.hxx>
66 #include <IntAna2d_Conic.hxx>
67 #include <IntAna2d_AnaIntersection.hxx>
68 #include <gp_Hypr2d.hxx>
69 #include <Adaptor3d_HCurve.hxx>
70 #include <Adaptor3d_HSurface.hxx>
71
72
73 #include <TColgp_Array2OfPnt.hxx>
74 #include <TColStd_HArray1OfReal.hxx>
75 #include <Adaptor3d_TopolTool.hxx>
76 #include <ElCLib.hxx>
77 //#define ICS_DEB
78 static 
79   void EstLimForInfExtr(const gp_Lin&   Line,
80                         const TheSurface& surface, 
81                         const Standard_Boolean IsOffSurf,
82                         const Standard_Integer nbsu, 
83                         const Standard_Boolean U1inf, 
84                         const Standard_Boolean U2inf, 
85                         const Standard_Boolean V1inf, 
86                         const Standard_Boolean V2inf, 
87                         Standard_Real& U1new, 
88                         Standard_Real& U2new, 
89                         Standard_Real& V1new, 
90                         Standard_Real& V2new, 
91                         Standard_Boolean& NoIntersection);
92
93 static 
94   void EstLimForInfRevl(const gp_Lin&   Line,
95                         const TheSurface& surface, 
96                         const Standard_Boolean U1inf, 
97                         const Standard_Boolean U2inf, 
98                         const Standard_Boolean V1inf, 
99                         const Standard_Boolean V2inf, 
100                         Standard_Real& U1new, 
101                         Standard_Real& U2new, 
102                         Standard_Real& V1new, 
103                         Standard_Real& V2new, 
104                         Standard_Boolean& NoIntersection);
105
106 static 
107   void EstLimForInfOffs(const gp_Lin&   Line,
108                         const TheSurface& surface, 
109                         const Standard_Integer nbsu, 
110                         const Standard_Boolean U1inf, 
111                         const Standard_Boolean U2inf, 
112                         const Standard_Boolean V1inf, 
113                         const Standard_Boolean V2inf, 
114                         Standard_Real& U1new, 
115                         Standard_Real& U2new, 
116                         Standard_Real& V1new, 
117                         Standard_Real& V2new, 
118                         Standard_Boolean& NoIntersection);
119
120 static
121   void EstLimForInfSurf(Standard_Real& U1new, 
122                         Standard_Real& U2new, 
123                         Standard_Real& V1new, 
124                         Standard_Real& V2new);
125
126 static
127   void SectionPointToParameters(const Intf_SectionPoint& Sp,
128                                 const IntCurveSurface_ThePolyhedron& Surf,
129                                 const IntCurveSurface_ThePolygon&    Curv,
130                                 Standard_Real& u,
131                                 Standard_Real& v,
132                                 Standard_Real& w);
133
134 static 
135   void IntCurveSurface_ComputeTransitions(const TheCurve& curve,
136                                           const Standard_Real w,
137                                           IntCurveSurface_TransitionOnCurve&   TransOnCurve,
138                                           const TheSurface& surface,
139                                           const Standard_Real u,
140                                           const Standard_Real v);
141
142 static 
143   void IntCurveSurface_ComputeParamsOnQuadric(const TheSurface& surface,
144                                               const gp_Pnt& P,
145                                               Standard_Real& u,
146                                               Standard_Real& v);
147
148 static 
149   void ProjectIntersectAndEstLim(const gp_Lin&        theLine,
150                                  const gp_Pln&        thePln,
151                                  const ProjLib_Plane& theBasCurvProj,
152                                  Standard_Real&       theVmin,
153                                  Standard_Real&       theVmax,
154                                  Standard_Boolean&    theNoIntersection);
155
156 //=======================================================================
157 //function : IntCurveSurface_Inter
158 //purpose  : 
159 //=======================================================================
160 IntCurveSurface_Inter::IntCurveSurface_Inter() 
161
162 }
163 //=======================================================================
164 //function : DoSurface
165 //purpose  : 
166 //=======================================================================
167 void IntCurveSurface_Inter::DoSurface(const TheSurface&   surface,
168                                       const Standard_Real u0,
169                                       const Standard_Real u1,
170                                       const Standard_Real v0,
171                                       const Standard_Real v1,
172                                       TColgp_Array2OfPnt& pntsOnSurface,
173                                       Bnd_Box&            boxSurface,
174                                       Standard_Real&      gap)
175 {
176   Standard_Integer iU = 0, iV = 0;
177   Standard_Real U = 0., V = 0;
178 // modified by NIZHNY-MKK  Mon Oct  3 17:38:45 2005
179 //   Standard_Real dU = fabs(u1-u0)/50., dV = fabs(v1-v0)/50.;
180   Standard_Real dU = (u1-u0)/50., dV = (v1-v0)/50.;
181   gp_Pnt aPnt;
182
183   for(iU = 0; iU < 50; iU++) {
184     
185     if(iU == 0)
186       U = u0;
187     else if(iU == 49)
188       U = u1;
189     else 
190       U = u0 + dU * ((Standard_Real)iU);
191     
192     for(iV = 0; iV < 50; iV++) {
193
194       if(iV == 0)
195         V = v0;
196       else if(iV == 49)
197         V = v1;
198       else
199         V = v0 + dV * ((Standard_Real)iV);
200
201       TheSurfaceTool::D0(surface,U,V,aPnt);
202       boxSurface.Add(aPnt);
203       pntsOnSurface.SetValue(iU+1,iV+1,aPnt);
204     }
205   }
206   Standard_Real Ures = TheSurfaceTool::UResolution(surface,dU);
207   Standard_Real Vres = TheSurfaceTool::VResolution(surface,dV);
208   gap = Max(Ures,Vres);
209 }
210
211 //=======================================================================
212 //function : DoNewBounds
213 //purpose  : 
214 //=======================================================================
215 void IntCurveSurface_Inter::DoNewBounds(
216                                         const TheSurface&           surface,
217                                         const Standard_Real         u0,
218                                         const Standard_Real         u1,
219                                         const Standard_Real         v0,
220                                         const Standard_Real         v1,
221                                         const TColgp_Array2OfPnt&   pntsOnSurface,
222                                         const TColStd_Array1OfReal& X,
223                                         const TColStd_Array1OfReal& Y,
224                                         const TColStd_Array1OfReal& Z,
225                                         TColStd_Array1OfReal&       Bounds)
226 {
227   Bounds.SetValue(1,u0);
228   Bounds.SetValue(2,u1);
229   Bounds.SetValue(3,v0);
230   Bounds.SetValue(4,v1);
231
232   Standard_Boolean isUClosed = (TheSurfaceTool::IsUClosed(surface) || TheSurfaceTool::IsUPeriodic(surface));
233   Standard_Boolean isVClosed = (TheSurfaceTool::IsVClosed(surface) || TheSurfaceTool::IsVPeriodic(surface));
234   Standard_Boolean checkU = (isUClosed) ? Standard_False : Standard_True;
235   Standard_Boolean checkV = (isVClosed) ? Standard_False : Standard_True;
236
237   Standard_Integer i = 0, j = 0, k = 0, iU = 0, iV = 0;
238   Standard_Integer iUmin = 50, iVmin = 50, iUmax = 1, iVmax = 1;
239
240   for(i = 1; i <= 2; i++) {
241     for(j = 1; j <= 2; j++) {
242       for(k = 1; k <= 2; k++) {
243         gp_Pnt aPoint(X(i),Y(j),Z(k));
244         Standard_Real DistMin = 1.e+100;
245         Standard_Integer diU = 0, diV = 0;
246         for(iU = 1; iU <= 50; iU++) {
247           for(iV = 1; iV <= 50; iV++) {
248             const gp_Pnt aP = pntsOnSurface.Value(iU,iV);
249             Standard_Real dist = aP.SquareDistance(aPoint);
250             if(dist < DistMin) {
251               DistMin = dist;
252               diU = iU;
253               diV = iV;
254             }
255           }
256         }
257         if(diU > 0 && diU < iUmin) iUmin = diU;
258         if(diU > 0 && diU > iUmax) iUmax = diU;
259         if(diV > 0 && diV < iVmin) iVmin = diV;
260         if(diV > 0 && diV > iVmax) iVmax = diV;
261       }
262     }
263   }
264
265 // modified by NIZHNY-MKK  Mon Oct  3 17:38:43 2005
266 //   Standard_Real dU = fabs(u1-u0)/50., dV = fabs(v1-v0)/50.;
267   Standard_Real dU = (u1-u0)/50., dV = (v1-v0)/50.;
268
269   Standard_Real USmin = u0 + dU * ((Standard_Real)(iUmin - 1));
270   Standard_Real USmax = u0 + dU * ((Standard_Real)(iUmax - 1));
271   Standard_Real VSmin = v0 + dV * ((Standard_Real)(iVmin - 1));
272   Standard_Real VSmax = v0 + dV * ((Standard_Real)(iVmax - 1));
273
274   if(USmin > USmax) {
275     Standard_Real tmp = USmax;
276     USmax = USmin;
277     USmin = tmp;
278   }
279   if(VSmin > VSmax) {
280     Standard_Real tmp = VSmax;
281     VSmax = VSmin;
282     VSmin = tmp;
283   }
284
285   USmin -= 1.5*dU;
286   if(USmin < u0) USmin = u0;
287   USmax += 1.5*dU;
288   if(USmax > u1) USmax = u1;
289   VSmin -= 1.5*dV;
290   if(VSmin < v0) VSmin = v0;
291   VSmax += 1.5*dV;
292   if(VSmax > v1) VSmax = v1;
293
294   if(checkU) {
295     Bounds.SetValue(1,USmin);
296     Bounds.SetValue(2,USmax);
297   }
298   if(checkV) {
299     Bounds.SetValue(3,VSmin);
300     Bounds.SetValue(4,VSmax);
301   }
302 }
303
304 //=======================================================================
305 //function : Perform
306 //purpose  : Decompose la surface si besoin est 
307 //=======================================================================
308 void IntCurveSurface_Inter::Perform(const TheCurve&   curve,
309                                     const TheSurface& surface)
310 {
311   ResetFields();
312   done = Standard_True;
313   Standard_Integer NbUOnS = TheSurfaceTool::NbUIntervals(surface,GeomAbs_C2);
314   Standard_Integer NbVOnS = TheSurfaceTool::NbVIntervals(surface,GeomAbs_C2);
315   Standard_Real U0,U1,V0,V1; 
316
317   if(NbUOnS > 1)
318   {
319     TColStd_Array1OfReal TabU(1,NbUOnS+1);
320     TheSurfaceTool::UIntervals(surface,TabU,GeomAbs_C2);
321     for(Standard_Integer iu = 1;iu <= NbUOnS; iu++)
322     {
323       U0 = TabU.Value(iu);
324       U1 = TabU.Value(iu+1);
325       if(NbVOnS > 1)
326       {
327         TColStd_Array1OfReal TabV(1,NbVOnS+1);
328         TheSurfaceTool::VIntervals(surface,TabV,GeomAbs_C2);
329         for(Standard_Integer iv = 1;iv <= NbVOnS; iv++)
330         {
331           // More than one interval on U and V param space.
332           V0 = TabV.Value(iv);
333           V1 = TabV.Value(iv+1);
334           Perform(curve,surface,U0,V0,U1,V1);
335         }
336       }
337       else
338       {
339         // More than one interval only on U param space.
340         V0 = TheSurfaceTool::FirstVParameter(surface);
341         V1 = TheSurfaceTool::LastVParameter(surface);
342         Perform(curve,surface,U0,V0,U1,V1);
343       }
344     }
345   }
346   else if(NbVOnS > 1)
347   {
348     // More than one interval only on V param space.
349     U0 = TheSurfaceTool::FirstUParameter(surface);
350     U1 = TheSurfaceTool::LastUParameter(surface);   
351     TColStd_Array1OfReal TabV(1,NbVOnS+1);
352     TheSurfaceTool::VIntervals(surface,TabV,GeomAbs_C2);
353     for(Standard_Integer iv = 1;iv <= NbVOnS; iv++)
354     {
355       V0 = TabV.Value(iv);
356       V1 = TabV.Value(iv+1);
357       Perform(curve,surface,U0,V0,U1,V1);
358     }
359   }
360   else
361   {
362     // One interval on U and V param space.
363     V0 = TheSurfaceTool::FirstVParameter(surface);
364     V1 = TheSurfaceTool::LastVParameter(surface);
365     U0 = TheSurfaceTool::FirstUParameter(surface);
366     U1 = TheSurfaceTool::LastUParameter(surface);
367
368     Perform(curve,surface,U0,V0,U1,V1);
369   }
370 }
371 //=======================================================================
372 //function : Perform
373 //purpose  : 
374 //=======================================================================
375 void IntCurveSurface_Inter::Perform(const TheCurve&   curve,
376                                     const TheSurface& surface,
377                                     const Standard_Real U1,const Standard_Real V1,
378                                     const Standard_Real U2,const Standard_Real V2)
379 {
380   // Protection from double type overflow.
381   // This may happen inside square magnitude computation based on normal,
382   // which was computed on bound parameteres (bug26525).
383   Standard_Real UU1 = U1, UU2 = U2, VV1 = V1, VV2 = V2;
384   if (U1 < -1.0e50)
385     UU1 = -1.0e50;
386   if (U2 >  1.0e50)
387     UU2 = 1.0e50;
388   if (V1 < -1.0e50)
389     VV1 = -1.0e50;
390   if (V2 >  1.0e50)
391     VV2 = 1.0e50;
392
393   GeomAbs_CurveType CurveType = TheCurveTool::GetType(curve);
394   
395   switch(CurveType) { 
396   case GeomAbs_Line:  
397     { 
398       PerformConicSurf(TheCurveTool::Line(curve),curve,surface,UU1,VV1,UU2,VV2);
399       break;
400     }
401   case GeomAbs_Circle: 
402     { 
403       PerformConicSurf(TheCurveTool::Circle(curve),curve,surface,UU1,VV1,UU2,VV2);
404       break;
405     }
406   case GeomAbs_Ellipse: 
407     { 
408       PerformConicSurf(TheCurveTool::Ellipse(curve),curve,surface,UU1,VV1,UU2,VV2);
409       break;
410     }    
411   case GeomAbs_Parabola:  
412     {
413       PerformConicSurf(TheCurveTool::Parabola(curve),curve,surface,UU1,VV1,UU2,VV2);
414       break;
415     }
416   case GeomAbs_Hyperbola: 
417     {
418       PerformConicSurf(TheCurveTool::Hyperbola(curve),curve,surface,UU1,VV1,UU2,VV2);
419       break;
420     }
421   default:    
422     {
423       Standard_Integer nbIntervalsOnCurve = TheCurveTool::NbIntervals(curve,GeomAbs_C2);
424       GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
425       if(   (SurfaceType != GeomAbs_Plane) 
426          && (SurfaceType != GeomAbs_Cylinder)
427          && (SurfaceType != GeomAbs_Cone) 
428          && (SurfaceType != GeomAbs_Sphere)) { 
429         
430         if(nbIntervalsOnCurve > 1) { 
431           TColStd_Array1OfReal TabW(1,nbIntervalsOnCurve+1);
432           TheCurveTool::Intervals(curve,TabW,GeomAbs_C2);
433           for(Standard_Integer i = 1; i<=nbIntervalsOnCurve; i++) { 
434             Standard_Real u1,u2;
435             u1 = TabW.Value(i);
436             u2 = TabW.Value(i+1);
437
438             Handle(TColStd_HArray1OfReal) aPars;
439             Standard_Real defl = 0.1;
440             Standard_Integer NbMin = 10;
441             TheCurveTool::SamplePars(curve, u1, u2, defl, NbMin, aPars);
442
443 //          IntCurveSurface_ThePolygon  polygon(curve,u1,u2,TheCurveTool::NbSamples(curve,u1,u2));
444             IntCurveSurface_ThePolygon  polygon(curve, aPars->Array1());
445             InternalPerform(curve,polygon,surface,UU1,VV1,UU2,VV2);
446           }
447         }
448         else {
449           Standard_Real u1,u2;
450           u1 = TheCurveTool::FirstParameter(curve);
451           u2 = TheCurveTool::LastParameter(curve);
452
453           Handle(TColStd_HArray1OfReal) aPars;
454           Standard_Real defl = 0.1;
455           Standard_Integer NbMin = 10;
456           TheCurveTool::SamplePars(curve, u1, u2, defl, NbMin, aPars);
457
458 //        IntCurveSurface_ThePolygon       polygon(curve,TheCurveTool::NbSamples(curve,u1,u2));
459           IntCurveSurface_ThePolygon  polygon(curve, aPars->Array1());
460           InternalPerform(curve,polygon,surface,UU1,VV1,UU2,VV2);
461         }
462       }
463       else {   //-- la surface est une quadrique 
464         InternalPerformCurveQuadric(curve,surface);
465       }
466     }
467   }
468 }
469 //=======================================================================
470 //function : Perform
471 //purpose  : 
472 //=======================================================================
473 void IntCurveSurface_Inter::Perform(const TheCurve&   curve,
474                                     const IntCurveSurface_ThePolygon& polygon,
475                                     const TheSurface& surface) {
476   ResetFields();
477   done = Standard_True;
478   Standard_Real u1,v1,u2,v2;
479   u1 = TheSurfaceTool::FirstUParameter(surface);
480   v1 = TheSurfaceTool::FirstVParameter(surface);
481   u2 = TheSurfaceTool::LastUParameter(surface);
482   v2 = TheSurfaceTool::LastVParameter(surface);
483   Standard_Integer nbsu,nbsv;
484   nbsu = TheSurfaceTool::NbSamplesU(surface,u1,u2);
485   nbsv = TheSurfaceTool::NbSamplesV(surface,v1,v2);
486   if(nbsu>40) nbsu=40;
487   if(nbsv>40) nbsv=40;
488   IntCurveSurface_ThePolyhedron    polyhedron(surface,nbsu,nbsv,u1,v1,u2,v2);
489   Perform(curve,polygon,surface,polyhedron);
490 }
491 //=======================================================================
492 //function : Perform
493 //purpose  : 
494 //=======================================================================
495 void IntCurveSurface_Inter::Perform(const TheCurve&   curve,
496                                     const TheSurface& surface,
497                                     const IntCurveSurface_ThePolyhedron& polyhedron) {
498   ResetFields();
499   done = Standard_True;
500   Standard_Real u1 = TheCurveTool::FirstParameter(curve);
501   Standard_Real u2 = TheCurveTool::LastParameter(curve);
502   IntCurveSurface_ThePolygon polygon(curve,TheCurveTool::NbSamples(curve,u1,u2));
503   Perform(curve,polygon,surface,polyhedron);
504 }
505 //=======================================================================
506 //function : Perform
507 //purpose  : 
508 //=======================================================================
509 void IntCurveSurface_Inter::Perform(const TheCurve&       curve,
510                                     const IntCurveSurface_ThePolygon&     polygon,
511                                     const TheSurface&     surface,
512                                     const IntCurveSurface_ThePolyhedron&  polyhedron) {
513   ResetFields();
514   done = Standard_True;
515   Standard_Real u1,v1,u2,v2;
516   u1 = TheSurfaceTool::FirstUParameter(surface);
517   v1 = TheSurfaceTool::FirstVParameter(surface);
518   u2 = TheSurfaceTool::LastUParameter(surface);
519   v2 = TheSurfaceTool::LastVParameter(surface);
520   InternalPerform(curve,polygon,surface,polyhedron,u1,v1,u2,v2);
521 }
522
523 //=======================================================================
524 //function : Perform
525 //purpose  : 
526 //=======================================================================
527 void IntCurveSurface_Inter::Perform(const TheCurve&       curve,
528                                     const IntCurveSurface_ThePolygon&     polygon,
529                                     const TheSurface&     surface,
530                                     const IntCurveSurface_ThePolyhedron&  polyhedron,
531                                     Bnd_BoundSortBox& BndBSB) {
532   ResetFields();
533   done = Standard_True;
534   Standard_Real u1,v1,u2,v2;
535   u1 = TheSurfaceTool::FirstUParameter(surface);
536   v1 = TheSurfaceTool::FirstVParameter(surface);
537   u2 = TheSurfaceTool::LastUParameter(surface);
538   v2 = TheSurfaceTool::LastVParameter(surface);
539   InternalPerform(curve,polygon,surface,polyhedron,u1,v1,u2,v2,BndBSB);
540 }
541 //=======================================================================
542 //function : InternalPerform
543 //purpose  : C a l c u l   d u   p o i n t   a p p r o c h e           
544 //==              p u i s   d u   p o i n t   E x a c t  
545 //=======================================================================
546 void IntCurveSurface_Inter::InternalPerform(const TheCurve&       curve,
547                                             const IntCurveSurface_ThePolygon&     polygon,
548                                             const TheSurface&     surface,
549                                             const IntCurveSurface_ThePolyhedron&  polyhedron,
550                                             const Standard_Real u0,
551                                             const Standard_Real v0,
552                                             const Standard_Real u1, 
553                                             const Standard_Real v1,
554                                             Bnd_BoundSortBox& BSB) {
555   IntCurveSurface_TheInterference  interference(polygon,polyhedron,BSB);
556   IntCurveSurface_TheCSFunction    theicsfunction(surface,curve);
557   IntCurveSurface_TheExactInter    intersectionExacte(theicsfunction,TOLTANGENCY);
558   math_FunctionSetRoot rsnld(intersectionExacte.Function());
559
560 //  Standard_Real    u,v,w,winit;
561   Standard_Real    u,v,w;
562   gp_Pnt P;
563   Standard_Real winf = polygon.InfParameter();
564   Standard_Real wsup = polygon.SupParameter();
565   Standard_Integer NbSectionPoints = interference.NbSectionPoints();
566   Standard_Integer NbTangentZones  = interference.NbTangentZones();
567
568
569   //-- Les interferences renvoient parfois de nombreuses fois (>20) les memes points 
570   
571   Standard_Integer i,NbStartPoints=NbSectionPoints;
572   for(i=1; i<= NbTangentZones; i++) {
573     const Intf_TangentZone& TZ = interference.ZoneValue(i);
574     Standard_Integer nbpnts = TZ.NumberOfPoints();
575     NbStartPoints+=nbpnts;
576   }
577   
578   if(NbStartPoints) { 
579     Standard_Real *TabU = new Standard_Real [NbStartPoints+1];
580     Standard_Real *TabV = new Standard_Real [NbStartPoints+1];
581     Standard_Real *TabW = new Standard_Real [NbStartPoints+1];
582     Standard_Integer IndexPoint=0;
583     
584     for(i=1; i<= NbSectionPoints; i++) {
585       const Intf_SectionPoint& SP = interference.PntValue(i);
586       SectionPointToParameters(SP,polyhedron,polygon,u,v,w);
587       TabU[IndexPoint]=u;
588       TabV[IndexPoint]=v;
589       TabW[IndexPoint]=w;
590       IndexPoint++;
591     }
592     for(i=1; i<= NbTangentZones; i++) {
593       const Intf_TangentZone& TZ = interference.ZoneValue(i);
594       Standard_Integer nbpnts = TZ.NumberOfPoints();
595       for(Standard_Integer j=1; j<=nbpnts; j++) { 
596         const Intf_SectionPoint& SP = TZ.GetPoint(j);
597         SectionPointToParameters(SP,polyhedron,polygon,u,v,w);  
598         TabU[IndexPoint]=u;
599         TabV[IndexPoint]=v;
600         TabW[IndexPoint]=w;
601         IndexPoint++;
602       }
603     }
604     
605     //-- Tri
606     Standard_Real su=0,sv=0,sw=0,ptol;
607     ptol = 10*Precision::PConfusion();
608     
609     //-- Tri suivant la variable W 
610     Standard_Boolean Triok;
611     do { 
612       Triok=Standard_True;
613       Standard_Integer im1;
614       for(i=1,im1=0;i<NbStartPoints;im1++,i++) {       
615         if(TabW[i] < TabW[im1]) { 
616           Standard_Real t=TabW[i]; TabW[i]=TabW[im1]; TabW[im1]=t;
617           t=TabU[i]; TabU[i]=TabU[im1]; TabU[im1]=t;
618           t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
619           Triok=Standard_False;
620         }
621       }
622     }
623     while(Triok==Standard_False);
624     
625     //-- On trie pour des meme W suivant U
626     do { 
627       Triok=Standard_True;
628       Standard_Integer im1;
629       for(i=1,im1=0;i<NbStartPoints;im1++,i++) {       
630 // modified by NIZHNY-MKK  Mon Oct  3 17:38:49 2005
631 //      if(Abs(TabW[i]-TabW[im1])<ptol) { 
632         if((TabW[i]-TabW[im1])<ptol) { 
633           TabW[i]=TabW[im1];
634           if(TabU[i] < TabU[im1]) { 
635             Standard_Real t=TabU[i]; TabU[i]=TabU[im1]; TabU[im1]=t;
636             t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
637             Triok=Standard_False;
638           }
639         }
640       }
641     }
642     while(Triok==Standard_False);
643     
644     //-- On trie pour des meme U et W suivant V
645     do { 
646       Triok=Standard_True;
647       Standard_Integer im1;
648       for(i=1,im1=0;i<NbStartPoints;im1++,i++) {       
649 // modified by NIZHNY-MKK  Mon Oct  3 17:38:52 2005
650 //      if((Abs(TabW[i]-TabW[im1])<ptol) && (Abs(TabU[i]-TabU[im1])<ptol)) { 
651         if(((TabW[i]-TabW[im1])<ptol) && ((TabU[i]-TabU[im1])<ptol)) { 
652           TabU[i]=TabU[im1];
653           if(TabV[i] < TabV[im1]) { 
654             Standard_Real t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
655             Triok=Standard_False;
656           }
657         }
658       }
659     }
660     while(Triok==Standard_False);
661     
662     
663     for(i=0;i<NbStartPoints; i++) { 
664       u=TabU[i]; v=TabV[i]; w=TabW[i];
665       if(i==0) {
666         su=u-1; 
667       }
668       if(Abs(u-su)>ptol || Abs(v-sv)>ptol || Abs(w-sw)>ptol) { 
669         intersectionExacte.Perform(u,v,w,rsnld,u0,u1,v0,v1,winf,wsup);
670         if(intersectionExacte.IsDone()) { 
671           if(!intersectionExacte.IsEmpty()) { 
672             P=intersectionExacte.Point();
673             w=intersectionExacte.ParameterOnCurve();
674             intersectionExacte.ParameterOnSurface(u,v);
675             AppendPoint(curve,w,surface,u,v);
676           }
677         }
678       }
679       su=TabU[i]; sv=TabV[i]; sw=TabW[i];
680     }
681     delete [] TabW;
682     delete [] TabV;
683     delete [] TabU;
684   }
685 }
686
687 //=======================================================================
688 //function : InternalPerform
689 //purpose  : 
690 //=======================================================================
691 void IntCurveSurface_Inter::InternalPerform(const TheCurve&       curve,
692                                             const IntCurveSurface_ThePolygon&     polygon,
693                                             const TheSurface&     surface,
694                                             const IntCurveSurface_ThePolyhedron&  polyhedron,
695                                             const Standard_Real u0,
696                                             const Standard_Real v0,
697                                             const Standard_Real u1, 
698                                             const Standard_Real v1) { 
699   IntCurveSurface_TheInterference  interference(polygon,polyhedron);
700   IntCurveSurface_TheCSFunction    theicsfunction(surface,curve);
701   IntCurveSurface_TheExactInter    intersectionExacte(theicsfunction,TOLTANGENCY);
702   math_FunctionSetRoot rsnld(intersectionExacte.Function());
703
704 //  Standard_Real    u,v,w,winit;
705   Standard_Real    u,v,w;
706   gp_Pnt P;
707   Standard_Real winf = polygon.InfParameter();
708   Standard_Real wsup = polygon.SupParameter();
709   Standard_Integer NbSectionPoints = interference.NbSectionPoints();
710   Standard_Integer NbTangentZones  = interference.NbTangentZones();
711
712
713   //-- Les interferences renvoient parfois de nombreuses fois (>20) les memes points 
714   
715   Standard_Integer i,NbStartPoints=NbSectionPoints;
716   for(i=1; i<= NbTangentZones; i++) {
717     const Intf_TangentZone& TZ = interference.ZoneValue(i);
718     Standard_Integer nbpnts = TZ.NumberOfPoints();
719     NbStartPoints+=nbpnts;
720   }
721   
722   if(NbStartPoints) { 
723     Standard_Real *TabU = new Standard_Real [NbStartPoints+1];
724     Standard_Real *TabV = new Standard_Real [NbStartPoints+1];
725     Standard_Real *TabW = new Standard_Real [NbStartPoints+1];
726     Standard_Integer IndexPoint=0;
727     
728     for(i=1; i<= NbSectionPoints; i++) {
729       const Intf_SectionPoint& SP = interference.PntValue(i);
730       SectionPointToParameters(SP,polyhedron,polygon,u,v,w);
731       TabU[IndexPoint]=u;
732       TabV[IndexPoint]=v;
733       TabW[IndexPoint]=w;
734       IndexPoint++;
735     }
736     for(i=1; i<= NbTangentZones; i++) {
737       const Intf_TangentZone& TZ = interference.ZoneValue(i);
738       Standard_Integer nbpnts = TZ.NumberOfPoints();
739       for(Standard_Integer j=1; j<=nbpnts; j++) { 
740         const Intf_SectionPoint& SP = TZ.GetPoint(j);
741         SectionPointToParameters(SP,polyhedron,polygon,u,v,w);  
742         TabU[IndexPoint]=u;
743         TabV[IndexPoint]=v;
744         TabW[IndexPoint]=w;
745         IndexPoint++;
746       }
747     }
748     
749     //-- Tri
750     Standard_Real su=0,sv=0,sw=0,ptol;
751     ptol = 10*Precision::PConfusion();
752     
753     //-- Tri suivant la variable W 
754     Standard_Boolean Triok;
755     do { 
756       Triok=Standard_True;
757       Standard_Integer im1;
758       for(i=1,im1=0;i<NbStartPoints;im1++,i++) {       
759         if(TabW[i] < TabW[im1]) { 
760           Standard_Real t=TabW[i]; TabW[i]=TabW[im1]; TabW[im1]=t;
761           t=TabU[i]; TabU[i]=TabU[im1]; TabU[im1]=t;
762           t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
763           Triok=Standard_False;
764         }
765       }
766     }
767     while(Triok==Standard_False);
768     
769     //-- On trie pour des meme W suivant U
770     do { 
771       Triok=Standard_True;
772       Standard_Integer im1;
773       for(i=1,im1=0;i<NbStartPoints;im1++,i++) {       
774 // modified by NIZHNY-MKK  Mon Oct  3 17:38:56 2005
775 //      if(Abs(TabW[i]-TabW[im1])<ptol) { 
776         if((TabW[i]-TabW[im1])<ptol) { 
777           TabW[i]=TabW[im1];
778           if(TabU[i] < TabU[im1]) { 
779             Standard_Real t=TabU[i]; TabU[i]=TabU[im1]; TabU[im1]=t;
780             t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
781             Triok=Standard_False;
782           }
783         }
784       }
785     }
786     while(Triok==Standard_False);
787     
788     //-- On trie pour des meme U et W suivant V
789     do { 
790       Triok=Standard_True;
791       Standard_Integer im1;
792       for(i=1,im1=0;i<NbStartPoints;im1++,i++) {       
793 // modified by NIZHNY-MKK  Mon Oct  3 17:38:58 2005
794 //      if((Abs(TabW[i]-TabW[im1])<ptol) && (Abs(TabU[i]-TabU[im1])<ptol)) { 
795         if(((TabW[i]-TabW[im1])<ptol) && ((TabU[i]-TabU[im1])<ptol)) { 
796           TabU[i]=TabU[im1];
797           if(TabV[i] < TabV[im1]) { 
798             Standard_Real t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
799             Triok=Standard_False;
800           }
801         }
802       }
803     }
804     while(Triok==Standard_False);
805     
806     
807     for(i=0;i<NbStartPoints; i++) { 
808       u=TabU[i]; v=TabV[i]; w=TabW[i];
809       if(i==0) {
810         su=u-1; 
811       }
812       if(Abs(u-su)>ptol || Abs(v-sv)>ptol || Abs(w-sw)>ptol) { 
813         intersectionExacte.Perform(u,v,w,rsnld,u0,u1,v0,v1,winf,wsup);
814         if(intersectionExacte.IsDone()) { 
815           if(!intersectionExacte.IsEmpty()) { 
816             P=intersectionExacte.Point();
817             w=intersectionExacte.ParameterOnCurve();
818             intersectionExacte.ParameterOnSurface(u,v);
819             AppendPoint(curve,w,surface,u,v);
820           }
821         }
822       }
823       su=TabU[i]; sv=TabV[i]; sw=TabW[i];
824     }
825     delete [] TabW;
826     delete [] TabV;
827     delete [] TabU;
828   }
829 }
830 //=======================================================================
831 //function : InternalPerformCurveQuadric
832 //purpose  : 
833 //=======================================================================
834 void IntCurveSurface_Inter::InternalPerformCurveQuadric(const TheCurve&                 curve,
835                                                         const TheSurface&               surface) { 
836   IntCurveSurface_TheQuadCurvExactInter QuadCurv(surface,curve);  
837   if(QuadCurv.IsDone()) { 
838     Standard_Integer NbRoots = QuadCurv.NbRoots();
839     Standard_Real u,v,w;
840     for(Standard_Integer i = 1; i<= NbRoots; i++) { 
841       w = QuadCurv.Root(i); 
842       IntCurveSurface_ComputeParamsOnQuadric(surface,TheCurveTool::Value(curve,w),u,v);
843       AppendPoint(curve,w,surface,u,v); 
844     }
845     //-- Intervals non traites .............................................
846   }
847 }
848
849 //=======================================================================
850 //function : InternalPerform
851 //purpose  : 
852 //=======================================================================
853 void IntCurveSurface_Inter::InternalPerform(const TheCurve&       curve,
854                                             const IntCurveSurface_ThePolygon&  polygon,
855                                             const TheSurface&     surface,
856                                             const Standard_Real U1,
857                                             const Standard_Real V1,
858                                             const Standard_Real U2,
859                                             const Standard_Real V2) {
860   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
861   if(   (SurfaceType != GeomAbs_Plane) 
862      && (SurfaceType != GeomAbs_Cylinder)
863      && (SurfaceType != GeomAbs_Cone) 
864      && (SurfaceType != GeomAbs_Sphere) ) {
865     if(SurfaceType != GeomAbs_BSplineSurface) {
866       Standard_Integer nbsu,nbsv;
867       nbsu = TheSurfaceTool::NbSamplesU(surface,U1,U2);
868       nbsv = TheSurfaceTool::NbSamplesV(surface,V1,V2);
869       if(nbsu>40) nbsu=40;
870       if(nbsv>40) nbsv=40;
871       IntCurveSurface_ThePolyhedron polyhedron(surface,nbsu,nbsv,U1,V1,U2,V2);
872       InternalPerform(curve,polygon,surface,polyhedron,U1,V1,U2,V2); 
873     }
874     else {
875       Handle(Adaptor3d_HSurface) aS = TheSurfaceTool::UTrim(surface, U1, U2, 1.e-9);
876       aS = aS->VTrim(V1, V2, 1.e-9);
877       Handle(Adaptor3d_TopolTool) aTopTool = new Adaptor3d_TopolTool(aS);
878       Standard_Real defl = 0.1;
879       aTopTool->SamplePnts(defl, 10, 10);
880
881       Standard_Integer nbpu = aTopTool->NbSamplesU();
882       Standard_Integer nbpv = aTopTool->NbSamplesV();
883       TColStd_Array1OfReal Upars(1, nbpu), Vpars(1, nbpv);
884       aTopTool->UParameters(Upars);
885       aTopTool->VParameters(Vpars);
886
887       IntCurveSurface_ThePolyhedron polyhedron(surface,Upars, Vpars);
888       InternalPerform(curve,polygon,surface,polyhedron,U1,V1,U2,V2); 
889     }
890   }
891   else {
892     IntCurveSurface_TheQuadCurvExactInter QuadCurv(surface,curve);  
893     if(QuadCurv.IsDone()) { 
894       Standard_Integer NbRoots = QuadCurv.NbRoots();
895       Standard_Real u,v,w;
896       for(Standard_Integer i = 1; i<= NbRoots; i++) { 
897         w = QuadCurv.Root(i); 
898         IntCurveSurface_ComputeParamsOnQuadric(surface,TheCurveTool::Value(curve,w),u,v);
899         AppendPoint(curve,w,surface,u,v); 
900       }
901       //-- Intervalles non traites .............................................
902     }
903   } //-- Fin : la Surface  est une quadrique
904 }
905 //=======================================================================
906 //function : PerformConicSurf
907 //purpose  : 
908 //=======================================================================
909 void IntCurveSurface_Inter::PerformConicSurf(const gp_Lin&   Line,
910                                              const TheCurve&       curve,
911                                              const TheSurface&     surface,
912                                              const Standard_Real U1,
913                                              const Standard_Real V1,
914                                              const Standard_Real U2,
915                                              const Standard_Real V2) {
916   
917
918   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
919   switch(SurfaceType) { 
920   case GeomAbs_Plane: 
921     {
922       IntAna_IntConicQuad LinPlane(Line,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE);
923       AppendIntAna(curve,surface,LinPlane);
924       break;
925     }
926   case GeomAbs_Cylinder:
927     {
928       IntAna_IntConicQuad LinCylinder(Line,TheSurfaceTool::Cylinder(surface));
929       AppendIntAna(curve,surface,LinCylinder);
930       break;
931     }
932   case GeomAbs_Sphere:
933     {
934       IntAna_IntConicQuad LinSphere(Line,TheSurfaceTool::Sphere(surface));
935       AppendIntAna(curve,surface,LinSphere);
936       break;
937     }
938   case GeomAbs_Torus: 
939     {
940       IntAna_IntLinTorus intlintorus(Line,TheSurfaceTool::Torus(surface));
941       if(intlintorus.IsDone()) { 
942         Standard_Integer nbp = intlintorus.NbPoints();
943         Standard_Real fi,theta,w;
944         for(Standard_Integer i = 1; i<= nbp; i++) { 
945           const gp_Pnt aDebPnt(intlintorus.Value(i));
946           (void )aDebPnt;
947           w = intlintorus.ParamOnLine(i);
948           intlintorus.ParamOnTorus(i,fi,theta);
949           AppendPoint(curve,w,surface,fi,theta);
950         }
951         break;
952       }
953     }  //-- Si Done retourne False, On passe dans Default !! 
954   case GeomAbs_Cone:
955     {
956       //OCC516(apo)->
957       const Standard_Real correction = 1.E+5*Precision::Angular();
958       gp_Cone cn = TheSurfaceTool::Cone(surface);
959       if(Abs(cn.SemiAngle()) < M_PI/2.0 - correction){
960         IntAna_IntConicQuad LinCone(Line,cn);
961         AppendIntAna(curve,surface,LinCone);
962         break;
963       }//<-OCC516(apo)
964     }   
965   default: 
966     {
967       Standard_Integer nbsu,nbsv;
968       nbsu = TheSurfaceTool::NbSamplesU(surface,U1,U2);
969       nbsv = TheSurfaceTool::NbSamplesV(surface,V1,V2);
970
971       Standard_Boolean U1inf = Precision::IsInfinite(U1);
972       Standard_Boolean U2inf = Precision::IsInfinite(U2);
973       Standard_Boolean V1inf = Precision::IsInfinite(V1);
974       Standard_Boolean V2inf = Precision::IsInfinite(V2);
975
976       Standard_Real U1new=U1, U2new=U2, V1new=V1, V2new=V2;
977
978       Standard_Boolean NoIntersection = Standard_False;
979
980       if(U1inf || U2inf || V1inf || V2inf ) {
981
982         if(SurfaceType == GeomAbs_SurfaceOfExtrusion) {
983
984           EstLimForInfExtr(Line, surface, Standard_False, nbsu, 
985                            U1inf, U2inf, V1inf, V2inf, 
986                            U1new, U2new, V1new, V2new, NoIntersection);
987
988         } 
989         else if(SurfaceType == GeomAbs_SurfaceOfRevolution) {
990
991           EstLimForInfRevl(Line, surface,  
992                            U1inf, U2inf, V1inf, V2inf, 
993                            U1new, U2new, V1new, V2new, NoIntersection);
994
995         } 
996         else if(SurfaceType == GeomAbs_OffsetSurface) {
997
998           EstLimForInfOffs(Line, surface, nbsu,  
999                            U1inf, U2inf, V1inf, V2inf, 
1000                            U1new, U2new, V1new, V2new, NoIntersection);
1001         }
1002         else {
1003             
1004           EstLimForInfSurf(U1new, U2new, V1new, V2new);
1005         }
1006
1007       }
1008
1009
1010       if(NoIntersection) return;
1011
1012
1013       // modified by NIZHNY-OFV  Mon Aug 20 14:56:47 2001 (60963 begin)
1014       if(nbsu<20) nbsu=20;
1015       if(nbsv<20) nbsv=20;
1016       // modified by NIZHNY-OFV  Mon Aug 20 14:57:06 2001 (60963 end)
1017       IntCurveSurface_ThePolyhedron polyhedron(surface,nbsu,nbsv,U1new,V1new,U2new,V2new);
1018       Intf_Tool                         bndTool;
1019       Bnd_Box                          boxLine;
1020       bndTool.LinBox(Line,polyhedron.Bounding(),boxLine);
1021       for(Standard_Integer nbseg=1; nbseg<= bndTool.NbSegments(); nbseg++) { 
1022         Standard_Real pinf = bndTool.BeginParam(nbseg);
1023         Standard_Real psup = bndTool.EndParam(nbseg);
1024         if((psup - pinf)<1e-10) { pinf-=1e-10; psup+=1e-10; } 
1025         IntCurveSurface_ThePolygon polygon(curve, pinf,psup,2);
1026         InternalPerform(curve,polygon,surface,polyhedron,U1new,V1new,U2new,V2new);
1027       }
1028     }
1029   }
1030 }
1031 //=======================================================================
1032 //function : PerformConicSurf
1033 //purpose  : 
1034 //=======================================================================
1035 void IntCurveSurface_Inter::PerformConicSurf(const gp_Circ&        Circle,
1036                                              const TheCurve&       curve,
1037                                              const TheSurface&     surface,
1038                                              const Standard_Real U1,
1039                                              const Standard_Real V1,
1040                                              const Standard_Real U2,
1041                                              const Standard_Real V2) { 
1042   
1043   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1044   switch(SurfaceType) { 
1045   case GeomAbs_Plane: 
1046     {
1047       IntAna_IntConicQuad CircPlane(Circle,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE,TOLERANCE);
1048       AppendIntAna(curve,surface,CircPlane);
1049       break;
1050     }
1051   case GeomAbs_Cylinder:
1052     {
1053       IntAna_IntConicQuad CircCylinder(Circle,TheSurfaceTool::Cylinder(surface));
1054       AppendIntAna(curve,surface,CircCylinder);
1055       break;
1056     }
1057   case GeomAbs_Cone:
1058     {
1059       IntAna_IntConicQuad CircCone(Circle,TheSurfaceTool::Cone(surface));
1060       AppendIntAna(curve,surface,CircCone);
1061       break;
1062     }   
1063   case GeomAbs_Sphere:
1064     {
1065       IntAna_IntConicQuad CircSphere(Circle,TheSurfaceTool::Sphere(surface));
1066       AppendIntAna(curve,surface,CircSphere);
1067       break;
1068     }
1069   default: 
1070     {
1071       IntCurveSurface_ThePolygon polygon(curve,NBSAMPLESONCIRCLE);
1072       InternalPerform(curve,polygon,surface,U1,V1,U2,V2);
1073     }
1074   }
1075 }
1076 //=======================================================================
1077 //function : PerformConicSurf
1078 //purpose  : 
1079 //=======================================================================
1080 void IntCurveSurface_Inter::PerformConicSurf(const gp_Elips&       Ellipse,
1081                                              const TheCurve&       curve,
1082                                              const TheSurface&     surface,
1083                                              const Standard_Real U1,
1084                                              const Standard_Real V1,
1085                                              const Standard_Real U2,
1086                                              const Standard_Real V2) { 
1087
1088   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1089   switch(SurfaceType) { 
1090   case GeomAbs_Plane: 
1091     {
1092       IntAna_IntConicQuad EllipsePlane(Ellipse,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE,TOLERANCE);
1093       AppendIntAna(curve,surface,EllipsePlane);
1094       break;
1095     }
1096   case GeomAbs_Cylinder:
1097     {
1098       IntAna_IntConicQuad EllipseCylinder(Ellipse,TheSurfaceTool::Cylinder(surface));
1099       AppendIntAna(curve,surface,EllipseCylinder);
1100       break;
1101     }
1102   case GeomAbs_Cone:
1103      {
1104       IntAna_IntConicQuad EllipseCone(Ellipse,TheSurfaceTool::Cone(surface));
1105       AppendIntAna(curve,surface,EllipseCone);
1106       break;
1107     }   
1108   case GeomAbs_Sphere:
1109     {
1110       IntAna_IntConicQuad EllipseSphere(Ellipse,TheSurfaceTool::Sphere(surface));
1111       AppendIntAna(curve,surface,EllipseSphere);
1112       break;
1113     }
1114   default: 
1115     {
1116       IntCurveSurface_ThePolygon polygon(curve,NBSAMPLESONELLIPSE);
1117       InternalPerform(curve,polygon,surface,U1,V1,U2,V2);
1118     }
1119   }
1120 }
1121 //=======================================================================
1122 //function : PerformConicSurf
1123 //purpose  : 
1124 //=======================================================================
1125 void IntCurveSurface_Inter::PerformConicSurf(const gp_Parab&       Parab,
1126                                              const TheCurve&       curve,
1127                                              const TheSurface&     surface,
1128                                              const Standard_Real U1,
1129                                              const Standard_Real V1,
1130                                              const Standard_Real U2,
1131                                              const Standard_Real V2) { 
1132
1133   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1134   switch(SurfaceType) { 
1135   case GeomAbs_Plane: 
1136     {
1137       IntAna_IntConicQuad ParabPlane(Parab,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE);
1138       AppendIntAna(curve,surface,ParabPlane);
1139       break;
1140     }
1141   case GeomAbs_Cylinder:
1142     {
1143       IntAna_IntConicQuad ParabCylinder(Parab,TheSurfaceTool::Cylinder(surface));
1144       AppendIntAna(curve,surface,ParabCylinder);
1145       break;
1146     }
1147   case GeomAbs_Cone:
1148      {
1149       IntAna_IntConicQuad ParabCone(Parab,TheSurfaceTool::Cone(surface));
1150       AppendIntAna(curve,surface,ParabCone);
1151       break;
1152     }   
1153   case GeomAbs_Sphere:
1154     {
1155       IntAna_IntConicQuad ParabSphere(Parab,TheSurfaceTool::Sphere(surface));
1156       AppendIntAna(curve,surface,ParabSphere);
1157       break;
1158     }
1159   default: 
1160     {
1161       Standard_Integer nbsu,nbsv;
1162       nbsu = TheSurfaceTool::NbSamplesU(surface,U1,U2);
1163       nbsv = TheSurfaceTool::NbSamplesV(surface,V1,V2);
1164       if(nbsu>40) nbsu=40;
1165       if(nbsv>40) nbsv=40;
1166       IntCurveSurface_ThePolyhedron polyhedron(surface,nbsu,nbsv,U1,V1,U2,V2);
1167       Intf_Tool                         bndTool;
1168       Bnd_Box                          boxParab;
1169       bndTool.ParabBox(Parab,polyhedron.Bounding(),boxParab);
1170       for(Standard_Integer nbseg=1; nbseg<= bndTool.NbSegments(); nbseg++) { 
1171         IntCurveSurface_ThePolygon polygon(curve,
1172                                            bndTool.BeginParam(nbseg),
1173                                            bndTool.EndParam(nbseg),
1174                                            NBSAMPLESONPARAB);
1175         InternalPerform(curve,polygon,surface,polyhedron,U1,V1,U2,V2);
1176       }
1177     }
1178   }
1179 }
1180 //=======================================================================
1181 //function : PerformConicSurf
1182 //purpose  : 
1183 //=======================================================================
1184 void IntCurveSurface_Inter::PerformConicSurf(const gp_Hypr&        Hypr,
1185                                              const TheCurve&       curve,
1186                                              const TheSurface&     surface,
1187                                              const Standard_Real U1,
1188                                              const Standard_Real V1,
1189                                              const Standard_Real U2,
1190                                              const Standard_Real V2) {
1191  
1192   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1193   switch(SurfaceType) { 
1194   case GeomAbs_Plane: 
1195     {
1196       IntAna_IntConicQuad HyprPlane(Hypr,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE);
1197       AppendIntAna(curve,surface,HyprPlane);
1198       break;
1199     }
1200   case GeomAbs_Cylinder:
1201     {
1202       IntAna_IntConicQuad HyprCylinder(Hypr,TheSurfaceTool::Cylinder(surface));
1203       AppendIntAna(curve,surface,HyprCylinder);
1204       break;
1205     }
1206   case GeomAbs_Cone:
1207      {
1208       IntAna_IntConicQuad HyprCone(Hypr,TheSurfaceTool::Cone(surface));
1209       AppendIntAna(curve,surface,HyprCone);
1210       break;
1211     }   
1212   case GeomAbs_Sphere:
1213     {
1214       IntAna_IntConicQuad HyprSphere(Hypr,TheSurfaceTool::Sphere(surface));
1215       AppendIntAna(curve,surface,HyprSphere);
1216       break;
1217     }
1218   default: 
1219     {
1220       Standard_Integer nbsu,nbsv;
1221       nbsu = TheSurfaceTool::NbSamplesU(surface,U1,U2);
1222       nbsv = TheSurfaceTool::NbSamplesV(surface,V1,V2);
1223       if(nbsu>40) nbsu=40;
1224       if(nbsv>40) nbsv=40;
1225       IntCurveSurface_ThePolyhedron polyhedron(surface,nbsu,nbsv,U1,V1,U2,V2);
1226       Intf_Tool                         bndTool;
1227       Bnd_Box                          boxHypr;
1228       bndTool.HyprBox(Hypr,polyhedron.Bounding(),boxHypr);
1229       for(Standard_Integer nbseg=1; nbseg<= bndTool.NbSegments(); nbseg++) { 
1230         IntCurveSurface_ThePolygon polygon(curve,
1231                                            bndTool.BeginParam(nbseg),
1232                                            bndTool.EndParam(nbseg),
1233                                            NBSAMPLESONHYPR);
1234         InternalPerform(curve,polygon,surface,polyhedron,U1,V1,U2,V2);
1235       }
1236     }
1237   }
1238 }
1239 //=======================================================================
1240 //function : AppendIntAna
1241 //purpose  : 
1242 //=======================================================================
1243 void IntCurveSurface_Inter::AppendIntAna(const TheCurve& curve,
1244                                          const TheSurface& surface,
1245                                          const IntAna_IntConicQuad& intana_ConicQuad) {
1246   if(intana_ConicQuad.IsDone()) { 
1247     if(intana_ConicQuad.IsInQuadric()) { 
1248       //-- cout<<" Courbe Dans la Quadrique !!! Non Traite !!!"<<endl;
1249     }
1250     else if(intana_ConicQuad.IsParallel()) { 
1251       //-- cout<<" Courbe // a la Quadrique !!! Non Traite !!!"<<endl;
1252     }
1253     else { 
1254       Standard_Integer nbp = intana_ConicQuad.NbPoints();
1255       Standard_Real u,v,w;
1256       for(Standard_Integer i = 1; i<= nbp; i++) { 
1257         gp_Pnt P(intana_ConicQuad.Point(i));
1258         w = intana_ConicQuad.ParamOnConic(i);
1259         IntCurveSurface_ComputeParamsOnQuadric(surface,P,u,v);
1260         AppendPoint(curve,w,surface,u,v);
1261       }
1262     }
1263   }
1264   else { 
1265     //-- cout<<" IntAna Conic Quad Not Done  "<<endl;
1266   }
1267 }
1268 //=======================================================================
1269 //function : AppendPoint
1270 //purpose  : 
1271 //=======================================================================
1272 void IntCurveSurface_Inter::AppendPoint(const TheCurve& curve,
1273                                         const Standard_Real lw,
1274                                         const TheSurface& surface,
1275                                         const Standard_Real su,
1276                                         const Standard_Real sv) {
1277  
1278   Standard_Real W0 = TheCurveTool::FirstParameter(curve);
1279   Standard_Real W1 = TheCurveTool::LastParameter(curve);
1280   Standard_Real U0 = TheSurfaceTool::FirstUParameter(surface);
1281   Standard_Real U1 = TheSurfaceTool::LastUParameter(surface);
1282   Standard_Real V0 = TheSurfaceTool::FirstVParameter(surface);
1283   Standard_Real V1 = TheSurfaceTool::LastVParameter(surface);
1284   //-- Test si la courbe est periodique
1285   Standard_Real w = lw, u = su, v = sv;
1286
1287   GeomAbs_CurveType aCType = TheCurveTool::GetType(curve);
1288
1289   if(TheCurveTool::IsPeriodic(curve) 
1290      || aCType == GeomAbs_Circle 
1291      || aCType == GeomAbs_Ellipse) {
1292     w = ElCLib::InPeriod(w, W0, W0 + TheCurveTool::Period(curve));
1293   }
1294
1295   if((W0 - w) >= TOLTANGENCY || (w - W1) >= TOLTANGENCY) return;
1296
1297   GeomAbs_SurfaceType aSType = TheSurfaceTool::GetType(surface);
1298   if (TheSurfaceTool::IsUPeriodic(surface)
1299       || aSType == GeomAbs_Cylinder
1300       || aSType == GeomAbs_Cone
1301       || aSType == GeomAbs_Sphere) {
1302     u = ElCLib::InPeriod(u, U0, U0 + TheSurfaceTool::UPeriod(surface));
1303   }
1304
1305   if (TheSurfaceTool::IsVPeriodic(surface)) {
1306     v = ElCLib::InPeriod(v, V0, V0 + TheSurfaceTool::VPeriod(surface));
1307   }
1308
1309   if((U0 - u) >= TOLTANGENCY || (u - U1) >= TOLTANGENCY) return;
1310   if((V0 - v) >= TOLTANGENCY || (v - V1) >= TOLTANGENCY) return;
1311
1312   IntCurveSurface_TransitionOnCurve   TransOnCurve;
1313   IntCurveSurface_ComputeTransitions(curve,w,TransOnCurve,
1314                                      surface,u,v);
1315   gp_Pnt P(TheCurveTool::Value(curve,w));
1316   IntCurveSurface_IntersectionPoint IP(P,u,v,w,TransOnCurve);
1317   Append(IP); //-- invoque la methode de IntCurveSurface_Intersection.
1318   
1319 }
1320
1321 //=======================================================================
1322 //function : AppendSegment
1323 //purpose  : 
1324 //=======================================================================
1325 void IntCurveSurface_Inter::AppendSegment(const TheCurve& ,
1326                                           const Standard_Real ,
1327                                           const Standard_Real ,
1328                                           const TheSurface& ) { 
1329   //cout<<" !!! Not Yet Implemented 
1330   //IntCurveSurface_Inter::Append(const IntCurveSurf ...)"<<endl;
1331 }
1332
1333 //=======================================================================
1334 //function : SectionPointToParameters
1335 //purpose  : P o i n t   d   i n t e r f e r e n c e   - - >  
1336 //    U , V       e t   W 
1337 //=======================================================================
1338 void SectionPointToParameters(const Intf_SectionPoint&              Sp,
1339                               const IntCurveSurface_ThePolyhedron&  Polyhedron,
1340                               const IntCurveSurface_ThePolygon&     Polygon,
1341                               Standard_Real& U,
1342                               Standard_Real& V,
1343                               Standard_Real& W) {
1344   
1345   Intf_PIType       typ;
1346   Standard_Integer  Adr1,Adr2;
1347   Standard_Real     Param,u,v;
1348   gp_Pnt P(Sp.Pnt());
1349   
1350   Standard_Integer Pt1,Pt2,Pt3;
1351   Standard_Real u1 = 0.,v1 = 0.,param;
1352   //----------------------------------------------------------------------
1353   //--          Calcul des parametres approches sur la surface          --
1354   //----------------------------------------------------------------------
1355   
1356   Sp.InfoSecond(typ,Adr1,Adr2,Param);
1357   switch(typ) { 
1358   case Intf_VERTEX:   //-- Adr1 est le numero du vertex
1359     {
1360       Polyhedron.Parameters(Adr1,u1,v1);
1361       break;
1362     }
1363   case Intf_EDGE:
1364     {
1365       Polyhedron.Parameters(Adr1,u1,v1);    
1366       Polyhedron.Parameters(Adr2,u,v);
1367       u1+= Param * (u-u1);
1368       v1+= Param * (v-v1);
1369       break;
1370     }
1371   case Intf_FACE:
1372     {
1373       Standard_Real ua,va,ub,vb,uc,vc,ca,cb,cc,cabc;
1374       Polyhedron.Triangle(Adr1,Pt1,Pt2,Pt3);
1375       gp_Pnt PA(Polyhedron.Point(Pt1));
1376       gp_Pnt PB(Polyhedron.Point(Pt2));
1377       gp_Pnt PC(Polyhedron.Point(Pt3));
1378       Polyhedron.Parameters(Pt1,ua,va);
1379       Polyhedron.Parameters(Pt2,ub,vb);
1380       Polyhedron.Parameters(Pt3,uc,vc);
1381       gp_Vec Normale(gp_Vec(PA,PB).Crossed(gp_Vec(PA,PC)));
1382       cc = (gp_Vec(PA,PB).Crossed(gp_Vec(PA,P))).Dot(Normale);
1383       ca = (gp_Vec(PB,PC).Crossed(gp_Vec(PB,P))).Dot(Normale);
1384       cb = (gp_Vec(PC,PA).Crossed(gp_Vec(PC,P))).Dot(Normale);
1385       cabc = ca + cb + cc;
1386       
1387       ca/=cabc;     cb/=cabc;       cc/=cabc;
1388       
1389       u1 = ca * ua + cb * ub + cc * uc;
1390       v1 = ca * va + cb * vb + cc * vc;
1391       break;
1392     }
1393   default: 
1394     {
1395       cout<<" Default dans SectionPointToParameters "<<endl;
1396       break;
1397     }
1398   }
1399   //----------------------------------------------------------------------
1400   //--                Calcul du point approche sur la Curve             --
1401   //----------------------------------------------------------------------
1402   Standard_Integer SegIndex;
1403
1404   Sp.InfoFirst(typ,SegIndex,param);
1405   W = Polygon.ApproxParamOnCurve(SegIndex,param);
1406   if(param>1.0 || param<0.0) { 
1407     //-- IntCurveSurface_ThePolyhedronTool::Dump(Polyhedron);
1408     //-- IntCurveSurface_ThePolygonTool::Dump(Polygon);
1409     }
1410   U = u1;
1411   V = v1;
1412 }
1413 //=======================================================================
1414 //function : IntCurveSurface_ComputeTransitions
1415 //purpose  : 
1416 //=======================================================================
1417 void IntCurveSurface_ComputeTransitions(const TheCurve& curve,
1418                                         const Standard_Real w,
1419                                         IntCurveSurface_TransitionOnCurve&   TransOnCurve,
1420                                         const TheSurface& surface,
1421                                         const Standard_Real u,
1422                                         const Standard_Real v) { 
1423   
1424   gp_Vec NSurf,D1U,D1V;//TgCurv;
1425   gp_Pnt Psurf;
1426   Standard_Real CosDir;
1427   
1428   
1429   TheSurfaceTool::D1(surface,u,v,Psurf,D1U,D1V);
1430   NSurf = D1U.Crossed(D1V);
1431   TheCurveTool::D1(curve,w,Psurf,D1U);
1432   Standard_Real Norm = NSurf.Magnitude();
1433   if(Norm>TOLERANCE_ANGULAIRE) { 
1434     D1U.Normalize();
1435     CosDir = NSurf.Dot(D1U);
1436     CosDir/=Norm;
1437     if( -CosDir > TOLERANCE_ANGULAIRE) { 
1438       //--  --Curve--->    <----Surface----
1439       TransOnCurve   = IntCurveSurface_In;
1440     }
1441     else if(CosDir > TOLERANCE_ANGULAIRE) { 
1442       //--  --Curve--->  ----Surface-->
1443       TransOnCurve   = IntCurveSurface_Out;
1444     }     
1445     else { 
1446       TransOnCurve   = IntCurveSurface_Tangent;
1447     }
1448   }
1449   else { 
1450     TransOnCurve   = IntCurveSurface_Tangent;
1451   }
1452 }
1453 //=======================================================================
1454 //function : IntCurveSurface_ComputeParamsOnQuadric
1455 //purpose  : 
1456 //=======================================================================
1457 void IntCurveSurface_ComputeParamsOnQuadric(const TheSurface& surface,
1458                                             const gp_Pnt& P,
1459                                             Standard_Real& u,
1460                                             Standard_Real& v) { 
1461   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1462   switch(SurfaceType) { 
1463   case GeomAbs_Plane: 
1464     {
1465       ElSLib::Parameters(TheSurfaceTool::Plane(surface),P,u,v);
1466         break;
1467       }
1468   case GeomAbs_Cylinder:
1469     {
1470       ElSLib::Parameters(TheSurfaceTool::Cylinder(surface),P,u,v);
1471         break;
1472       }
1473   case GeomAbs_Cone:
1474     {
1475       ElSLib::Parameters(TheSurfaceTool::Cone(surface),P,u,v);
1476         break;
1477       }   
1478   case GeomAbs_Sphere:
1479     {
1480       ElSLib::Parameters(TheSurfaceTool::Sphere(surface),P,u,v);
1481         break;
1482       }
1483   default:
1484     break;
1485   }
1486 }
1487 //=======================================================================
1488 //function : EstLimForInfExtr
1489 //purpose  : Estimation of limits for infinite surfaces
1490 //=======================================================================
1491 void EstLimForInfExtr(const gp_Lin&   Line,
1492                       const TheSurface& surface, 
1493                       const Standard_Boolean IsOffSurf,
1494                       const Standard_Integer nbsu, 
1495                       const Standard_Boolean U1inf, 
1496                       const Standard_Boolean U2inf, 
1497                       const Standard_Boolean V1inf, 
1498                       const Standard_Boolean V2inf, 
1499                       Standard_Real& U1new, 
1500                       Standard_Real& U2new, 
1501                       Standard_Real& V1new, 
1502                       Standard_Real& V2new, 
1503                       Standard_Boolean& NoIntersection)
1504
1505 {
1506
1507   NoIntersection = Standard_False;
1508
1509   Handle(Adaptor3d_HSurface) aBasSurf; 
1510
1511   if(IsOffSurf) aBasSurf = TheSurfaceTool::BasisSurface(surface);
1512
1513   gp_Dir aDirOfExt;
1514
1515   if(IsOffSurf) aDirOfExt = aBasSurf->Direction();
1516   else aDirOfExt = TheSurfaceTool::Direction(surface);
1517
1518   Standard_Real tolang = TOLERANCE_ANGULAIRE;
1519
1520   if(aDirOfExt.IsParallel(Line.Direction(), tolang)) {
1521     NoIntersection = Standard_True;
1522     return;
1523   }
1524
1525   if((V1inf || V2inf) && !(U1inf || U2inf)) {
1526
1527     Standard_Real vmin = RealLast(), vmax = -vmin;
1528     gp_Lin aL;
1529     Standard_Real step = (U2new - U1new) / nbsu;
1530     Standard_Real u = U1new, v;
1531     gp_Pnt aP;
1532     Extrema_POnCurv aP1, aP2;
1533     Standard_Integer i;
1534
1535     for(i = 0; i <= nbsu; i++) {
1536
1537       TheSurfaceTool::D0(surface, u, 0., aP);
1538       aL.SetLocation(aP);
1539       aL.SetDirection(aDirOfExt);
1540
1541       Extrema_ExtElC aExtr(aL, Line, tolang);
1542       
1543       if(!aExtr.IsDone()) return;
1544
1545       if(aExtr.IsParallel()) {
1546         NoIntersection = Standard_True;
1547         return;
1548       }
1549
1550       aExtr.Points(1, aP1, aP2);
1551       v = aP1.Parameter();
1552       vmin = Min(vmin, v);
1553       vmax = Max(vmax, v);
1554
1555       u += step;
1556
1557     }
1558
1559     vmin = vmin - Abs(vmin) - 10.;
1560     vmax = vmax + Abs(vmax) + 10.;
1561
1562     V1new = Max(V1new, vmin);
1563     V2new = Min(V2new, vmax);
1564
1565   }
1566   else if(U1inf || U2inf) {
1567
1568     Standard_Real umin = RealLast(), umax = -umin;
1569     Standard_Real u0 = Min(Max(0., U1new), U2new);
1570     Standard_Real v0 = Min(Max(0., V1new), V2new);
1571     gp_Pnt aP;
1572     TheSurfaceTool::D0(surface, u0, v0, aP);
1573     gp_Pln aRefPln(aP, aDirOfExt);
1574
1575     Handle(Adaptor3d_HCurve) aBasCurv;
1576
1577     if(IsOffSurf) aBasCurv = aBasSurf->BasisCurve();
1578     else aBasCurv = TheSurfaceTool::BasisCurve(surface);
1579
1580     ProjLib_Plane Projector(aRefPln);
1581
1582     Projector.Project(Line);
1583
1584     if(!Projector.IsDone()) return;
1585
1586     gp_Lin2d Line2d = Projector.Line();
1587
1588     GeomAbs_CurveType aCurvTyp = aBasCurv->GetType();
1589
1590     if(aCurvTyp == GeomAbs_Line) {
1591       
1592       Projector.Project(aBasCurv->Line());
1593
1594       if(!Projector.IsDone()) return;
1595
1596       gp_Lin2d aL2d = Projector.Line();
1597
1598       IntAna2d_AnaIntersection anInter(Line2d, aL2d);
1599
1600       if(!anInter.IsDone()) return;
1601
1602       if(anInter.IsEmpty() || anInter.IdenticalElements() || 
1603                               anInter.ParallelElements()     ) {
1604         NoIntersection = Standard_True;
1605         return;
1606       }
1607
1608       const IntAna2d_IntPoint& anIntPnt = anInter.Point(1);
1609       umin = umax = anIntPnt.ParamOnSecond();
1610     }
1611     else if(aCurvTyp == GeomAbs_Parabola || aCurvTyp == GeomAbs_Hyperbola) {
1612
1613       IntAna2d_Conic aCon(Line2d);
1614       IntAna2d_AnaIntersection anInter;
1615       
1616       if(aCurvTyp == GeomAbs_Parabola) {
1617         Projector.Project(aBasCurv->Parabola());
1618         if(!Projector.IsDone()) return;
1619
1620         const gp_Parab2d& aP2d = Projector.Parabola();
1621
1622         anInter.Perform(aP2d, aCon);
1623       }
1624       else {
1625         Projector.Project(aBasCurv->Hyperbola());
1626         if(!Projector.IsDone()) return;
1627
1628         const gp_Hypr2d& aH2d = Projector.Hyperbola();
1629         anInter.Perform(aH2d, aCon);
1630       }
1631
1632       if(!anInter.IsDone()) return;
1633
1634       if(anInter.IsEmpty()) {
1635         NoIntersection = Standard_True;
1636         return;
1637       }
1638
1639       Standard_Integer i, nbint = anInter.NbPoints();
1640       for(i = 1; i <= nbint; i++) {
1641
1642         const IntAna2d_IntPoint& anIntPnt = anInter.Point(i);
1643
1644         umin = Min(anIntPnt.ParamOnFirst(), umin);
1645         umax = Max(anIntPnt.ParamOnFirst(), umax);
1646  
1647       }
1648
1649     }
1650     else {
1651       return;
1652     }
1653       
1654     umin = umin - Abs(umin) - 10;
1655     umax = umax + Abs(umax) + 10;
1656
1657     U1new = Max(U1new, umin);
1658     U2new = Min(U2new, umax);
1659     
1660     if(V1inf || V2inf) {
1661       EstLimForInfExtr(Line, surface, IsOffSurf, nbsu, 
1662                        Standard_False, Standard_False, V1inf, V2inf,
1663                        U1new, U2new, V1new, V2new, NoIntersection);
1664     }
1665   }
1666     
1667   return;
1668   
1669 }
1670 //=======================================================================
1671 //function : ProjectIntersectAndEstLim
1672 //purpose  : project <theLine> and it's X-axe symmetric line to <thePln> and
1673 //           intersect resulting curve with <theBasCurvProj>.
1674 //           Then estimate max and min parameters of intersection on
1675 //           <theBasCurvProj>.
1676 //           Is called from EstLimForInfRevl()
1677 //=======================================================================
1678 void ProjectIntersectAndEstLim(const gp_Lin&        theLine,
1679                                const gp_Pln&        thePln,
1680                                const ProjLib_Plane& theBasCurvProj,
1681                                Standard_Real&       theVmin,
1682                                Standard_Real&       theVmax,
1683                                Standard_Boolean&    theNoIntersection)
1684 {
1685   ProjLib_Plane aLineProj( thePln, theLine );
1686   if (!aLineProj.IsDone()) {
1687 #ifdef OCCT_DEBUG
1688   cout
1689   << "Info: IntCurveSurface_Inter::ProjectIntersectAndEstLim(), !aLineProj.IsDone()"
1690   << endl;
1691 #endif
1692     return;
1693   }
1694   gp_Lin2d aLin2d = aLineProj.Line();
1695
1696   // make a second line X-axe symmetric to the first one 
1697   gp_Pnt2d aP1 = aLin2d.Location();
1698   gp_Pnt2d aP2 (aP1.XY() + aLin2d.Direction().XY());
1699   gp_Pnt2d aP1sym (aP1.X(), -aP1.Y());
1700   gp_Pnt2d aP2sym (aP2.X(), -aP2.Y());
1701   gp_Lin2d aLin2dsym (aP1sym, gp_Vec2d(aP1sym,aP2sym));
1702   
1703   // intersect projections
1704   IntAna2d_Conic aCon    (aLin2d);
1705   IntAna2d_Conic aConSym (aLin2dsym);
1706   IntAna2d_AnaIntersection anIntersect;
1707   IntAna2d_AnaIntersection anIntersectSym;
1708   
1709   switch (theBasCurvProj.GetType()) {
1710   case GeomAbs_Line:
1711     anIntersectSym.Perform(theBasCurvProj.Line(), aConSym); 
1712     anIntersect.Perform(theBasCurvProj.Line(), aCon); break;
1713   case GeomAbs_Hyperbola:                         
1714     anIntersectSym.Perform(theBasCurvProj.Hyperbola(), aConSym); 
1715     anIntersect.Perform(theBasCurvProj.Hyperbola(), aCon); break;
1716   case GeomAbs_Parabola:                          
1717     anIntersectSym.Perform(theBasCurvProj.Parabola(), aConSym); 
1718     anIntersect.Perform(theBasCurvProj.Parabola(), aCon); break;
1719   default:
1720     return; // not infinite curve
1721   }
1722
1723   // retrieve params of intersections
1724   Standard_Integer aNbIntPnt = anIntersect.IsDone() ? anIntersect.NbPoints() : 0;
1725   Standard_Integer aNbIntPntSym = anIntersectSym.IsDone() ? anIntersectSym.NbPoints() : 0;
1726   Standard_Integer iPnt, aNbPnt = Max (aNbIntPnt, aNbIntPntSym);
1727   
1728   if (aNbPnt == 0) {
1729     theNoIntersection = Standard_True;
1730     return;
1731   }
1732   Standard_Real aParam;
1733   for (iPnt = 1; iPnt <= aNbPnt; iPnt++)
1734   {
1735     if (iPnt <= aNbIntPnt) {
1736       const IntAna2d_IntPoint& aIntPnt = anIntersect.Point(iPnt);
1737       aParam = aIntPnt.ParamOnFirst();
1738       theVmin = Min (theVmin, aParam );
1739       theVmax = Max (theVmax, aParam );
1740     }
1741     if (iPnt <= aNbIntPntSym) {
1742       const IntAna2d_IntPoint& aIntPnt = anIntersectSym.Point(iPnt);
1743       aParam = aIntPnt.ParamOnFirst();
1744       theVmin = Min (theVmin, aParam );
1745       theVmax = Max (theVmax, aParam );
1746     }
1747   }
1748
1749   return;
1750 }
1751 //=======================================================================
1752 //function : EstLimForInfRevl
1753 //purpose  : Estimate V1 and V2 to pass to InternalPerform() if they are
1754 //           infinite for Surface of Revolution
1755 //           Algo: intersect projections of Line and basis curve on the
1756 //           plane passing through revolution axe
1757 //=======================================================================
1758 void EstLimForInfRevl(const gp_Lin&   Line,
1759                       const TheSurface& surface, 
1760                       const Standard_Boolean U1inf, 
1761                       const Standard_Boolean U2inf, 
1762                       const Standard_Boolean V1inf, 
1763                       const Standard_Boolean V2inf, 
1764                       Standard_Real& U1new, 
1765                       Standard_Real& U2new, 
1766                       Standard_Real& V1new, 
1767                       Standard_Real& V2new, 
1768                       Standard_Boolean& NoIntersection)
1769 {
1770
1771   NoIntersection = Standard_False;
1772
1773   if (U1inf || U2inf) {
1774     if (U1inf)
1775       U1new = Max (0., U1new);
1776     else
1777       U2new = Min (2 * M_PI, U2new);
1778     if (! V1inf && !V2inf) return;
1779   }
1780
1781   Handle(Adaptor3d_HCurve) aBasisCurve = TheSurfaceTool::BasisCurve(surface);
1782   gp_Ax1 aRevAx = TheSurfaceTool::AxeOfRevolution(surface);
1783   gp_Vec aXVec = aRevAx.Direction();
1784   Standard_Real aTolAng = Precision::Angular();
1785
1786   // make plane to project a basis curve
1787   gp_Pnt O = aRevAx.Location();
1788   Standard_Real aU = 0.;
1789   gp_Pnt P = aBasisCurve->Value(aU);
1790   while (O.SquareDistance(P) <= Precision::PConfusion() ||
1791          aXVec.IsParallel( gp_Vec(O,P), aTolAng)) {
1792     aU += 1.;
1793     P = aBasisCurve->Value(aU);
1794     if (aU > 3)
1795       // basis curve is a line coinciding with aXVec, P is any not on aXVec
1796       P = gp_Pnt(aU, aU+1, aU+2);
1797   }
1798   gp_Vec aNVec = aXVec ^ gp_Vec(O,P);
1799   gp_Pln aPln (gp_Ax3 (O, aNVec ,aXVec));
1800
1801   // project basic curve
1802   ProjLib_Plane aBasCurvProj(aPln);
1803   switch (aBasisCurve->GetType()) {
1804   case GeomAbs_Line:
1805     aBasCurvProj.Project(aBasisCurve->Line     ()); break;
1806   case GeomAbs_Hyperbola:                       
1807     aBasCurvProj.Project(aBasisCurve->Hyperbola()); break;
1808   case GeomAbs_Parabola:                        
1809     aBasCurvProj.Project(aBasisCurve->Parabola ()); break;
1810   default:
1811     return; // not infinite curve
1812   }
1813   if (!aBasCurvProj.IsDone()) {
1814 #ifdef OCCT_DEBUG
1815     cout << "Info: IntCurveSurface_Inter::EstLimForInfRevl(), !aBasCurvProj.IsDone()" << endl;
1816 #endif
1817     return;
1818   }
1819   // make plane to project Line
1820   if (aXVec.IsParallel( Line.Direction(), aTolAng)) {
1821     P = Line.Location();
1822     while (O.SquareDistance(P) <= Precision::PConfusion()) {
1823       aU += 1.;
1824       P = gp_Pnt(aU, aU+1, aU+2); // any not on aXVec
1825     }
1826     aNVec = aXVec ^ gp_Vec(O,P);
1827   } else
1828     aNVec = aXVec.Crossed( Line.Direction() );
1829   
1830   aPln  = gp_Pln (gp_Ax3 (O, aNVec ,aXVec));
1831
1832   // make a second plane perpendicular to the first one, rotated around aXVec
1833   gp_Pln aPlnPrp = aPln.Rotated (gp_Ax1 (O,aXVec), M_PI/2.);
1834   
1835   // project Line and it's X-axe symmetric one to plane and intersect
1836   // resulting curve with projection of Basic Curev
1837   Standard_Real aVmin = RealLast(), aVmax = -aVmin;
1838   Standard_Boolean aNoInt1 = Standard_False, aNoInt2 = Standard_False;
1839   ProjectIntersectAndEstLim (Line, aPln,    aBasCurvProj, aVmin, aVmax, aNoInt1);
1840   ProjectIntersectAndEstLim (Line, aPlnPrp, aBasCurvProj, aVmin, aVmax, aNoInt2);
1841
1842   if (aNoInt1 && aNoInt2) {
1843     NoIntersection = Standard_True;
1844     return;
1845   }
1846   
1847   aVmin = aVmin - Abs(aVmin) - 10;
1848   aVmax = aVmax + Abs(aVmax) + 10;
1849
1850   if (V1inf) V1new = aVmin;
1851   if (V2inf) V2new = aVmax;
1852
1853   //cout << "EstLimForInfRevl: Vmin " << V1new << " Vmax " << V2new << endl;
1854   
1855   return;
1856 }
1857
1858 //=======================================================================
1859 //function : EstLimForInfOffs
1860 //purpose  : 
1861 //=======================================================================
1862 void EstLimForInfOffs(const gp_Lin&   Line,
1863                       const TheSurface& surface, 
1864                       const Standard_Integer nbsu, 
1865                       const Standard_Boolean U1inf, 
1866                       const Standard_Boolean U2inf, 
1867                       const Standard_Boolean V1inf, 
1868                       const Standard_Boolean V2inf, 
1869                       Standard_Real& U1new, 
1870                       Standard_Real& U2new, 
1871                       Standard_Real& V1new, 
1872                       Standard_Real& V2new, 
1873                       Standard_Boolean& NoIntersection)
1874 {
1875
1876   NoIntersection = Standard_False;
1877
1878   const Handle(Adaptor3d_HSurface)& aBasSurf = TheSurfaceTool::BasisSurface(surface);
1879   Standard_Real anOffVal = TheSurfaceTool::OffsetValue(surface);
1880
1881   GeomAbs_SurfaceType aTypeOfBasSurf = aBasSurf->GetType();
1882
1883   //  case for plane, cylinder and cone - make equivalent surface;
1884   if(aTypeOfBasSurf == GeomAbs_Plane) {
1885
1886     gp_Pln aPln = aBasSurf->Plane();
1887     gp_Vec aT = aPln.Position().XDirection()^aPln.Position().YDirection();
1888     aT *= anOffVal;
1889     aPln.Translate(aT);
1890     IntAna_IntConicQuad LinPlane(Line,aPln,TOLERANCE_ANGULAIRE);
1891
1892     if(!LinPlane.IsDone()) return;
1893
1894     if(LinPlane.IsParallel() || LinPlane.IsInQuadric()) {
1895
1896       NoIntersection = Standard_True;
1897       return;
1898
1899     }
1900     
1901     Standard_Real u, v;
1902     ElSLib::Parameters(aPln, LinPlane.Point(1), u, v);
1903     U1new = Max(U1new, u - 10.);
1904     U2new = Min(U2new, u + 10.);
1905     V1new = Max(V1new, v - 10.);
1906     V2new = Min(V2new, v + 10.);
1907
1908   }
1909   else if(aTypeOfBasSurf == GeomAbs_Cylinder) {
1910
1911     gp_Cylinder aCyl = aBasSurf->Cylinder();
1912
1913     Standard_Real aR = aCyl.Radius();
1914     gp_Ax3 anA = aCyl.Position();
1915
1916     if (anA.Direct()) 
1917       aR += anOffVal;
1918     else 
1919       aR -= anOffVal;
1920
1921     if ( aR >= TOLTANGENCY ) {
1922       aCyl.SetRadius(aR);
1923     }
1924     else if ( aR <= -TOLTANGENCY ){
1925       anA.Rotate(gp_Ax1(anA.Location(), anA.Direction()), M_PI);
1926       aCyl.SetPosition(anA);
1927 // modified by NIZHNY-MKK  Mon Oct  3 17:37:54 2005
1928 //       aCyl.SetRadius(Abs(aR));
1929       aCyl.SetRadius(-aR);
1930     }
1931     else {
1932
1933       NoIntersection = Standard_True;
1934       return;
1935       
1936     }
1937
1938     IntAna_IntConicQuad LinCylinder(Line, aCyl);
1939
1940     if(!LinCylinder.IsDone()) return;
1941
1942     if(LinCylinder.IsParallel() || LinCylinder.IsInQuadric()) {
1943
1944       NoIntersection = Standard_True;
1945       return;
1946
1947     }
1948     
1949     Standard_Integer i, nbp = LinCylinder.NbPoints();
1950     Standard_Real vmin = RealLast(), vmax = -vmin, u, v;
1951
1952     for(i = 1; i <= nbp; i++) {
1953
1954       ElSLib::Parameters(aCyl, LinCylinder.Point(i), u, v);
1955       vmin = Min(vmin, v);
1956       vmax = Max(vmax, v);
1957
1958     } 
1959
1960     V1new = Max(V1new, vmin - Abs(vmin) - 10.);
1961     V2new = Min(V2new, vmax + Abs(vmax) + 10.);
1962
1963   }
1964   else if(aTypeOfBasSurf == GeomAbs_Cone) {
1965
1966     gp_Cone aCon = aBasSurf->Cone();
1967     Standard_Real anAng = aCon.SemiAngle();
1968     Standard_Real aR = aCon.RefRadius() + anOffVal * Cos(anAng);
1969     gp_Ax3 anA = aCon.Position();
1970     if ( aR >= 0.) {
1971
1972       gp_Vec aZ( anA.Direction());
1973       aZ *= - anOffVal * Sin(anAng);
1974       anA.Translate(aZ);
1975       aCon.SetPosition(anA);
1976       aCon.SetRadius(aR);
1977       aCon.SetSemiAngle(anAng);
1978
1979     }
1980     else {
1981
1982       return;
1983
1984     }
1985
1986     IntAna_IntConicQuad LinCone(Line, aCon);
1987
1988     if(!LinCone.IsDone()) return;
1989
1990     if(LinCone.IsParallel() || LinCone.IsInQuadric()) {
1991
1992       NoIntersection = Standard_True;
1993       return;
1994
1995     }
1996     
1997     Standard_Integer i, nbp = LinCone.NbPoints();
1998     Standard_Real vmin = RealLast(), vmax = -vmin, u, v;
1999
2000     for(i = 1; i <= nbp; i++) {
2001
2002       ElSLib::Parameters(aCon, LinCone.Point(i), u, v);
2003       vmin = Min(vmin, v);
2004       vmax = Max(vmax, v);
2005
2006     } 
2007
2008     V1new = Max(V1new, vmin - Abs(vmin) - 10.);
2009     V2new = Min(V2new, vmax + Abs(vmax) + 10.);
2010
2011   }
2012   else if(aTypeOfBasSurf == GeomAbs_SurfaceOfExtrusion) {
2013
2014     Standard_Real anU1 = U1new, anU2 = U2new;
2015
2016     EstLimForInfExtr(Line, surface, Standard_True, nbsu, 
2017                      U1inf, U2inf, V1inf, V2inf, 
2018                      U1new, U2new, V1new, V2new, NoIntersection);
2019
2020     if(NoIntersection) return;
2021
2022     if(U1inf || U2inf) {
2023
2024       GeomAbs_CurveType aBasCurvType = aBasSurf->BasisCurve()->GetType();
2025       if(aBasCurvType == GeomAbs_Line) {
2026         U1new = Max(anU1, -1.e10);
2027         U2new = Min(anU2,  1.e10);
2028       } 
2029       else if(aBasCurvType == GeomAbs_Parabola) {
2030         gp_Parab aPrb  = aBasSurf->BasisCurve()->Parabola();
2031         Standard_Real aF = aPrb.Focal();
2032         Standard_Real dU = 2.e5 * Sqrt(aF);
2033         U1new = Max(anU1, -dU);
2034         U2new = Min(anU2,  dU);
2035       }
2036       else if(aBasCurvType == GeomAbs_Hyperbola) {
2037         U1new = Max(anU1, -30.);
2038         U2new = Min(anU2,  30.);
2039       }
2040       else {
2041         U1new = Max(anU1, -1.e10);
2042         U2new = Min(anU2,  1.e10);
2043       } 
2044
2045     }
2046
2047   }
2048   else if(aTypeOfBasSurf == GeomAbs_SurfaceOfRevolution) {
2049
2050     GeomAbs_CurveType aBasCurvType = aBasSurf->BasisCurve()->GetType();
2051     if(aBasCurvType == GeomAbs_Line) {
2052       V1new = Max(V1new, -1.e10);
2053       V2new = Min(V2new,  1.e10);
2054     }   
2055     else if(aBasCurvType == GeomAbs_Parabola) {
2056       gp_Parab aPrb  = aBasSurf->BasisCurve()->Parabola();
2057       Standard_Real aF = aPrb.Focal();
2058       Standard_Real dV = 2.e5 * Sqrt(aF);
2059       V1new = Max(V1new, -dV);
2060       V2new = Min(V2new,  dV);
2061     }
2062     else if(aBasCurvType == GeomAbs_Hyperbola) {
2063      V1new  = Max(V1new, -30.);
2064      V2new  = Min(V2new,  30.);
2065     }
2066     else {
2067      V1new  = Max(V1new, -1.e10);
2068      V2new  = Min(V2new,  1.e10);
2069     }   
2070
2071   }
2072   else {
2073
2074     V1new = Max(V1new, -1.e10);
2075     V2new = Min(V2new,  1.e10);
2076   }
2077 }
2078 //=======================================================================
2079 //function : EstLimForInfSurf
2080 //purpose  : 
2081 //=======================================================================
2082 void EstLimForInfSurf(Standard_Real& U1new, 
2083                       Standard_Real& U2new, 
2084                       Standard_Real& V1new, 
2085                       Standard_Real& V2new)
2086 {
2087     U1new = Max(U1new, -1.e10);
2088     U2new = Min(U2new,  1.e10);
2089     V1new = Max(V1new, -1.e10);
2090     V2new = Min(V2new,  1.e10);
2091 }