0025104: Prism from BSpline curve can not be chamfered
[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   Standard_Boolean isAnaProcessed = Standard_True;
920   switch(SurfaceType) { 
921   case GeomAbs_Plane: 
922     {
923       IntAna_IntConicQuad LinPlane(Line,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE);
924       AppendIntAna(curve,surface,LinPlane);
925       break;
926     }
927   case GeomAbs_Cylinder:
928     {
929       IntAna_IntConicQuad LinCylinder(Line,TheSurfaceTool::Cylinder(surface));
930       AppendIntAna(curve,surface,LinCylinder);
931       break;
932     }
933   case GeomAbs_Sphere:
934     {
935       IntAna_IntConicQuad LinSphere(Line,TheSurfaceTool::Sphere(surface));
936       AppendIntAna(curve,surface,LinSphere);
937       break;
938     }
939   case GeomAbs_Torus: 
940     {
941       IntAna_IntLinTorus intlintorus(Line, TheSurfaceTool::Torus(surface));
942       if (intlintorus.IsDone()) {
943         Standard_Integer nbp = intlintorus.NbPoints();
944         Standard_Real fi, theta, w;
945         for (Standard_Integer i = 1; i <= nbp; i++) {
946           const gp_Pnt aDebPnt(intlintorus.Value(i));
947           (void)aDebPnt;
948           w = intlintorus.ParamOnLine(i);
949           intlintorus.ParamOnTorus(i, fi, theta);
950           AppendPoint(curve, w, surface, fi, theta);
951         }
952       }
953       else
954         isAnaProcessed = Standard_False;
955       break;
956     }
957   case GeomAbs_Cone:
958     {
959       const Standard_Real correction = 1.E+5*Precision::Angular();
960       gp_Cone cn = TheSurfaceTool::Cone(surface);
961       if(Abs(cn.SemiAngle()) < M_PI/2.0 - correction) {
962         IntAna_IntConicQuad LinCone(Line, cn);
963         AppendIntAna(curve, surface, LinCone);
964       }
965       else
966         isAnaProcessed = Standard_False;
967       break;
968     }
969   default:
970     isAnaProcessed = Standard_False;
971   }
972   if (!isAnaProcessed)
973   {
974     Standard_Integer nbsu,nbsv;
975     nbsu = TheSurfaceTool::NbSamplesU(surface,U1,U2);
976     nbsv = TheSurfaceTool::NbSamplesV(surface,V1,V2);
977
978     Standard_Boolean U1inf = Precision::IsInfinite(U1);
979     Standard_Boolean U2inf = Precision::IsInfinite(U2);
980     Standard_Boolean V1inf = Precision::IsInfinite(V1);
981     Standard_Boolean V2inf = Precision::IsInfinite(V2);
982
983     Standard_Real U1new=U1, U2new=U2, V1new=V1, V2new=V2;
984
985     Standard_Boolean NoIntersection = Standard_False;
986
987     if(U1inf || U2inf || V1inf || V2inf ) {
988
989       if(SurfaceType == GeomAbs_SurfaceOfExtrusion) {
990
991         EstLimForInfExtr(Line, surface, Standard_False, nbsu,
992                          U1inf, U2inf, V1inf, V2inf,
993                          U1new, U2new, V1new, V2new, NoIntersection);
994
995       }
996       else if (SurfaceType == GeomAbs_SurfaceOfRevolution) {
997
998         EstLimForInfRevl(Line, surface,
999                          U1inf, U2inf, V1inf, V2inf,
1000                          U1new, U2new, V1new, V2new, NoIntersection);
1001
1002       }
1003       else if (SurfaceType == GeomAbs_OffsetSurface) {
1004
1005         EstLimForInfOffs(Line, surface, nbsu,
1006                          U1inf, U2inf, V1inf, V2inf,
1007                          U1new, U2new, V1new, V2new, NoIntersection);
1008       }
1009       else {
1010
1011         EstLimForInfSurf(U1new, U2new, V1new, V2new);
1012       }
1013
1014     }
1015
1016
1017     if(NoIntersection) return;
1018
1019
1020     // modified by NIZHNY-OFV  Mon Aug 20 14:56:47 2001 (60963 begin)
1021     if(nbsu<20) nbsu=20;
1022     if(nbsv<20) nbsv=20;
1023     // modified by NIZHNY-OFV  Mon Aug 20 14:57:06 2001 (60963 end)
1024     IntCurveSurface_ThePolyhedron polyhedron(surface,nbsu,nbsv,U1new,V1new,U2new,V2new);
1025     Intf_Tool                     bndTool;
1026     Bnd_Box                       boxLine;
1027     bndTool.LinBox(Line,polyhedron.Bounding(),boxLine);
1028     for(Standard_Integer nbseg=1; nbseg<= bndTool.NbSegments(); nbseg++) { 
1029       Standard_Real pinf = bndTool.BeginParam(nbseg);
1030       Standard_Real psup = bndTool.EndParam(nbseg);
1031       if((psup - pinf)<1e-10) { pinf-=1e-10; psup+=1e-10; } 
1032       IntCurveSurface_ThePolygon polygon(curve, pinf,psup,2);
1033       InternalPerform(curve,polygon,surface,polyhedron,U1new,V1new,U2new,V2new);
1034     }
1035   }
1036 }
1037 //=======================================================================
1038 //function : PerformConicSurf
1039 //purpose  : 
1040 //=======================================================================
1041 void IntCurveSurface_Inter::PerformConicSurf(const gp_Circ&        Circle,
1042                                              const TheCurve&       curve,
1043                                              const TheSurface&     surface,
1044                                              const Standard_Real U1,
1045                                              const Standard_Real V1,
1046                                              const Standard_Real U2,
1047                                              const Standard_Real V2) { 
1048   
1049   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1050   switch(SurfaceType) { 
1051   case GeomAbs_Plane: 
1052     {
1053       IntAna_IntConicQuad CircPlane(Circle,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE,TOLERANCE);
1054       AppendIntAna(curve,surface,CircPlane);
1055       break;
1056     }
1057   case GeomAbs_Cylinder:
1058     {
1059       IntAna_IntConicQuad CircCylinder(Circle,TheSurfaceTool::Cylinder(surface));
1060       AppendIntAna(curve,surface,CircCylinder);
1061       break;
1062     }
1063   case GeomAbs_Cone:
1064     {
1065       IntAna_IntConicQuad CircCone(Circle,TheSurfaceTool::Cone(surface));
1066       AppendIntAna(curve,surface,CircCone);
1067       break;
1068     }   
1069   case GeomAbs_Sphere:
1070     {
1071       IntAna_IntConicQuad CircSphere(Circle,TheSurfaceTool::Sphere(surface));
1072       AppendIntAna(curve,surface,CircSphere);
1073       break;
1074     }
1075   default: 
1076     {
1077       IntCurveSurface_ThePolygon polygon(curve,NBSAMPLESONCIRCLE);
1078       InternalPerform(curve,polygon,surface,U1,V1,U2,V2);
1079     }
1080   }
1081 }
1082 //=======================================================================
1083 //function : PerformConicSurf
1084 //purpose  : 
1085 //=======================================================================
1086 void IntCurveSurface_Inter::PerformConicSurf(const gp_Elips&       Ellipse,
1087                                              const TheCurve&       curve,
1088                                              const TheSurface&     surface,
1089                                              const Standard_Real U1,
1090                                              const Standard_Real V1,
1091                                              const Standard_Real U2,
1092                                              const Standard_Real V2) { 
1093
1094   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1095   switch(SurfaceType) { 
1096   case GeomAbs_Plane: 
1097     {
1098       IntAna_IntConicQuad EllipsePlane(Ellipse,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE,TOLERANCE);
1099       AppendIntAna(curve,surface,EllipsePlane);
1100       break;
1101     }
1102   case GeomAbs_Cylinder:
1103     {
1104       IntAna_IntConicQuad EllipseCylinder(Ellipse,TheSurfaceTool::Cylinder(surface));
1105       AppendIntAna(curve,surface,EllipseCylinder);
1106       break;
1107     }
1108   case GeomAbs_Cone:
1109      {
1110       IntAna_IntConicQuad EllipseCone(Ellipse,TheSurfaceTool::Cone(surface));
1111       AppendIntAna(curve,surface,EllipseCone);
1112       break;
1113     }   
1114   case GeomAbs_Sphere:
1115     {
1116       IntAna_IntConicQuad EllipseSphere(Ellipse,TheSurfaceTool::Sphere(surface));
1117       AppendIntAna(curve,surface,EllipseSphere);
1118       break;
1119     }
1120   default: 
1121     {
1122       IntCurveSurface_ThePolygon polygon(curve,NBSAMPLESONELLIPSE);
1123       InternalPerform(curve,polygon,surface,U1,V1,U2,V2);
1124     }
1125   }
1126 }
1127 //=======================================================================
1128 //function : PerformConicSurf
1129 //purpose  : 
1130 //=======================================================================
1131 void IntCurveSurface_Inter::PerformConicSurf(const gp_Parab&       Parab,
1132                                              const TheCurve&       curve,
1133                                              const TheSurface&     surface,
1134                                              const Standard_Real U1,
1135                                              const Standard_Real V1,
1136                                              const Standard_Real U2,
1137                                              const Standard_Real V2) { 
1138
1139   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1140   switch(SurfaceType) { 
1141   case GeomAbs_Plane: 
1142     {
1143       IntAna_IntConicQuad ParabPlane(Parab,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE);
1144       AppendIntAna(curve,surface,ParabPlane);
1145       break;
1146     }
1147   case GeomAbs_Cylinder:
1148     {
1149       IntAna_IntConicQuad ParabCylinder(Parab,TheSurfaceTool::Cylinder(surface));
1150       AppendIntAna(curve,surface,ParabCylinder);
1151       break;
1152     }
1153   case GeomAbs_Cone:
1154      {
1155       IntAna_IntConicQuad ParabCone(Parab,TheSurfaceTool::Cone(surface));
1156       AppendIntAna(curve,surface,ParabCone);
1157       break;
1158     }   
1159   case GeomAbs_Sphere:
1160     {
1161       IntAna_IntConicQuad ParabSphere(Parab,TheSurfaceTool::Sphere(surface));
1162       AppendIntAna(curve,surface,ParabSphere);
1163       break;
1164     }
1165   default: 
1166     {
1167       Standard_Integer nbsu,nbsv;
1168       nbsu = TheSurfaceTool::NbSamplesU(surface,U1,U2);
1169       nbsv = TheSurfaceTool::NbSamplesV(surface,V1,V2);
1170       if(nbsu>40) nbsu=40;
1171       if(nbsv>40) nbsv=40;
1172       IntCurveSurface_ThePolyhedron polyhedron(surface,nbsu,nbsv,U1,V1,U2,V2);
1173       Intf_Tool                         bndTool;
1174       Bnd_Box                          boxParab;
1175       bndTool.ParabBox(Parab,polyhedron.Bounding(),boxParab);
1176       for(Standard_Integer nbseg=1; nbseg<= bndTool.NbSegments(); nbseg++) { 
1177         IntCurveSurface_ThePolygon polygon(curve,
1178                                            bndTool.BeginParam(nbseg),
1179                                            bndTool.EndParam(nbseg),
1180                                            NBSAMPLESONPARAB);
1181         InternalPerform(curve,polygon,surface,polyhedron,U1,V1,U2,V2);
1182       }
1183     }
1184   }
1185 }
1186 //=======================================================================
1187 //function : PerformConicSurf
1188 //purpose  : 
1189 //=======================================================================
1190 void IntCurveSurface_Inter::PerformConicSurf(const gp_Hypr&        Hypr,
1191                                              const TheCurve&       curve,
1192                                              const TheSurface&     surface,
1193                                              const Standard_Real U1,
1194                                              const Standard_Real V1,
1195                                              const Standard_Real U2,
1196                                              const Standard_Real V2) {
1197  
1198   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1199   switch(SurfaceType) { 
1200   case GeomAbs_Plane: 
1201     {
1202       IntAna_IntConicQuad HyprPlane(Hypr,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE);
1203       AppendIntAna(curve,surface,HyprPlane);
1204       break;
1205     }
1206   case GeomAbs_Cylinder:
1207     {
1208       IntAna_IntConicQuad HyprCylinder(Hypr,TheSurfaceTool::Cylinder(surface));
1209       AppendIntAna(curve,surface,HyprCylinder);
1210       break;
1211     }
1212   case GeomAbs_Cone:
1213      {
1214       IntAna_IntConicQuad HyprCone(Hypr,TheSurfaceTool::Cone(surface));
1215       AppendIntAna(curve,surface,HyprCone);
1216       break;
1217     }   
1218   case GeomAbs_Sphere:
1219     {
1220       IntAna_IntConicQuad HyprSphere(Hypr,TheSurfaceTool::Sphere(surface));
1221       AppendIntAna(curve,surface,HyprSphere);
1222       break;
1223     }
1224   default: 
1225     {
1226       Standard_Integer nbsu,nbsv;
1227       nbsu = TheSurfaceTool::NbSamplesU(surface,U1,U2);
1228       nbsv = TheSurfaceTool::NbSamplesV(surface,V1,V2);
1229       if(nbsu>40) nbsu=40;
1230       if(nbsv>40) nbsv=40;
1231       IntCurveSurface_ThePolyhedron polyhedron(surface,nbsu,nbsv,U1,V1,U2,V2);
1232       Intf_Tool                         bndTool;
1233       Bnd_Box                          boxHypr;
1234       bndTool.HyprBox(Hypr,polyhedron.Bounding(),boxHypr);
1235       for(Standard_Integer nbseg=1; nbseg<= bndTool.NbSegments(); nbseg++) { 
1236         IntCurveSurface_ThePolygon polygon(curve,
1237                                            bndTool.BeginParam(nbseg),
1238                                            bndTool.EndParam(nbseg),
1239                                            NBSAMPLESONHYPR);
1240         InternalPerform(curve,polygon,surface,polyhedron,U1,V1,U2,V2);
1241       }
1242     }
1243   }
1244 }
1245 //=======================================================================
1246 //function : AppendIntAna
1247 //purpose  : 
1248 //=======================================================================
1249 void IntCurveSurface_Inter::AppendIntAna(const TheCurve& curve,
1250                                          const TheSurface& surface,
1251                                          const IntAna_IntConicQuad& intana_ConicQuad) {
1252   if(intana_ConicQuad.IsDone()) { 
1253     if(intana_ConicQuad.IsInQuadric()) { 
1254       //-- cout<<" Courbe Dans la Quadrique !!! Non Traite !!!"<<endl;
1255     }
1256     else if(intana_ConicQuad.IsParallel()) { 
1257       //-- cout<<" Courbe // a la Quadrique !!! Non Traite !!!"<<endl;
1258     }
1259     else { 
1260       Standard_Integer nbp = intana_ConicQuad.NbPoints();
1261       Standard_Real u,v,w;
1262       for(Standard_Integer i = 1; i<= nbp; i++) { 
1263         gp_Pnt P(intana_ConicQuad.Point(i));
1264         w = intana_ConicQuad.ParamOnConic(i);
1265         IntCurveSurface_ComputeParamsOnQuadric(surface,P,u,v);
1266         AppendPoint(curve,w,surface,u,v);
1267       }
1268     }
1269   }
1270   else { 
1271     //-- cout<<" IntAna Conic Quad Not Done  "<<endl;
1272   }
1273 }
1274 //=======================================================================
1275 //function : AppendPoint
1276 //purpose  : 
1277 //=======================================================================
1278 void IntCurveSurface_Inter::AppendPoint(const TheCurve& curve,
1279                                         const Standard_Real lw,
1280                                         const TheSurface& surface,
1281                                         const Standard_Real su,
1282                                         const Standard_Real sv) {
1283  
1284   Standard_Real W0 = TheCurveTool::FirstParameter(curve);
1285   Standard_Real W1 = TheCurveTool::LastParameter(curve);
1286   Standard_Real U0 = TheSurfaceTool::FirstUParameter(surface);
1287   Standard_Real U1 = TheSurfaceTool::LastUParameter(surface);
1288   Standard_Real V0 = TheSurfaceTool::FirstVParameter(surface);
1289   Standard_Real V1 = TheSurfaceTool::LastVParameter(surface);
1290   //-- Test si la courbe est periodique
1291   Standard_Real w = lw, u = su, v = sv;
1292
1293   GeomAbs_CurveType aCType = TheCurveTool::GetType(curve);
1294
1295   if(TheCurveTool::IsPeriodic(curve) 
1296      || aCType == GeomAbs_Circle 
1297      || aCType == GeomAbs_Ellipse) {
1298     w = ElCLib::InPeriod(w, W0, W0 + TheCurveTool::Period(curve));
1299   }
1300
1301   if((W0 - w) >= TOLTANGENCY || (w - W1) >= TOLTANGENCY) return;
1302
1303   GeomAbs_SurfaceType aSType = TheSurfaceTool::GetType(surface);
1304   if (TheSurfaceTool::IsUPeriodic(surface)
1305       || aSType == GeomAbs_Cylinder
1306       || aSType == GeomAbs_Cone
1307       || aSType == GeomAbs_Sphere) {
1308     u = ElCLib::InPeriod(u, U0, U0 + TheSurfaceTool::UPeriod(surface));
1309   }
1310
1311   if (TheSurfaceTool::IsVPeriodic(surface)) {
1312     v = ElCLib::InPeriod(v, V0, V0 + TheSurfaceTool::VPeriod(surface));
1313   }
1314
1315   if((U0 - u) >= TOLTANGENCY || (u - U1) >= TOLTANGENCY) return;
1316   if((V0 - v) >= TOLTANGENCY || (v - V1) >= TOLTANGENCY) return;
1317
1318   IntCurveSurface_TransitionOnCurve   TransOnCurve;
1319   IntCurveSurface_ComputeTransitions(curve,w,TransOnCurve,
1320                                      surface,u,v);
1321   gp_Pnt P(TheCurveTool::Value(curve,w));
1322   IntCurveSurface_IntersectionPoint IP(P,u,v,w,TransOnCurve);
1323   Append(IP); //-- invoque la methode de IntCurveSurface_Intersection.
1324   
1325 }
1326
1327 //=======================================================================
1328 //function : AppendSegment
1329 //purpose  : 
1330 //=======================================================================
1331 void IntCurveSurface_Inter::AppendSegment(const TheCurve& ,
1332                                           const Standard_Real ,
1333                                           const Standard_Real ,
1334                                           const TheSurface& ) { 
1335   //cout<<" !!! Not Yet Implemented 
1336   //IntCurveSurface_Inter::Append(const IntCurveSurf ...)"<<endl;
1337 }
1338
1339 //=======================================================================
1340 //function : SectionPointToParameters
1341 //purpose  : P o i n t   d   i n t e r f e r e n c e   - - >  
1342 //    U , V       e t   W 
1343 //=======================================================================
1344 void SectionPointToParameters(const Intf_SectionPoint&              Sp,
1345                               const IntCurveSurface_ThePolyhedron&  Polyhedron,
1346                               const IntCurveSurface_ThePolygon&     Polygon,
1347                               Standard_Real& U,
1348                               Standard_Real& V,
1349                               Standard_Real& W) {
1350   
1351   Intf_PIType       typ;
1352   Standard_Integer  Adr1,Adr2;
1353   Standard_Real     Param,u,v;
1354   gp_Pnt P(Sp.Pnt());
1355   
1356   Standard_Integer Pt1,Pt2,Pt3;
1357   Standard_Real u1 = 0.,v1 = 0.,param;
1358   //----------------------------------------------------------------------
1359   //--          Calcul des parametres approches sur la surface          --
1360   //----------------------------------------------------------------------
1361   
1362   Sp.InfoSecond(typ,Adr1,Adr2,Param);
1363   switch(typ) { 
1364   case Intf_VERTEX:   //-- Adr1 est le numero du vertex
1365     {
1366       Polyhedron.Parameters(Adr1,u1,v1);
1367       break;
1368     }
1369   case Intf_EDGE:
1370     {
1371       Polyhedron.Parameters(Adr1,u1,v1);    
1372       Polyhedron.Parameters(Adr2,u,v);
1373       u1+= Param * (u-u1);
1374       v1+= Param * (v-v1);
1375       break;
1376     }
1377   case Intf_FACE:
1378     {
1379       Standard_Real ua,va,ub,vb,uc,vc,ca,cb,cc,cabc;
1380       Polyhedron.Triangle(Adr1,Pt1,Pt2,Pt3);
1381       gp_Pnt PA(Polyhedron.Point(Pt1));
1382       gp_Pnt PB(Polyhedron.Point(Pt2));
1383       gp_Pnt PC(Polyhedron.Point(Pt3));
1384       Polyhedron.Parameters(Pt1,ua,va);
1385       Polyhedron.Parameters(Pt2,ub,vb);
1386       Polyhedron.Parameters(Pt3,uc,vc);
1387       gp_Vec Normale(gp_Vec(PA,PB).Crossed(gp_Vec(PA,PC)));
1388       cc = (gp_Vec(PA,PB).Crossed(gp_Vec(PA,P))).Dot(Normale);
1389       ca = (gp_Vec(PB,PC).Crossed(gp_Vec(PB,P))).Dot(Normale);
1390       cb = (gp_Vec(PC,PA).Crossed(gp_Vec(PC,P))).Dot(Normale);
1391       cabc = ca + cb + cc;
1392       
1393       ca/=cabc;     cb/=cabc;       cc/=cabc;
1394       
1395       u1 = ca * ua + cb * ub + cc * uc;
1396       v1 = ca * va + cb * vb + cc * vc;
1397       break;
1398     }
1399   default: 
1400     {
1401       cout<<" Default dans SectionPointToParameters "<<endl;
1402       break;
1403     }
1404   }
1405   //----------------------------------------------------------------------
1406   //--                Calcul du point approche sur la Curve             --
1407   //----------------------------------------------------------------------
1408   Standard_Integer SegIndex;
1409
1410   Sp.InfoFirst(typ,SegIndex,param);
1411   W = Polygon.ApproxParamOnCurve(SegIndex,param);
1412   if(param>1.0 || param<0.0) { 
1413     //-- IntCurveSurface_ThePolyhedronTool::Dump(Polyhedron);
1414     //-- IntCurveSurface_ThePolygonTool::Dump(Polygon);
1415     }
1416   U = u1;
1417   V = v1;
1418 }
1419 //=======================================================================
1420 //function : IntCurveSurface_ComputeTransitions
1421 //purpose  : 
1422 //=======================================================================
1423 void IntCurveSurface_ComputeTransitions(const TheCurve& curve,
1424                                         const Standard_Real w,
1425                                         IntCurveSurface_TransitionOnCurve&   TransOnCurve,
1426                                         const TheSurface& surface,
1427                                         const Standard_Real u,
1428                                         const Standard_Real v) { 
1429   
1430   gp_Vec NSurf,D1U,D1V;//TgCurv;
1431   gp_Pnt Psurf;
1432   Standard_Real CosDir;
1433   
1434   
1435   TheSurfaceTool::D1(surface,u,v,Psurf,D1U,D1V);
1436   NSurf = D1U.Crossed(D1V);
1437   TheCurveTool::D1(curve,w,Psurf,D1U);
1438   Standard_Real Norm = NSurf.Magnitude();
1439   if(Norm>TOLERANCE_ANGULAIRE &&
1440      D1U.SquareMagnitude() > TOLERANCE_ANGULAIRE) { 
1441     D1U.Normalize();
1442     CosDir = NSurf.Dot(D1U);
1443     CosDir/=Norm;
1444     if( -CosDir > TOLERANCE_ANGULAIRE) { 
1445       //--  --Curve--->    <----Surface----
1446       TransOnCurve   = IntCurveSurface_In;
1447     }
1448     else if(CosDir > TOLERANCE_ANGULAIRE) { 
1449       //--  --Curve--->  ----Surface-->
1450       TransOnCurve   = IntCurveSurface_Out;
1451     }     
1452     else { 
1453       TransOnCurve   = IntCurveSurface_Tangent;
1454     }
1455   }
1456   else { 
1457     TransOnCurve   = IntCurveSurface_Tangent;
1458   }
1459 }
1460 //=======================================================================
1461 //function : IntCurveSurface_ComputeParamsOnQuadric
1462 //purpose  : 
1463 //=======================================================================
1464 void IntCurveSurface_ComputeParamsOnQuadric(const TheSurface& surface,
1465                                             const gp_Pnt& P,
1466                                             Standard_Real& u,
1467                                             Standard_Real& v) { 
1468   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1469   switch(SurfaceType) { 
1470   case GeomAbs_Plane: 
1471     {
1472       ElSLib::Parameters(TheSurfaceTool::Plane(surface),P,u,v);
1473         break;
1474       }
1475   case GeomAbs_Cylinder:
1476     {
1477       ElSLib::Parameters(TheSurfaceTool::Cylinder(surface),P,u,v);
1478         break;
1479       }
1480   case GeomAbs_Cone:
1481     {
1482       ElSLib::Parameters(TheSurfaceTool::Cone(surface),P,u,v);
1483         break;
1484       }   
1485   case GeomAbs_Sphere:
1486     {
1487       ElSLib::Parameters(TheSurfaceTool::Sphere(surface),P,u,v);
1488         break;
1489       }
1490   default:
1491     break;
1492   }
1493 }
1494 //=======================================================================
1495 //function : EstLimForInfExtr
1496 //purpose  : Estimation of limits for infinite surfaces
1497 //=======================================================================
1498 void EstLimForInfExtr(const gp_Lin&   Line,
1499                       const TheSurface& surface, 
1500                       const Standard_Boolean IsOffSurf,
1501                       const Standard_Integer nbsu, 
1502                       const Standard_Boolean U1inf, 
1503                       const Standard_Boolean U2inf, 
1504                       const Standard_Boolean V1inf, 
1505                       const Standard_Boolean V2inf, 
1506                       Standard_Real& U1new, 
1507                       Standard_Real& U2new, 
1508                       Standard_Real& V1new, 
1509                       Standard_Real& V2new, 
1510                       Standard_Boolean& NoIntersection)
1511
1512 {
1513
1514   NoIntersection = Standard_False;
1515
1516   Handle(Adaptor3d_HSurface) aBasSurf; 
1517
1518   if(IsOffSurf) aBasSurf = TheSurfaceTool::BasisSurface(surface);
1519
1520   gp_Dir aDirOfExt;
1521
1522   if(IsOffSurf) aDirOfExt = aBasSurf->Direction();
1523   else aDirOfExt = TheSurfaceTool::Direction(surface);
1524
1525   Standard_Real tolang = TOLERANCE_ANGULAIRE;
1526
1527   if(aDirOfExt.IsParallel(Line.Direction(), tolang)) {
1528     NoIntersection = Standard_True;
1529     return;
1530   }
1531
1532   if((V1inf || V2inf) && !(U1inf || U2inf)) {
1533
1534     Standard_Real vmin = RealLast(), vmax = -vmin;
1535     gp_Lin aL;
1536     Standard_Real step = (U2new - U1new) / nbsu;
1537     Standard_Real u = U1new, v;
1538     gp_Pnt aP;
1539     Extrema_POnCurv aP1, aP2;
1540     Standard_Integer i;
1541
1542     for(i = 0; i <= nbsu; i++) {
1543
1544       TheSurfaceTool::D0(surface, u, 0., aP);
1545       aL.SetLocation(aP);
1546       aL.SetDirection(aDirOfExt);
1547
1548       Extrema_ExtElC aExtr(aL, Line, tolang);
1549       
1550       if(!aExtr.IsDone()) return;
1551
1552       if(aExtr.IsParallel()) {
1553         NoIntersection = Standard_True;
1554         return;
1555       }
1556
1557       aExtr.Points(1, aP1, aP2);
1558       v = aP1.Parameter();
1559       vmin = Min(vmin, v);
1560       vmax = Max(vmax, v);
1561
1562       u += step;
1563
1564     }
1565
1566     vmin = vmin - Abs(vmin) - 10.;
1567     vmax = vmax + Abs(vmax) + 10.;
1568
1569     V1new = Max(V1new, vmin);
1570     V2new = Min(V2new, vmax);
1571
1572   }
1573   else if(U1inf || U2inf) {
1574
1575     Standard_Real umin = RealLast(), umax = -umin;
1576     Standard_Real u0 = Min(Max(0., U1new), U2new);
1577     Standard_Real v0 = Min(Max(0., V1new), V2new);
1578     gp_Pnt aP;
1579     TheSurfaceTool::D0(surface, u0, v0, aP);
1580     gp_Pln aRefPln(aP, aDirOfExt);
1581
1582     Handle(Adaptor3d_HCurve) aBasCurv;
1583
1584     if(IsOffSurf) aBasCurv = aBasSurf->BasisCurve();
1585     else aBasCurv = TheSurfaceTool::BasisCurve(surface);
1586
1587     ProjLib_Plane Projector(aRefPln);
1588
1589     Projector.Project(Line);
1590
1591     if(!Projector.IsDone()) return;
1592
1593     gp_Lin2d Line2d = Projector.Line();
1594
1595     GeomAbs_CurveType aCurvTyp = aBasCurv->GetType();
1596
1597     if(aCurvTyp == GeomAbs_Line) {
1598       
1599       Projector.Project(aBasCurv->Line());
1600
1601       if(!Projector.IsDone()) return;
1602
1603       gp_Lin2d aL2d = Projector.Line();
1604
1605       IntAna2d_AnaIntersection anInter(Line2d, aL2d);
1606
1607       if(!anInter.IsDone()) return;
1608
1609       if(anInter.IsEmpty() || anInter.IdenticalElements() || 
1610                               anInter.ParallelElements()     ) {
1611         NoIntersection = Standard_True;
1612         return;
1613       }
1614
1615       const IntAna2d_IntPoint& anIntPnt = anInter.Point(1);
1616       umin = umax = anIntPnt.ParamOnSecond();
1617     }
1618     else if(aCurvTyp == GeomAbs_Parabola || aCurvTyp == GeomAbs_Hyperbola) {
1619
1620       IntAna2d_Conic aCon(Line2d);
1621       IntAna2d_AnaIntersection anInter;
1622       
1623       if(aCurvTyp == GeomAbs_Parabola) {
1624         Projector.Project(aBasCurv->Parabola());
1625         if(!Projector.IsDone()) return;
1626
1627         const gp_Parab2d& aP2d = Projector.Parabola();
1628
1629         anInter.Perform(aP2d, aCon);
1630       }
1631       else {
1632         Projector.Project(aBasCurv->Hyperbola());
1633         if(!Projector.IsDone()) return;
1634
1635         const gp_Hypr2d& aH2d = Projector.Hyperbola();
1636         anInter.Perform(aH2d, aCon);
1637       }
1638
1639       if(!anInter.IsDone()) return;
1640
1641       if(anInter.IsEmpty()) {
1642         NoIntersection = Standard_True;
1643         return;
1644       }
1645
1646       Standard_Integer i, nbint = anInter.NbPoints();
1647       for(i = 1; i <= nbint; i++) {
1648
1649         const IntAna2d_IntPoint& anIntPnt = anInter.Point(i);
1650
1651         umin = Min(anIntPnt.ParamOnFirst(), umin);
1652         umax = Max(anIntPnt.ParamOnFirst(), umax);
1653  
1654       }
1655
1656     }
1657     else {
1658       return;
1659     }
1660       
1661     umin = umin - Abs(umin) - 10;
1662     umax = umax + Abs(umax) + 10;
1663
1664     U1new = Max(U1new, umin);
1665     U2new = Min(U2new, umax);
1666     
1667     if(V1inf || V2inf) {
1668       EstLimForInfExtr(Line, surface, IsOffSurf, nbsu, 
1669                        Standard_False, Standard_False, V1inf, V2inf,
1670                        U1new, U2new, V1new, V2new, NoIntersection);
1671     }
1672   }
1673     
1674   return;
1675   
1676 }
1677 //=======================================================================
1678 //function : ProjectIntersectAndEstLim
1679 //purpose  : project <theLine> and it's X-axe symmetric line to <thePln> and
1680 //           intersect resulting curve with <theBasCurvProj>.
1681 //           Then estimate max and min parameters of intersection on
1682 //           <theBasCurvProj>.
1683 //           Is called from EstLimForInfRevl()
1684 //=======================================================================
1685 void ProjectIntersectAndEstLim(const gp_Lin&        theLine,
1686                                const gp_Pln&        thePln,
1687                                const ProjLib_Plane& theBasCurvProj,
1688                                Standard_Real&       theVmin,
1689                                Standard_Real&       theVmax,
1690                                Standard_Boolean&    theNoIntersection)
1691 {
1692   ProjLib_Plane aLineProj( thePln, theLine );
1693   if (!aLineProj.IsDone()) {
1694 #ifdef OCCT_DEBUG
1695   cout
1696   << "Info: IntCurveSurface_Inter::ProjectIntersectAndEstLim(), !aLineProj.IsDone()"
1697   << endl;
1698 #endif
1699     return;
1700   }
1701   gp_Lin2d aLin2d = aLineProj.Line();
1702
1703   // make a second line X-axe symmetric to the first one 
1704   gp_Pnt2d aP1 = aLin2d.Location();
1705   gp_Pnt2d aP2 (aP1.XY() + aLin2d.Direction().XY());
1706   gp_Pnt2d aP1sym (aP1.X(), -aP1.Y());
1707   gp_Pnt2d aP2sym (aP2.X(), -aP2.Y());
1708   gp_Lin2d aLin2dsym (aP1sym, gp_Vec2d(aP1sym,aP2sym));
1709   
1710   // intersect projections
1711   IntAna2d_Conic aCon    (aLin2d);
1712   IntAna2d_Conic aConSym (aLin2dsym);
1713   IntAna2d_AnaIntersection anIntersect;
1714   IntAna2d_AnaIntersection anIntersectSym;
1715   
1716   switch (theBasCurvProj.GetType()) {
1717   case GeomAbs_Line:
1718     anIntersectSym.Perform(theBasCurvProj.Line(), aConSym); 
1719     anIntersect.Perform(theBasCurvProj.Line(), aCon); break;
1720   case GeomAbs_Hyperbola:                         
1721     anIntersectSym.Perform(theBasCurvProj.Hyperbola(), aConSym); 
1722     anIntersect.Perform(theBasCurvProj.Hyperbola(), aCon); break;
1723   case GeomAbs_Parabola:                          
1724     anIntersectSym.Perform(theBasCurvProj.Parabola(), aConSym); 
1725     anIntersect.Perform(theBasCurvProj.Parabola(), aCon); break;
1726   default:
1727     return; // not infinite curve
1728   }
1729
1730   // retrieve params of intersections
1731   Standard_Integer aNbIntPnt = anIntersect.IsDone() ? anIntersect.NbPoints() : 0;
1732   Standard_Integer aNbIntPntSym = anIntersectSym.IsDone() ? anIntersectSym.NbPoints() : 0;
1733   Standard_Integer iPnt, aNbPnt = Max (aNbIntPnt, aNbIntPntSym);
1734   
1735   if (aNbPnt == 0) {
1736     theNoIntersection = Standard_True;
1737     return;
1738   }
1739   Standard_Real aParam;
1740   for (iPnt = 1; iPnt <= aNbPnt; iPnt++)
1741   {
1742     if (iPnt <= aNbIntPnt) {
1743       const IntAna2d_IntPoint& aIntPnt = anIntersect.Point(iPnt);
1744       aParam = aIntPnt.ParamOnFirst();
1745       theVmin = Min (theVmin, aParam );
1746       theVmax = Max (theVmax, aParam );
1747     }
1748     if (iPnt <= aNbIntPntSym) {
1749       const IntAna2d_IntPoint& aIntPnt = anIntersectSym.Point(iPnt);
1750       aParam = aIntPnt.ParamOnFirst();
1751       theVmin = Min (theVmin, aParam );
1752       theVmax = Max (theVmax, aParam );
1753     }
1754   }
1755
1756   return;
1757 }
1758 //=======================================================================
1759 //function : EstLimForInfRevl
1760 //purpose  : Estimate V1 and V2 to pass to InternalPerform() if they are
1761 //           infinite for Surface of Revolution
1762 //           Algo: intersect projections of Line and basis curve on the
1763 //           plane passing through revolution axe
1764 //=======================================================================
1765 void EstLimForInfRevl(const gp_Lin&   Line,
1766                       const TheSurface& surface, 
1767                       const Standard_Boolean U1inf, 
1768                       const Standard_Boolean U2inf, 
1769                       const Standard_Boolean V1inf, 
1770                       const Standard_Boolean V2inf, 
1771                       Standard_Real& U1new, 
1772                       Standard_Real& U2new, 
1773                       Standard_Real& V1new, 
1774                       Standard_Real& V2new, 
1775                       Standard_Boolean& NoIntersection)
1776 {
1777
1778   NoIntersection = Standard_False;
1779
1780   if (U1inf || U2inf) {
1781     if (U1inf)
1782       U1new = Max (0., U1new);
1783     else
1784       U2new = Min (2 * M_PI, U2new);
1785     if (! V1inf && !V2inf) return;
1786   }
1787
1788   Handle(Adaptor3d_HCurve) aBasisCurve = TheSurfaceTool::BasisCurve(surface);
1789   gp_Ax1 aRevAx = TheSurfaceTool::AxeOfRevolution(surface);
1790   gp_Vec aXVec = aRevAx.Direction();
1791   Standard_Real aTolAng = Precision::Angular();
1792
1793   // make plane to project a basis curve
1794   gp_Pnt O = aRevAx.Location();
1795   Standard_Real aU = 0.;
1796   gp_Pnt P = aBasisCurve->Value(aU);
1797   while (O.SquareDistance(P) <= Precision::PConfusion() ||
1798          aXVec.IsParallel( gp_Vec(O,P), aTolAng)) {
1799     aU += 1.;
1800     P = aBasisCurve->Value(aU);
1801     if (aU > 3)
1802       // basis curve is a line coinciding with aXVec, P is any not on aXVec
1803       P = gp_Pnt(aU, aU+1, aU+2);
1804   }
1805   gp_Vec aNVec = aXVec ^ gp_Vec(O,P);
1806   gp_Pln aPln (gp_Ax3 (O, aNVec ,aXVec));
1807
1808   // project basic curve
1809   ProjLib_Plane aBasCurvProj(aPln);
1810   switch (aBasisCurve->GetType()) {
1811   case GeomAbs_Line:
1812     aBasCurvProj.Project(aBasisCurve->Line     ()); break;
1813   case GeomAbs_Hyperbola:                       
1814     aBasCurvProj.Project(aBasisCurve->Hyperbola()); break;
1815   case GeomAbs_Parabola:                        
1816     aBasCurvProj.Project(aBasisCurve->Parabola ()); break;
1817   default:
1818     return; // not infinite curve
1819   }
1820   if (!aBasCurvProj.IsDone()) {
1821 #ifdef OCCT_DEBUG
1822     cout << "Info: IntCurveSurface_Inter::EstLimForInfRevl(), !aBasCurvProj.IsDone()" << endl;
1823 #endif
1824     return;
1825   }
1826   // make plane to project Line
1827   if (aXVec.IsParallel( Line.Direction(), aTolAng)) {
1828     P = Line.Location();
1829     while (O.SquareDistance(P) <= Precision::PConfusion()) {
1830       aU += 1.;
1831       P = gp_Pnt(aU, aU+1, aU+2); // any not on aXVec
1832     }
1833     aNVec = aXVec ^ gp_Vec(O,P);
1834   } else
1835     aNVec = aXVec.Crossed( Line.Direction() );
1836   
1837   aPln  = gp_Pln (gp_Ax3 (O, aNVec ,aXVec));
1838
1839   // make a second plane perpendicular to the first one, rotated around aXVec
1840   gp_Pln aPlnPrp = aPln.Rotated (gp_Ax1 (O,aXVec), M_PI/2.);
1841   
1842   // project Line and it's X-axe symmetric one to plane and intersect
1843   // resulting curve with projection of Basic Curev
1844   Standard_Real aVmin = RealLast(), aVmax = -aVmin;
1845   Standard_Boolean aNoInt1 = Standard_False, aNoInt2 = Standard_False;
1846   ProjectIntersectAndEstLim (Line, aPln,    aBasCurvProj, aVmin, aVmax, aNoInt1);
1847   ProjectIntersectAndEstLim (Line, aPlnPrp, aBasCurvProj, aVmin, aVmax, aNoInt2);
1848
1849   if (aNoInt1 && aNoInt2) {
1850     NoIntersection = Standard_True;
1851     return;
1852   }
1853   
1854   aVmin = aVmin - Abs(aVmin) - 10;
1855   aVmax = aVmax + Abs(aVmax) + 10;
1856
1857   if (V1inf) V1new = aVmin;
1858   if (V2inf) V2new = aVmax;
1859
1860   //cout << "EstLimForInfRevl: Vmin " << V1new << " Vmax " << V2new << endl;
1861   
1862   return;
1863 }
1864
1865 //=======================================================================
1866 //function : EstLimForInfOffs
1867 //purpose  : 
1868 //=======================================================================
1869 void EstLimForInfOffs(const gp_Lin&   Line,
1870                       const TheSurface& surface, 
1871                       const Standard_Integer nbsu, 
1872                       const Standard_Boolean U1inf, 
1873                       const Standard_Boolean U2inf, 
1874                       const Standard_Boolean V1inf, 
1875                       const Standard_Boolean V2inf, 
1876                       Standard_Real& U1new, 
1877                       Standard_Real& U2new, 
1878                       Standard_Real& V1new, 
1879                       Standard_Real& V2new, 
1880                       Standard_Boolean& NoIntersection)
1881 {
1882
1883   NoIntersection = Standard_False;
1884
1885   const Handle(Adaptor3d_HSurface)& aBasSurf = TheSurfaceTool::BasisSurface(surface);
1886   Standard_Real anOffVal = TheSurfaceTool::OffsetValue(surface);
1887
1888   GeomAbs_SurfaceType aTypeOfBasSurf = aBasSurf->GetType();
1889
1890   //  case for plane, cylinder and cone - make equivalent surface;
1891   if(aTypeOfBasSurf == GeomAbs_Plane) {
1892
1893     gp_Pln aPln = aBasSurf->Plane();
1894     gp_Vec aT = aPln.Position().XDirection()^aPln.Position().YDirection();
1895     aT *= anOffVal;
1896     aPln.Translate(aT);
1897     IntAna_IntConicQuad LinPlane(Line,aPln,TOLERANCE_ANGULAIRE);
1898
1899     if(!LinPlane.IsDone()) return;
1900
1901     if(LinPlane.IsParallel() || LinPlane.IsInQuadric()) {
1902
1903       NoIntersection = Standard_True;
1904       return;
1905
1906     }
1907     
1908     Standard_Real u, v;
1909     ElSLib::Parameters(aPln, LinPlane.Point(1), u, v);
1910     U1new = Max(U1new, u - 10.);
1911     U2new = Min(U2new, u + 10.);
1912     V1new = Max(V1new, v - 10.);
1913     V2new = Min(V2new, v + 10.);
1914
1915   }
1916   else if(aTypeOfBasSurf == GeomAbs_Cylinder) {
1917
1918     gp_Cylinder aCyl = aBasSurf->Cylinder();
1919
1920     Standard_Real aR = aCyl.Radius();
1921     gp_Ax3 anA = aCyl.Position();
1922
1923     if (anA.Direct()) 
1924       aR += anOffVal;
1925     else 
1926       aR -= anOffVal;
1927
1928     if ( aR >= TOLTANGENCY ) {
1929       aCyl.SetRadius(aR);
1930     }
1931     else if ( aR <= -TOLTANGENCY ){
1932       anA.Rotate(gp_Ax1(anA.Location(), anA.Direction()), M_PI);
1933       aCyl.SetPosition(anA);
1934 // modified by NIZHNY-MKK  Mon Oct  3 17:37:54 2005
1935 //       aCyl.SetRadius(Abs(aR));
1936       aCyl.SetRadius(-aR);
1937     }
1938     else {
1939
1940       NoIntersection = Standard_True;
1941       return;
1942       
1943     }
1944
1945     IntAna_IntConicQuad LinCylinder(Line, aCyl);
1946
1947     if(!LinCylinder.IsDone()) return;
1948
1949     if(LinCylinder.IsParallel() || LinCylinder.IsInQuadric()) {
1950
1951       NoIntersection = Standard_True;
1952       return;
1953
1954     }
1955     
1956     Standard_Integer i, nbp = LinCylinder.NbPoints();
1957     Standard_Real vmin = RealLast(), vmax = -vmin, u, v;
1958
1959     for(i = 1; i <= nbp; i++) {
1960
1961       ElSLib::Parameters(aCyl, LinCylinder.Point(i), u, v);
1962       vmin = Min(vmin, v);
1963       vmax = Max(vmax, v);
1964
1965     } 
1966
1967     V1new = Max(V1new, vmin - Abs(vmin) - 10.);
1968     V2new = Min(V2new, vmax + Abs(vmax) + 10.);
1969
1970   }
1971   else if(aTypeOfBasSurf == GeomAbs_Cone) {
1972
1973     gp_Cone aCon = aBasSurf->Cone();
1974     Standard_Real anAng = aCon.SemiAngle();
1975     Standard_Real aR = aCon.RefRadius() + anOffVal * Cos(anAng);
1976     gp_Ax3 anA = aCon.Position();
1977     if ( aR >= 0.) {
1978
1979       gp_Vec aZ( anA.Direction());
1980       aZ *= - anOffVal * Sin(anAng);
1981       anA.Translate(aZ);
1982       aCon.SetPosition(anA);
1983       aCon.SetRadius(aR);
1984       aCon.SetSemiAngle(anAng);
1985
1986     }
1987     else {
1988
1989       return;
1990
1991     }
1992
1993     IntAna_IntConicQuad LinCone(Line, aCon);
1994
1995     if(!LinCone.IsDone()) return;
1996
1997     if(LinCone.IsParallel() || LinCone.IsInQuadric()) {
1998
1999       NoIntersection = Standard_True;
2000       return;
2001
2002     }
2003     
2004     Standard_Integer i, nbp = LinCone.NbPoints();
2005     Standard_Real vmin = RealLast(), vmax = -vmin, u, v;
2006
2007     for(i = 1; i <= nbp; i++) {
2008
2009       ElSLib::Parameters(aCon, LinCone.Point(i), u, v);
2010       vmin = Min(vmin, v);
2011       vmax = Max(vmax, v);
2012
2013     } 
2014
2015     V1new = Max(V1new, vmin - Abs(vmin) - 10.);
2016     V2new = Min(V2new, vmax + Abs(vmax) + 10.);
2017
2018   }
2019   else if(aTypeOfBasSurf == GeomAbs_SurfaceOfExtrusion) {
2020
2021     Standard_Real anU1 = U1new, anU2 = U2new;
2022
2023     EstLimForInfExtr(Line, surface, Standard_True, nbsu, 
2024                      U1inf, U2inf, V1inf, V2inf, 
2025                      U1new, U2new, V1new, V2new, NoIntersection);
2026
2027     if(NoIntersection) return;
2028
2029     if(U1inf || U2inf) {
2030
2031       GeomAbs_CurveType aBasCurvType = aBasSurf->BasisCurve()->GetType();
2032       if(aBasCurvType == GeomAbs_Line) {
2033         U1new = Max(anU1, -1.e10);
2034         U2new = Min(anU2,  1.e10);
2035       } 
2036       else if(aBasCurvType == GeomAbs_Parabola) {
2037         gp_Parab aPrb  = aBasSurf->BasisCurve()->Parabola();
2038         Standard_Real aF = aPrb.Focal();
2039         Standard_Real dU = 2.e5 * Sqrt(aF);
2040         U1new = Max(anU1, -dU);
2041         U2new = Min(anU2,  dU);
2042       }
2043       else if(aBasCurvType == GeomAbs_Hyperbola) {
2044         U1new = Max(anU1, -30.);
2045         U2new = Min(anU2,  30.);
2046       }
2047       else {
2048         U1new = Max(anU1, -1.e10);
2049         U2new = Min(anU2,  1.e10);
2050       } 
2051
2052     }
2053
2054   }
2055   else if(aTypeOfBasSurf == GeomAbs_SurfaceOfRevolution) {
2056
2057     GeomAbs_CurveType aBasCurvType = aBasSurf->BasisCurve()->GetType();
2058     if(aBasCurvType == GeomAbs_Line) {
2059       V1new = Max(V1new, -1.e10);
2060       V2new = Min(V2new,  1.e10);
2061     }   
2062     else if(aBasCurvType == GeomAbs_Parabola) {
2063       gp_Parab aPrb  = aBasSurf->BasisCurve()->Parabola();
2064       Standard_Real aF = aPrb.Focal();
2065       Standard_Real dV = 2.e5 * Sqrt(aF);
2066       V1new = Max(V1new, -dV);
2067       V2new = Min(V2new,  dV);
2068     }
2069     else if(aBasCurvType == GeomAbs_Hyperbola) {
2070      V1new  = Max(V1new, -30.);
2071      V2new  = Min(V2new,  30.);
2072     }
2073     else {
2074      V1new  = Max(V1new, -1.e10);
2075      V2new  = Min(V2new,  1.e10);
2076     }   
2077
2078   }
2079   else {
2080
2081     V1new = Max(V1new, -1.e10);
2082     V2new = Min(V2new,  1.e10);
2083   }
2084 }
2085 //=======================================================================
2086 //function : EstLimForInfSurf
2087 //purpose  : 
2088 //=======================================================================
2089 void EstLimForInfSurf(Standard_Real& U1new, 
2090                       Standard_Real& U2new, 
2091                       Standard_Real& V1new, 
2092                       Standard_Real& V2new)
2093 {
2094     U1new = Max(U1new, -1.e10);
2095     U2new = Min(U2new,  1.e10);
2096     V1new = Max(V1new, -1.e10);
2097     V2new = Min(V2new,  1.e10);
2098 }