Warning, /include/opencascade/IntCurveSurface_Inter.gxx is written in an unsupported language. File is not indexed.
0001 // Created on: 1993-04-09
0002 // Created by: Laurent BUCHARD
0003 // Copyright (c) 1993-1999 Matra Datavision
0004 // Copyright (c) 1999-2014 OPEN CASCADE SAS
0005 //
0006 // This file is part of Open CASCADE Technology software library.
0007 //
0008 // This library is free software; you can redistribute it and/or modify it under
0009 // the terms of the GNU Lesser General Public License version 2.1 as published
0010 // by the Free Software Foundation, with special exception defined in the file
0011 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0012 // distribution for complete text of the license and disclaimer of any warranty.
0013 //
0014 // Alternatively, this file may be used under the terms of Open CASCADE
0015 // commercial license or contractual agreement.
0016
0017 //#ifndef OCCT_DEBUG
0018 //#define No_Standard_RangeError
0019 //#define No_Standard_OutOfRange
0020 //#endif
0021
0022
0023 #define TOLTANGENCY 0.00000001
0024 #define TOLERANCE_ANGULAIRE 1.e-12//0.00000001
0025 #define TOLERANCE 0.00000001
0026
0027 #define NBSAMPLESONCIRCLE 32
0028 #define NBSAMPLESONELLIPSE 32
0029 #define NBSAMPLESONPARAB 16
0030 #define NBSAMPLESONHYPR 32
0031
0032 #include <ElSLib.hxx>
0033 #include <Intf_SectionPoint.hxx>
0034 #include <Intf_TangentZone.hxx>
0035 #include <Intf_Tool.hxx>
0036 #include <math_FunctionSetRoot.hxx>
0037 #include <IntCurveSurface_IntersectionPoint.hxx>
0038 #include <IntCurveSurface_IntersectionSegment.hxx>
0039 #include <IntAna_IntConicQuad.hxx>
0040 #include <IntAna_Quadric.hxx>
0041 #include <gp_Vec.hxx>
0042 #include <gp_Dir.hxx>
0043 #include <Precision.hxx>
0044 #include <IntAna_IntLinTorus.hxx>
0045
0046 #include <GeomAbs_Shape.hxx>
0047 #include <GeomAbs_CurveType.hxx>
0048 #include <TColStd_Array1OfReal.hxx>
0049 #include <TColgp_Array1OfPnt.hxx>
0050 #include <gp_Pnt.hxx>
0051 #include <gp_Pln.hxx>
0052 #include <gp_Lin.hxx>
0053 #include <GProp_PEquation.hxx>
0054 #include <GProp_PGProps.hxx>
0055 #include <GProp_PrincipalProps.hxx>
0056 #include <Extrema_ExtElC.hxx>
0057 #include <Extrema_POnCurv.hxx>
0058
0059
0060
0061 #include <ProjLib_Plane.hxx>
0062 #include <IntAna2d_AnaIntersection.hxx>
0063 #include <gp_Lin2d.hxx>
0064 #include <IntAna2d_IntPoint.hxx>
0065 #include <gp_Parab2d.hxx>
0066 #include <IntAna2d_Conic.hxx>
0067 #include <IntAna2d_AnaIntersection.hxx>
0068 #include <gp_Hypr2d.hxx>
0069 #include <Adaptor3d_Curve.hxx>
0070 #include <Adaptor3d_Surface.hxx>
0071
0072
0073 #include <TColgp_Array2OfPnt.hxx>
0074 #include <TColStd_HArray1OfReal.hxx>
0075 #include <Adaptor3d_TopolTool.hxx>
0076 #include <ElCLib.hxx>
0077 //#define ICS_DEB
0078 static
0079 void EstLimForInfExtr(const gp_Lin& Line,
0080 const TheSurface& surface,
0081 const Standard_Boolean IsOffSurf,
0082 const Standard_Integer nbsu,
0083 const Standard_Boolean U1inf,
0084 const Standard_Boolean U2inf,
0085 const Standard_Boolean V1inf,
0086 const Standard_Boolean V2inf,
0087 Standard_Real& U1new,
0088 Standard_Real& U2new,
0089 Standard_Real& V1new,
0090 Standard_Real& V2new,
0091 Standard_Boolean& NoIntersection);
0092
0093 static
0094 void EstLimForInfRevl(const gp_Lin& Line,
0095 const TheSurface& surface,
0096 const Standard_Boolean U1inf,
0097 const Standard_Boolean U2inf,
0098 const Standard_Boolean V1inf,
0099 const Standard_Boolean V2inf,
0100 Standard_Real& U1new,
0101 Standard_Real& U2new,
0102 Standard_Real& V1new,
0103 Standard_Real& V2new,
0104 Standard_Boolean& NoIntersection);
0105
0106 static
0107 void EstLimForInfOffs(const gp_Lin& Line,
0108 const TheSurface& surface,
0109 const Standard_Integer nbsu,
0110 const Standard_Boolean U1inf,
0111 const Standard_Boolean U2inf,
0112 const Standard_Boolean V1inf,
0113 const Standard_Boolean V2inf,
0114 Standard_Real& U1new,
0115 Standard_Real& U2new,
0116 Standard_Real& V1new,
0117 Standard_Real& V2new,
0118 Standard_Boolean& NoIntersection);
0119
0120 static
0121 void EstLimForInfSurf(Standard_Real& U1new,
0122 Standard_Real& U2new,
0123 Standard_Real& V1new,
0124 Standard_Real& V2new);
0125
0126 static
0127 void SectionPointToParameters(const Intf_SectionPoint& Sp,
0128 const IntCurveSurface_ThePolyhedron& Surf,
0129 const IntCurveSurface_ThePolygon& Curv,
0130 Standard_Real& u,
0131 Standard_Real& v,
0132 Standard_Real& w);
0133
0134 static
0135 void IntCurveSurface_ComputeTransitions(const TheCurve& curve,
0136 const Standard_Real w,
0137 IntCurveSurface_TransitionOnCurve& TransOnCurve,
0138 const TheSurface& surface,
0139 const Standard_Real u,
0140 const Standard_Real v);
0141
0142 static
0143 void IntCurveSurface_ComputeParamsOnQuadric(const TheSurface& surface,
0144 const gp_Pnt& P,
0145 Standard_Real& u,
0146 Standard_Real& v);
0147
0148 static
0149 void ProjectIntersectAndEstLim(const gp_Lin& theLine,
0150 const gp_Pln& thePln,
0151 const ProjLib_Plane& theBasCurvProj,
0152 Standard_Real& theVmin,
0153 Standard_Real& theVmax,
0154 Standard_Boolean& theNoIntersection);
0155
0156 //=======================================================================
0157 //function : IntCurveSurface_Inter
0158 //purpose :
0159 //=======================================================================
0160 IntCurveSurface_Inter::IntCurveSurface_Inter()
0161 {
0162 }
0163 //=======================================================================
0164 //function : DoSurface
0165 //purpose :
0166 //=======================================================================
0167 void IntCurveSurface_Inter::DoSurface(const TheSurface& surface,
0168 const Standard_Real u0,
0169 const Standard_Real u1,
0170 const Standard_Real v0,
0171 const Standard_Real v1,
0172 TColgp_Array2OfPnt& pntsOnSurface,
0173 Bnd_Box& boxSurface,
0174 Standard_Real& gap)
0175 {
0176 Standard_Integer iU = 0, iV = 0;
0177 Standard_Real U = 0., V = 0;
0178 // modified by NIZHNY-MKK Mon Oct 3 17:38:45 2005
0179 // Standard_Real dU = fabs(u1-u0)/50., dV = fabs(v1-v0)/50.;
0180 Standard_Real dU = (u1-u0)/50., dV = (v1-v0)/50.;
0181 gp_Pnt aPnt;
0182
0183 for(iU = 0; iU < 50; iU++) {
0184
0185 if(iU == 0)
0186 U = u0;
0187 else if(iU == 49)
0188 U = u1;
0189 else
0190 U = u0 + dU * ((Standard_Real)iU);
0191
0192 for(iV = 0; iV < 50; iV++) {
0193
0194 if(iV == 0)
0195 V = v0;
0196 else if(iV == 49)
0197 V = v1;
0198 else
0199 V = v0 + dV * ((Standard_Real)iV);
0200
0201 TheSurfaceTool::D0(surface,U,V,aPnt);
0202 boxSurface.Add(aPnt);
0203 pntsOnSurface.SetValue(iU+1,iV+1,aPnt);
0204 }
0205 }
0206 Standard_Real Ures = TheSurfaceTool::UResolution(surface,dU);
0207 Standard_Real Vres = TheSurfaceTool::VResolution(surface,dV);
0208 gap = Max(Ures,Vres);
0209 }
0210
0211 //=======================================================================
0212 //function : DoNewBounds
0213 //purpose :
0214 //=======================================================================
0215 void IntCurveSurface_Inter::DoNewBounds(
0216 const TheSurface& surface,
0217 const Standard_Real u0,
0218 const Standard_Real u1,
0219 const Standard_Real v0,
0220 const Standard_Real v1,
0221 const TColgp_Array2OfPnt& pntsOnSurface,
0222 const TColStd_Array1OfReal& X,
0223 const TColStd_Array1OfReal& Y,
0224 const TColStd_Array1OfReal& Z,
0225 TColStd_Array1OfReal& Bounds)
0226 {
0227 Bounds.SetValue(1,u0);
0228 Bounds.SetValue(2,u1);
0229 Bounds.SetValue(3,v0);
0230 Bounds.SetValue(4,v1);
0231
0232 Standard_Boolean isUClosed = (TheSurfaceTool::IsUClosed(surface) || TheSurfaceTool::IsUPeriodic(surface));
0233 Standard_Boolean isVClosed = (TheSurfaceTool::IsVClosed(surface) || TheSurfaceTool::IsVPeriodic(surface));
0234 Standard_Boolean checkU = (isUClosed) ? Standard_False : Standard_True;
0235 Standard_Boolean checkV = (isVClosed) ? Standard_False : Standard_True;
0236
0237 Standard_Integer i = 0, j = 0, k = 0, iU = 0, iV = 0;
0238 Standard_Integer iUmin = 50, iVmin = 50, iUmax = 1, iVmax = 1;
0239
0240 for(i = 1; i <= 2; i++) {
0241 for(j = 1; j <= 2; j++) {
0242 for(k = 1; k <= 2; k++) {
0243 gp_Pnt aPoint(X(i),Y(j),Z(k));
0244 Standard_Real DistMin = 1.e+100;
0245 Standard_Integer diU = 0, diV = 0;
0246 for(iU = 1; iU <= 50; iU++) {
0247 for(iV = 1; iV <= 50; iV++) {
0248 const gp_Pnt aP = pntsOnSurface.Value(iU,iV);
0249 Standard_Real dist = aP.SquareDistance(aPoint);
0250 if(dist < DistMin) {
0251 DistMin = dist;
0252 diU = iU;
0253 diV = iV;
0254 }
0255 }
0256 }
0257 if(diU > 0 && diU < iUmin) iUmin = diU;
0258 if(diU > 0 && diU > iUmax) iUmax = diU;
0259 if(diV > 0 && diV < iVmin) iVmin = diV;
0260 if(diV > 0 && diV > iVmax) iVmax = diV;
0261 }
0262 }
0263 }
0264
0265 // modified by NIZHNY-MKK Mon Oct 3 17:38:43 2005
0266 // Standard_Real dU = fabs(u1-u0)/50., dV = fabs(v1-v0)/50.;
0267 Standard_Real dU = (u1-u0)/50., dV = (v1-v0)/50.;
0268
0269 Standard_Real USmin = u0 + dU * ((Standard_Real)(iUmin - 1));
0270 Standard_Real USmax = u0 + dU * ((Standard_Real)(iUmax - 1));
0271 Standard_Real VSmin = v0 + dV * ((Standard_Real)(iVmin - 1));
0272 Standard_Real VSmax = v0 + dV * ((Standard_Real)(iVmax - 1));
0273
0274 if(USmin > USmax) {
0275 Standard_Real tmp = USmax;
0276 USmax = USmin;
0277 USmin = tmp;
0278 }
0279 if(VSmin > VSmax) {
0280 Standard_Real tmp = VSmax;
0281 VSmax = VSmin;
0282 VSmin = tmp;
0283 }
0284
0285 USmin -= 1.5*dU;
0286 if(USmin < u0) USmin = u0;
0287 USmax += 1.5*dU;
0288 if(USmax > u1) USmax = u1;
0289 VSmin -= 1.5*dV;
0290 if(VSmin < v0) VSmin = v0;
0291 VSmax += 1.5*dV;
0292 if(VSmax > v1) VSmax = v1;
0293
0294 if(checkU) {
0295 Bounds.SetValue(1,USmin);
0296 Bounds.SetValue(2,USmax);
0297 }
0298 if(checkV) {
0299 Bounds.SetValue(3,VSmin);
0300 Bounds.SetValue(4,VSmax);
0301 }
0302 }
0303
0304 //=======================================================================
0305 //function : Perform
0306 //purpose : Decompose la surface si besoin est
0307 //=======================================================================
0308 void IntCurveSurface_Inter::Perform(const TheCurve& curve,
0309 const TheSurface& surface)
0310 {
0311 ResetFields();
0312 done = Standard_True;
0313 Standard_Integer NbUOnS = TheSurfaceTool::NbUIntervals(surface,GeomAbs_C2);
0314 Standard_Integer NbVOnS = TheSurfaceTool::NbVIntervals(surface,GeomAbs_C2);
0315 Standard_Real U0,U1,V0,V1;
0316
0317 if(NbUOnS > 1)
0318 {
0319 TColStd_Array1OfReal TabU(1,NbUOnS+1);
0320 TheSurfaceTool::UIntervals(surface,TabU,GeomAbs_C2);
0321 for(Standard_Integer iu = 1;iu <= NbUOnS; iu++)
0322 {
0323 U0 = TabU.Value(iu);
0324 U1 = TabU.Value(iu+1);
0325 if(NbVOnS > 1)
0326 {
0327 TColStd_Array1OfReal TabV(1,NbVOnS+1);
0328 TheSurfaceTool::VIntervals(surface,TabV,GeomAbs_C2);
0329 for(Standard_Integer iv = 1;iv <= NbVOnS; iv++)
0330 {
0331 // More than one interval on U and V param space.
0332 V0 = TabV.Value(iv);
0333 V1 = TabV.Value(iv+1);
0334 Perform(curve,surface,U0,V0,U1,V1);
0335 }
0336 }
0337 else
0338 {
0339 // More than one interval only on U param space.
0340 V0 = TheSurfaceTool::FirstVParameter(surface);
0341 V1 = TheSurfaceTool::LastVParameter(surface);
0342 Perform(curve,surface,U0,V0,U1,V1);
0343 }
0344 }
0345 }
0346 else if(NbVOnS > 1)
0347 {
0348 // More than one interval only on V param space.
0349 U0 = TheSurfaceTool::FirstUParameter(surface);
0350 U1 = TheSurfaceTool::LastUParameter(surface);
0351 TColStd_Array1OfReal TabV(1,NbVOnS+1);
0352 TheSurfaceTool::VIntervals(surface,TabV,GeomAbs_C2);
0353 for(Standard_Integer iv = 1;iv <= NbVOnS; iv++)
0354 {
0355 V0 = TabV.Value(iv);
0356 V1 = TabV.Value(iv+1);
0357 Perform(curve,surface,U0,V0,U1,V1);
0358 }
0359 }
0360 else
0361 {
0362 // One interval on U and V param space.
0363 V0 = TheSurfaceTool::FirstVParameter(surface);
0364 V1 = TheSurfaceTool::LastVParameter(surface);
0365 U0 = TheSurfaceTool::FirstUParameter(surface);
0366 U1 = TheSurfaceTool::LastUParameter(surface);
0367
0368 Perform(curve,surface,U0,V0,U1,V1);
0369 }
0370 }
0371 //=======================================================================
0372 //function : Perform
0373 //purpose :
0374 //=======================================================================
0375 void IntCurveSurface_Inter::Perform(const TheCurve& curve,
0376 const TheSurface& surface,
0377 const Standard_Real U1,const Standard_Real V1,
0378 const Standard_Real U2,const Standard_Real V2)
0379 {
0380 // Protection from double type overflow.
0381 // This may happen inside square magnitude computation based on normal,
0382 // which was computed on bound parameteres (bug26525).
0383 Standard_Real UU1 = U1, UU2 = U2, VV1 = V1, VV2 = V2;
0384 if (U1 < -1.0e50)
0385 UU1 = -1.0e50;
0386 if (U2 > 1.0e50)
0387 UU2 = 1.0e50;
0388 if (V1 < -1.0e50)
0389 VV1 = -1.0e50;
0390 if (V2 > 1.0e50)
0391 VV2 = 1.0e50;
0392
0393 GeomAbs_CurveType CurveType = TheCurveTool::GetType(curve);
0394
0395 switch(CurveType) {
0396 case GeomAbs_Line:
0397 {
0398 PerformConicSurf(TheCurveTool::Line(curve),curve,surface,UU1,VV1,UU2,VV2);
0399 break;
0400 }
0401 case GeomAbs_Circle:
0402 {
0403 PerformConicSurf(TheCurveTool::Circle(curve),curve,surface,UU1,VV1,UU2,VV2);
0404 break;
0405 }
0406 case GeomAbs_Ellipse:
0407 {
0408 PerformConicSurf(TheCurveTool::Ellipse(curve),curve,surface,UU1,VV1,UU2,VV2);
0409 break;
0410 }
0411 case GeomAbs_Parabola:
0412 {
0413 PerformConicSurf(TheCurveTool::Parabola(curve),curve,surface,UU1,VV1,UU2,VV2);
0414 break;
0415 }
0416 case GeomAbs_Hyperbola:
0417 {
0418 PerformConicSurf(TheCurveTool::Hyperbola(curve),curve,surface,UU1,VV1,UU2,VV2);
0419 break;
0420 }
0421 default:
0422 {
0423 Standard_Integer nbIntervalsOnCurve = TheCurveTool::NbIntervals(curve,GeomAbs_C2);
0424 GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
0425 if( (SurfaceType != GeomAbs_Plane)
0426 && (SurfaceType != GeomAbs_Cylinder)
0427 && (SurfaceType != GeomAbs_Cone)
0428 && (SurfaceType != GeomAbs_Sphere)) {
0429
0430 if(nbIntervalsOnCurve > 1) {
0431 TColStd_Array1OfReal TabW(1,nbIntervalsOnCurve+1);
0432 TheCurveTool::Intervals(curve,TabW,GeomAbs_C2);
0433 for(Standard_Integer i = 1; i<=nbIntervalsOnCurve; i++) {
0434 Standard_Real u1,u2;
0435 u1 = TabW.Value(i);
0436 u2 = TabW.Value(i+1);
0437
0438 Handle(TColStd_HArray1OfReal) aPars;
0439 Standard_Real defl = 0.1;
0440 Standard_Integer NbMin = 10;
0441 TheCurveTool::SamplePars(curve, u1, u2, defl, NbMin, aPars);
0442
0443 // IntCurveSurface_ThePolygon polygon(curve,u1,u2,TheCurveTool::NbSamples(curve,u1,u2));
0444 IntCurveSurface_ThePolygon polygon(curve, aPars->Array1());
0445 InternalPerform(curve,polygon,surface,UU1,VV1,UU2,VV2);
0446 }
0447 }
0448 else {
0449 Standard_Real u1,u2;
0450 u1 = TheCurveTool::FirstParameter(curve);
0451 u2 = TheCurveTool::LastParameter(curve);
0452
0453 Handle(TColStd_HArray1OfReal) aPars;
0454 Standard_Real defl = 0.1;
0455 Standard_Integer NbMin = 10;
0456 TheCurveTool::SamplePars(curve, u1, u2, defl, NbMin, aPars);
0457
0458 // IntCurveSurface_ThePolygon polygon(curve,TheCurveTool::NbSamples(curve,u1,u2));
0459 IntCurveSurface_ThePolygon polygon(curve, aPars->Array1());
0460 InternalPerform(curve,polygon,surface,UU1,VV1,UU2,VV2);
0461 }
0462 }
0463 else { //-- la surface est une quadrique
0464 InternalPerformCurveQuadric(curve,surface);
0465 }
0466 }
0467 }
0468 }
0469 //=======================================================================
0470 //function : Perform
0471 //purpose :
0472 //=======================================================================
0473 void IntCurveSurface_Inter::Perform(const TheCurve& curve,
0474 const IntCurveSurface_ThePolygon& polygon,
0475 const TheSurface& surface) {
0476 ResetFields();
0477 done = Standard_True;
0478 Standard_Real u1,v1,u2,v2;
0479 u1 = TheSurfaceTool::FirstUParameter(surface);
0480 v1 = TheSurfaceTool::FirstVParameter(surface);
0481 u2 = TheSurfaceTool::LastUParameter(surface);
0482 v2 = TheSurfaceTool::LastVParameter(surface);
0483 Standard_Integer nbsu,nbsv;
0484 nbsu = TheSurfaceTool::NbSamplesU(surface,u1,u2);
0485 nbsv = TheSurfaceTool::NbSamplesV(surface,v1,v2);
0486 if(nbsu>40) nbsu=40;
0487 if(nbsv>40) nbsv=40;
0488 IntCurveSurface_ThePolyhedron polyhedron(surface,nbsu,nbsv,u1,v1,u2,v2);
0489 Perform(curve,polygon,surface,polyhedron);
0490 }
0491 //=======================================================================
0492 //function : Perform
0493 //purpose :
0494 //=======================================================================
0495 void IntCurveSurface_Inter::Perform(const TheCurve& curve,
0496 const TheSurface& surface,
0497 const IntCurveSurface_ThePolyhedron& polyhedron) {
0498 ResetFields();
0499 done = Standard_True;
0500 Standard_Real u1 = TheCurveTool::FirstParameter(curve);
0501 Standard_Real u2 = TheCurveTool::LastParameter(curve);
0502 IntCurveSurface_ThePolygon polygon(curve,TheCurveTool::NbSamples(curve,u1,u2));
0503 Perform(curve,polygon,surface,polyhedron);
0504 }
0505 //=======================================================================
0506 //function : Perform
0507 //purpose :
0508 //=======================================================================
0509 void IntCurveSurface_Inter::Perform(const TheCurve& curve,
0510 const IntCurveSurface_ThePolygon& polygon,
0511 const TheSurface& surface,
0512 const IntCurveSurface_ThePolyhedron& polyhedron) {
0513 ResetFields();
0514 done = Standard_True;
0515 Standard_Real u1,v1,u2,v2;
0516 u1 = TheSurfaceTool::FirstUParameter(surface);
0517 v1 = TheSurfaceTool::FirstVParameter(surface);
0518 u2 = TheSurfaceTool::LastUParameter(surface);
0519 v2 = TheSurfaceTool::LastVParameter(surface);
0520 InternalPerform(curve,polygon,surface,polyhedron,u1,v1,u2,v2);
0521 }
0522
0523 //=======================================================================
0524 //function : Perform
0525 //purpose :
0526 //=======================================================================
0527 void IntCurveSurface_Inter::Perform(const TheCurve& curve,
0528 const IntCurveSurface_ThePolygon& polygon,
0529 const TheSurface& surface,
0530 const IntCurveSurface_ThePolyhedron& polyhedron,
0531 Bnd_BoundSortBox& BndBSB) {
0532 ResetFields();
0533 done = Standard_True;
0534 Standard_Real u1,v1,u2,v2;
0535 u1 = TheSurfaceTool::FirstUParameter(surface);
0536 v1 = TheSurfaceTool::FirstVParameter(surface);
0537 u2 = TheSurfaceTool::LastUParameter(surface);
0538 v2 = TheSurfaceTool::LastVParameter(surface);
0539 InternalPerform(curve,polygon,surface,polyhedron,u1,v1,u2,v2,BndBSB);
0540 }
0541 //=======================================================================
0542 //function : InternalPerform
0543 //purpose : C a l c u l d u p o i n t a p p r o c h e
0544 //== p u i s d u p o i n t E x a c t
0545 //=======================================================================
0546 void IntCurveSurface_Inter::InternalPerform(const TheCurve& curve,
0547 const IntCurveSurface_ThePolygon& polygon,
0548 const TheSurface& surface,
0549 const IntCurveSurface_ThePolyhedron& polyhedron,
0550 const Standard_Real u0,
0551 const Standard_Real v0,
0552 const Standard_Real u1,
0553 const Standard_Real v1,
0554 Bnd_BoundSortBox& BSB) {
0555 IntCurveSurface_TheInterference interference(polygon,polyhedron,BSB);
0556 IntCurveSurface_TheCSFunction theicsfunction(surface,curve);
0557 IntCurveSurface_TheExactInter intersectionExacte(theicsfunction,TOLTANGENCY);
0558 math_FunctionSetRoot rsnld(intersectionExacte.Function());
0559
0560 // Standard_Real u,v,w,winit;
0561 Standard_Real u,v,w;
0562 gp_Pnt P;
0563 Standard_Real winf = polygon.InfParameter();
0564 Standard_Real wsup = polygon.SupParameter();
0565 Standard_Integer NbSectionPoints = interference.NbSectionPoints();
0566 Standard_Integer NbTangentZones = interference.NbTangentZones();
0567
0568
0569 //-- Les interferences renvoient parfois de nombreuses fois (>20) les memes points
0570
0571 Standard_Integer i,NbStartPoints=NbSectionPoints;
0572 for(i=1; i<= NbTangentZones; i++) {
0573 const Intf_TangentZone& TZ = interference.ZoneValue(i);
0574 Standard_Integer nbpnts = TZ.NumberOfPoints();
0575 NbStartPoints+=nbpnts;
0576 }
0577
0578 if(NbStartPoints) {
0579 Standard_Real *TabU = new Standard_Real [NbStartPoints+1];
0580 Standard_Real *TabV = new Standard_Real [NbStartPoints+1];
0581 Standard_Real *TabW = new Standard_Real [NbStartPoints+1];
0582 Standard_Integer IndexPoint=0;
0583
0584 for(i=1; i<= NbSectionPoints; i++) {
0585 const Intf_SectionPoint& SP = interference.PntValue(i);
0586 SectionPointToParameters(SP,polyhedron,polygon,u,v,w);
0587 TabU[IndexPoint]=u;
0588 TabV[IndexPoint]=v;
0589 TabW[IndexPoint]=w;
0590 IndexPoint++;
0591 }
0592 for(i=1; i<= NbTangentZones; i++) {
0593 const Intf_TangentZone& TZ = interference.ZoneValue(i);
0594 Standard_Integer nbpnts = TZ.NumberOfPoints();
0595 for(Standard_Integer j=1; j<=nbpnts; j++) {
0596 const Intf_SectionPoint& SP = TZ.GetPoint(j);
0597 SectionPointToParameters(SP,polyhedron,polygon,u,v,w);
0598 TabU[IndexPoint]=u;
0599 TabV[IndexPoint]=v;
0600 TabW[IndexPoint]=w;
0601 IndexPoint++;
0602 }
0603 }
0604
0605 //-- Tri
0606 Standard_Real su=0,sv=0,sw=0,ptol;
0607 ptol = 10*Precision::PConfusion();
0608
0609 //-- Tri suivant la variable W
0610 Standard_Boolean Triok;
0611 do {
0612 Triok=Standard_True;
0613 Standard_Integer im1;
0614 for(i=1,im1=0;i<NbStartPoints;im1++,i++) {
0615 if(TabW[i] < TabW[im1]) {
0616 Standard_Real t=TabW[i]; TabW[i]=TabW[im1]; TabW[im1]=t;
0617 t=TabU[i]; TabU[i]=TabU[im1]; TabU[im1]=t;
0618 t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
0619 Triok=Standard_False;
0620 }
0621 }
0622 }
0623 while(Triok==Standard_False);
0624
0625 //-- On trie pour des meme W suivant U
0626 do {
0627 Triok=Standard_True;
0628 Standard_Integer im1;
0629 for(i=1,im1=0;i<NbStartPoints;im1++,i++) {
0630 // modified by NIZHNY-MKK Mon Oct 3 17:38:49 2005
0631 // if(Abs(TabW[i]-TabW[im1])<ptol) {
0632 if((TabW[i]-TabW[im1])<ptol) {
0633 TabW[i]=TabW[im1];
0634 if(TabU[i] < TabU[im1]) {
0635 Standard_Real t=TabU[i]; TabU[i]=TabU[im1]; TabU[im1]=t;
0636 t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
0637 Triok=Standard_False;
0638 }
0639 }
0640 }
0641 }
0642 while(Triok==Standard_False);
0643
0644 //-- On trie pour des meme U et W suivant V
0645 do {
0646 Triok=Standard_True;
0647 Standard_Integer im1;
0648 for(i=1,im1=0;i<NbStartPoints;im1++,i++) {
0649 // modified by NIZHNY-MKK Mon Oct 3 17:38:52 2005
0650 // if((Abs(TabW[i]-TabW[im1])<ptol) && (Abs(TabU[i]-TabU[im1])<ptol)) {
0651 if(((TabW[i]-TabW[im1])<ptol) && ((TabU[i]-TabU[im1])<ptol)) {
0652 TabU[i]=TabU[im1];
0653 if(TabV[i] < TabV[im1]) {
0654 Standard_Real t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
0655 Triok=Standard_False;
0656 }
0657 }
0658 }
0659 }
0660 while(Triok==Standard_False);
0661
0662
0663 for(i=0;i<NbStartPoints; i++) {
0664 u=TabU[i]; v=TabV[i]; w=TabW[i];
0665 if(i==0) {
0666 su=u-1;
0667 }
0668 if(Abs(u-su)>ptol || Abs(v-sv)>ptol || Abs(w-sw)>ptol) {
0669 intersectionExacte.Perform(u,v,w,rsnld,u0,u1,v0,v1,winf,wsup);
0670 if(intersectionExacte.IsDone()) {
0671 if(!intersectionExacte.IsEmpty()) {
0672 P=intersectionExacte.Point();
0673 w=intersectionExacte.ParameterOnCurve();
0674 intersectionExacte.ParameterOnSurface(u,v);
0675 AppendPoint(curve,w,surface,u,v);
0676 }
0677 }
0678 }
0679 su=TabU[i]; sv=TabV[i]; sw=TabW[i];
0680 }
0681 delete [] TabW;
0682 delete [] TabV;
0683 delete [] TabU;
0684 }
0685 }
0686
0687 //=======================================================================
0688 //function : InternalPerform
0689 //purpose :
0690 //=======================================================================
0691 void IntCurveSurface_Inter::InternalPerform(const TheCurve& curve,
0692 const IntCurveSurface_ThePolygon& polygon,
0693 const TheSurface& surface,
0694 const IntCurveSurface_ThePolyhedron& polyhedron,
0695 const Standard_Real u0,
0696 const Standard_Real v0,
0697 const Standard_Real u1,
0698 const Standard_Real v1) {
0699 IntCurveSurface_TheInterference interference(polygon,polyhedron);
0700 IntCurveSurface_TheCSFunction theicsfunction(surface,curve);
0701 IntCurveSurface_TheExactInter intersectionExacte(theicsfunction,TOLTANGENCY);
0702 math_FunctionSetRoot rsnld(intersectionExacte.Function());
0703
0704 // Standard_Real u,v,w,winit;
0705 Standard_Real u,v,w;
0706 gp_Pnt P;
0707 Standard_Real winf = polygon.InfParameter();
0708 Standard_Real wsup = polygon.SupParameter();
0709 Standard_Integer NbSectionPoints = interference.NbSectionPoints();
0710 Standard_Integer NbTangentZones = interference.NbTangentZones();
0711
0712
0713 //-- Les interferences renvoient parfois de nombreuses fois (>20) les memes points
0714
0715 Standard_Integer i,NbStartPoints=NbSectionPoints;
0716 for(i=1; i<= NbTangentZones; i++) {
0717 const Intf_TangentZone& TZ = interference.ZoneValue(i);
0718 Standard_Integer nbpnts = TZ.NumberOfPoints();
0719 NbStartPoints+=nbpnts;
0720 }
0721
0722 if(NbStartPoints) {
0723 Standard_Real *TabU = new Standard_Real [NbStartPoints+1];
0724 Standard_Real *TabV = new Standard_Real [NbStartPoints+1];
0725 Standard_Real *TabW = new Standard_Real [NbStartPoints+1];
0726 Standard_Integer IndexPoint=0;
0727
0728 for(i=1; i<= NbSectionPoints; i++) {
0729 const Intf_SectionPoint& SP = interference.PntValue(i);
0730 SectionPointToParameters(SP,polyhedron,polygon,u,v,w);
0731 TabU[IndexPoint]=u;
0732 TabV[IndexPoint]=v;
0733 TabW[IndexPoint]=w;
0734 IndexPoint++;
0735 }
0736 for(i=1; i<= NbTangentZones; i++) {
0737 const Intf_TangentZone& TZ = interference.ZoneValue(i);
0738 Standard_Integer nbpnts = TZ.NumberOfPoints();
0739 for(Standard_Integer j=1; j<=nbpnts; j++) {
0740 const Intf_SectionPoint& SP = TZ.GetPoint(j);
0741 SectionPointToParameters(SP,polyhedron,polygon,u,v,w);
0742 TabU[IndexPoint]=u;
0743 TabV[IndexPoint]=v;
0744 TabW[IndexPoint]=w;
0745 IndexPoint++;
0746 }
0747 }
0748
0749 //-- Tri
0750 Standard_Real su=0,sv=0,sw=0,ptol;
0751 ptol = 10*Precision::PConfusion();
0752
0753 //-- Tri suivant la variable W
0754 Standard_Boolean Triok;
0755 do {
0756 Triok=Standard_True;
0757 Standard_Integer im1;
0758 for(i=1,im1=0;i<NbStartPoints;im1++,i++) {
0759 if(TabW[i] < TabW[im1]) {
0760 Standard_Real t=TabW[i]; TabW[i]=TabW[im1]; TabW[im1]=t;
0761 t=TabU[i]; TabU[i]=TabU[im1]; TabU[im1]=t;
0762 t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
0763 Triok=Standard_False;
0764 }
0765 }
0766 }
0767 while(Triok==Standard_False);
0768
0769 //-- On trie pour des meme W suivant U
0770 do {
0771 Triok=Standard_True;
0772 Standard_Integer im1;
0773 for(i=1,im1=0;i<NbStartPoints;im1++,i++) {
0774 // modified by NIZHNY-MKK Mon Oct 3 17:38:56 2005
0775 // if(Abs(TabW[i]-TabW[im1])<ptol) {
0776 if((TabW[i]-TabW[im1])<ptol) {
0777 TabW[i]=TabW[im1];
0778 if(TabU[i] < TabU[im1]) {
0779 Standard_Real t=TabU[i]; TabU[i]=TabU[im1]; TabU[im1]=t;
0780 t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
0781 Triok=Standard_False;
0782 }
0783 }
0784 }
0785 }
0786 while(Triok==Standard_False);
0787
0788 //-- On trie pour des meme U et W suivant V
0789 do {
0790 Triok=Standard_True;
0791 Standard_Integer im1;
0792 for(i=1,im1=0;i<NbStartPoints;im1++,i++) {
0793 // modified by NIZHNY-MKK Mon Oct 3 17:38:58 2005
0794 // if((Abs(TabW[i]-TabW[im1])<ptol) && (Abs(TabU[i]-TabU[im1])<ptol)) {
0795 if(((TabW[i]-TabW[im1])<ptol) && ((TabU[i]-TabU[im1])<ptol)) {
0796 TabU[i]=TabU[im1];
0797 if(TabV[i] < TabV[im1]) {
0798 Standard_Real t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
0799 Triok=Standard_False;
0800 }
0801 }
0802 }
0803 }
0804 while(Triok==Standard_False);
0805
0806
0807 for(i=0;i<NbStartPoints; i++) {
0808 u=TabU[i]; v=TabV[i]; w=TabW[i];
0809 if(i==0) {
0810 su=u-1;
0811 }
0812 if(Abs(u-su)>ptol || Abs(v-sv)>ptol || Abs(w-sw)>ptol) {
0813 intersectionExacte.Perform(u,v,w,rsnld,u0,u1,v0,v1,winf,wsup);
0814 if(intersectionExacte.IsDone()) {
0815 if(!intersectionExacte.IsEmpty()) {
0816 P=intersectionExacte.Point();
0817 w=intersectionExacte.ParameterOnCurve();
0818 intersectionExacte.ParameterOnSurface(u,v);
0819 AppendPoint(curve,w,surface,u,v);
0820 }
0821 }
0822 }
0823 su=TabU[i]; sv=TabV[i]; sw=TabW[i];
0824 }
0825 delete [] TabW;
0826 delete [] TabV;
0827 delete [] TabU;
0828 }
0829 }
0830 //=======================================================================
0831 //function : InternalPerformCurveQuadric
0832 //purpose :
0833 //=======================================================================
0834 void IntCurveSurface_Inter::InternalPerformCurveQuadric(const TheCurve& curve,
0835 const TheSurface& surface) {
0836 IntCurveSurface_TheQuadCurvExactInter QuadCurv(surface,curve);
0837 if(QuadCurv.IsDone()) {
0838 Standard_Integer NbRoots = QuadCurv.NbRoots();
0839 Standard_Real u,v,w;
0840 for(Standard_Integer i = 1; i<= NbRoots; i++) {
0841 w = QuadCurv.Root(i);
0842 IntCurveSurface_ComputeParamsOnQuadric(surface,TheCurveTool::Value(curve,w),u,v);
0843 AppendPoint(curve,w,surface,u,v);
0844 }
0845 //-- Intervals non traites .............................................
0846 }
0847 }
0848
0849 //=======================================================================
0850 //function : InternalPerform
0851 //purpose :
0852 //=======================================================================
0853 void IntCurveSurface_Inter::InternalPerform(const TheCurve& curve,
0854 const IntCurveSurface_ThePolygon& polygon,
0855 const TheSurface& surface,
0856 const Standard_Real U1,
0857 const Standard_Real V1,
0858 const Standard_Real U2,
0859 const Standard_Real V2) {
0860 GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
0861 if( (SurfaceType != GeomAbs_Plane)
0862 && (SurfaceType != GeomAbs_Cylinder)
0863 && (SurfaceType != GeomAbs_Cone)
0864 && (SurfaceType != GeomAbs_Sphere) ) {
0865 if(SurfaceType != GeomAbs_BSplineSurface) {
0866 Standard_Integer nbsu,nbsv;
0867 nbsu = TheSurfaceTool::NbSamplesU(surface,U1,U2);
0868 nbsv = TheSurfaceTool::NbSamplesV(surface,V1,V2);
0869 if(nbsu>40) nbsu=40;
0870 if(nbsv>40) nbsv=40;
0871 IntCurveSurface_ThePolyhedron polyhedron(surface,nbsu,nbsv,U1,V1,U2,V2);
0872 InternalPerform(curve,polygon,surface,polyhedron,U1,V1,U2,V2);
0873 }
0874 else {
0875 Handle(Adaptor3d_Surface) aS = TheSurfaceTool::UTrim(surface, U1, U2, 1.e-9);
0876 aS = aS->VTrim(V1, V2, 1.e-9);
0877 Handle(Adaptor3d_TopolTool) aTopTool = new Adaptor3d_TopolTool(aS);
0878 Standard_Real defl = 0.1;
0879 aTopTool->SamplePnts(defl, 10, 10);
0880
0881 Standard_Integer nbpu = aTopTool->NbSamplesU();
0882 Standard_Integer nbpv = aTopTool->NbSamplesV();
0883 TColStd_Array1OfReal Upars(1, nbpu), Vpars(1, nbpv);
0884 aTopTool->UParameters(Upars);
0885 aTopTool->VParameters(Vpars);
0886
0887 IntCurveSurface_ThePolyhedron polyhedron(surface,Upars, Vpars);
0888 InternalPerform(curve,polygon,surface,polyhedron,U1,V1,U2,V2);
0889 }
0890 }
0891 else {
0892 IntCurveSurface_TheQuadCurvExactInter QuadCurv(surface,curve);
0893 if(QuadCurv.IsDone()) {
0894 Standard_Integer NbRoots = QuadCurv.NbRoots();
0895 Standard_Real u,v,w;
0896 for(Standard_Integer i = 1; i<= NbRoots; i++) {
0897 w = QuadCurv.Root(i);
0898 IntCurveSurface_ComputeParamsOnQuadric(surface,TheCurveTool::Value(curve,w),u,v);
0899 AppendPoint(curve,w,surface,u,v);
0900 }
0901 //-- Intervalles non traites .............................................
0902 }
0903 } //-- Fin : la Surface est une quadrique
0904 }
0905 //=======================================================================
0906 //function : PerformConicSurf
0907 //purpose :
0908 //=======================================================================
0909 void IntCurveSurface_Inter::PerformConicSurf(const gp_Lin& Line,
0910 const TheCurve& curve,
0911 const TheSurface& surface,
0912 const Standard_Real U1,
0913 const Standard_Real V1,
0914 const Standard_Real U2,
0915 const Standard_Real V2) {
0916
0917
0918 GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
0919 Standard_Boolean isAnaProcessed = Standard_True;
0920 switch(SurfaceType) {
0921 case GeomAbs_Plane:
0922 {
0923 IntAna_IntConicQuad LinPlane(Line,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE);
0924 AppendIntAna(curve,surface,LinPlane);
0925 break;
0926 }
0927 case GeomAbs_Cylinder:
0928 {
0929 IntAna_IntConicQuad LinCylinder(Line,TheSurfaceTool::Cylinder(surface));
0930 AppendIntAna(curve,surface,LinCylinder);
0931 break;
0932 }
0933 case GeomAbs_Sphere:
0934 {
0935 IntAna_IntConicQuad LinSphere(Line,TheSurfaceTool::Sphere(surface));
0936 AppendIntAna(curve,surface,LinSphere);
0937 break;
0938 }
0939 case GeomAbs_Torus:
0940 {
0941 IntAna_IntLinTorus intlintorus(Line, TheSurfaceTool::Torus(surface));
0942 if (intlintorus.IsDone()) {
0943 Standard_Integer nbp = intlintorus.NbPoints();
0944 Standard_Real fi, theta, w;
0945 for (Standard_Integer i = 1; i <= nbp; i++) {
0946 const gp_Pnt aDebPnt(intlintorus.Value(i));
0947 (void)aDebPnt;
0948 w = intlintorus.ParamOnLine(i);
0949 intlintorus.ParamOnTorus(i, fi, theta);
0950 AppendPoint(curve, w, surface, fi, theta);
0951 }
0952 }
0953 else
0954 isAnaProcessed = Standard_False;
0955 break;
0956 }
0957 case GeomAbs_Cone:
0958 {
0959 const Standard_Real correction = 1.E+5*Precision::Angular();
0960 gp_Cone cn = TheSurfaceTool::Cone(surface);
0961 if(Abs(cn.SemiAngle()) < M_PI/2.0 - correction) {
0962 IntAna_IntConicQuad LinCone(Line, cn);
0963 AppendIntAna(curve, surface, LinCone);
0964 }
0965 else
0966 isAnaProcessed = Standard_False;
0967 break;
0968 }
0969 default:
0970 isAnaProcessed = Standard_False;
0971 }
0972 if (!isAnaProcessed)
0973 {
0974 Standard_Integer nbsu,nbsv;
0975 nbsu = TheSurfaceTool::NbSamplesU(surface,U1,U2);
0976 nbsv = TheSurfaceTool::NbSamplesV(surface,V1,V2);
0977
0978 Standard_Boolean U1inf = Precision::IsInfinite(U1);
0979 Standard_Boolean U2inf = Precision::IsInfinite(U2);
0980 Standard_Boolean V1inf = Precision::IsInfinite(V1);
0981 Standard_Boolean V2inf = Precision::IsInfinite(V2);
0982
0983 Standard_Real U1new=U1, U2new=U2, V1new=V1, V2new=V2;
0984
0985 Standard_Boolean NoIntersection = Standard_False;
0986
0987 if(U1inf || U2inf || V1inf || V2inf ) {
0988
0989 if(SurfaceType == GeomAbs_SurfaceOfExtrusion) {
0990
0991 EstLimForInfExtr(Line, surface, Standard_False, nbsu,
0992 U1inf, U2inf, V1inf, V2inf,
0993 U1new, U2new, V1new, V2new, NoIntersection);
0994
0995 }
0996 else if (SurfaceType == GeomAbs_SurfaceOfRevolution) {
0997
0998 EstLimForInfRevl(Line, surface,
0999 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 //-- std::cout<<" Courbe Dans la Quadrique !!! Non Traite !!!"<<std::endl;
1255 myIsParallel = Standard_True;
1256 }
1257 else if(intana_ConicQuad.IsParallel()) {
1258 //-- std::cout<<" Courbe // a la Quadrique !!! Non Traite !!!"<<std::endl;
1259 myIsParallel = Standard_True;
1260 }
1261 else {
1262 Standard_Integer nbp = intana_ConicQuad.NbPoints();
1263 Standard_Real u,v,w;
1264 for(Standard_Integer i = 1; i<= nbp; i++) {
1265 gp_Pnt P(intana_ConicQuad.Point(i));
1266 w = intana_ConicQuad.ParamOnConic(i);
1267 IntCurveSurface_ComputeParamsOnQuadric(surface,P,u,v);
1268 AppendPoint(curve,w,surface,u,v);
1269 }
1270 }
1271 }
1272 else {
1273 //-- std::cout<<" IntAna Conic Quad Not Done "<<std::endl;
1274 }
1275 }
1276 //=======================================================================
1277 //function : AppendPoint
1278 //purpose :
1279 //=======================================================================
1280 void IntCurveSurface_Inter::AppendPoint(const TheCurve& curve,
1281 const Standard_Real lw,
1282 const TheSurface& surface,
1283 const Standard_Real su,
1284 const Standard_Real sv) {
1285
1286 Standard_Real W0 = TheCurveTool::FirstParameter(curve);
1287 Standard_Real W1 = TheCurveTool::LastParameter(curve);
1288 Standard_Real U0 = TheSurfaceTool::FirstUParameter(surface);
1289 Standard_Real U1 = TheSurfaceTool::LastUParameter(surface);
1290 Standard_Real V0 = TheSurfaceTool::FirstVParameter(surface);
1291 Standard_Real V1 = TheSurfaceTool::LastVParameter(surface);
1292 //-- Test si la courbe est periodique
1293 Standard_Real w = lw, u = su, v = sv;
1294
1295 GeomAbs_CurveType aCType = TheCurveTool::GetType(curve);
1296
1297 if(TheCurveTool::IsPeriodic(curve)
1298 || aCType == GeomAbs_Circle
1299 || aCType == GeomAbs_Ellipse) {
1300 w = ElCLib::InPeriod(w, W0, W0 + TheCurveTool::Period(curve));
1301 }
1302
1303 if((W0 - w) >= TOLTANGENCY || (w - W1) >= TOLTANGENCY) return;
1304
1305 GeomAbs_SurfaceType aSType = TheSurfaceTool::GetType(surface);
1306 if (TheSurfaceTool::IsUPeriodic(surface)
1307 || aSType == GeomAbs_Cylinder
1308 || aSType == GeomAbs_Cone
1309 || aSType == GeomAbs_Sphere) {
1310 u = ElCLib::InPeriod(u, U0, U0 + TheSurfaceTool::UPeriod(surface));
1311 }
1312
1313 if (TheSurfaceTool::IsVPeriodic(surface)) {
1314 v = ElCLib::InPeriod(v, V0, V0 + TheSurfaceTool::VPeriod(surface));
1315 }
1316
1317 if((U0 - u) >= TOLTANGENCY || (u - U1) >= TOLTANGENCY) return;
1318 if((V0 - v) >= TOLTANGENCY || (v - V1) >= TOLTANGENCY) return;
1319
1320 IntCurveSurface_TransitionOnCurve TransOnCurve;
1321 IntCurveSurface_ComputeTransitions(curve,w,TransOnCurve,
1322 surface,u,v);
1323 gp_Pnt P(TheCurveTool::Value(curve,w));
1324 IntCurveSurface_IntersectionPoint IP(P,u,v,w,TransOnCurve);
1325 Append(IP); //-- invoque la methode de IntCurveSurface_Intersection.
1326
1327 }
1328
1329 //=======================================================================
1330 //function : AppendSegment
1331 //purpose :
1332 //=======================================================================
1333 void IntCurveSurface_Inter::AppendSegment(const TheCurve& ,
1334 const Standard_Real ,
1335 const Standard_Real ,
1336 const TheSurface& ) {
1337 //std::cout<<" !!! Not Yet Implemented
1338 //IntCurveSurface_Inter::Append(const IntCurveSurf ...)"<<std::endl;
1339 }
1340
1341 //=======================================================================
1342 //function : SectionPointToParameters
1343 //purpose : P o i n t d i n t e r f e r e n c e - - >
1344 // U , V e t W
1345 //=======================================================================
1346 void SectionPointToParameters(const Intf_SectionPoint& Sp,
1347 const IntCurveSurface_ThePolyhedron& Polyhedron,
1348 const IntCurveSurface_ThePolygon& Polygon,
1349 Standard_Real& U,
1350 Standard_Real& V,
1351 Standard_Real& W) {
1352
1353 Intf_PIType typ;
1354 Standard_Integer Adr1,Adr2;
1355 Standard_Real Param,u,v;
1356 gp_Pnt P(Sp.Pnt());
1357
1358 Standard_Integer Pt1,Pt2,Pt3;
1359 Standard_Real u1 = 0.,v1 = 0.,param;
1360 //----------------------------------------------------------------------
1361 //-- Calcul des parametres approches sur la surface --
1362 //----------------------------------------------------------------------
1363
1364 Sp.InfoSecond(typ,Adr1,Adr2,Param);
1365 switch(typ) {
1366 case Intf_VERTEX: //-- Adr1 est le numero du vertex
1367 {
1368 Polyhedron.Parameters(Adr1,u1,v1);
1369 break;
1370 }
1371 case Intf_EDGE:
1372 {
1373 Polyhedron.Parameters(Adr1,u1,v1);
1374 Polyhedron.Parameters(Adr2,u,v);
1375 u1+= Param * (u-u1);
1376 v1+= Param * (v-v1);
1377 break;
1378 }
1379 case Intf_FACE:
1380 {
1381 Standard_Real ua,va,ub,vb,uc,vc,ca,cb,cc,cabc;
1382 Polyhedron.Triangle(Adr1,Pt1,Pt2,Pt3);
1383 gp_Pnt PA(Polyhedron.Point(Pt1));
1384 gp_Pnt PB(Polyhedron.Point(Pt2));
1385 gp_Pnt PC(Polyhedron.Point(Pt3));
1386 Polyhedron.Parameters(Pt1,ua,va);
1387 Polyhedron.Parameters(Pt2,ub,vb);
1388 Polyhedron.Parameters(Pt3,uc,vc);
1389 gp_Vec Normale(gp_Vec(PA,PB).Crossed(gp_Vec(PA,PC)));
1390 cc = (gp_Vec(PA,PB).Crossed(gp_Vec(PA,P))).Dot(Normale);
1391 ca = (gp_Vec(PB,PC).Crossed(gp_Vec(PB,P))).Dot(Normale);
1392 cb = (gp_Vec(PC,PA).Crossed(gp_Vec(PC,P))).Dot(Normale);
1393 cabc = ca + cb + cc;
1394
1395 ca/=cabc; cb/=cabc; cc/=cabc;
1396
1397 u1 = ca * ua + cb * ub + cc * uc;
1398 v1 = ca * va + cb * vb + cc * vc;
1399 break;
1400 }
1401 default:
1402 {
1403 std::cout<<" Default dans SectionPointToParameters "<<std::endl;
1404 break;
1405 }
1406 }
1407 //----------------------------------------------------------------------
1408 //-- Calcul du point approche sur la Curve --
1409 //----------------------------------------------------------------------
1410 Standard_Integer SegIndex;
1411
1412 Sp.InfoFirst(typ,SegIndex,param);
1413 W = Polygon.ApproxParamOnCurve(SegIndex,param);
1414 if(param>1.0 || param<0.0) {
1415 //-- IntCurveSurface_ThePolyhedronTool::Dump(Polyhedron);
1416 //-- IntCurveSurface_ThePolygonTool::Dump(Polygon);
1417 }
1418 U = u1;
1419 V = v1;
1420 }
1421 //=======================================================================
1422 //function : IntCurveSurface_ComputeTransitions
1423 //purpose :
1424 //=======================================================================
1425 void IntCurveSurface_ComputeTransitions(const TheCurve& curve,
1426 const Standard_Real w,
1427 IntCurveSurface_TransitionOnCurve& TransOnCurve,
1428 const TheSurface& surface,
1429 const Standard_Real u,
1430 const Standard_Real v) {
1431
1432 gp_Vec NSurf,D1U,D1V;//TgCurv;
1433 gp_Pnt Psurf;
1434 Standard_Real CosDir;
1435
1436
1437 TheSurfaceTool::D1(surface,u,v,Psurf,D1U,D1V);
1438 NSurf = D1U.Crossed(D1V);
1439 TheCurveTool::D1(curve,w,Psurf,D1U);
1440 Standard_Real Norm = NSurf.Magnitude();
1441 if(Norm>TOLERANCE_ANGULAIRE &&
1442 D1U.SquareMagnitude() > TOLERANCE_ANGULAIRE) {
1443 D1U.Normalize();
1444 CosDir = NSurf.Dot(D1U);
1445 CosDir/=Norm;
1446 if( -CosDir > TOLERANCE_ANGULAIRE) {
1447 //-- --Curve---> <----Surface----
1448 TransOnCurve = IntCurveSurface_In;
1449 }
1450 else if(CosDir > TOLERANCE_ANGULAIRE) {
1451 //-- --Curve---> ----Surface-->
1452 TransOnCurve = IntCurveSurface_Out;
1453 }
1454 else {
1455 TransOnCurve = IntCurveSurface_Tangent;
1456 }
1457 }
1458 else {
1459 TransOnCurve = IntCurveSurface_Tangent;
1460 }
1461 }
1462 //=======================================================================
1463 //function : IntCurveSurface_ComputeParamsOnQuadric
1464 //purpose :
1465 //=======================================================================
1466 void IntCurveSurface_ComputeParamsOnQuadric(const TheSurface& surface,
1467 const gp_Pnt& P,
1468 Standard_Real& u,
1469 Standard_Real& v) {
1470 GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1471 switch(SurfaceType) {
1472 case GeomAbs_Plane:
1473 {
1474 ElSLib::Parameters(TheSurfaceTool::Plane(surface),P,u,v);
1475 break;
1476 }
1477 case GeomAbs_Cylinder:
1478 {
1479 ElSLib::Parameters(TheSurfaceTool::Cylinder(surface),P,u,v);
1480 break;
1481 }
1482 case GeomAbs_Cone:
1483 {
1484 ElSLib::Parameters(TheSurfaceTool::Cone(surface),P,u,v);
1485 break;
1486 }
1487 case GeomAbs_Sphere:
1488 {
1489 ElSLib::Parameters(TheSurfaceTool::Sphere(surface),P,u,v);
1490 break;
1491 }
1492 default:
1493 break;
1494 }
1495 }
1496 //=======================================================================
1497 //function : EstLimForInfExtr
1498 //purpose : Estimation of limits for infinite surfaces
1499 //=======================================================================
1500 void EstLimForInfExtr(const gp_Lin& Line,
1501 const TheSurface& surface,
1502 const Standard_Boolean IsOffSurf,
1503 const Standard_Integer nbsu,
1504 const Standard_Boolean U1inf,
1505 const Standard_Boolean U2inf,
1506 const Standard_Boolean V1inf,
1507 const Standard_Boolean V2inf,
1508 Standard_Real& U1new,
1509 Standard_Real& U2new,
1510 Standard_Real& V1new,
1511 Standard_Real& V2new,
1512 Standard_Boolean& NoIntersection)
1513
1514 {
1515
1516 NoIntersection = Standard_False;
1517
1518 Handle(Adaptor3d_Surface) aBasSurf;
1519
1520 if(IsOffSurf) aBasSurf = TheSurfaceTool::BasisSurface(surface);
1521
1522 gp_Dir aDirOfExt;
1523
1524 if(IsOffSurf) aDirOfExt = aBasSurf->Direction();
1525 else aDirOfExt = TheSurfaceTool::Direction(surface);
1526
1527 Standard_Real tolang = TOLERANCE_ANGULAIRE;
1528
1529 if(aDirOfExt.IsParallel(Line.Direction(), tolang)) {
1530 NoIntersection = Standard_True;
1531 return;
1532 }
1533
1534 if((V1inf || V2inf) && !(U1inf || U2inf)) {
1535
1536 Standard_Real vmin = RealLast(), vmax = -vmin;
1537 gp_Lin aL;
1538 Standard_Real step = (U2new - U1new) / nbsu;
1539 Standard_Real u = U1new, v;
1540 gp_Pnt aP;
1541 Extrema_POnCurv aP1, aP2;
1542 Standard_Integer i;
1543
1544 for(i = 0; i <= nbsu; i++) {
1545
1546 TheSurfaceTool::D0(surface, u, 0., aP);
1547 aL.SetLocation(aP);
1548 aL.SetDirection(aDirOfExt);
1549
1550 Extrema_ExtElC aExtr(aL, Line, tolang);
1551
1552 if(!aExtr.IsDone()) return;
1553
1554 if(aExtr.IsParallel()) {
1555 NoIntersection = Standard_True;
1556 return;
1557 }
1558
1559 aExtr.Points(1, aP1, aP2);
1560 v = aP1.Parameter();
1561 vmin = Min(vmin, v);
1562 vmax = Max(vmax, v);
1563
1564 u += step;
1565
1566 }
1567
1568 vmin = vmin - Abs(vmin) - 10.;
1569 vmax = vmax + Abs(vmax) + 10.;
1570
1571 V1new = Max(V1new, vmin);
1572 V2new = Min(V2new, vmax);
1573
1574 }
1575 else if(U1inf || U2inf) {
1576
1577 Standard_Real umin = RealLast(), umax = -umin;
1578 Standard_Real u0 = Min(Max(0., U1new), U2new);
1579 Standard_Real v0 = Min(Max(0., V1new), V2new);
1580 gp_Pnt aP;
1581 TheSurfaceTool::D0(surface, u0, v0, aP);
1582 gp_Pln aRefPln(aP, aDirOfExt);
1583
1584 Handle(Adaptor3d_Curve) aBasCurv;
1585
1586 if(IsOffSurf) aBasCurv = aBasSurf->BasisCurve();
1587 else aBasCurv = TheSurfaceTool::BasisCurve(surface);
1588
1589 ProjLib_Plane Projector(aRefPln);
1590
1591 Projector.Project(Line);
1592
1593 if(!Projector.IsDone()) return;
1594
1595 gp_Lin2d Line2d = Projector.Line();
1596
1597 GeomAbs_CurveType aCurvTyp = aBasCurv->GetType();
1598
1599 if(aCurvTyp == GeomAbs_Line) {
1600
1601 Projector.Project(aBasCurv->Line());
1602
1603 if(!Projector.IsDone()) return;
1604
1605 gp_Lin2d aL2d = Projector.Line();
1606
1607 IntAna2d_AnaIntersection anInter(Line2d, aL2d);
1608
1609 if(!anInter.IsDone()) return;
1610
1611 if(anInter.IsEmpty() || anInter.IdenticalElements() ||
1612 anInter.ParallelElements() ) {
1613 NoIntersection = Standard_True;
1614 return;
1615 }
1616
1617 const IntAna2d_IntPoint& anIntPnt = anInter.Point(1);
1618 umin = umax = anIntPnt.ParamOnSecond();
1619 }
1620 else if(aCurvTyp == GeomAbs_Parabola || aCurvTyp == GeomAbs_Hyperbola) {
1621
1622 IntAna2d_Conic aCon(Line2d);
1623 IntAna2d_AnaIntersection anInter;
1624
1625 if(aCurvTyp == GeomAbs_Parabola) {
1626 Projector.Project(aBasCurv->Parabola());
1627 if(!Projector.IsDone()) return;
1628
1629 const gp_Parab2d& aP2d = Projector.Parabola();
1630
1631 anInter.Perform(aP2d, aCon);
1632 }
1633 else {
1634 Projector.Project(aBasCurv->Hyperbola());
1635 if(!Projector.IsDone()) return;
1636
1637 const gp_Hypr2d& aH2d = Projector.Hyperbola();
1638 anInter.Perform(aH2d, aCon);
1639 }
1640
1641 if(!anInter.IsDone()) return;
1642
1643 if(anInter.IsEmpty()) {
1644 NoIntersection = Standard_True;
1645 return;
1646 }
1647
1648 Standard_Integer i, nbint = anInter.NbPoints();
1649 for(i = 1; i <= nbint; i++) {
1650
1651 const IntAna2d_IntPoint& anIntPnt = anInter.Point(i);
1652
1653 umin = Min(anIntPnt.ParamOnFirst(), umin);
1654 umax = Max(anIntPnt.ParamOnFirst(), umax);
1655
1656 }
1657
1658 }
1659 else {
1660 return;
1661 }
1662
1663 umin = umin - Abs(umin) - 10;
1664 umax = umax + Abs(umax) + 10;
1665
1666 U1new = Max(U1new, umin);
1667 U2new = Min(U2new, umax);
1668
1669 if(V1inf || V2inf) {
1670 EstLimForInfExtr(Line, surface, IsOffSurf, nbsu,
1671 Standard_False, Standard_False, V1inf, V2inf,
1672 U1new, U2new, V1new, V2new, NoIntersection);
1673 }
1674 }
1675
1676 return;
1677
1678 }
1679 //=======================================================================
1680 //function : ProjectIntersectAndEstLim
1681 //purpose : project <theLine> and it's X-axe symmetric line to <thePln> and
1682 // intersect resulting curve with <theBasCurvProj>.
1683 // Then estimate max and min parameters of intersection on
1684 // <theBasCurvProj>.
1685 // Is called from EstLimForInfRevl()
1686 //=======================================================================
1687 void ProjectIntersectAndEstLim(const gp_Lin& theLine,
1688 const gp_Pln& thePln,
1689 const ProjLib_Plane& theBasCurvProj,
1690 Standard_Real& theVmin,
1691 Standard_Real& theVmax,
1692 Standard_Boolean& theNoIntersection)
1693 {
1694 ProjLib_Plane aLineProj( thePln, theLine );
1695 if (!aLineProj.IsDone()) {
1696 #ifdef OCCT_DEBUG
1697 std::cout
1698 << "Info: IntCurveSurface_Inter::ProjectIntersectAndEstLim(), !aLineProj.IsDone()"
1699 << std::endl;
1700 #endif
1701 return;
1702 }
1703 gp_Lin2d aLin2d = aLineProj.Line();
1704
1705 // make a second line X-axe symmetric to the first one
1706 gp_Pnt2d aP1 = aLin2d.Location();
1707 gp_Pnt2d aP2 (aP1.XY() + aLin2d.Direction().XY());
1708 gp_Pnt2d aP1sym (aP1.X(), -aP1.Y());
1709 gp_Pnt2d aP2sym (aP2.X(), -aP2.Y());
1710 gp_Lin2d aLin2dsym (aP1sym, gp_Vec2d(aP1sym,aP2sym));
1711
1712 // intersect projections
1713 IntAna2d_Conic aCon (aLin2d);
1714 IntAna2d_Conic aConSym (aLin2dsym);
1715 IntAna2d_AnaIntersection anIntersect;
1716 IntAna2d_AnaIntersection anIntersectSym;
1717
1718 switch (theBasCurvProj.GetType()) {
1719 case GeomAbs_Line:
1720 anIntersectSym.Perform(theBasCurvProj.Line(), aConSym);
1721 anIntersect.Perform(theBasCurvProj.Line(), aCon); break;
1722 case GeomAbs_Hyperbola:
1723 anIntersectSym.Perform(theBasCurvProj.Hyperbola(), aConSym);
1724 anIntersect.Perform(theBasCurvProj.Hyperbola(), aCon); break;
1725 case GeomAbs_Parabola:
1726 anIntersectSym.Perform(theBasCurvProj.Parabola(), aConSym);
1727 anIntersect.Perform(theBasCurvProj.Parabola(), aCon); break;
1728 default:
1729 return; // not infinite curve
1730 }
1731
1732 // retrieve params of intersections
1733 Standard_Integer aNbIntPnt = anIntersect.IsDone() ? anIntersect.NbPoints() : 0;
1734 Standard_Integer aNbIntPntSym = anIntersectSym.IsDone() ? anIntersectSym.NbPoints() : 0;
1735 Standard_Integer iPnt, aNbPnt = Max (aNbIntPnt, aNbIntPntSym);
1736
1737 if (aNbPnt == 0) {
1738 theNoIntersection = Standard_True;
1739 return;
1740 }
1741 Standard_Real aParam;
1742 for (iPnt = 1; iPnt <= aNbPnt; iPnt++)
1743 {
1744 if (iPnt <= aNbIntPnt) {
1745 const IntAna2d_IntPoint& aIntPnt = anIntersect.Point(iPnt);
1746 aParam = aIntPnt.ParamOnFirst();
1747 theVmin = Min (theVmin, aParam );
1748 theVmax = Max (theVmax, aParam );
1749 }
1750 if (iPnt <= aNbIntPntSym) {
1751 const IntAna2d_IntPoint& aIntPnt = anIntersectSym.Point(iPnt);
1752 aParam = aIntPnt.ParamOnFirst();
1753 theVmin = Min (theVmin, aParam );
1754 theVmax = Max (theVmax, aParam );
1755 }
1756 }
1757
1758 return;
1759 }
1760 //=======================================================================
1761 //function : EstLimForInfRevl
1762 //purpose : Estimate V1 and V2 to pass to InternalPerform() if they are
1763 // infinite for Surface of Revolution
1764 // Algo: intersect projections of Line and basis curve on the
1765 // plane passing through revolution axe
1766 //=======================================================================
1767 void EstLimForInfRevl(const gp_Lin& Line,
1768 const TheSurface& surface,
1769 const Standard_Boolean U1inf,
1770 const Standard_Boolean U2inf,
1771 const Standard_Boolean V1inf,
1772 const Standard_Boolean V2inf,
1773 Standard_Real& U1new,
1774 Standard_Real& U2new,
1775 Standard_Real& V1new,
1776 Standard_Real& V2new,
1777 Standard_Boolean& NoIntersection)
1778 {
1779
1780 NoIntersection = Standard_False;
1781
1782 if (U1inf || U2inf) {
1783 if (U1inf)
1784 U1new = Max (0., U1new);
1785 else
1786 U2new = Min (2 * M_PI, U2new);
1787 if (! V1inf && !V2inf) return;
1788 }
1789
1790 Handle(Adaptor3d_Curve) aBasisCurve = TheSurfaceTool::BasisCurve(surface);
1791 gp_Ax1 aRevAx = TheSurfaceTool::AxeOfRevolution(surface);
1792 gp_Vec aXVec = aRevAx.Direction();
1793 Standard_Real aTolAng = Precision::Angular();
1794
1795 // make plane to project a basis curve
1796 gp_Pnt O = aRevAx.Location();
1797 Standard_Real aU = 0.;
1798 gp_Pnt P = aBasisCurve->Value(aU);
1799 while (O.SquareDistance(P) <= Precision::PConfusion() ||
1800 aXVec.IsParallel( gp_Vec(O,P), aTolAng)) {
1801 aU += 1.;
1802 P = aBasisCurve->Value(aU);
1803 if (aU > 3)
1804 // basis curve is a line coinciding with aXVec, P is any not on aXVec
1805 P = gp_Pnt(aU, aU+1, aU+2);
1806 }
1807 gp_Vec aNVec = aXVec ^ gp_Vec(O,P);
1808 gp_Pln aPln (gp_Ax3 (O, aNVec ,aXVec));
1809
1810 // project basic curve
1811 ProjLib_Plane aBasCurvProj(aPln);
1812 switch (aBasisCurve->GetType()) {
1813 case GeomAbs_Line:
1814 aBasCurvProj.Project(aBasisCurve->Line ()); break;
1815 case GeomAbs_Hyperbola:
1816 aBasCurvProj.Project(aBasisCurve->Hyperbola()); break;
1817 case GeomAbs_Parabola:
1818 aBasCurvProj.Project(aBasisCurve->Parabola ()); break;
1819 default:
1820 return; // not infinite curve
1821 }
1822 if (!aBasCurvProj.IsDone()) {
1823 #ifdef OCCT_DEBUG
1824 std::cout << "Info: IntCurveSurface_Inter::EstLimForInfRevl(), !aBasCurvProj.IsDone()" << std::endl;
1825 #endif
1826 return;
1827 }
1828 // make plane to project Line
1829 if (aXVec.IsParallel( Line.Direction(), aTolAng)) {
1830 P = Line.Location();
1831 while (O.SquareDistance(P) <= Precision::PConfusion()) {
1832 aU += 1.;
1833 P = gp_Pnt(aU, aU+1, aU+2); // any not on aXVec
1834 }
1835 aNVec = aXVec ^ gp_Vec(O,P);
1836 } else
1837 aNVec = aXVec.Crossed( Line.Direction() );
1838
1839 aPln = gp_Pln (gp_Ax3 (O, aNVec ,aXVec));
1840
1841 // make a second plane perpendicular to the first one, rotated around aXVec
1842 gp_Pln aPlnPrp = aPln.Rotated (gp_Ax1 (O,aXVec), M_PI/2.);
1843
1844 // project Line and it's X-axe symmetric one to plane and intersect
1845 // resulting curve with projection of Basic Curev
1846 Standard_Real aVmin = RealLast(), aVmax = -aVmin;
1847 Standard_Boolean aNoInt1 = Standard_False, aNoInt2 = Standard_False;
1848 ProjectIntersectAndEstLim (Line, aPln, aBasCurvProj, aVmin, aVmax, aNoInt1);
1849 ProjectIntersectAndEstLim (Line, aPlnPrp, aBasCurvProj, aVmin, aVmax, aNoInt2);
1850
1851 if (aNoInt1 && aNoInt2) {
1852 NoIntersection = Standard_True;
1853 return;
1854 }
1855
1856 aVmin = aVmin - Abs(aVmin) - 10;
1857 aVmax = aVmax + Abs(aVmax) + 10;
1858
1859 if (V1inf) V1new = aVmin;
1860 if (V2inf) V2new = aVmax;
1861
1862 //std::cout << "EstLimForInfRevl: Vmin " << V1new << " Vmax " << V2new << std::endl;
1863
1864 return;
1865 }
1866
1867 //=======================================================================
1868 //function : EstLimForInfOffs
1869 //purpose :
1870 //=======================================================================
1871 void EstLimForInfOffs(const gp_Lin& Line,
1872 const TheSurface& surface,
1873 const Standard_Integer nbsu,
1874 const Standard_Boolean U1inf,
1875 const Standard_Boolean U2inf,
1876 const Standard_Boolean V1inf,
1877 const Standard_Boolean V2inf,
1878 Standard_Real& U1new,
1879 Standard_Real& U2new,
1880 Standard_Real& V1new,
1881 Standard_Real& V2new,
1882 Standard_Boolean& NoIntersection)
1883 {
1884
1885 NoIntersection = Standard_False;
1886
1887 const Handle(Adaptor3d_Surface)& aBasSurf = TheSurfaceTool::BasisSurface(surface);
1888 Standard_Real anOffVal = TheSurfaceTool::OffsetValue(surface);
1889
1890 GeomAbs_SurfaceType aTypeOfBasSurf = aBasSurf->GetType();
1891
1892 // case for plane, cylinder and cone - make equivalent surface;
1893 if(aTypeOfBasSurf == GeomAbs_Plane) {
1894
1895 gp_Pln aPln = aBasSurf->Plane();
1896 gp_Vec aT = aPln.Position().XDirection()^aPln.Position().YDirection();
1897 aT *= anOffVal;
1898 aPln.Translate(aT);
1899 IntAna_IntConicQuad LinPlane(Line,aPln,TOLERANCE_ANGULAIRE);
1900
1901 if(!LinPlane.IsDone()) return;
1902
1903 if(LinPlane.IsParallel() || LinPlane.IsInQuadric()) {
1904
1905 NoIntersection = Standard_True;
1906 return;
1907
1908 }
1909
1910 Standard_Real u, v;
1911 ElSLib::Parameters(aPln, LinPlane.Point(1), u, v);
1912 U1new = Max(U1new, u - 10.);
1913 U2new = Min(U2new, u + 10.);
1914 V1new = Max(V1new, v - 10.);
1915 V2new = Min(V2new, v + 10.);
1916
1917 }
1918 else if(aTypeOfBasSurf == GeomAbs_Cylinder) {
1919
1920 gp_Cylinder aCyl = aBasSurf->Cylinder();
1921
1922 Standard_Real aR = aCyl.Radius();
1923 gp_Ax3 anA = aCyl.Position();
1924
1925 if (anA.Direct())
1926 aR += anOffVal;
1927 else
1928 aR -= anOffVal;
1929
1930 if ( aR >= TOLTANGENCY ) {
1931 aCyl.SetRadius(aR);
1932 }
1933 else if ( aR <= -TOLTANGENCY ){
1934 anA.Rotate(gp_Ax1(anA.Location(), anA.Direction()), M_PI);
1935 aCyl.SetPosition(anA);
1936 // modified by NIZHNY-MKK Mon Oct 3 17:37:54 2005
1937 // aCyl.SetRadius(Abs(aR));
1938 aCyl.SetRadius(-aR);
1939 }
1940 else {
1941
1942 NoIntersection = Standard_True;
1943 return;
1944
1945 }
1946
1947 IntAna_IntConicQuad LinCylinder(Line, aCyl);
1948
1949 if(!LinCylinder.IsDone()) return;
1950
1951 if(LinCylinder.IsParallel() || LinCylinder.IsInQuadric()) {
1952
1953 NoIntersection = Standard_True;
1954 return;
1955
1956 }
1957
1958 Standard_Integer i, nbp = LinCylinder.NbPoints();
1959 Standard_Real vmin = RealLast(), vmax = -vmin, u, v;
1960
1961 for(i = 1; i <= nbp; i++) {
1962
1963 ElSLib::Parameters(aCyl, LinCylinder.Point(i), u, v);
1964 vmin = Min(vmin, v);
1965 vmax = Max(vmax, v);
1966
1967 }
1968
1969 V1new = Max(V1new, vmin - Abs(vmin) - 10.);
1970 V2new = Min(V2new, vmax + Abs(vmax) + 10.);
1971
1972 }
1973 else if(aTypeOfBasSurf == GeomAbs_Cone) {
1974
1975 gp_Cone aCon = aBasSurf->Cone();
1976 Standard_Real anAng = aCon.SemiAngle();
1977 Standard_Real aR = aCon.RefRadius() + anOffVal * Cos(anAng);
1978 gp_Ax3 anA = aCon.Position();
1979 if ( aR >= 0.) {
1980
1981 gp_Vec aZ( anA.Direction());
1982 aZ *= - anOffVal * Sin(anAng);
1983 anA.Translate(aZ);
1984 aCon.SetPosition(anA);
1985 aCon.SetRadius(aR);
1986 aCon.SetSemiAngle(anAng);
1987
1988 }
1989 else {
1990
1991 return;
1992
1993 }
1994
1995 IntAna_IntConicQuad LinCone(Line, aCon);
1996
1997 if(!LinCone.IsDone()) return;
1998
1999 if(LinCone.IsParallel() || LinCone.IsInQuadric()) {
2000
2001 NoIntersection = Standard_True;
2002 return;
2003
2004 }
2005
2006 Standard_Integer i, nbp = LinCone.NbPoints();
2007 Standard_Real vmin = RealLast(), vmax = -vmin, u, v;
2008
2009 for(i = 1; i <= nbp; i++) {
2010
2011 ElSLib::Parameters(aCon, LinCone.Point(i), u, v);
2012 vmin = Min(vmin, v);
2013 vmax = Max(vmax, v);
2014
2015 }
2016
2017 V1new = Max(V1new, vmin - Abs(vmin) - 10.);
2018 V2new = Min(V2new, vmax + Abs(vmax) + 10.);
2019
2020 }
2021 else if(aTypeOfBasSurf == GeomAbs_SurfaceOfExtrusion) {
2022
2023 Standard_Real anU1 = U1new, anU2 = U2new;
2024
2025 EstLimForInfExtr(Line, surface, Standard_True, nbsu,
2026 U1inf, U2inf, V1inf, V2inf,
2027 U1new, U2new, V1new, V2new, NoIntersection);
2028
2029 if(NoIntersection) return;
2030
2031 if(U1inf || U2inf) {
2032
2033 GeomAbs_CurveType aBasCurvType = aBasSurf->BasisCurve()->GetType();
2034 if(aBasCurvType == GeomAbs_Line) {
2035 U1new = Max(anU1, -1.e10);
2036 U2new = Min(anU2, 1.e10);
2037 }
2038 else if(aBasCurvType == GeomAbs_Parabola) {
2039 gp_Parab aPrb = aBasSurf->BasisCurve()->Parabola();
2040 Standard_Real aF = aPrb.Focal();
2041 Standard_Real dU = 2.e5 * Sqrt(aF);
2042 U1new = Max(anU1, -dU);
2043 U2new = Min(anU2, dU);
2044 }
2045 else if(aBasCurvType == GeomAbs_Hyperbola) {
2046 U1new = Max(anU1, -30.);
2047 U2new = Min(anU2, 30.);
2048 }
2049 else {
2050 U1new = Max(anU1, -1.e10);
2051 U2new = Min(anU2, 1.e10);
2052 }
2053
2054 }
2055
2056 }
2057 else if(aTypeOfBasSurf == GeomAbs_SurfaceOfRevolution) {
2058
2059 GeomAbs_CurveType aBasCurvType = aBasSurf->BasisCurve()->GetType();
2060 if(aBasCurvType == GeomAbs_Line) {
2061 V1new = Max(V1new, -1.e10);
2062 V2new = Min(V2new, 1.e10);
2063 }
2064 else if(aBasCurvType == GeomAbs_Parabola) {
2065 gp_Parab aPrb = aBasSurf->BasisCurve()->Parabola();
2066 Standard_Real aF = aPrb.Focal();
2067 Standard_Real dV = 2.e5 * Sqrt(aF);
2068 V1new = Max(V1new, -dV);
2069 V2new = Min(V2new, dV);
2070 }
2071 else if(aBasCurvType == GeomAbs_Hyperbola) {
2072 V1new = Max(V1new, -30.);
2073 V2new = Min(V2new, 30.);
2074 }
2075 else {
2076 V1new = Max(V1new, -1.e10);
2077 V2new = Min(V2new, 1.e10);
2078 }
2079
2080 }
2081 else {
2082
2083 V1new = Max(V1new, -1.e10);
2084 V2new = Min(V2new, 1.e10);
2085 }
2086 }
2087 //=======================================================================
2088 //function : EstLimForInfSurf
2089 //purpose :
2090 //=======================================================================
2091 void EstLimForInfSurf(Standard_Real& U1new,
2092 Standard_Real& U2new,
2093 Standard_Real& V1new,
2094 Standard_Real& V2new)
2095 {
2096 U1new = Max(U1new, -1.e10);
2097 U2new = Min(U2new, 1.e10);
2098 V1new = Max(V1new, -1.e10);
2099 V2new = Min(V2new, 1.e10);
2100 }
