projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0031668: Visualization - WebGL sample doesn't work on Emscripten 1.39
[occt.git] / src / Extrema / Extrema_ExtPExtS.cxx
CommitLineData
b311480e 1// Created on: 1999-09-16
2// Created by: Edward AGAPOV
3// Copyright (c) 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
17#include <Standard_NotImplemented.hxx>
18#include <Standard_OutOfRange.hxx>
19#include <StdFail_NotDone.hxx>
20#include <Adaptor3d_HCurve.hxx>
21#include <ElCLib.hxx>
22#include <Extrema_ExtPElC.hxx>
23#include <Extrema_ExtPExtS.hxx>
24#include <Extrema_POnCurv.hxx>
25#include <Extrema_POnSurf.hxx>
26#include <Precision.hxx>
27#include <gp.hxx>
28#include <gp_Ax2.hxx>
29#include <gp_Dir.hxx>
30#include <gp_Lin.hxx>
31#include <gp_Pln.hxx>
32#include <gp_Pnt.hxx>
33#include <gp_Vec.hxx>
34#include <math_FunctionSetRoot.hxx>
35#include <math_Vector.hxx>
6b84c3f7 36#include <GeomAdaptor_HSurfaceOfLinearExtrusion.hxx>
0734c53d 37
7fd59977 38
92efcf78 39IMPLEMENT_STANDARD_RTTIEXT(Extrema_ExtPExtS,Standard_Transient)
40
7fd59977 41static gp_Ax2 GetPosition (const Handle(Adaptor3d_HCurve)& C);
42
43static void PerformExtPElC (Extrema_ExtPElC& E,
44 const gp_Pnt& P,
45 const Handle(Adaptor3d_HCurve)& C,
46 const Standard_Real Tol);
47
48static Standard_Boolean
49 IsCaseAnalyticallyComputable (const GeomAbs_CurveType& theType,
50 const gp_Ax2& theCurvePos,
51 const gp_Dir& theSurfaceDirection) ;
52
53static gp_Pnt GetValue(const Standard_Real U,
54 const Handle(Adaptor3d_HCurve)& C);
55//=======================================================================
56//function : Project
57//purpose : Returns the projection of a point <Point> on a plane
58// <ThePlane> along a direction <TheDir>.
59//=======================================================================
60static gp_Pnt ProjectPnt(const gp_Ax2& ThePlane,
61 const gp_Dir& TheDir,
62 const gp_Pnt& Point)
63{
64 gp_Vec PO(Point,ThePlane.Location());
65 Standard_Real Alpha = PO * gp_Vec(ThePlane.Direction());
66 Alpha /= TheDir * ThePlane.Direction();
67 gp_Pnt P;
68 P.SetXYZ(Point.XYZ() + Alpha * TheDir.XYZ());
69 return P;
70}
71
72//=======================================================================
73//function : IsOriginalPnt
74//purpose :
75//=======================================================================
76
77static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
78 const Extrema_POnSurf* Points,
79 const Standard_Integer NbPoints)
80{
81 for (Standard_Integer i=1; i<=NbPoints; i++) {
82 if (Points[i-1].Value().IsEqual(P, Precision::Confusion())) {
83 return Standard_False;
84 }
85 }
86 return Standard_True;
87}
88
89//=======================================================================
90//function : MakePreciser
91//purpose :
92//=======================================================================
93
94void Extrema_ExtPExtS::MakePreciser (Standard_Real& U,
95 const gp_Pnt& P,
96 const Standard_Boolean isMin,
97 const gp_Ax2& OrtogSection) const
98{
99 if (U>myusup) {
100 U = myusup;
101 } else if (U<myuinf) {
102 U = myuinf;
103 } else {
104
105 Standard_Real step = (myusup - myuinf) / 30, D2e, D2next ,D2prev;
106 gp_Pnt
107 Pe = ProjectPnt (OrtogSection, myDirection, GetValue(U,myC)),
108 Pprev = ProjectPnt (OrtogSection, myDirection, GetValue(U-step, myC)),
109 Pnext = ProjectPnt (OrtogSection, myDirection, GetValue(U+step, myC));
110 D2e = P.SquareDistance(Pe),
111 D2next = P.SquareDistance(Pnext),
112 D2prev = P.SquareDistance(Pprev);
113 Standard_Boolean notFound;
114 if (isMin)
115 notFound = (D2e > D2prev || D2e > D2next);
116 else
117 notFound = (D2e < D2prev || D2e < D2next);
118
119 if (notFound && (D2e < D2next && isMin)) {
120 step = -step;
121 D2next = D2prev;
122 Pnext = Pprev;
123 }
124 while (notFound) {
125 U = U + step;
126 if (U > myusup) {
aabe3a17 127 U = myusup;
128 break;
7fd59977 129 }
130 if (U < myuinf) {
aabe3a17 131 U = myuinf;
132 break;
7fd59977 133 }
134 D2e = D2next;
135 Pe = Pnext;
136 Pnext = ProjectPnt (OrtogSection, myDirection, GetValue(U+step, myC));
137 D2next = P.SquareDistance(Pnext);
138 if (isMin)
aabe3a17 139 notFound = D2e > D2next;
7fd59977 140 else
aabe3a17 141 notFound = D2e < D2next;
7fd59977 142 }
143 }
144}
145//=============================================================================
146
cbff1e55 147Extrema_ExtPExtS::Extrema_ExtPExtS()
148: myuinf(0.0),
149 myusup(0.0),
150 mytolu(0.0),
151 myvinf(0.0),
152 myvsup(0.0),
153 mytolv(0.0),
154 myIsAnalyticallyComputable(Standard_False),
155 myDone(Standard_False),
638ad7f3 156 myNbExt(0)
7fd59977 157{
7fd59977 158}
159
160//=============================================================================
161
0734c53d 162Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt& theP,
6b84c3f7 163 const Handle(GeomAdaptor_HSurfaceOfLinearExtrusion)& theS,
0734c53d 164 const Standard_Real theUmin,
165 const Standard_Real theUsup,
166 const Standard_Real theVmin,
167 const Standard_Real theVsup,
168 const Standard_Real theTolU,
169 const Standard_Real theTolV)
cbff1e55 170: myuinf(theUmin),
171 myusup(theUsup),
172 mytolu(theTolU),
173 myvinf(theVmin),
174 myvsup(theVsup),
175 mytolv(theTolV),
176 myS (theS),
177 myIsAnalyticallyComputable(Standard_False),
178 myDone(Standard_False),
638ad7f3 179 myNbExt(0)
7fd59977 180{
0734c53d 181 Initialize (theS,
182 theUmin,
183 theUsup,
184 theVmin,
185 theVsup,
186 theTolU,
187 theTolV);
188
189 Perform (theP);
7fd59977 190}
191//=============================================================================
192
0734c53d 193Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt& theP,
6b84c3f7 194 const Handle(GeomAdaptor_HSurfaceOfLinearExtrusion)& theS,
0734c53d 195 const Standard_Real theTolU,
196 const Standard_Real theTolV)
cbff1e55 197: myuinf(theS->FirstUParameter()),
198 myusup(theS->LastUParameter()),
199 mytolu(theTolU),
200 myvinf(theS->FirstVParameter()),
201 myvsup(theS->LastVParameter()),
202 mytolv(theTolV),
203 myS (theS),
204 myIsAnalyticallyComputable(Standard_False),
205 myDone(Standard_False),
638ad7f3 206 myNbExt(0)
7fd59977 207{
0734c53d 208 Initialize (theS,
209 theS->FirstUParameter(),
210 theS->LastUParameter(),
211 theS->FirstVParameter(),
212 theS->LastVParameter(),
213 theTolU,
214 theTolV);
215
216 Perform (theP);
7fd59977 217}
cbff1e55 218
7fd59977 219//=======================================================================
220//function : Initialize
221//purpose :
222//=======================================================================
223
6b84c3f7 224void Extrema_ExtPExtS::Initialize (const Handle(GeomAdaptor_HSurfaceOfLinearExtrusion)& theS,
0734c53d 225 const Standard_Real theUinf,
226 const Standard_Real theUsup,
227 const Standard_Real theVinf,
228 const Standard_Real theVsup,
229 const Standard_Real theTolU,
230 const Standard_Real theTolV)
7fd59977 231{
0734c53d 232 myuinf = theUinf;
233 myusup = theUsup;
234 mytolu = theTolU;
7fd59977 235
0734c53d 236 myvinf = theVinf;
237 myvsup = theVsup;
238 mytolv = theTolV;
7fd59977 239
638ad7f3 240 myIsAnalyticallyComputable = Standard_False;
241 myDone = Standard_False;
242 myNbExt = 0;
243
0734c53d 244 Handle(Adaptor3d_HCurve) anACurve = theS->BasisCurve();
7fd59977 245
0734c53d 246 myF.Initialize (theS->ChangeSurface());
7fd59977 247 myC = anACurve;
0734c53d 248 myS = theS;
7fd59977 249 myPosition = GetPosition(myC);
0734c53d 250 myDirection = theS->Direction();
7fd59977 251 myIsAnalyticallyComputable = //Standard_False;
0734c53d 252 IsCaseAnalyticallyComputable (myC->GetType(), myPosition, myDirection);
7fd59977 253
254 if (!myIsAnalyticallyComputable)
0734c53d 255 {
256 myExtPS.Initialize (theS->ChangeSurface(),
257 32,
258 32,
259 theUinf,
260 theUsup,
261 theVinf,
262 theVsup,
263 theTolU,
264 theTolV);
265 }
7fd59977 266}
267
268
269//=======================================================================
270//function : Perform
271//purpose :
272//=======================================================================
273
274void Extrema_ExtPExtS::Perform (const gp_Pnt& P)
275{
aabe3a17 276 const Standard_Integer NbExtMax = 4; //dimension of arrays
277 //myPoint[] and mySqDist[]
278 //For "analytical" case
7fd59977 279 myDone = Standard_False;
280 myNbExt = 0;
281
282 if (!myIsAnalyticallyComputable) {
283 myExtPS.Perform(P);
284 myDone = myExtPS.IsDone();
285// modified by NIZHNY-EAP Wed Nov 17 12:59:08 1999 ___BEGIN___
286 myNbExt = myExtPS.NbExt();
287// modified by NIZHNY-EAP Wed Nov 17 12:59:09 1999 ___END___
288 return;
289 }
290
291 gp_Pnt Pe, Pp = ProjectPnt(myPosition,myDirection,P);
292 Extrema_ExtPElC anExt;
293 PerformExtPElC(anExt, Pp, myC, mytolu);
294 if (!anExt.IsDone()) return;
295
296 gp_Ax2 anOrtogSection (P, myDirection);
297 Standard_Real U,V;
298 Standard_Boolean
299 isMin,
300 isSimpleCase =
301 myDirection.IsParallel(myPosition.Direction(),Precision::Angular());
302 Standard_Integer i, aNbExt = anExt.NbExt();
303 math_Vector UV(1,2), Tol(1,2), UVinf(1,2), UVsup(1,2);
304 Tol(1) = mytolu; Tol(2) = mytolv;
305 UVinf(1) = myuinf; UVinf(2) = myvinf;
306 UVsup(1) = myusup; UVsup(2) = myvsup;
307
308
309 for (i=1; i<=aNbExt; i++) {
310 Extrema_POnCurv POC=anExt.Point(i);
311 U = POC.Parameter();
312 //// modified by jgv, 23.12.2008 for OCC17194 ////
313 if (myC->IsPeriodic())
aabe3a17 314 {
315 Standard_Real U2 = U;
316 ElCLib::AdjustPeriodic(myuinf, myuinf + 2.*M_PI, Precision::PConfusion(), U, U2);
317 }
7fd59977 318 //////////////////////////////////////////////////
319 gp_Pnt E = POC.Value();
320 Pe = ProjectPnt(anOrtogSection, myDirection, E);
321
322 if (isSimpleCase) {
323 V = gp_Vec(E,Pe) * gp_Vec(myDirection);
324 // modified by NIZHNY-MKK Thu Sep 18 14:46:14 2003.BEGIN
325 // myPoint[++myNbExt] = Extrema_POnSurf(U, V, Pe);
326 // myValue[myNbExt] = anExt.Value(i);
327 myPoint[myNbExt] = Extrema_POnSurf(U, V, Pe);
328 mySqDist[myNbExt] = anExt.SquareDistance(i);
329 myNbExt++;
aabe3a17 330 if(myNbExt == NbExtMax)
331 {
332 break;
333 }
7fd59977 334 // modified by NIZHNY-MKK Thu Sep 18 14:46:18 2003.END
335 }
336 else {
337 myF.SetPoint(P);
338 isMin = anExt.IsMin(i);//( Pp.Distance(GetValue(U+10,myC)) > anExt.Value(i) );
339
340 MakePreciser(U, P, isMin, anOrtogSection);
341 E = GetValue(U, myC);
342 Pe = ProjectPnt (anOrtogSection, myDirection, E),
343 V = gp_Vec(E,Pe) * gp_Vec(myDirection);
344 UV(1) = U; UV(2) = V;
859a47c3 345 math_FunctionSetRoot aFSR (myF, Tol);
346 aFSR.Perform(myF, UV, UVinf, UVsup);
7fd59977 347// for (Standard_Integer k=1 ; k <= myF.NbExt();
348 Standard_Integer k;
349 for ( k=1 ; k <= myF.NbExt(); k++) {
aabe3a17 350 if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
351 // modified by NIZHNY-MKK Thu Sep 18 14:46:41 2003.BEGIN
352 // myPoint[++myNbExt] = myF.Point(k);
353 // myValue[myNbExt] = myF.Value(k);
354 myPoint[myNbExt] = myF.Point(k);
355 mySqDist[myNbExt] = myF.SquareDistance(k);
356 myNbExt++;
357 if(myNbExt == NbExtMax)
358 {
359 break;
360 }
361 // modified by NIZHNY-MKK Thu Sep 18 14:46:43 2003.END
362 }
363 }
364 if(myNbExt == NbExtMax)
365 {
366 break;
7fd59977 367 }
368 // try symmetric point
aabe3a17 369 myF.SetPoint(P); //To clear previous solutions
7fd59977 370 U *= -1;
371 MakePreciser(U, P, isMin, anOrtogSection);
372 E = GetValue(U, myC);
373 Pe = ProjectPnt (anOrtogSection, myDirection, E),
374 V = gp_Vec(E,Pe) * gp_Vec(myDirection);
375 UV(1) = U; UV(2) = V;
376
377 aFSR.Perform (myF,UV,UVinf,UVsup);
378
379 for (k=1 ; k <= myF.NbExt(); k++) {
aabe3a17 380 if(myF.SquareDistance(k) > Precision::Confusion()*Precision::Confusion())
381 {
382 //Additional checking solution: FSR sometimes is wrong
383 //when starting point is far from solution.
384 Standard_Real dist = Sqrt(myF.SquareDistance(k));
385 math_Vector Vals(1, 2);
386 const Extrema_POnSurf& PonS=myF.Point(k);
387 Standard_Real u, v;
388 PonS.Parameter(u, v);
389 UV(1) = u;
390 UV(2) = v;
391 myF.Value(UV, Vals);
392 gp_Vec du, dv;
393 myS->D1(u, v, Pe, du, dv);
394 Standard_Real mdu = du.Magnitude();
395 Standard_Real mdv = dv.Magnitude();
396 u = Abs(Vals(1));
397 v = Abs(Vals(2));
398 if(mdu > Precision::PConfusion())
399 {
400 if(u / dist / mdu > Precision::PConfusion())
401 {
402 continue;
403 }
404 }
405 if(mdv > Precision::PConfusion())
406 {
407 if(v / dist / mdv > Precision::PConfusion())
408 {
409 continue;
410 }
411 }
412
413 }
414 if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
415 // modified by NIZHNY-MKK Thu Sep 18 14:46:59 2003.BEGIN
416 // myPoint[++myNbExt] = myF.Point(k);
417 // myValue[myNbExt] = myF.Value(k);
418 myPoint[myNbExt] = myF.Point(k);
419 mySqDist[myNbExt] = myF.SquareDistance(k);
420 myNbExt++;
421 if(myNbExt == NbExtMax)
422 {
423 break;
424 }
425 // modified by NIZHNY-MKK Thu Sep 18 14:47:04 2003.END
426 }
427 }
428 if(myNbExt == NbExtMax)
429 {
430 break;
7fd59977 431 }
432 }
433 }
434 myDone = Standard_True;
435 return;
436}
437
438//=============================================================================
439
440Standard_Boolean Extrema_ExtPExtS::IsDone () const { return myDone; }
441//=============================================================================
442
443Standard_Integer Extrema_ExtPExtS::NbExt () const
444{
9775fa61 445 if (!IsDone()) { throw StdFail_NotDone(); }
7fd59977 446 if (myIsAnalyticallyComputable)
447 return myNbExt;
448 else
449 return myExtPS.NbExt();
450}
451//=============================================================================
452
453Standard_Real Extrema_ExtPExtS::SquareDistance (const Standard_Integer N) const
454{
638ad7f3 455 if ((N < 1) || (N > NbExt()))
456 {
457 throw Standard_OutOfRange();
458 }
7fd59977 459 if (myIsAnalyticallyComputable)
460 // modified by NIZHNY-MKK Thu Sep 18 14:48:39 2003.BEGIN
461 // return myValue[N];
462 return mySqDist[N-1];
463 // modified by NIZHNY-MKK Thu Sep 18 14:48:42 2003.END
464 else
465 return myExtPS.SquareDistance(N);
466}
467//=============================================================================
468
5d99f2c8 469const Extrema_POnSurf& Extrema_ExtPExtS::Point (const Standard_Integer N) const
7fd59977 470{
638ad7f3 471 if ((N < 1) || (N > NbExt()))
472 {
473 throw Standard_OutOfRange();
474 }
7fd59977 475 if (myIsAnalyticallyComputable) {
476 // modified by NIZHNY-MKK Thu Sep 18 14:47:40 2003.BEGIN
477 // return myPoint[N];
478 return myPoint[N-1];
479 }
480 // modified by NIZHNY-MKK Thu Sep 18 14:47:43 2003.END
481 else
482 return myExtPS.Point(N);
483}
484//=============================================================================
485
486
487static gp_Ax2 GetPosition (const Handle(Adaptor3d_HCurve)& C)
488{
489 switch (C->GetType()) {
490 case GeomAbs_Line: {
491 gp_Lin L = C->Line();
492 gp_Pln Pln = gp_Pln(L.Location(), L.Direction());
493 //:abv 30.05.02: OCC - use constructor instead of Set...s() to avoid exception
494 gp_Ax2 Pos ( Pln.Location(), Pln.Position().Direction(), Pln.Position().XDirection() );
495// Pos.SetAxis(Pln.XAxis());
496// Pos.SetXDirection(Pln.YAxis().Direction());
497// Pos.SetYDirection(Pln.Position().Direction());
498 return Pos;
499 }
500 case GeomAbs_Circle:
501 return C->Circle().Position();
502 case GeomAbs_Ellipse:
503 return C->Ellipse().Position();
504 case GeomAbs_Hyperbola:
505 return C->Hyperbola().Position();
506 case GeomAbs_Parabola:
507 return C->Parabola().Position();
508 default:
509 return gp_Ax2 ();
510 }
511}
512//=============================================================================
513
514static void PerformExtPElC (Extrema_ExtPElC& E,
515 const gp_Pnt& P,
516 const Handle(Adaptor3d_HCurve)& C,
517 const Standard_Real Tol)
518{
519 switch (C->GetType()) {
520 case GeomAbs_Hyperbola:
521 E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
522 return;
523 case GeomAbs_Line:
524 E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
525 return;
526 case GeomAbs_Circle:
c6541a0c 527 E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
7fd59977 528 return;
529 case GeomAbs_Ellipse:
c6541a0c 530 E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
7fd59977 531 return;
532 case GeomAbs_Parabola:
533 E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
534 return;
7fd59977 535 default:
536 return;
7fd59977 537 }
538}
539
540//=======================================================================
541//function : IsCaseAnalyticallyComputable
542//purpose :
543//=======================================================================
544
545static Standard_Boolean IsCaseAnalyticallyComputable
546 (const GeomAbs_CurveType& theType,
547 const gp_Ax2& theCurvePos,
548 const gp_Dir& theSurfaceDirection)
549{
550 // check type
551 switch (theType) {
552 case GeomAbs_Line:
553 case GeomAbs_Circle:
554 case GeomAbs_Ellipse:
555 case GeomAbs_Hyperbola:
556 case GeomAbs_Parabola:
557 break;
558 default:
559 return Standard_False;
560 }
561 // check if it is a plane
562 if (Abs(theCurvePos.Direction() * theSurfaceDirection) <= gp::Resolution())
563 return Standard_False;
564 else
565 return Standard_True;
566}
567//=======================================================================
568//function : GetValue
569//purpose :
570//=======================================================================
571
572static gp_Pnt GetValue(const Standard_Real U,
573 const Handle(Adaptor3d_HCurve)& C)
574{
575 switch (C->GetType()) {
576 case GeomAbs_Line:
577 return ElCLib::Value(U, C->Line());
578 case GeomAbs_Circle:
579 return ElCLib::Value(U, C->Circle());
580 case GeomAbs_Ellipse:
581 return ElCLib::Value(U, C->Ellipse());
582 case GeomAbs_Hyperbola:
583 return ElCLib::Value(U, C->Hyperbola());
584 case GeomAbs_Parabola:
585 return ElCLib::Value(U, C->Parabola());
586 default:
587 return gp_Pnt();
588 }
589}
590//=======================================================================
591//function : GetU
592//purpose :
593//=======================================================================
0797d9d3 594//#ifdef OCCT_DEBUG
7fd59977 595//static Standard_Real GetU(const gp_Vec& vec,
596// const gp_Pnt& P,
597// const Handle(Adaptor3d_HCurve)& C)
598//{
599// switch (C->GetType()) {
600// case GeomAbs_Line:
601// return ElCLib::Parameter(C->Line().Translated(vec), P);
602// case GeomAbs_Circle:
603// return ElCLib::Parameter(C->Circle().Translated(vec), P);
604// case GeomAbs_Ellipse:
605// return ElCLib::Parameter(C->Ellipse().Translated(vec), P);
606// case GeomAbs_Hyperbola:
607// return ElCLib::Parameter(C->Hyperbola().Translated(vec), P);
608// case GeomAbs_Parabola:
609// return ElCLib::Parameter(C->Parabola().Translated(vec), P);
610// default:
611// return 0;
612// }
613//}
614//#endif
Open CASCADE Technology repository