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