0031035: Coding - uninitialized class fields reported by Visual Studio Code Analysis
[occt.git] / src / IntAna / IntAna_Curve.cxx
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
227 Standard_Boolean IntAna_Curve::IsOpen() const
228 {
229   return(RestrictedInf && RestrictedSup);
230 }
231
232 //=======================================================================
233 //function : Domain
234 //purpose  : 
235 //=======================================================================
236 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 //=======================================================================
253 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 //=======================================================================
263 Standard_Boolean IntAna_Curve::IsFirstOpen() const
264 {
265   return(firstbounded);
266 }
267
268 //=======================================================================
269 //function : IsLastOpen
270 //purpose  : 
271 //=======================================================================
272 Standard_Boolean IntAna_Curve::IsLastOpen() const
273 {
274   return(lastbounded);
275 }
276 //=======================================================================
277 //function : SetIsFirstOpen
278 //purpose  : 
279 //=======================================================================
280 void IntAna_Curve::SetIsFirstOpen(const Standard_Boolean Flag)
281 {
282   firstbounded = Flag;
283 }
284
285 //=======================================================================
286 //function : SetIsLastOpen
287 //purpose  : 
288 //=======================================================================
289 void IntAna_Curve::SetIsLastOpen(const Standard_Boolean Flag)
290 {
291   lastbounded = Flag;
292 }
293
294 //=======================================================================
295 //function : InternalUVValue
296 //purpose  : 
297 //=======================================================================
298 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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
456 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 //=======================================================================
546 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 //=======================================================================
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 }