1 // Created on: 1999-09-16
2 // Created by: Edward AGAPOV
3 // Copyright (c) 1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
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.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
17 #include <Standard_NotImplemented.hxx>
18 #include <Standard_OutOfRange.hxx>
19 #include <StdFail_NotDone.hxx>
20 #include <Adaptor3d_Curve.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>
34 #include <math_FunctionSetRoot.hxx>
35 #include <math_Vector.hxx>
36 #include <GeomAdaptor_SurfaceOfLinearExtrusion.hxx>
39 IMPLEMENT_STANDARD_RTTIEXT(Extrema_ExtPExtS,Standard_Transient)
41 static gp_Ax2 GetPosition (const Handle(Adaptor3d_Curve)& C);
43 static void PerformExtPElC (Extrema_ExtPElC& E,
45 const Handle(Adaptor3d_Curve)& C,
46 const Standard_Real Tol);
48 static Standard_Boolean
49 IsCaseAnalyticallyComputable (const GeomAbs_CurveType& theType,
50 const gp_Ax2& theCurvePos,
51 const gp_Dir& theSurfaceDirection) ;
53 static gp_Pnt GetValue(const Standard_Real U,
54 const Handle(Adaptor3d_Curve)& C);
55 //=======================================================================
57 //purpose : Returns the projection of a point <Point> on a plane
58 // <ThePlane> along a direction <TheDir>.
59 //=======================================================================
60 static gp_Pnt ProjectPnt(const gp_Ax2& ThePlane,
64 gp_Vec PO(Point,ThePlane.Location());
65 Standard_Real Alpha = PO * gp_Vec(ThePlane.Direction());
66 Alpha /= TheDir * ThePlane.Direction();
68 P.SetXYZ(Point.XYZ() + Alpha * TheDir.XYZ());
72 //=======================================================================
73 //function : IsOriginalPnt
75 //=======================================================================
77 static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
78 const Extrema_POnSurf* Points,
79 const Standard_Integer NbPoints)
81 for (Standard_Integer i=1; i<=NbPoints; i++) {
82 if (Points[i-1].Value().IsEqual(P, Precision::Confusion())) {
83 return Standard_False;
89 //=======================================================================
90 //function : MakePreciser
92 //=======================================================================
94 void Extrema_ExtPExtS::MakePreciser (Standard_Real& U,
96 const Standard_Boolean isMin,
97 const gp_Ax2& OrtogSection) const
101 } else if (U<myuinf) {
105 Standard_Real step = (myusup - myuinf) / 30, D2e, D2next ,D2prev;
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;
115 notFound = (D2e > D2prev || D2e > D2next);
117 notFound = (D2e < D2prev || D2e < D2next);
119 if (notFound && (D2e < D2next && isMin)) {
136 Pnext = ProjectPnt (OrtogSection, myDirection, GetValue(U+step, myC));
137 D2next = P.SquareDistance(Pnext);
139 notFound = D2e > D2next;
141 notFound = D2e < D2next;
145 //=============================================================================
147 Extrema_ExtPExtS::Extrema_ExtPExtS()
154 myIsAnalyticallyComputable(Standard_False),
155 myDone(Standard_False),
158 for (size_t anIdx = 0; anIdx < sizeof (mySqDist) / sizeof (mySqDist[0]); anIdx++)
160 mySqDist[anIdx] = RealLast();
164 //=============================================================================
166 Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt& theP,
167 const Handle(GeomAdaptor_SurfaceOfLinearExtrusion)& theS,
168 const Standard_Real theUmin,
169 const Standard_Real theUsup,
170 const Standard_Real theVmin,
171 const Standard_Real theVsup,
172 const Standard_Real theTolU,
173 const Standard_Real theTolV)
181 myIsAnalyticallyComputable(Standard_False),
182 myDone(Standard_False),
185 for (size_t anIdx = 0; anIdx < sizeof (mySqDist) / sizeof (mySqDist[0]); anIdx++)
187 mySqDist[anIdx] = RealLast();
199 //=============================================================================
201 Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt& theP,
202 const Handle(GeomAdaptor_SurfaceOfLinearExtrusion)& theS,
203 const Standard_Real theTolU,
204 const Standard_Real theTolV)
205 : myuinf(theS->FirstUParameter()),
206 myusup(theS->LastUParameter()),
208 myvinf(theS->FirstVParameter()),
209 myvsup(theS->LastVParameter()),
212 myIsAnalyticallyComputable(Standard_False),
213 myDone(Standard_False),
216 for (size_t anIdx = 0; anIdx < sizeof (mySqDist) / sizeof (mySqDist[0]); anIdx++)
218 mySqDist[anIdx] = RealLast();
221 theS->FirstUParameter(),
222 theS->LastUParameter(),
223 theS->FirstVParameter(),
224 theS->LastVParameter(),
231 //=======================================================================
232 //function : Initialize
234 //=======================================================================
236 void Extrema_ExtPExtS::Initialize (const Handle(GeomAdaptor_SurfaceOfLinearExtrusion)& theS,
237 const Standard_Real theUinf,
238 const Standard_Real theUsup,
239 const Standard_Real theVinf,
240 const Standard_Real theVsup,
241 const Standard_Real theTolU,
242 const Standard_Real theTolV)
252 myIsAnalyticallyComputable = Standard_False;
253 myDone = Standard_False;
256 Handle(Adaptor3d_Curve) anACurve = theS->BasisCurve();
258 myF.Initialize (*theS);
261 myPosition = GetPosition(myC);
262 myDirection = theS->Direction();
263 myIsAnalyticallyComputable = //Standard_False;
264 IsCaseAnalyticallyComputable (myC->GetType(), myPosition, myDirection);
266 if (!myIsAnalyticallyComputable)
268 myExtPS.Initialize (*theS,
281 //=======================================================================
284 //=======================================================================
286 void Extrema_ExtPExtS::Perform (const gp_Pnt& P)
288 const Standard_Integer NbExtMax = 4; //dimension of arrays
289 //myPoint[] and mySqDist[]
290 //For "analytical" case
291 myDone = Standard_False;
294 if (!myIsAnalyticallyComputable) {
296 myDone = myExtPS.IsDone();
297 // modified by NIZHNY-EAP Wed Nov 17 12:59:08 1999 ___BEGIN___
298 myNbExt = myExtPS.NbExt();
299 // modified by NIZHNY-EAP Wed Nov 17 12:59:09 1999 ___END___
303 gp_Pnt Pe, Pp = ProjectPnt(myPosition,myDirection,P);
304 Extrema_ExtPElC anExt;
305 PerformExtPElC(anExt, Pp, myC, mytolu);
306 if (!anExt.IsDone()) return;
308 gp_Ax2 anOrtogSection (P, myDirection);
313 myDirection.IsParallel(myPosition.Direction(),Precision::Angular());
314 Standard_Integer i, aNbExt = anExt.NbExt();
315 math_Vector UV(1,2), Tol(1,2), UVinf(1,2), UVsup(1,2);
316 Tol(1) = mytolu; Tol(2) = mytolv;
317 UVinf(1) = myuinf; UVinf(2) = myvinf;
318 UVsup(1) = myusup; UVsup(2) = myvsup;
321 for (i=1; i<=aNbExt; i++) {
322 Extrema_POnCurv POC=anExt.Point(i);
324 //// modified by jgv, 23.12.2008 for OCC17194 ////
325 if (myC->IsPeriodic())
327 Standard_Real U2 = U;
328 ElCLib::AdjustPeriodic(myuinf, myuinf + 2.*M_PI, Precision::PConfusion(), U, U2);
330 //////////////////////////////////////////////////
331 gp_Pnt E = POC.Value();
332 Pe = ProjectPnt(anOrtogSection, myDirection, E);
335 V = gp_Vec(E,Pe) * gp_Vec(myDirection);
336 // modified by NIZHNY-MKK Thu Sep 18 14:46:14 2003.BEGIN
337 // myPoint[++myNbExt] = Extrema_POnSurf(U, V, Pe);
338 // myValue[myNbExt] = anExt.Value(i);
339 myPoint[myNbExt] = Extrema_POnSurf(U, V, Pe);
340 mySqDist[myNbExt] = anExt.SquareDistance(i);
342 if(myNbExt == NbExtMax)
346 // modified by NIZHNY-MKK Thu Sep 18 14:46:18 2003.END
350 isMin = anExt.IsMin(i);//( Pp.Distance(GetValue(U+10,myC)) > anExt.Value(i) );
352 MakePreciser(U, P, isMin, anOrtogSection);
353 E = GetValue(U, myC);
354 Pe = ProjectPnt (anOrtogSection, myDirection, E),
355 V = gp_Vec(E,Pe) * gp_Vec(myDirection);
356 UV(1) = U; UV(2) = V;
357 math_FunctionSetRoot aFSR (myF, Tol);
358 aFSR.Perform(myF, UV, UVinf, UVsup);
359 // for (Standard_Integer k=1 ; k <= myF.NbExt();
361 for ( k=1 ; k <= myF.NbExt(); k++) {
362 if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
363 // modified by NIZHNY-MKK Thu Sep 18 14:46:41 2003.BEGIN
364 // myPoint[++myNbExt] = myF.Point(k);
365 // myValue[myNbExt] = myF.Value(k);
366 myPoint[myNbExt] = myF.Point(k);
367 mySqDist[myNbExt] = myF.SquareDistance(k);
369 if(myNbExt == NbExtMax)
373 // modified by NIZHNY-MKK Thu Sep 18 14:46:43 2003.END
376 if(myNbExt == NbExtMax)
380 // try symmetric point
381 myF.SetPoint(P); //To clear previous solutions
383 MakePreciser(U, P, isMin, anOrtogSection);
384 E = GetValue(U, myC);
385 Pe = ProjectPnt (anOrtogSection, myDirection, E),
386 V = gp_Vec(E,Pe) * gp_Vec(myDirection);
387 UV(1) = U; UV(2) = V;
389 aFSR.Perform (myF,UV,UVinf,UVsup);
391 for (k=1 ; k <= myF.NbExt(); k++) {
392 if(myF.SquareDistance(k) > Precision::Confusion()*Precision::Confusion())
394 //Additional checking solution: FSR sometimes is wrong
395 //when starting point is far from solution.
396 Standard_Real dist = Sqrt(myF.SquareDistance(k));
397 math_Vector Vals(1, 2);
398 const Extrema_POnSurf& PonS=myF.Point(k);
400 PonS.Parameter(u, v);
405 myS->D1(u, v, Pe, du, dv);
406 Standard_Real mdu = du.Magnitude();
407 Standard_Real mdv = dv.Magnitude();
410 if(mdu > Precision::PConfusion())
412 if(u / dist / mdu > Precision::PConfusion())
417 if(mdv > Precision::PConfusion())
419 if(v / dist / mdv > Precision::PConfusion())
426 if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
427 // modified by NIZHNY-MKK Thu Sep 18 14:46:59 2003.BEGIN
428 // myPoint[++myNbExt] = myF.Point(k);
429 // myValue[myNbExt] = myF.Value(k);
430 myPoint[myNbExt] = myF.Point(k);
431 mySqDist[myNbExt] = myF.SquareDistance(k);
433 if(myNbExt == NbExtMax)
437 // modified by NIZHNY-MKK Thu Sep 18 14:47:04 2003.END
440 if(myNbExt == NbExtMax)
446 myDone = Standard_True;
450 //=============================================================================
452 Standard_Boolean Extrema_ExtPExtS::IsDone () const { return myDone; }
453 //=============================================================================
455 Standard_Integer Extrema_ExtPExtS::NbExt () const
457 if (!IsDone()) { throw StdFail_NotDone(); }
458 if (myIsAnalyticallyComputable)
461 return myExtPS.NbExt();
463 //=============================================================================
465 Standard_Real Extrema_ExtPExtS::SquareDistance (const Standard_Integer N) const
467 if ((N < 1) || (N > NbExt()))
469 throw Standard_OutOfRange();
471 if (myIsAnalyticallyComputable)
472 // modified by NIZHNY-MKK Thu Sep 18 14:48:39 2003.BEGIN
473 // return myValue[N];
474 return mySqDist[N-1];
475 // modified by NIZHNY-MKK Thu Sep 18 14:48:42 2003.END
477 return myExtPS.SquareDistance(N);
479 //=============================================================================
481 const Extrema_POnSurf& Extrema_ExtPExtS::Point (const Standard_Integer N) const
483 if ((N < 1) || (N > NbExt()))
485 throw Standard_OutOfRange();
487 if (myIsAnalyticallyComputable) {
488 // modified by NIZHNY-MKK Thu Sep 18 14:47:40 2003.BEGIN
489 // return myPoint[N];
492 // modified by NIZHNY-MKK Thu Sep 18 14:47:43 2003.END
494 return myExtPS.Point(N);
496 //=============================================================================
499 static gp_Ax2 GetPosition (const Handle(Adaptor3d_Curve)& C)
501 switch (C->GetType()) {
503 gp_Lin L = C->Line();
504 gp_Pln Pln = gp_Pln(L.Location(), L.Direction());
505 //:abv 30.05.02: OCC - use constructor instead of Set...s() to avoid exception
506 gp_Ax2 Pos ( Pln.Location(), Pln.Position().Direction(), Pln.Position().XDirection() );
507 // Pos.SetAxis(Pln.XAxis());
508 // Pos.SetXDirection(Pln.YAxis().Direction());
509 // Pos.SetYDirection(Pln.Position().Direction());
513 return C->Circle().Position();
514 case GeomAbs_Ellipse:
515 return C->Ellipse().Position();
516 case GeomAbs_Hyperbola:
517 return C->Hyperbola().Position();
518 case GeomAbs_Parabola:
519 return C->Parabola().Position();
524 //=============================================================================
526 static void PerformExtPElC (Extrema_ExtPElC& E,
528 const Handle(Adaptor3d_Curve)& C,
529 const Standard_Real Tol)
531 switch (C->GetType()) {
532 case GeomAbs_Hyperbola:
533 E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
536 E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
539 E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
541 case GeomAbs_Ellipse:
542 E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
544 case GeomAbs_Parabola:
545 E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
552 //=======================================================================
553 //function : IsCaseAnalyticallyComputable
555 //=======================================================================
557 static Standard_Boolean IsCaseAnalyticallyComputable
558 (const GeomAbs_CurveType& theType,
559 const gp_Ax2& theCurvePos,
560 const gp_Dir& theSurfaceDirection)
566 case GeomAbs_Ellipse:
567 case GeomAbs_Hyperbola:
568 case GeomAbs_Parabola:
571 return Standard_False;
573 // check if it is a plane
574 if (Abs(theCurvePos.Direction() * theSurfaceDirection) <= gp::Resolution())
575 return Standard_False;
577 return Standard_True;
579 //=======================================================================
580 //function : GetValue
582 //=======================================================================
584 static gp_Pnt GetValue(const Standard_Real U,
585 const Handle(Adaptor3d_Curve)& C)
587 switch (C->GetType()) {
589 return ElCLib::Value(U, C->Line());
591 return ElCLib::Value(U, C->Circle());
592 case GeomAbs_Ellipse:
593 return ElCLib::Value(U, C->Ellipse());
594 case GeomAbs_Hyperbola:
595 return ElCLib::Value(U, C->Hyperbola());
596 case GeomAbs_Parabola:
597 return ElCLib::Value(U, C->Parabola());
602 //=======================================================================
605 //=======================================================================
607 //static Standard_Real GetU(const gp_Vec& vec,
609 // const Handle(Adaptor3d_Curve)& C)
611 // switch (C->GetType()) {
612 // case GeomAbs_Line:
613 // return ElCLib::Parameter(C->Line().Translated(vec), P);
614 // case GeomAbs_Circle:
615 // return ElCLib::Parameter(C->Circle().Translated(vec), P);
616 // case GeomAbs_Ellipse:
617 // return ElCLib::Parameter(C->Ellipse().Translated(vec), P);
618 // case GeomAbs_Hyperbola:
619 // return ElCLib::Parameter(C->Hyperbola().Translated(vec), P);
620 // case GeomAbs_Parabola:
621 // return ElCLib::Parameter(C->Parabola().Translated(vec), P);