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