1 // Created on: 1992-06-30
2 // Created by: Laurent BUCHARD
3 // Copyright (c) 1992-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 #ifndef OCCT_DEBUG
18 #define No_Standard_RangeError
19 #define No_Standard_OutOfRange
20 #endif
21
22 //----------------------------------------------------------------------
23 //-- Differents constructeurs sont proposes qui correspondent aux
24 //-- polynomes en Z :
25 //-- A(Sin(Theta),Cos(Theta)) Z**2
26 //-- + B(Sin(Theta),Cos(Theta)) Z
27 //-- + C(Sin(Theta),Cos(Theta))
28 //--
29 //-- Une Courbe est definie sur un domaine
30 //--
31 //-- Value retourne le point de parametre U(Theta),V(Theta)
32 //-- ou V est la solution du polynome A V**2 + B V + C
33 //-- (Selon les cas, on prend V+ ou V-)
34 //--
35 //-- D1u calcule le vecteur tangent a la courbe
36 //-- et retourne le booleen Standard_False si ce calcul ne peut
37 //-- pas etre mene a bien.
38 //----------------------------------------------------------------------
39
40 #include <algorithm>
41
42 #include <ElSLib.hxx>
43 #include <gp_Cone.hxx>
44 #include <gp_Cylinder.hxx>
45 #include <gp_Pnt.hxx>
46 #include <gp_Vec.hxx>
47 #include <gp_XYZ.hxx>
48 #include <IntAna_Curve.hxx>
49 #include <math_DirectPolynomialRoots.hxx>
50 #include <Precision.hxx>
51 #include <Standard_DomainError.hxx>
52
53 //=======================================================================
54 //function : IntAna_Curve
55 //purpose :
56 //=======================================================================
IntAna_Curve()57 IntAna_Curve::IntAna_Curve()
58 : Z0Cte(0.0),
59 Z0Sin(0.0),
60 Z0Cos(0.0),
61 Z0SinSin(0.0),
62 Z0CosCos(0.0),
63 Z0CosSin(0.0),
64 Z1Cte(0.0),
65 Z1Sin(0.0),
66 Z1Cos(0.0),
67 Z1SinSin(0.0),
68 Z1CosCos(0.0),
69 Z1CosSin(0.0),
70 Z2Cte(0.0),
71 Z2Sin(0.0),
72 Z2Cos(0.0),
73 Z2SinSin(0.0),
74 Z2CosCos(0.0),
75 Z2CosSin(0.0),
76 TwoCurves(Standard_False),
77 TakeZPositive(Standard_False),
78 Tolerance(0.0),
79 DomainInf(0.0),
80 DomainSup(0.0),
81 RestrictedInf(Standard_False),
82 RestrictedSup(Standard_False),
83 firstbounded(Standard_False),
84 lastbounded(Standard_False),
85 typequadric(GeomAbs_OtherSurface),
86 RCyl(0.0),
87 Angle(0.0),
88 myFirstParameter(0.0),
89 myLastParameter(0.0)
90 {
91 }
92 //=======================================================================
93 //function : SetConeQuadValues
94 //purpose : Description de l intersection Cone Quadrique
95 //=======================================================================
SetConeQuadValues(const gp_Cone & Cone,const Standard_Real Qxx,const Standard_Real Qyy,const Standard_Real Qzz,const Standard_Real Qxy,const Standard_Real Qxz,const Standard_Real Qyz,const Standard_Real Qx,const Standard_Real Qy,const Standard_Real Qz,const Standard_Real Q1,const Standard_Real TOL,const Standard_Real DomInf,const Standard_Real DomSup,const Standard_Boolean twocurves,const Standard_Boolean takezpositive)96 void IntAna_Curve::SetConeQuadValues(const gp_Cone& Cone,
97 const Standard_Real Qxx,
98 const Standard_Real Qyy,
99 const Standard_Real Qzz,
100 const Standard_Real Qxy,
101 const Standard_Real Qxz,
102 const Standard_Real Qyz,
103 const Standard_Real Qx,
104 const Standard_Real Qy,
105 const Standard_Real Qz,
106 const Standard_Real Q1,
107 const Standard_Real TOL,
108 const Standard_Real DomInf,
109 const Standard_Real DomSup,
110 const Standard_Boolean twocurves,
111 const Standard_Boolean takezpositive)
112 {
113
114 Ax3 = Cone.Position();
115 RCyl = Cone.RefRadius();
116
117 Angle = Cone.SemiAngle();
118 Standard_Real UnSurTgAngle = 1.0/(Tan(Cone.SemiAngle()));
119
120 typequadric= GeomAbs_Cone;
121
122 TwoCurves = twocurves; //-- deux Z pour un meme parametre
123 TakeZPositive = takezpositive; //-- Prendre sur la courbe le Z Positif
124 //-- ( -B + Sqrt()) et non (-B - Sqrt())
125
126
127 Z0Cte = Q1; //-- Attention On a Z?Cos Cos(t)
128 Z0Sin = 0.0; //-- et Non 2 Z?Cos Cos(t) !!!
129 Z0Cos = 0.0; //-- Ce pour tous les Parametres
130 Z0CosCos = 0.0; //-- ie pas de Coefficient 2
131 Z0SinSin = 0.0; //-- devant les termes CS C S
132 Z0CosSin = 0.0;
133
134 Z1Cte = 2.0*(UnSurTgAngle)*Qz;
135 Z1Sin = Qy+Qy;
136 Z1Cos = Qx+Qx;
137 Z1CosCos = 0.0;
138 Z1SinSin = 0.0;
139 Z1CosSin = 0.0;
140
141 Z2Cte = Qzz * UnSurTgAngle*UnSurTgAngle;
142 Z2Sin = (UnSurTgAngle+UnSurTgAngle)*Qyz;
143 Z2Cos = (UnSurTgAngle+UnSurTgAngle)*Qxz;
144 Z2CosCos = Qxx;
145 Z2SinSin = Qyy;
146 Z2CosSin = Qxy;
147
148 Tolerance = TOL;
149 DomainInf = DomInf;
150 DomainSup = DomSup;
151
152 RestrictedInf = RestrictedSup = Standard_True; //-- Le Domaine est Borne
153 firstbounded = lastbounded = Standard_False;
154
155 myFirstParameter = DomainInf;
156 myLastParameter = (TwoCurves) ? DomainSup + DomainSup - DomainInf :
157 DomainSup;
158 }
159
160 //=======================================================================
161 //function : SetCylinderQuadValues
162 //purpose : Description de l intersection Cylindre Quadrique
163 //=======================================================================
SetCylinderQuadValues(const gp_Cylinder & Cyl,const Standard_Real Qxx,const Standard_Real Qyy,const Standard_Real Qzz,const Standard_Real Qxy,const Standard_Real Qxz,const Standard_Real Qyz,const Standard_Real Qx,const Standard_Real Qy,const Standard_Real Qz,const Standard_Real Q1,const Standard_Real TOL,const Standard_Real DomInf,const Standard_Real DomSup,const Standard_Boolean twocurves,const Standard_Boolean takezpositive)164 void IntAna_Curve::SetCylinderQuadValues(const gp_Cylinder& Cyl,
165 const Standard_Real Qxx,
166 const Standard_Real Qyy,
167 const Standard_Real Qzz,
168 const Standard_Real Qxy,
169 const Standard_Real Qxz,
170 const Standard_Real Qyz,
171 const Standard_Real Qx,
172 const Standard_Real Qy,
173 const Standard_Real Qz,
174 const Standard_Real Q1,
175 const Standard_Real TOL,
176 const Standard_Real DomInf,
177 const Standard_Real DomSup,
178 const Standard_Boolean twocurves,
179 const Standard_Boolean takezpositive)
180 {
181
182 Ax3 = Cyl.Position();
183 RCyl = Cyl.Radius();
184 typequadric= GeomAbs_Cylinder;
185
186 TwoCurves = twocurves; //-- deux Z pour un meme parametre
187 TakeZPositive = takezpositive; //-- Prendre sur la courbe le Z Positif
188 Standard_Real RCylmul2 = RCyl+RCyl; //-- ( -B + Sqrt())
189
190 Z0Cte = Q1;
191 Z0Sin = RCylmul2*Qy;
192 Z0Cos = RCylmul2*Qx;
193 Z0CosCos = Qxx*RCyl*RCyl;
194 Z0SinSin = Qyy*RCyl*RCyl;
195 Z0CosSin = RCyl*RCyl*Qxy;
196
197 Z1Cte = Qz+Qz;
198 Z1Sin = RCylmul2*Qyz;
199 Z1Cos = RCylmul2*Qxz;
200 Z1CosCos = 0.0;
201 Z1SinSin = 0.0;
202 Z1CosSin = 0.0;
203
204 Z2Cte = Qzz;
205 Z2Sin = 0.0;
206 Z2Cos = 0.0;
207 Z2CosCos = 0.0;
208 Z2SinSin = 0.0;
209 Z2CosSin = 0.0;
210
211 Tolerance = TOL;
212 DomainInf = DomInf;
213 DomainSup = DomSup;
214
215 RestrictedInf = RestrictedSup = Standard_True;
216 firstbounded = lastbounded = Standard_False;
217
218 myFirstParameter = DomainInf;
219 myLastParameter = (TwoCurves) ? DomainSup + DomainSup - DomainInf :
220 DomainSup;
221 }
222
223 //=======================================================================
224 //function : IsOpen
225 //purpose :
226 //=======================================================================
IsOpen() const227 Standard_Boolean IntAna_Curve::IsOpen() const
228 {
229 return(RestrictedInf && RestrictedSup);
230 }
231
232 //=======================================================================
233 //function : Domain
234 //purpose :
235 //=======================================================================
Domain(Standard_Real & theFirst,Standard_Real & theLast) const236 void IntAna_Curve::Domain(Standard_Real& theFirst,
237 Standard_Real& theLast) const
238 {
239 if (RestrictedInf && RestrictedSup)
240 {
241 theFirst = myFirstParameter;
242 theLast = myLastParameter;
243 }
244 else
245 {
246 throw Standard_DomainError("IntAna_Curve::Domain");
247 }
248 }
249 //=======================================================================
250 //function : IsConstant
251 //purpose :
252 //=======================================================================
IsConstant() const253 Standard_Boolean IntAna_Curve::IsConstant() const
254 {
255 //-- ??? Pas facile de decider a la seule vue des Param.
256 return(Standard_False);
257 }
258
259 //=======================================================================
260 //function : IsFirstOpen
261 //purpose :
262 //=======================================================================
IsFirstOpen() const263 Standard_Boolean IntAna_Curve::IsFirstOpen() const
264 {
265 return(firstbounded);
266 }
267
268 //=======================================================================
269 //function : IsLastOpen
270 //purpose :
271 //=======================================================================
IsLastOpen() const272 Standard_Boolean IntAna_Curve::IsLastOpen() const
273 {
274 return(lastbounded);
275 }
276 //=======================================================================
277 //function : SetIsFirstOpen
278 //purpose :
279 //=======================================================================
SetIsFirstOpen(const Standard_Boolean Flag)280 void IntAna_Curve::SetIsFirstOpen(const Standard_Boolean Flag)
281 {
282 firstbounded = Flag;
283 }
284
285 //=======================================================================
286 //function : SetIsLastOpen
287 //purpose :
288 //=======================================================================
SetIsLastOpen(const Standard_Boolean Flag)289 void IntAna_Curve::SetIsLastOpen(const Standard_Boolean Flag)
290 {
291 lastbounded = Flag;
292 }
293
294 //=======================================================================
295 //function : InternalUVValue
296 //purpose :
297 //=======================================================================
InternalUVValue(const Standard_Real theta,Standard_Real & Param1,Standard_Real & Param2,Standard_Real & A,Standard_Real & B,Standard_Real & C,Standard_Real & cost,Standard_Real & sint,Standard_Real & SigneSqrtDis) const298 void IntAna_Curve::InternalUVValue(const Standard_Real theta,
299 Standard_Real& Param1,
300 Standard_Real& Param2,
301 Standard_Real& A,
302 Standard_Real& B,
303 Standard_Real& C,
304 Standard_Real& cost,
305 Standard_Real& sint,
306 Standard_Real& SigneSqrtDis) const
307 {
308 const Standard_Real aRelTolp = 1.0+Epsilon(1.0), aRelTolm = 1.0-Epsilon(1.0);
309
310 // Infinitesimal step of increasing curve parameter. See comment below.
311 const Standard_Real aDT = 100.0*Epsilon(DomainSup + DomainSup - DomainInf);
312
313 Standard_Real Theta=theta;
314 Standard_Boolean SecondSolution=Standard_False;
315
316 if ((Theta<DomainInf*aRelTolm) ||
317 ((Theta>DomainSup*aRelTolp) && (!TwoCurves)) ||
318 (Theta>(DomainSup + DomainSup - DomainInf)*aRelTolp))
319 {
320 SigneSqrtDis = 0.;
321 throw Standard_DomainError("IntAna_Curve::Domain");
322 }
323
324 if (Abs(Theta - DomainSup) < aDT)
325 {
326 // Point of Null-discriminant.
327 Theta = DomainSup;
328 }
329 else if (Theta>DomainSup)
330 {
331 Theta = DomainSup + DomainSup - Theta;
332 SecondSolution=Standard_True;
333 }
334
335 Param1=Theta;
336
337 if(!TwoCurves) {
338 SecondSolution=TakeZPositive;
339 }
340 //
341 cost = Cos(Theta);
342 sint = Sin(Theta);
343 const Standard_Real aSin2t = Sin(Theta + Theta);
344 const Standard_Real aCos2t = Cos(Theta + Theta);
345
346 A=Z2Cte+sint*(Z2Sin+sint*Z2SinSin)+cost*(Z2Cos+cost*Z2CosCos)
347 + Z2CosSin*aSin2t;
348
349 const Standard_Real aDA = cost*Z2Sin - sint*Z2Cos +
350 aSin2t*(Z2SinSin - Z2CosCos) +
351 aCos2t*(Z2CosSin * Z2CosSin);
352
353 B=Z1Cte+sint*(Z1Sin+sint*Z1SinSin)+cost*(Z1Cos+cost*Z1CosCos)
354 + Z1CosSin*aSin2t;
355
356 const Standard_Real aDB = Z1Sin*cost - Z1Cos*sint +
357 aSin2t*(Z1SinSin - Z1CosCos) +
358 aCos2t*(Z1CosSin + Z1CosSin);
359
360 C=Z0Cte+sint*(Z0Sin+sint*Z0SinSin)+cost*(Z0Cos+cost*Z0CosCos)
361 + Z0CosSin*aSin2t;
362
363 const Standard_Real aDC = Z0Sin*cost - Z0Cos*sint +
364 aSin2t*(Z0SinSin - Z0CosCos) +
365 aCos2t*(Z0CosSin + Z0CosSin);
366
367 Standard_Real aDiscriminant = B*B-4.0*A*C;
368
369 // We consider that infinitesimal dt = aDT.
370 // Error of discriminant computation is equal to
371 // (d(Disc)/dt)*dt, where 1st derivative d(Disc)/dt = 2*B*aDB - 4*(A*aDC + C*aDA).
372
373 const Standard_Real aTolD = 2.0*aDT*Abs(B*aDB - 2.0*(A*aDC + C*aDA));
374
375 if (aDiscriminant < aTolD)
376 aDiscriminant = 0.0;
377
378 if (Abs(A) <= Precision::PConfusion())
379 {
380 if (Abs(B) <= Precision::PConfusion())
381 {
382 Param2 = 0.0;
383 }
384 else
385 {
386 Param2 = -C / B;
387 }
388 }
389 else
390 {
391 SigneSqrtDis = (SecondSolution) ? Sqrt(aDiscriminant) : -Sqrt(aDiscriminant);
392 Param2 = (-B + SigneSqrtDis) / (A + A);
393 }
394 }
395
396 //=======================================================================
397 //function : Value
398 //purpose :
399 //=======================================================================
Value(const Standard_Real theta)400 gp_Pnt IntAna_Curve::Value(const Standard_Real theta)
401 {
402 Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
403 //
404 A=0.0; B=0.0; C=0.0;
405 U=0.0; V=0.0;
406 sint=0.0; cost=0.0;
407 SigneSqrtDis=0.0;
408 InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
409 //-- checked the parameter U and Raises Domain Error if Error
410 return(InternalValue(U,V));
411 }
412 //=======================================================================
413 //function : D1u
414 //purpose :
415 //=======================================================================
D1u(const Standard_Real theta,gp_Pnt & Pt,gp_Vec & Vec)416 Standard_Boolean IntAna_Curve::D1u(const Standard_Real theta,
417 gp_Pnt& Pt,
418 gp_Vec& Vec)
419 {
420 //-- Pour detecter le cas ou le calcul est impossible
421 Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
422 A=0.0; B=0.0; C=0.0;
423 U=0.0; V=0.0;
424 sint=0.0; cost=0.0;
425 //
426 InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
427 //
428 Pt = Value(theta);
429 if (Abs(A)<1.0e-7 || Abs(SigneSqrtDis)<1.0e-10) return(Standard_False);
430
431
432 //-- Approximation de la derivee (mieux que le calcul mathematique!)
433 Standard_Real dtheta = (DomainSup - DomainInf)*1.0e-6;
434 Standard_Real theta2 = theta+dtheta;
435 if ((theta2<DomainInf) || ((theta2>DomainSup) && (!TwoCurves))
436 || (theta2>(DomainSup + DomainSup - DomainInf + 1.0e-14)))
437 {
438 dtheta = -dtheta;
439 theta2 = theta+dtheta;
440 }
441 gp_Pnt P2 = Value(theta2);
442 dtheta = 1.0/dtheta;
443 Vec.SetCoord((P2.X()-Pt.X())*dtheta,
444 (P2.Y()-Pt.Y())*dtheta,
445 (P2.Z()-Pt.Z())*dtheta);
446
447 return(Standard_True);
448 }
449 //=======================================================================
450 //function : FindParameter
451 //purpose : Projects P to the ALine. Returns the list of parameters as a results
452 // of projection.
453 // Sometimes aline can be self-intersected line (see bug #29807 where
454 // ALine goes through the cone apex).
455 //=======================================================================
FindParameter(const gp_Pnt & theP,TColStd_ListOfReal & theParams) const456 void IntAna_Curve::FindParameter(const gp_Pnt& theP,
457 TColStd_ListOfReal& theParams) const
458 {
459 const Standard_Real aPIpPI = M_PI + M_PI,
460 anEpsAng = 1.e-8,
461 InternalPrecision = 1.e-8, //precision of internal algorithm of values computation
462 aSqTolPrecision = Precision::SquareConfusion(); //for boundary points to check their coincidence with others
463
464 Standard_Real aTheta = 0.0;
465 //
466 switch (typequadric)
467 {
468 case GeomAbs_Cylinder:
469 {
470 Standard_Real aZ;
471 ElSLib::CylinderParameters(Ax3, RCyl, theP, aTheta, aZ);
472 }
473 break;
474
475 case GeomAbs_Cone:
476 {
477 Standard_Real aZ;
478 ElSLib::ConeParameters(Ax3, RCyl, Angle, theP, aTheta, aZ);
479 }
480 break;
481
482 default:
483 return;
484 }
485 //
486 if (!firstbounded && (DomainInf > aTheta) && ((DomainInf - aTheta) <= anEpsAng))
487 {
488 aTheta = DomainInf;
489 }
490 else if (!lastbounded && (aTheta > DomainSup) && ((aTheta - DomainSup) <= anEpsAng))
491 {
492 aTheta = DomainSup;
493 }
494 //
495 if (aTheta < DomainInf)
496 {
497 aTheta = aTheta + aPIpPI;
498 }
499 else if (aTheta > DomainSup)
500 {
501 aTheta = aTheta - aPIpPI;
502 }
503
504 const Standard_Integer aMaxPar = 5;
505 Standard_Real aParams[aMaxPar] = {DomainInf, DomainSup, aTheta,
506 (TwoCurves)? DomainSup + DomainSup - aTheta : RealLast(),
507 (TwoCurves) ? DomainSup + DomainSup - DomainInf : RealLast()};
508
509 std::sort(aParams, aParams + aMaxPar - 1);
510
511 for (Standard_Integer i = 0; i < aMaxPar; i++)
512 {
513 if (aParams[i] > myLastParameter)
514 break;
515
516 if (aParams[i] < myFirstParameter)
517 continue;
518
519 if (i && (aParams[i] - aParams[i - 1]) < Precision::PConfusion())
520 continue;
521
522 Standard_Real U = 0.0, V= 0.0,
523 A = 0.0, B = 0.0, C = 0.0,
524 sint = 0.0, cost = 0.0, SigneSqrtDis = 0.0;
525 InternalUVValue(aParams[i], U, V, A, B, C,
526 cost, sint, SigneSqrtDis);
527 const gp_Pnt aP(InternalValue(U, V));
528
529 Standard_Real aSqTol;
530 if (aParams[i] == aTheta ||
531 (TwoCurves && aParams[i] == DomainSup + DomainSup - aTheta))
532 aSqTol = InternalPrecision;
533 else
534 aSqTol = aSqTolPrecision;
535
536 if (aP.SquareDistance(theP) < aSqTol)
537 {
538 theParams.Append(aParams[i]);
539 }
540 }
541 }
542 //=======================================================================
543 //function : InternalValue
544 //purpose :
545 //=======================================================================
InternalValue(const Standard_Real U,const Standard_Real _V) const546 gp_Pnt IntAna_Curve::InternalValue(const Standard_Real U,
547 const Standard_Real _V) const
548 {
549 //-- std::cout<<" ["<<U<<","<<V<<"]";
550 Standard_Real V = _V;
551 if(V > 100000.0 ) { V= 100000.0; }
552 if(V < -100000.0 ) { V=-100000.0; }
553
554 switch(typequadric) {
555 case GeomAbs_Cone:
556 {
557 //------------------------------------------------
558 //-- Parametrage : X = V * Cos(U) ---
559 //-- Y = V * Sin(U) ---
560 //-- Z = (V-RCyl) / Tan(SemiAngle)--
561 //------------------------------------------------
562 //-- Angle Vaut Cone.SemiAngle()
563 return(ElSLib::ConeValue(U,(V-RCyl)/Sin(Angle),Ax3,RCyl,Angle));
564 }
565 break;
566
567 case GeomAbs_Cylinder:
568 return(ElSLib::CylinderValue(U,V,Ax3,RCyl));
569 case GeomAbs_Sphere:
570 return(ElSLib::SphereValue(U,V,Ax3,RCyl));
571 default:
572 return(gp_Pnt(0.0,0.0,0.0));
573 }
574 }
575
576 //
577 //=======================================================================
578 //function : SetDomain
579 //purpose :
580 //=======================================================================
SetDomain(const Standard_Real theFirst,const Standard_Real theLast)581 void IntAna_Curve::SetDomain(const Standard_Real theFirst,
582 const Standard_Real theLast)
583 {
584 if (theLast <= theFirst)
585 {
586 throw Standard_DomainError("IntAna_Curve::Domain");
587 }
588 //
589 myFirstParameter = theFirst;
590 myLastParameter = theLast;
591 }
592