0031035: Coding - uninitialized class fields reported by Visual Studio Code Analysis
[occt.git] / src / IntAna / IntAna_IntConicQuad.cxx
1 // Created on: 1992-07-27
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 #define CREATE  IntAna_IntConicQuad::IntAna_IntConicQuad
23 #define PERFORM void IntAna_IntConicQuad::Perform
24
25
26 #include <ElCLib.hxx>
27 #include <gp_Ax3.hxx>
28 #include <gp_Circ.hxx>
29 #include <gp_Circ2d.hxx>
30 #include <gp_Elips.hxx>
31 #include <gp_Hypr.hxx>
32 #include <gp_Lin.hxx>
33 #include <gp_Lin2d.hxx>
34 #include <gp_Parab.hxx>
35 #include <gp_Pln.hxx>
36 #include <gp_Pnt.hxx>
37 #include <gp_Vec.hxx>
38 #include <IntAna2d_AnaIntersection.hxx>
39 #include <IntAna2d_IntPoint.hxx>
40 #include <IntAna_IntConicQuad.hxx>
41 #include <IntAna_QuadQuadGeo.hxx>
42 #include <IntAna_Quadric.hxx>
43 #include <IntAna_ResultType.hxx>
44 #include <math_DirectPolynomialRoots.hxx>
45 #include <math_TrigonometricFunctionRoots.hxx>
46 #include <Standard_DomainError.hxx>
47 #include <Standard_OutOfRange.hxx>
48 #include <StdFail_NotDone.hxx>
49
50 static Standard_Real PIpPI = M_PI + M_PI;
51 //=============================================================================
52 //==                                          E m p t y   C o n s t r u c t o r
53 //== 
54 CREATE(void) {
55   done=Standard_False;
56   parallel = Standard_False;
57   inquadric = Standard_False;
58   nbpts = 0;
59   memset (paramonc, 0, sizeof (paramonc));
60 }
61 //=============================================================================
62 //==                                                 L i n e  -   Q u a d r i c  
63 //==
64 CREATE(const gp_Lin& L,const IntAna_Quadric& Quad) {
65   Perform(L,Quad);
66 }
67
68 PERFORM(const gp_Lin& L,const IntAna_Quadric& Quad) {
69
70   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte;
71   done=inquadric=parallel=Standard_False;
72
73   //----------------------------------------------------------------------
74   //-- Substitution de x=t Lx + Lx0       ( exprime dans                 )
75   //--                 y=t Ly + Ly0       (  le systeme de coordonnees   )
76   //--                 z=t Lz + Lz0       (  canonique                   )
77   //--
78   //-- Dans     Qxx x**2 + Qyy y**2 + Qzz z**2
79   //--          + 2 ( Qxy x y  + Qxz x z  + Qyz y z  )
80   //--          + 2 ( Qx x + Qy y + Qz z )
81   //--          + QCte
82   //--
83   //-- Done un polynome en t : A2 t**2 + A1 t + A0 avec :
84   //----------------------------------------------------------------------
85
86
87   Standard_Real Lx0,Ly0,Lz0,Lx,Ly,Lz;
88
89
90   nbpts=0;
91
92   L.Direction().Coord(Lx,Ly,Lz);
93   L.Location().Coord(Lx0,Ly0,Lz0);
94  
95   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte);
96   
97   Standard_Real A0=(QCte + Qxx*Lx0*Lx0 + Qyy*Ly0*Ly0  + Qzz*Lz0*Lz0
98           + 2.0* (  Lx0*( Qx + Qxy*Ly0 + Qxz*Lz0) 
99                   + Ly0*( Qy + Qyz*Lz0 )
100                   + Qz*Lz0 )
101           );
102           
103   
104   Standard_Real A1=2.0*( Lx*( Qx + Qxx*Lx0 + Qxy*Ly0 + Qxz*Lz0 )
105               +Ly*( Qy + Qxy*Lx0 + Qyy*Ly0 + Qyz*Lz0 )
106               +Lz*( Qz + Qxz*Lx0 + Qyz*Ly0 + Qzz*Lz0 ));
107
108   Standard_Real A2=(Qxx*Lx*Lx + Qyy*Ly*Ly + Qzz*Lz*Lz
109            + 2.0*( Lx*( Qxy*Ly + Qxz*Lz )
110                   + Qyz*Ly*Lz));
111
112   math_DirectPolynomialRoots LinQuadPol(A2,A1,A0);
113   
114   if( LinQuadPol.IsDone()) {
115     done=Standard_True;
116     if(LinQuadPol.InfiniteRoots()) {
117       inquadric=Standard_True;
118     }
119     else {
120       nbpts= LinQuadPol.NbSolutions();
121       
122       for(Standard_Integer i=1 ; i<=nbpts; i++) {
123         Standard_Real t= LinQuadPol.Value(i);
124         paramonc[i-1] = t;
125         pnts[i-1]=gp_Pnt( Lx0+Lx*t
126                          ,Ly0+Ly*t
127                          ,Lz0+Lz*t);
128       }
129     }
130   }
131 }
132
133 //=============================================================================
134 //==                                            C i r c l e   -   Q u a d r i c  
135 //==
136 CREATE(const gp_Circ& C,const IntAna_Quadric& Quad) {
137   Perform(C,Quad); 
138 }
139
140 PERFORM(const gp_Circ& C,const IntAna_Quadric& Quad) {
141   
142   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte;
143   
144   //----------------------------------------------------------------------
145   //-- Dans le repere liee a C.Position() : 
146   //-- xC = R * Cos[t]
147   //-- yC = R * Sin[t]
148   //-- zC = 0
149   //--
150   //-- On exprime la quadrique dans ce repere et on substitue 
151   //-- xC,yC et zC    a    x,y et z 
152   //-- 
153   //-- On Obtient un polynome en Cos[t] et Sin[t] de degre 2
154   //----------------------------------------------------------------------
155   done=inquadric=parallel=Standard_False;  
156
157   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte);
158   Quad.NewCoefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte,C.Position());
159   
160   Standard_Real R=C.Radius();
161   Standard_Real RR=R*R;
162   
163   Standard_Real P_CosCos = RR * Qxx;    //-- Cos Cos
164   Standard_Real P_SinSin = RR * Qyy;    //-- Sin Sin
165   Standard_Real P_Sin    = R  * Qy;     //-- 2 Sin
166   Standard_Real P_Cos    = R  * Qx;     //-- 2 Cos
167   Standard_Real P_CosSin = RR * Qxy;    //-- 2 Cos Sin
168   Standard_Real P_Cte    = QCte;        //-- 1
169   
170   math_TrigonometricFunctionRoots CircQuadPol( P_CosCos-P_SinSin
171                                               ,P_CosSin
172                                               ,P_Cos+P_Cos
173                                               ,P_Sin+P_Sin
174                                               ,P_Cte+P_SinSin
175                                               ,0.0,PIpPI);
176   
177   if(CircQuadPol.IsDone()) {
178     done=Standard_True;
179     if(CircQuadPol.InfiniteRoots()) {
180       inquadric=Standard_True;
181     }
182     else {
183       nbpts= CircQuadPol.NbSolutions();
184       
185       for(Standard_Integer i=1 ; i<=nbpts; i++) {
186         Standard_Real t= CircQuadPol.Value(i);
187         paramonc[i-1] = t;
188         pnts[i-1]     = ElCLib::CircleValue(t,C.Position(),R);
189       }
190     }
191   }
192 }
193
194
195 //=============================================================================
196 //==                                                  E l i p s - Q u a d r i c 
197 //==
198 CREATE(const gp_Elips& E,const IntAna_Quadric& Quad) {
199   Perform(E,Quad);
200 }
201
202 PERFORM(const gp_Elips& E,const IntAna_Quadric& Quad) {
203   
204   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte;
205   
206   done=inquadric=parallel=Standard_False;  
207   
208   //----------------------------------------------------------------------
209   //-- Dans le repere liee a E.Position() : 
210   //-- xE = R * Cos[t]
211   //-- yE = r * Sin[t]
212   //-- zE = 0
213   //--
214   //-- On exprime la quadrique dans ce repere et on substitue 
215   //-- xE,yE et zE    a    x,y et z 
216   //-- 
217   //-- On Obtient un polynome en Cos[t] et Sin[t] de degre 2
218   //----------------------------------------------------------------------
219   
220   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte);
221   Quad.NewCoefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte,E.Position());
222   
223   Standard_Real R=E.MajorRadius();  
224   Standard_Real r=E.MinorRadius();   
225
226   
227   Standard_Real P_CosCos = R*R * Qxx;    //-- Cos Cos
228   Standard_Real P_SinSin = r*r * Qyy;    //-- Sin Sin
229   Standard_Real P_Sin    = r   * Qy;     //-- 2 Sin
230   Standard_Real P_Cos    = R   * Qx;     //-- 2 Cos
231   Standard_Real P_CosSin = R*r * Qxy;    //-- 2 Cos Sin
232   Standard_Real P_Cte    = QCte;         //-- 1
233   
234   math_TrigonometricFunctionRoots ElipsQuadPol( P_CosCos-P_SinSin
235                                               ,P_CosSin
236                                               ,P_Cos+P_Cos
237                                               ,P_Sin+P_Sin
238                                               ,P_Cte+P_SinSin
239                                               ,0.0,PIpPI);
240   
241   if(ElipsQuadPol.IsDone()) {
242     done=Standard_True;
243     if(ElipsQuadPol.InfiniteRoots()) {
244       inquadric=Standard_True;
245     }
246     else {
247       nbpts= ElipsQuadPol.NbSolutions();
248       for(Standard_Integer i=1 ; i<=nbpts; i++) {
249         Standard_Real t= ElipsQuadPol.Value(i);
250         paramonc[i-1] = t;
251         pnts[i-1]     = ElCLib::EllipseValue(t,E.Position(),R,r);
252       }
253     }
254   }
255 }
256
257 //=============================================================================
258 //==                                                  P a r a b - Q u a d r i c 
259 //==
260 CREATE(const gp_Parab& P,const IntAna_Quadric& Quad) {
261   Perform(P,Quad);
262 }
263
264 PERFORM(const gp_Parab& P,const IntAna_Quadric& Quad) {
265   
266   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte;
267   
268   done=inquadric=parallel=Standard_False;
269
270   //----------------------------------------------------------------------
271   //-- Dans le repere liee a P.Position() : 
272   //-- xP = y*y / (2 p)
273   //-- yP = y
274   //-- zP = 0
275   //--
276   //-- On exprime la quadrique dans ce repere et on substitue 
277   //-- xP,yP et zP    a    x,y et z 
278   //-- 
279   //-- On Obtient un polynome en y de degre 4
280   //----------------------------------------------------------------------
281   
282   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte);
283   Quad.NewCoefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte,P.Position());
284   
285   Standard_Real f=P.Focal();
286   Standard_Real Un_Sur_2p = 0.25 / f;
287
288   Standard_Real A4 = Qxx * Un_Sur_2p * Un_Sur_2p;
289   Standard_Real A3 = (Qxy+Qxy) * Un_Sur_2p;
290   Standard_Real A2 = Qyy + (Qx+Qx) * Un_Sur_2p;
291   Standard_Real A1 = Qy+Qy;
292   Standard_Real A0 = QCte;
293   
294   math_DirectPolynomialRoots ParabQuadPol(A4,A3,A2,A1,A0);
295   
296   if( ParabQuadPol.IsDone()) {
297     done=Standard_True;
298     if(ParabQuadPol.InfiniteRoots()) {
299       inquadric=Standard_True;
300     }
301     else {
302       nbpts= ParabQuadPol.NbSolutions();
303       
304       for(Standard_Integer i=1 ; i<=nbpts; i++) {
305         Standard_Real t= ParabQuadPol.Value(i);
306         paramonc[i-1] = t;
307         pnts[i-1]     = ElCLib::ParabolaValue(t,P.Position(),f);
308       }
309     }
310   }  
311 }
312
313 //=============================================================================
314 //==                                                    H y p r - Q u a d r i c 
315 //==
316 CREATE(const gp_Hypr& H,const IntAna_Quadric& Quad) {
317   Perform(H,Quad);
318 }
319
320 PERFORM(const gp_Hypr& H,const IntAna_Quadric& Quad) {
321   
322   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte;
323
324   done=inquadric=parallel=Standard_False;  
325
326   //----------------------------------------------------------------------
327   //-- Dans le repere liee a P.Position() : 
328   //-- xH = R Ch[t]
329   //-- yH = r Sh[t]
330   //-- zH = 0
331   //--
332   //-- On exprime la quadrique dans ce repere et on substitue 
333   //-- xP,yP et zP    a    x,y et z 
334   //-- 
335   //-- On Obtient un polynome en Exp[t] de degre 4
336   //----------------------------------------------------------------------
337   
338   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte);
339   Quad.NewCoefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte,H.Position());
340   
341   Standard_Real R=H.MajorRadius();
342   Standard_Real r=H.MinorRadius();
343
344   Standard_Real RR=R*R;
345   Standard_Real rr=r*r;
346   Standard_Real Rr=R*r;
347
348   Standard_Real A4 = RR * Qxx + Rr* ( Qxy + Qxy ) + rr * Qyy;
349   Standard_Real A3 = 4.0* ( R*Qx + r*Qy );
350   Standard_Real A2 = 2.0* ( (QCte+QCte) + Qxx*RR - Qyy*rr );
351   Standard_Real A1 = 4.0* ( R*Qx - r*Qy);
352   Standard_Real A0 = Qxx*RR - Rr*(Qxy+Qxy) + Qyy*rr; 
353   
354   math_DirectPolynomialRoots HyperQuadPol(A4,A3,A2,A1,A0);
355   
356   if( HyperQuadPol.IsDone()) {
357     done=Standard_True;
358     if(HyperQuadPol.InfiniteRoots()) {
359       inquadric=Standard_True;
360     }
361     else {
362       nbpts= HyperQuadPol.NbSolutions();
363       Standard_Integer bonnessolutions = 0;
364       for(Standard_Integer i=1 ; i<=nbpts; i++) {
365         Standard_Real t= HyperQuadPol.Value(i);
366         if(t>=RealEpsilon()) {
367           Standard_Real Lnt = Log(t);
368           paramonc[bonnessolutions] = Lnt;
369           pnts[bonnessolutions]     = ElCLib::HyperbolaValue(Lnt,H.Position(),R,r);
370           bonnessolutions++;
371         }
372       }
373       nbpts=bonnessolutions;
374     }
375   } 
376 }
377 //=============================================================================
378
379
380
381
382 IntAna_IntConicQuad::IntAna_IntConicQuad (const gp_Lin& L, const gp_Pln& P,
383                                           const Standard_Real Tolang,
384                                           const Standard_Real Tol,
385                                           const Standard_Real Len) {
386   Perform(L,P,Tolang,Tol,Len);
387 }
388
389
390 IntAna_IntConicQuad::IntAna_IntConicQuad (const gp_Circ& C, const gp_Pln& P,
391                                           const Standard_Real Tolang,
392                                           const Standard_Real Tol) {
393   Perform(C,P,Tolang,Tol);
394 }
395
396
397 IntAna_IntConicQuad::IntAna_IntConicQuad (const gp_Elips& E, const gp_Pln& P,
398                                           const Standard_Real Tolang,
399                                           const Standard_Real Tol) {
400   Perform(E,P,Tolang,Tol);
401 }
402
403
404 IntAna_IntConicQuad::IntAna_IntConicQuad (const gp_Parab& Pb, const gp_Pln& P,
405                                           const Standard_Real Tolang){
406   Perform(Pb,P,Tolang);
407 }
408
409
410 IntAna_IntConicQuad::IntAna_IntConicQuad (const gp_Hypr& H, const gp_Pln& P,
411                                           const Standard_Real Tolang){
412   Perform(H,P,Tolang);
413 }
414
415
416 void IntAna_IntConicQuad::Perform (const gp_Lin& L, const gp_Pln& P,
417                                    const Standard_Real Tolang,
418                                    const Standard_Real Tol,
419                                    const Standard_Real Len) {
420   
421   // Tolang represente la tolerance angulaire a partir de laquelle on considere
422   // que l angle entre 2 vecteurs est nul. On raisonnera sur le cosinus de cet
423   // angle, (on a Cos(t) equivalent a t au voisinage de Pi/2).
424   
425   done=Standard_False;
426   
427   Standard_Real A,B,C,D;
428   Standard_Real Al,Bl,Cl;
429   Standard_Real Dis,Direc;
430   
431   P.Coefficients(A,B,C,D);
432   gp_Pnt Orig(L.Location());
433   L.Direction().Coord(Al,Bl,Cl);
434
435   Direc=A*Al+B*Bl+C*Cl;
436   Dis = A*Orig.X() + B*Orig.Y() + C*Orig.Z() + D;
437   //
438   parallel=Standard_False;
439   if (Abs(Direc) < Tolang) {
440     parallel=Standard_True;
441     if (Len!=0 && Direc!=0) {
442       //check the distance from bounding point of the line to the plane
443       gp_Pnt aP1, aP2;
444       //
445       aP1.SetCoord(Orig.X()-Dis*A, Orig.Y()-Dis*B, Orig.Z()-Dis*C);
446       aP2.SetCoord(aP1.X()+Len*Al, aP1.Y()+Len*Bl, aP1.Z()+Len*Cl);
447       if (P.Distance(aP2) > Tol) {
448         parallel=Standard_False;
449       } 
450     }
451   }
452   if (parallel) {
453     if (Abs(Dis) < Tolang) {
454       inquadric=Standard_True;
455     }
456     else {
457       inquadric=Standard_False;
458     }
459   }
460   else {
461     parallel=Standard_False;
462     inquadric=Standard_False;
463     nbpts = 1;
464     paramonc [0] = - Dis/Direc;
465     pnts[0].SetCoord(Orig.X()+paramonc[0]*Al,
466                      Orig.Y()+paramonc[0]*Bl,
467                      Orig.Z()+paramonc[0]*Cl);
468   }
469   done=Standard_True;
470 }
471
472
473 void IntAna_IntConicQuad::Perform (const gp_Circ& C, const gp_Pln& P,
474                                   const Standard_Real Tolang,
475                                   const Standard_Real Tol)
476 {
477   
478   done=Standard_False;
479   
480   gp_Pln Plconic(gp_Ax3(C.Position()));
481   IntAna_QuadQuadGeo IntP(Plconic,P,Tolang,Tol);
482   if (!IntP.IsDone()) {return;}
483   if (IntP.TypeInter() == IntAna_Empty) {
484     parallel=Standard_True;
485     Standard_Real distmax = P.Distance(C.Location()) + C.Radius()*Tolang;
486     if (distmax < Tol) {
487       inquadric = Standard_True;
488     }
489     else {
490       inquadric = Standard_False;
491     }
492     done=Standard_True;
493   }
494   else     if(IntP.TypeInter() == IntAna_Same) { 
495     inquadric = Standard_True;
496     done = Standard_True;
497   }
498   else {
499     inquadric=Standard_False;
500     parallel=Standard_False;
501     gp_Lin Ligsol(IntP.Line(1));
502     
503     gp_Vec V0(Plconic.Location(),Ligsol.Location());
504     gp_Vec Axex(Plconic.Position().XDirection());
505     gp_Vec Axey(Plconic.Position().YDirection());
506     
507     gp_Pnt2d Orig(Axex.Dot(V0),Axey.Dot(V0));
508     gp_Vec2d Dire(Axex.Dot(Ligsol.Direction()),
509                   Axey.Dot(Ligsol.Direction()));
510     
511     gp_Lin2d Ligs(Orig,Dire);
512     gp_Pnt2d Pnt2dBid(0.0,0.0);
513     gp_Dir2d Dir2dBid(1.0,0.0);
514     gp_Ax2d Ax2dBid(Pnt2dBid,Dir2dBid);
515     gp_Circ2d Cir(Ax2dBid,C.Radius());
516     
517     IntAna2d_AnaIntersection Int2d(Ligs,Cir);
518     
519     if (!Int2d.IsDone()) {return;}
520     
521     nbpts=Int2d.NbPoints();
522     for (Standard_Integer i=1; i<=nbpts; i++) {
523       
524       gp_Pnt2d resul(Int2d.Point(i).Value());
525       Standard_Real X= resul.X();
526       Standard_Real Y= resul.Y();
527       pnts[i-1].SetCoord(Plconic.Location().X() + X*Axex.X() + Y*Axey.X(),
528                          Plconic.Location().Y() + X*Axex.Y() + Y*Axey.Y(),
529                          Plconic.Location().Z() + X*Axex.Z() + Y*Axey.Z());
530       paramonc[i-1]=Int2d.Point(i).ParamOnSecond();
531     }
532     done=Standard_True;
533   }
534 }
535
536
537 void IntAna_IntConicQuad::Perform (const gp_Elips& E, const gp_Pln& Pln,
538                                    const Standard_Real,
539                                    const Standard_Real){
540   Perform(E,Pln);
541 }
542
543
544 void IntAna_IntConicQuad::Perform (const gp_Parab& P, const gp_Pln& Pln,
545                                    const Standard_Real){
546   Perform(P,Pln);
547 }
548
549
550 void IntAna_IntConicQuad::Perform (const gp_Hypr& H, const gp_Pln& Pln,
551                                    const Standard_Real){
552   Perform(H,Pln);
553 }
554
555