1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15
16 #include <ElCLib.hxx>
17 #include <ElSLib.hxx>
18 #include <gp.hxx>
19 #include <gp_Cone.hxx>
20 #include <gp_Cylinder.hxx>
21 #include <gp_Pln.hxx>
22 #include <gp_Pnt.hxx>
23 #include <gp_Sphere.hxx>
24 #include <gp_Torus.hxx>
25 #include <gp_Vec.hxx>
26 #include <IntSurf_Quadric.hxx>
27 #include <StdFail_NotDone.hxx>
28
29 // ============================================================
IntSurf_Quadric()30 IntSurf_Quadric::IntSurf_Quadric ():typ(GeomAbs_OtherSurface),
31 prm1(0.), prm2(0.), prm3(0.), prm4(0.), ax3direc(Standard_False)
32 {
33 }
34 // ============================================================
IntSurf_Quadric(const gp_Pln & P)35 IntSurf_Quadric::IntSurf_Quadric (const gp_Pln& P):
36 ax3(P.Position()),typ(GeomAbs_Plane)
37 {
38 ax3direc = ax3.Direct();
39 P.Coefficients(prm1,prm2,prm3,prm4);
40 }
41 // ============================================================
IntSurf_Quadric(const gp_Cylinder & C)42 IntSurf_Quadric::IntSurf_Quadric (const gp_Cylinder& C):
43
44 ax3(C.Position()),lin(ax3.Axis()),typ(GeomAbs_Cylinder)
45 {
46 prm2=prm3=prm4=0.0;
47 ax3direc=ax3.Direct();
48 prm1=C.Radius();
49 }
50 // ============================================================
IntSurf_Quadric(const gp_Sphere & S)51 IntSurf_Quadric::IntSurf_Quadric (const gp_Sphere& S):
52
53 ax3(S.Position()),lin(ax3.Axis()),typ(GeomAbs_Sphere)
54 {
55 prm2=prm3=prm4=0.0;
56 ax3direc = ax3.Direct();
57 prm1=S.Radius();
58 }
59 // ============================================================
IntSurf_Quadric(const gp_Cone & C)60 IntSurf_Quadric::IntSurf_Quadric (const gp_Cone& C):
61
62 ax3(C.Position()),typ(GeomAbs_Cone)
63 {
64 ax3direc = ax3.Direct();
65 lin.SetPosition(ax3.Axis());
66 prm1 = C.RefRadius();
67 prm2 = C.SemiAngle();
68 prm3 = Cos(prm2);
69 prm4 = 0.0;
70 }
71 // ============================================================
IntSurf_Quadric(const gp_Torus & T)72 IntSurf_Quadric::IntSurf_Quadric (const gp_Torus& T):
73
74 ax3(T.Position()),typ(GeomAbs_Torus)
75 {
76 ax3direc = ax3.Direct();
77 lin.SetPosition(ax3.Axis());
78 prm1 = T.MajorRadius();
79 prm2 = T.MinorRadius();
80 prm3 = 0.0;
81 prm4 = 0.0;
82 }
83 // ============================================================
SetValue(const gp_Pln & P)84 void IntSurf_Quadric::SetValue (const gp_Pln& P)
85 {
86 typ = GeomAbs_Plane;
87 ax3 = P.Position();
88 ax3direc = ax3.Direct();
89 P.Coefficients(prm1,prm2,prm3,prm4);
90 }
91 // ============================================================
SetValue(const gp_Cylinder & C)92 void IntSurf_Quadric::SetValue (const gp_Cylinder& C)
93 {
94 typ = GeomAbs_Cylinder;
95 ax3 = C.Position();
96 ax3direc = ax3.Direct();
97 lin.SetPosition(ax3.Axis());
98 prm1 = C.Radius();
99 prm2=prm3=prm4=0.0;
100 }
101 // ============================================================
SetValue(const gp_Sphere & S)102 void IntSurf_Quadric::SetValue (const gp_Sphere& S)
103 {
104 typ = GeomAbs_Sphere;
105 ax3 = S.Position();
106 ax3direc = ax3.Direct();
107 lin.SetPosition(ax3.Axis());
108 prm1 = S.Radius();
109 prm2=prm3=prm4=0.0;
110 }
111 // ============================================================
SetValue(const gp_Cone & C)112 void IntSurf_Quadric::SetValue (const gp_Cone& C)
113 {
114 typ = GeomAbs_Cone;
115 ax3 = C.Position();
116 ax3direc = ax3.Direct();
117 lin.SetPosition(ax3.Axis());
118 prm1 = C.RefRadius();
119 prm2 = C.SemiAngle();
120 prm3 = Cos(prm2);
121 prm4 = 0.0;
122 }
123 // ============================================================
SetValue(const gp_Torus & T)124 void IntSurf_Quadric::SetValue (const gp_Torus& T)
125 {
126 typ = GeomAbs_Torus;
127 ax3 = T.Position();
128 ax3direc = ax3.Direct();
129 lin.SetPosition(ax3.Axis());
130 prm1 = T.MajorRadius();
131 prm2 = T.MinorRadius();
132 prm3 = 0.0;
133 prm4 = 0.0;
134 }
135 // ============================================================
Distance(const gp_Pnt & P) const136 Standard_Real IntSurf_Quadric::Distance (const gp_Pnt& P) const {
137 switch (typ) {
138 case GeomAbs_Plane: // plan
139 return prm1*P.X() + prm2*P.Y() + prm3*P.Z() + prm4;
140 case GeomAbs_Cylinder: // cylindre
141 return (lin.Distance(P) - prm1);
142 case GeomAbs_Sphere: // sphere
143 return (lin.Location().Distance(P) - prm1);
144 case GeomAbs_Cone: // cone
145 {
146 Standard_Real dist = lin.Distance(P);
147 Standard_Real U,V;
148 ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
149 gp_Pnt Pp = ElSLib::ConeValue(U,V,ax3,prm1,prm2);
150 Standard_Real distp = lin.Distance(Pp);
151 dist = (dist-distp)/prm3;
152 return(dist);
153 }
154 case GeomAbs_Torus: // torus
155 {
156 gp_Pnt O, Pp, PT;
157 //
158 O = ax3.Location();
159 gp_Vec OZ (ax3.Direction());
160 Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
161 //
162 gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
163 ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
164 PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
165 //
166 Standard_Real dist = P.Distance(PT) - prm2;
167 return dist;
168 }
169 default:
170 {
171 }
172 break;
173 }
174 return(0.0);
175 }
176 // ============================================================
Gradient(const gp_Pnt & P) const177 gp_Vec IntSurf_Quadric::Gradient (const gp_Pnt& P) const {
178 gp_Vec grad;
179 switch (typ) {
180 case GeomAbs_Plane: // plan
181 grad.SetCoord(prm1,prm2,prm3);
182 break;
183 case GeomAbs_Cylinder: // cylindre
184 {
185 gp_XYZ PP(lin.Location().XYZ());
186 PP.Add(ElCLib::Parameter(lin,P)*lin.Direction().XYZ());
187 grad.SetXYZ(P.XYZ()-PP);
188 Standard_Real N = grad.Magnitude();
189 if(N>1e-14) { grad.Divide(N); }
190 else { grad.SetCoord(0.0,0.0,0.0); }
191 }
192 break;
193 case GeomAbs_Sphere: // sphere
194 {
195 gp_XYZ PP(P.XYZ());
196 grad.SetXYZ((PP-lin.Location().XYZ()));
197 Standard_Real N = grad.Magnitude();
198 if(N>1e-14) { grad.Divide(N); }
199 else { grad.SetCoord(0.0,0.0,0.0); }
200 }
201 break;
202 case GeomAbs_Cone: // cone
203 {
204 Standard_Real U,V;
205 ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
206 gp_Pnt Pp = ElSLib::ConeValue(U,V,ax3,prm1,prm2);
207 gp_Vec D1u,D1v;
208 ElSLib::ConeD1(U,V,ax3,prm1,prm2,Pp,D1u,D1v);
209 grad=D1u.Crossed(D1v);
210 if(ax3direc==Standard_False) {
211 grad.Reverse();
212 }
213 grad.Normalize();
214 }
215 break;
216 case GeomAbs_Torus: // torus
217 {
218 gp_Pnt O, Pp, PT;
219 //
220 O = ax3.Location();
221 gp_Vec OZ (ax3.Direction());
222 Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
223 //
224 gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
225 ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
226 PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
227 //
228 grad.SetXYZ(P.XYZ() - PT.XYZ());
229 Standard_Real N = grad.Magnitude();
230 if(N>1e-14) { grad.Divide(N); }
231 else { grad.SetCoord(0., 0., 0.); }
232 }
233 break;
234 default:
235 {}
236 break;
237 }
238 return grad;
239 }
240 // ============================================================
ValAndGrad(const gp_Pnt & P,Standard_Real & Dist,gp_Vec & Grad) const241 void IntSurf_Quadric::ValAndGrad (const gp_Pnt& P,
242 Standard_Real& Dist,
243 gp_Vec& Grad) const
244 {
245
246 switch (typ) {
247 case GeomAbs_Plane:
248 {
249 Dist = prm1*P.X() + prm2*P.Y() + prm3*P.Z() + prm4;
250 Grad.SetCoord(prm1,prm2,prm3);
251 }
252 break;
253 case GeomAbs_Cylinder:
254 {
255 Dist = lin.Distance(P) - prm1;
256 gp_XYZ PP(lin.Location().XYZ());
257 PP.Add(ElCLib::Parameter(lin,P)*lin.Direction().XYZ());
258 Grad.SetXYZ((P.XYZ()-PP));
259 Standard_Real N = Grad.Magnitude();
260 if(N>1e-14) { Grad.Divide(N); }
261 else { Grad.SetCoord(0.0,0.0,0.0); }
262 }
263 break;
264 case GeomAbs_Sphere:
265 {
266 Dist = lin.Location().Distance(P) - prm1;
267 gp_XYZ PP(P.XYZ());
268 Grad.SetXYZ((PP-lin.Location().XYZ()));
269 Standard_Real N = Grad.Magnitude();
270 if(N>1e-14) { Grad.Divide(N); }
271 else { Grad.SetCoord(0.0,0.0,0.0); }
272 }
273 break;
274 case GeomAbs_Cone:
275 {
276 Standard_Real dist = lin.Distance(P);
277 Standard_Real U,V;
278 gp_Vec D1u,D1v;
279 gp_Pnt Pp;
280 ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
281 ElSLib::ConeD1(U,V,ax3,prm1,prm2,Pp,D1u,D1v);
282 Standard_Real distp = lin.Distance(Pp);
283 dist = (dist-distp)/prm3;
284 Dist = dist;
285 Grad=D1u.Crossed(D1v);
286 if(ax3direc==Standard_False) {
287 Grad.Reverse();
288 }
289 //-- lbr le 7 mars 96
290 //-- Si le gardient est nul, on est sur l axe
291 //-- et dans ce cas dist vaut 0
292 //-- On peut donc renvoyer une valeur quelconque.
293 if( Grad.X() > 1e-13
294 || Grad.Y() > 1e-13
295 || Grad.Z() > 1e-13) {
296 Grad.Normalize();
297 }
298 }
299 break;
300 case GeomAbs_Torus:
301 {
302 gp_Pnt O, Pp, PT;
303 //
304 O = ax3.Location();
305 gp_Vec OZ (ax3.Direction());
306 Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
307 //
308 gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
309 ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
310 PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
311 //
312 Dist = P.Distance(PT) - prm2;
313 //
314 Grad.SetXYZ(P.XYZ()-PT.XYZ());
315 Standard_Real N = Grad.Magnitude();
316 if(N>1e-14) { Grad.Divide(N); }
317 else { Grad.SetCoord(0., 0., 0.); }
318 }
319 break;
320 default:
321 {}
322 break;
323 }
324 }
325 // ============================================================
Value(const Standard_Real U,const Standard_Real V) const326 gp_Pnt IntSurf_Quadric::Value(const Standard_Real U,
327 const Standard_Real V) const
328 {
329 switch (typ) {
330
331 case GeomAbs_Plane:
332 return ElSLib::PlaneValue(U,V,ax3);
333 case GeomAbs_Cylinder:
334 return ElSLib::CylinderValue(U,V,ax3,prm1);
335 case GeomAbs_Sphere:
336 return ElSLib::SphereValue(U,V,ax3,prm1);
337 case GeomAbs_Cone:
338 return ElSLib::ConeValue(U,V,ax3,prm1,prm2);
339 case GeomAbs_Torus:
340 return ElSLib::TorusValue(U,V,ax3,prm1,prm2);
341 default:
342 {
343 gp_Pnt p(0,0,0);
344 return(p);
345 }
346 //break;
347 }
348 // pop : pour NT
349 // return gp_Pnt(0,0,0);
350 }
351 // ============================================================
D1(const Standard_Real U,const Standard_Real V,gp_Pnt & P,gp_Vec & D1U,gp_Vec & D1V) const352 void IntSurf_Quadric::D1(const Standard_Real U,
353 const Standard_Real V,
354 gp_Pnt& P,
355 gp_Vec& D1U,
356 gp_Vec& D1V) const
357 {
358 switch (typ) {
359 case GeomAbs_Plane:
360 ElSLib::PlaneD1(U,V,ax3,P,D1U,D1V);
361 break;
362 case GeomAbs_Cylinder:
363 ElSLib::CylinderD1(U,V,ax3,prm1,P,D1U,D1V);
364 break;
365 case GeomAbs_Sphere:
366 ElSLib::SphereD1(U,V,ax3,prm1,P,D1U,D1V);
367 break;
368 case GeomAbs_Cone:
369 ElSLib::ConeD1(U,V,ax3,prm1,prm2,P,D1U,D1V);
370 break;
371 case GeomAbs_Torus:
372 ElSLib::TorusD1(U,V,ax3,prm1,prm2,P,D1U,D1V);
373 break;
374 default:
375 {
376 }
377 break;
378 }
379 }
380 // ============================================================
DN(const Standard_Real U,const Standard_Real V,const Standard_Integer Nu,const Standard_Integer Nv) const381 gp_Vec IntSurf_Quadric::DN(const Standard_Real U,
382 const Standard_Real V,
383 const Standard_Integer Nu,
384 const Standard_Integer Nv) const
385 {
386 switch (typ) {
387 case GeomAbs_Plane:
388 return ElSLib::PlaneDN(U,V,ax3,Nu,Nv);
389 case GeomAbs_Cylinder:
390 return ElSLib::CylinderDN(U,V,ax3,prm1,Nu,Nv);
391 case GeomAbs_Sphere:
392 return ElSLib::SphereDN(U,V,ax3,prm1,Nu,Nv);
393 case GeomAbs_Cone:
394 return ElSLib::ConeDN(U,V,ax3,prm1,prm2,Nu,Nv);
395 case GeomAbs_Torus:
396 return ElSLib::TorusDN(U,V,ax3,prm1,prm2,Nu,Nv);
397 default:
398 {
399 gp_Vec v(0,0,0);
400 return(v);
401 }
402 //break;
403 }
404 // pop : pour NT
405 // return gp_Vec(0,0,0);
406 }
407 // ============================================================
Normale(const Standard_Real U,const Standard_Real V) const408 gp_Vec IntSurf_Quadric::Normale(const Standard_Real U,
409 const Standard_Real V) const
410 {
411 switch (typ) {
412 case GeomAbs_Plane:
413 if(ax3direc)
414 return ax3.Direction();
415 else
416 return ax3.Direction().Reversed();
417 case GeomAbs_Cylinder:
418 return Normale(Value(U,V));
419 case GeomAbs_Sphere:
420 return Normale(Value(U,V));
421 case GeomAbs_Cone:
422 {
423 gp_Pnt P;
424 gp_Vec D1u,D1v;
425 ElSLib::ConeD1(U,V,ax3,prm1,prm2,P,D1u,D1v);
426 if(D1u.Magnitude()<0.0000001) {
427 gp_Vec Vn(0.0,0.0,0.0);
428 return(Vn);
429 }
430 return(D1u.Crossed(D1v));
431 }
432 case GeomAbs_Torus:
433 return Normale(Value(U,V));
434 default:
435 {
436 gp_Vec v(0,0,0);
437 return(v);
438 }
439 // break;
440 }
441 // pop : pour NT
442 // return gp_Vec(0,0,0);
443 }
444 // ============================================================
Normale(const gp_Pnt & P) const445 gp_Vec IntSurf_Quadric::Normale (const gp_Pnt& P) const
446 {
447 switch (typ) {
448 case GeomAbs_Plane:
449 if(ax3direc)
450 return ax3.Direction();
451 else
452 return ax3.Direction().Reversed();
453 case GeomAbs_Cylinder:
454 {
455 if(ax3direc) {
456 return lin.Normal(P).Direction();
457 }
458 else {
459 gp_Dir D(lin.Normal(P).Direction());
460 D.Reverse();
461 return(D);
462 }
463 }
464 case GeomAbs_Sphere:
465 {
466 if(ax3direc) {
467 gp_Vec ax3P(ax3.Location(),P);
468 return gp_Dir(ax3P);
469 }
470 else {
471 gp_Vec Pax3(P,ax3.Location());
472 return gp_Dir(Pax3);
473 }
474 }
475 case GeomAbs_Cone:
476 {
477 Standard_Real U,V;
478 ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
479 return Normale(U,V);
480 }
481 case GeomAbs_Torus:
482 {
483 gp_Pnt O, Pp, PT;
484 //
485 O = ax3.Location();
486 gp_Vec OZ (ax3.Direction());
487 Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
488 //
489 gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
490 ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
491 PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
492 if (PT.SquareDistance(P) < 1e-14) {
493 return gp_Dir(OZ);
494 }
495 gp_Dir aD(ax3direc ? gp_Vec(PT, P) : gp_Vec(P, PT));
496 return aD;
497 }
498 default:
499 {
500 gp_Vec v(0,0,0);
501 return(v);
502 } // break;
503 }
504 }
505 // ============================================================
Parameters(const gp_Pnt & P,Standard_Real & U,Standard_Real & V) const506 void IntSurf_Quadric::Parameters (const gp_Pnt& P,
507 Standard_Real& U,
508 Standard_Real& V) const
509 {
510 switch (typ) {
511 case GeomAbs_Plane:
512 ElSLib::PlaneParameters(ax3,P,U,V);
513 break;
514 case GeomAbs_Cylinder:
515 ElSLib::CylinderParameters(ax3,prm1,P,U,V);
516 break;
517 case GeomAbs_Sphere:
518 ElSLib::SphereParameters(ax3,prm1,P,U,V);
519 break;
520 case GeomAbs_Cone:
521 ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
522 break;
523 case GeomAbs_Torus:
524 ElSLib::TorusParameters(ax3,prm1,prm2,P,U,V);
525 break;
526 default:
527 break;
528 }
529 }
530 // ============================================================
531
532
533
534
535
536