projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0024081: Memory corruption when projecting point on surface of revolution
[occt.git] / src / Extrema / Extrema_ExtPRevS.cxx
CommitLineData
b311480e 1// Created on: 1999-09-21
2// Created by: Edward AGAPOV
3// Copyright (c) 1999-1999 Matra Datavision
4// Copyright (c) 1999-2012 OPEN CASCADE SAS
5//
6// The content of this file is subject to the Open CASCADE Technology Public
7// License Version 6.5 (the "License"). You may not use the content of this file
8// except in compliance with the License. Please obtain a copy of the License
9// at http://www.opencascade.org and read it completely before using this file.
10//
11// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13//
14// The Original Code and all software distributed under the License is
15// distributed on an "AS IS" basis, without warranty of any kind, and the
16// Initial Developer hereby disclaims all such warranties, including without
17// limitation, any warranties of merchantability, fitness for a particular
18// purpose or non-infringement. Please see the License for the specific terms
19// and conditions governing the rights and limitations under the License.
20
7fd59977 21
22
23#include <Extrema_ExtPRevS.ixx>
24#include <Adaptor3d_HCurve.hxx>
25#include <Extrema_ExtPElC.hxx>
26#include <Extrema_GenExtPS.hxx>
27#include <Extrema_POnCurv.hxx>
28#include <Extrema_POnSurf.hxx>
29#include <gp_Ax1.hxx>
30#include <gp_Ax2.hxx>
31#include <gp_Lin.hxx>
32#include <gp_Pln.hxx>
33#include <gp_Pnt.hxx>
34#include <gp_Trsf.hxx>
35#include <gp_Vec.hxx>
36#include <Precision.hxx>
37#include <ElCLib.hxx>
38
39static gp_Ax2 GetPosition (const Adaptor3d_SurfaceOfRevolution& S)//const Handle(Adaptor_HCurve)& C)
40{
41 Handle(Adaptor3d_HCurve) C = S.BasisCurve();
42
43 switch (C->GetType()) {
44
45 case GeomAbs_Line: {
46 gp_Lin L = C->Line();
47 gp_Dir N = S.AxeOfRevolution().Direction();
48 if (N.IsParallel(L.Direction(), Precision::Angular())) {
49 gp_Vec OO (L.Location(), S.AxeOfRevolution().Location());
50 if (OO.Magnitude() <= gp::Resolution()) {
51 OO = gp_Vec(L.Location(), ElCLib::Value(100,L));
52 if (N.IsParallel(OO, Precision::Angular()))
53 return gp_Ax2(); // Line and axe of revolution coinside
54 }
55 N ^= OO;
56 }
57 else {
58 N ^= L.Direction();
59 }
60 return gp_Ax2 (L.Location(), N, L.Direction());
61 }
62 case GeomAbs_Circle:
63 return C->Circle().Position();
64 case GeomAbs_Ellipse:
65 return C->Ellipse().Position();
66 case GeomAbs_Hyperbola:
67 return C->Hyperbola().Position();
68 case GeomAbs_Parabola:
69 return C->Parabola().Position();
70 default:
71 return gp_Ax2();
72 }
73}
74
75//=======================================================================
76//function : HasSingularity
77//purpose :
78//=======================================================================
79
80static Standard_Boolean HasSingularity(const Adaptor3d_SurfaceOfRevolution& S)
81{
82
83 const Handle(Adaptor3d_HCurve) C = S.BasisCurve();
84 gp_Dir N = S.AxeOfRevolution().Direction();
85 gp_Pnt P = S.AxeOfRevolution().Location();
86
87 gp_Lin L(P, N);
88
89 P = C->Value(C->FirstParameter());
90
08cd2f6b 91 if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;
7fd59977 92
93 P = C->Value(C->LastParameter());
94
08cd2f6b 95 if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;
7fd59977 96
97 return Standard_False;
98}
99
100//=============================================================================
101
102static void PerformExtPElC (Extrema_ExtPElC& E,
103 const gp_Pnt& P,
104 const Handle(Adaptor3d_HCurve)& C,
105 const Standard_Real Tol)
106{
107 switch (C->GetType()) {
108 case GeomAbs_Hyperbola:
109 E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
110 return;
111 case GeomAbs_Line:
112 E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
113 return;
114 case GeomAbs_Circle:
c6541a0c 115 E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
7fd59977 116 return;
117 case GeomAbs_Ellipse:
c6541a0c 118 E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
7fd59977 119 return;
120 case GeomAbs_Parabola:
121 E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
122 return;
7fd59977 123 default:
124 return ;
7fd59977 125 }
126}
127
128//=======================================================================
129//function : IsCaseAnalyticallyComputable
130//purpose :
131//=======================================================================
132
133static Standard_Boolean IsCaseAnalyticallyComputable
134 (const GeomAbs_CurveType& theType,
135 const gp_Ax2& theCurvePos,
136 const gp_Ax1& AxeOfRevolution)
137{
138 // check type
139 switch (theType) {
140 case GeomAbs_Line:
141 case GeomAbs_Circle:
142 case GeomAbs_Ellipse:
143 case GeomAbs_Hyperbola:
144 case GeomAbs_Parabola:
145 break;
146 default:
147 return Standard_False;
148 }
149// the axe of revolution must be in the plane of the curve.
150 gp_Pln pl(theCurvePos.Location(), theCurvePos.Direction());
151 gp_Pnt p1 = AxeOfRevolution.Location();
152 Standard_Real dist = 100., dist2 = dist * dist;
153 Standard_Real aThreshold = Precision::Angular() * Precision::Angular() * dist2;
154 gp_Pnt p2 = AxeOfRevolution.Location().XYZ() + dist * AxeOfRevolution.Direction().XYZ();
155
156 if((pl.SquareDistance(p1) < aThreshold) &&
157 (pl.SquareDistance(p2) < aThreshold))
158 return Standard_True;
159 return Standard_False;
160 // gp_Vec V (AxeOfRevolution.Location(),theCurvePos.Location());
161 // if (Abs( V * theCurvePos.Direction()) <= gp::Resolution())
162 // return Standard_True;
163 // else
164 // return Standard_False;
165}
166
167//=======================================================================
168//function : IsOriginalPnt
169//purpose :
170//=======================================================================
171
172static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
173 const Extrema_POnSurf* Points,
174 const Standard_Integer NbPoints)
175{
176 for (Standard_Integer i=1; i<=NbPoints; i++) {
6f94f1df 177 if (Points[i-1].Value().IsEqual(P, Precision::Confusion())) {
7fd59977 178 return Standard_False;
179 }
180 }
181 return Standard_True;
182}
183//=======================================================================
184//function : IsExtremum
185//purpose :
186//=======================================================================
187
188static Standard_Boolean IsExtremum (const Standard_Real U, const Standard_Real V,
189 const gp_Pnt& P, const Adaptor3d_SurfacePtr& S,
190 gp_Pnt& E, Standard_Real& Dist2,
191 const Standard_Boolean IsVSup,
192 const Standard_Boolean IsMin)
193{
194 E = S->Value(U,V);
195 Dist2 = P.SquareDistance(E);
196 if (IsMin)
197 return (Dist2 < P.SquareDistance(S->Value(U+1,V)) &&
198 Dist2 < P.SquareDistance(S->Value(U-1,V)) &&
199 Dist2 < P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
200 else
201 return (Dist2 > P.SquareDistance(S->Value(U+1,V)) &&
202 Dist2 > P.SquareDistance(S->Value(U-1,V)) &&
203 Dist2 > P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
204}
205//=======================================================================
206//function : Extrema_ExtPRevS
207//purpose :
208//=======================================================================
209
210Extrema_ExtPRevS::Extrema_ExtPRevS()
211{
212 myS=NULL;
213 myDone = Standard_False;
214}
215//=======================================================================
216//function : Extrema_ExtPRevS
217//purpose :
218//=======================================================================
219
220Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
221 const Adaptor3d_SurfaceOfRevolution& S,
222 const Standard_Real Umin,
223 const Standard_Real Usup,
224 const Standard_Real Vmin,
225 const Standard_Real Vsup,
226 const Standard_Real TolU,
227 const Standard_Real TolV)
228{
229 myS=NULL;
230 Initialize (S,
231 Umin, Usup, Vmin, Vsup,
232 TolU, TolV);
233 Perform(P);
234}
235//=======================================================================
236//function : Extrema_ExtPRevS
237//purpose :
238//=======================================================================
239
240Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
241 const Adaptor3d_SurfaceOfRevolution& S,
242 const Standard_Real TolU,
243 const Standard_Real TolV)
244{
245 myS=NULL;
246 Initialize (S,
247 S.FirstUParameter(), S.LastUParameter(),
248 S.FirstVParameter(), S.LastVParameter(),
249 TolU, TolV);
250 Perform(P);
251}
252//=======================================================================
253//function : Initialize
254//purpose :
255//=======================================================================
256
257void Extrema_ExtPRevS::Initialize(const Adaptor3d_SurfaceOfRevolution& S,
258 const Standard_Real Umin,
259 const Standard_Real Usup,
260 const Standard_Real Vmin,
261 const Standard_Real Vsup,
262 const Standard_Real TolU,
263 const Standard_Real TolV)
264{
265 myvinf=Vmin;
266 myvsup=Vsup;
267 mytolv=TolV;
268
269 Handle(Adaptor3d_HCurve) anACurve = S.BasisCurve();
270 if (!myS || myS != (Adaptor3d_SurfacePtr)&S) {
271 myS = Adaptor3d_SurfacePtr(&S);
272 myPosition = GetPosition(S);
273 myIsAnalyticallyComputable =
274 IsCaseAnalyticallyComputable (anACurve->GetType(),myPosition,S.AxeOfRevolution());
275 }
276 if (!myIsAnalyticallyComputable) {
277
278 Standard_Integer nbu = 32, nbv = 32;
279
280 if(HasSingularity(S)) {nbv = 100;}
281
282 myExtPS.Initialize(S, nbu, nbv,
283 Umin, Usup, Vmin, Vsup,
284 TolU, TolV);
285 }
286}
287//=======================================================================
288//function : Perform
289//purpose :
290//=======================================================================
291
292void Extrema_ExtPRevS::Perform(const gp_Pnt& P)
293{
294 myDone = Standard_False;
295 myNbExt = 0;
296
297 if (!myIsAnalyticallyComputable) {
298
299 myExtPS.Perform(P);
300 myDone = myExtPS.IsDone();
301 myNbExt = myExtPS.NbExt();
302 return;
303 }
304
305 Handle(Adaptor3d_HCurve) anACurve = myS->BasisCurve();
306
307 gp_Ax1 Ax = myS->AxeOfRevolution();
308 gp_Vec Dir = Ax.Direction(), Z = myPosition.Direction();
309 gp_Pnt O = Ax.Location();
310
311 Standard_Real OPdir = gp_Vec(O,P).Dot(Dir);
312 gp_Pnt Pp = P.Translated(Dir.Multiplied(-OPdir));
313 if (O.IsEqual(Pp,Precision::Confusion())) // P is on the AxeOfRevolution
314 return;
315
316 Standard_Real U,V;
317 gp_Pnt P1, Ppp;
318 Standard_Real OPpz = gp_Vec(O,Pp).Dot(Z);
319 if (Abs(OPpz) <= gp::Resolution()) {
320 Ppp = Pp;
321 U = 0;
322 }
323 else {
324 Ppp = Pp.Translated(Z.Multiplied(-OPpz));
325 if (O.IsEqual(Ppp,Precision::Confusion()))
c6541a0c 326 U = M_PI/2;
7fd59977 327 else {
328 U = gp_Vec(O,Ppp).AngleWithRef(gp_Vec(O,Pp),Dir);
329 }
330 }
331
c6541a0c
D		 332 gp_Vec OPpp (O,Ppp), OPq (O, myS->Value(M_PI/2,0));
333 if (U != M_PI/2) {
7fd59977 334 if (Abs(OPq.Magnitude()) <= gp::Resolution())
c6541a0c 335 OPq = gp_Vec(O, myS->Value(M_PI/2,anACurve->LastParameter()/10));
7fd59977 336 if (OPpp.AngleWithRef(OPq,Dir) < 0)
c6541a0c 337 U += M_PI;
7fd59977 338 }
339
340 gp_Trsf T;
341 T.SetRotation(Ax, -U);
342 P1 = P.Transformed(T);
343
344 gp_Pnt E;
345 Standard_Real Dist2;
346 Standard_Integer i;
347
348 Extrema_ExtPElC anExt;
349 PerformExtPElC(anExt, P1, anACurve, mytolv);
350
351 if (anExt.IsDone()) {
352 myDone = Standard_True;
353 for (i=1; i<=anExt.NbExt(); i++) {
354 Extrema_POnCurv POC=anExt.Point(i);
355 V = POC.Parameter();
356 if (V > myvsup) {
357 // if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
358 // Standard_True, anExt.IsMin(i))) continue;
359 Standard_Real newV = myvsup;
360
361 if((anACurve->GetType() == GeomAbs_Circle) ||
362 (anACurve->GetType() == GeomAbs_Ellipse)) {
c6541a0c 363 newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
7fd59977 364
365 if (newV > myvsup) {
a20400b3 366 newV -= 2. * M_PI;
367
368 if (newV + mytolv < myvinf) {
369 newV = myvsup;
370 } else if (newV < myvinf) {
371 newV = myvinf;
372 }
7fd59977 373 }
374 }
375 V = newV;
376
377 if ( !IsExtremum (U, V, P, myS, E, Dist2,
378 Standard_True, anExt.IsMin(i))) {
379 continue;
380 }
381 } else if (V < myvinf) {
382 // if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
383 // Standard_False, anExt.IsMin(i))) continue;
384
385 Standard_Real newV = myvinf;
386
387 if((anACurve->GetType() == GeomAbs_Circle) ||
388 (anACurve->GetType() == GeomAbs_Ellipse)) {
c6541a0c 389 newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
7fd59977 390
a20400b3 391 if(newV < myvinf) {
392 newV += 2. * M_PI;
393
394 if (newV - mytolv > myvsup) {
395 newV = myvinf;
396 } else if (newV > myvsup) {
397 newV = myvsup;
398 }
399 }
7fd59977 400 }
401 V = newV;
402
403 if ( !IsExtremum (U, V, P, myS, E, Dist2,
404 Standard_False, anExt.IsMin(i))) continue;
405 } else {
406 E = myS->Value(U,V);
407 Dist2 = P.SquareDistance(E);
408 }
409 if (IsOriginalPnt(E, myPoint, myNbExt)) {
6f94f1df 410 myPoint[myNbExt] = Extrema_POnSurf(U,V,E);
7fd59977 411 mySqDist[myNbExt] = Dist2;
6f94f1df 412 myNbExt++;
7fd59977 413 }
414 }
415 }
c6541a0c 416 T.SetRotation(Ax, M_PI);
7fd59977 417 P1.Transform(T);
418
419 PerformExtPElC(anExt, P1, anACurve, mytolv);
420 if (anExt.IsDone()) {
421 myDone = Standard_True;
422
c6541a0c 423 U += M_PI;
7fd59977 424
425 for (i=1; i<=anExt.NbExt(); i++) {
426 Extrema_POnCurv POC=anExt.Point(i);
427 V = POC.Parameter();
428 if (V > myvsup) {
429 // if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
430 // Standard_True, anExt.IsMin(i))) continue;
431
432 Standard_Real newV = myvsup;
433
434 if((anACurve->GetType() == GeomAbs_Circle) ||
435 (anACurve->GetType() == GeomAbs_Ellipse)) {
c6541a0c 436 newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
7fd59977 437
438 if (newV > myvsup) {
a20400b3 439 newV -= 2. * M_PI;
440
441 if (newV + mytolv < myvinf) {
442 newV = myvsup;
443 } else if (newV < myvinf) {
444 newV = myvinf;
445 }
7fd59977 446 }
447 }
448 V = newV;
449
450 if ( !IsExtremum (U, V, P, myS, E, Dist2,
451 Standard_True, anExt.IsMin(i))) continue;
452 } else if (V < myvinf) {
453 // if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
454 // Standard_False, anExt.IsMin(i))) continue;
455 Standard_Real newV = myvinf;
456
457 if((anACurve->GetType() == GeomAbs_Circle) ||
458 (anACurve->GetType() == GeomAbs_Ellipse)) {
c6541a0c 459 newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
7fd59977 460
a20400b3 461 if(newV < myvinf) {
462 newV += 2. * M_PI;
463
464 if (newV - mytolv > myvsup) {
465 newV = myvinf;
466 } else if (newV > myvsup) {
467 newV = myvsup;
468 }
469 }
7fd59977 470 }
471 V = newV;
472
473 if ( !IsExtremum (U, V, P, myS, E, Dist2,
474 Standard_False, anExt.IsMin(i))) continue;
475 } else {
476 E = myS->Value(U,V);
477 Dist2 = P.SquareDistance(E);
478 }
479 if (IsOriginalPnt(E, myPoint, myNbExt)) {
6f94f1df 480 myPoint[myNbExt] = Extrema_POnSurf(U,V,E);
7fd59977 481 mySqDist[myNbExt] = Dist2;
6f94f1df 482 myNbExt++;
7fd59977 483 }
484 }
485 }
486}
487
488
489//=======================================================================
490//function : IsDone
491//purpose :
492//=======================================================================
493
494Standard_Boolean Extrema_ExtPRevS::IsDone() const
495{
496 return myDone;
497}
498
499
500//=======================================================================
501//function : NbExt
502//purpose :
503//=======================================================================
504
505Standard_Integer Extrema_ExtPRevS::NbExt() const
506{
507 if (!IsDone()) { StdFail_NotDone::Raise(); }
508 return myNbExt;
509}
510
511
512//=======================================================================
513//function : Value
514//purpose :
515//=======================================================================
516
517Standard_Real Extrema_ExtPRevS::SquareDistance(const Standard_Integer N) const
518{
519 if (!IsDone()) { StdFail_NotDone::Raise(); }
520 if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
521 if (myIsAnalyticallyComputable)
6f94f1df 522 return mySqDist[N-1];
7fd59977 523 else
524 return myExtPS.SquareDistance(N);
525}
526//=======================================================================
527//function : Point
528//purpose :
529//=======================================================================
530
531Extrema_POnSurf Extrema_ExtPRevS::Point(const Standard_Integer N) const
532{
533 if (!IsDone()) { StdFail_NotDone::Raise(); }
534 if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
535 if (myIsAnalyticallyComputable)
6f94f1df 536 return myPoint[N-1];
7fd59977 537 else
538 return myExtPS.Point(N);
539}
540
541
Open CASCADE Technology repository