1 // Created on: 1995-07-18
2 // Created by: Modelistation
3 // Copyright (c) 1995-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
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.
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.
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.
22 // Modified by skv - Thu Sep 30 15:21:07 2004 OCC593
25 #include <Extrema_GenExtPS.ixx>
26 #include <StdFail_NotDone.hxx>
27 #include <Standard_OutOfRange.hxx>
28 #include <TColStd_Array2OfInteger.hxx>
29 #include <TColStd_Array2OfReal.hxx>
30 #include <TColgp_Array2OfPnt.hxx>
31 #include <math_FunctionSetRoot.hxx>
32 #include <math_Vector.hxx>
33 #include <math_NewtonFunctionSetRoot.hxx>
34 #include <GeomAbs_IsoType.hxx>
35 #include <Bnd_Sphere.hxx>
36 #include <Extrema_HUBTreeOfSphere.hxx>
37 #include <Extrema_ExtFlag.hxx>
38 #include <Bnd_Array1OfSphere.hxx>
39 #include <Bnd_HArray1OfSphere.hxx>
40 #include <Precision.hxx>
41 #include <Geom_OffsetSurface.hxx>
42 #include <Geom_RectangularTrimmedSurface.hxx>
43 #include <Geom_BSplineSurface.hxx>
44 #include <Geom_BezierSurface.hxx>
45 #include <Adaptor3d_HSurface.hxx>
46 #include <Adaptor3d_HCurve.hxx>
47 #include <Adaptor3d_Curve.hxx>
48 #include <Geom_BSplineCurve.hxx>
49 #include <Geom_BezierCurve.hxx>
50 //IMPLEMENT_HARRAY1(Extrema_HArray1OfSphere)
53 class Bnd_SphereUBTreeSelector : public Extrema_UBTreeOfSphere::Selector
57 Bnd_SphereUBTreeSelector (const Handle(Bnd_HArray1OfSphere)& theSphereArray,
60 mySphereArray(theSphereArray),
63 //myXYZ = gp_Pnt(0, 0, 0);
66 void DefineCheckPoint( const gp_Pnt& theXYZ )
69 Bnd_Sphere& Sphere() const
72 virtual Standard_Boolean Reject( const Bnd_Sphere &theBnd ) const = 0;
74 virtual Standard_Boolean Accept(const Standard_Integer& theObj) = 0;
78 const Handle(Bnd_HArray1OfSphere)& mySphereArray;
83 class Bnd_SphereUBTreeSelectorMin : public Bnd_SphereUBTreeSelector
86 Bnd_SphereUBTreeSelectorMin (const Handle(Bnd_HArray1OfSphere)& theSphereArray,
88 : Bnd_SphereUBTreeSelector(theSphereArray, theSol),
92 void SetMinDist( const Standard_Real theMinDist )
93 { myMinDist = theMinDist; }
95 Standard_Real MinDist() const
98 Standard_Boolean Reject( const Bnd_Sphere &theBnd ) const
100 Bnd_SphereUBTreeSelectorMin* me =
101 const_cast<Bnd_SphereUBTreeSelectorMin*>(this);
102 // myMinDist is decreased each time a nearer object is found
103 return theBnd.IsOut( myXYZ.XYZ(), me->myMinDist );
106 Standard_Boolean Accept(const Standard_Integer&);
109 Standard_Real myMinDist;
112 Standard_Boolean Bnd_SphereUBTreeSelectorMin::Accept(const Standard_Integer& theInd)
114 const Bnd_Sphere& aSph = mySphereArray->Value(theInd);
115 Standard_Real aCurDist;
117 if ( (aCurDist = aSph.SquareDistance(myXYZ.XYZ())) < mySol.SquareDistance(myXYZ.XYZ()) )
120 if ( aCurDist < myMinDist )
121 myMinDist = aCurDist;
123 return Standard_True;
126 return Standard_False;
129 class Bnd_SphereUBTreeSelectorMax : public Bnd_SphereUBTreeSelector
132 Bnd_SphereUBTreeSelectorMax (const Handle(Bnd_HArray1OfSphere)& theSphereArray,
134 : Bnd_SphereUBTreeSelector(theSphereArray, theSol),
138 void SetMaxDist( const Standard_Real theMaxDist )
139 { myMaxDist = theMaxDist; }
141 Standard_Real MaxDist() const
142 { return myMaxDist; }
144 Standard_Boolean Reject( const Bnd_Sphere &theBnd ) const
146 Bnd_SphereUBTreeSelectorMax* me =
147 const_cast<Bnd_SphereUBTreeSelectorMax*>(this);
148 // myMaxDist is decreased each time a nearer object is found
149 return theBnd.IsOut( myXYZ.XYZ(), me->myMaxDist );
152 Standard_Boolean Accept(const Standard_Integer&);
155 Standard_Real myMaxDist;
158 Standard_Boolean Bnd_SphereUBTreeSelectorMax::Accept(const Standard_Integer& theInd)
160 const Bnd_Sphere& aSph = mySphereArray->Value(theInd);
161 Standard_Real aCurDist;
163 if ( (aCurDist = aSph.SquareDistance(myXYZ.XYZ())) > mySol.SquareDistance(myXYZ.XYZ()) )
166 if ( aCurDist > myMaxDist )
167 myMaxDist = aCurDist;
169 return Standard_True;
172 return Standard_False;
176 * This function computes the point on surface parameters on edge.
177 * if it coincides with theParam0 or theParam1, it is returned.
179 static Extrema_POnSurfParams ComputeEdgeParameters
180 (const Standard_Boolean IsUEdge,
181 const Extrema_POnSurfParams &theParam0,
182 const Extrema_POnSurfParams &theParam1,
183 const Adaptor3d_SurfacePtr &theSurf,
184 const gp_Pnt &thePoint,
185 const Standard_Real theDiffTol)
187 const Standard_Real aSqrDist01 =
188 theParam0.Value().SquareDistance(theParam1.Value());
190 if (aSqrDist01 <= theDiffTol) {
191 // The points are confused. Get the first point and change its type.
194 const Standard_Real aDiffDist =
195 Abs(theParam0.GetSqrDistance() - theParam1.GetSqrDistance());
197 if (aDiffDist >= aSqrDist01 - theDiffTol) {
198 // The shortest distance is one of the nodes.
199 if (theParam0.GetSqrDistance() > theParam1.GetSqrDistance()) {
200 // The shortest distance is the point 1.
203 // The shortest distance is the point 0.
207 // The shortest distance is inside the edge.
208 gp_XYZ aPoP(thePoint.XYZ().Subtracted(theParam0.Value().XYZ()));
209 gp_XYZ aPoP1(theParam1.Value().XYZ().Subtracted(theParam0.Value().XYZ()));
210 Standard_Real aRatio = aPoP.Dot(aPoP1)/aSqrDist01;
214 theParam0.Parameter(aU[0], aV[0]);
215 theParam1.Parameter(aU[1], aV[1]);
217 Standard_Real aUPar = aU[0];
218 Standard_Real aVPar = aV[0];
221 aUPar += aRatio*(aU[1] - aU[0]);
223 aVPar += aRatio*(aV[1] - aV[0]);
226 Extrema_POnSurfParams aParam(aUPar, aVPar, theSurf->Value(aUPar, aVPar));
227 Standard_Integer anIndices[2];
229 theParam0.GetIndices(anIndices[0], anIndices[1]);
230 aParam.SetElementType(IsUEdge ? Extrema_UIsoEdge : Extrema_VIsoEdge);
231 aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value()));
232 aParam.SetIndices(anIndices[0], anIndices[1]);
239 //=============================================================================
241 /*-----------------------------------------------------------------------------
243 Find all extremum distances between point P and surface
244 S using sampling (NbU,NbV).
247 The algorithm bases on the hypothesis that sampling is precise enough,
248 if there exist N extreme distances between the point and the surface,
249 so there also exist N extrema between the point and the grid.
250 So, the algorithm consists in starting from extrema of the grid to find the
251 extrema of the surface.
252 The extrema are calculated by the algorithm math_FunctionSetRoot with the
254 - F: Extrema_FuncExtPS created from P and S,
255 - UV: math_Vector the components which of are parameters of the extremum on the
257 - Tol: Min(TolU,TolV), (Prov.:math_FunctionSetRoot does not autorize a vector)
258 - UVinf: math_Vector the components which of are lower limits of u and v,
259 - UVsup: math_Vector the components which of are upper limits of u and v.
262 a- Creation of the table of distances (TbDist(0,NbU+1,0,NbV+1)):
263 The table is expanded at will; lines 0 and NbU+1 and
264 columns 0 and NbV+1 are initialized at RealFirst() or RealLast()
265 to simplify the tests carried out at stage b
266 (there is no need to test if the point is on border of the grid).
267 b- Calculation of extrema:
268 First the minimums and then the maximums are found. These 2 procedured
269 pass in a similar way:
271 - 'borders' of table TbDist (RealLast() in case of minimums
272 and RealLast() in case of maximums),
273 - table TbSel(0,NbU+1,0,NbV+1) of selection of points for
274 calculation of local extremum (0). When a point will selected,
275 it will not be selectable, as well as the ajacent points
276 (8 at least). The corresponding addresses will be set to 1.
277 b.b- Calculation of minimums (or maximums):
278 All distances from table TbDist are parsed in a loop:
279 - search minimum (or maximum) in the grid,
280 - calculate extremum on the surface,
281 - update table TbSel.
282 -----------------------------------------------------------------------------*/
284 static Standard_Boolean IsoIsDeg (const Adaptor3d_Surface& S,
285 const Standard_Real Param,
286 const GeomAbs_IsoType IT,
287 const Standard_Real TolMin,
288 const Standard_Real TolMax)
290 Standard_Real U1=0.,U2=0.,V1=0.,V2=0.,T;
291 Standard_Boolean Along = Standard_True;
292 U1 = S.FirstUParameter();
293 U2 = S.LastUParameter();
294 V1 = S.FirstVParameter();
295 V2 = S.LastVParameter();
298 Standard_Real Step,D1NormMax;
299 if (IT == GeomAbs_IsoV)
303 for (T=U1;T<=U2;T=T+Step)
305 S.D1(T,Param,P,D1U,D1V);
306 D1NormMax=Max(D1NormMax,D1U.Magnitude());
309 if (D1NormMax >TolMax || D1NormMax < TolMin )
310 Along = Standard_False;
316 for (T=V1;T<=V2;T=T+Step)
318 S.D1(Param,T,P,D1U,D1V);
319 D1NormMax=Max(D1NormMax,D1V.Magnitude());
322 if (D1NormMax >TolMax || D1NormMax < TolMin )
323 Along = Standard_False;
329 //----------------------------------------------------------
330 Extrema_GenExtPS::Extrema_GenExtPS()
332 myDone = Standard_False;
333 myInit = Standard_False;
334 myFlag = Extrema_ExtFlag_MINMAX;
335 myAlgo = Extrema_ExtAlgo_Grad;
339 Extrema_GenExtPS::Extrema_GenExtPS (const gp_Pnt& P,
340 const Adaptor3d_Surface& S,
341 const Standard_Integer NbU,
342 const Standard_Integer NbV,
343 const Standard_Real TolU,
344 const Standard_Real TolV,
345 const Extrema_ExtFlag F,
346 const Extrema_ExtAlgo A)
347 : myF (P,S), myFlag(F), myAlgo(A)
349 Initialize(S, NbU, NbV, TolU, TolV);
353 Extrema_GenExtPS::Extrema_GenExtPS (const gp_Pnt& P,
354 const Adaptor3d_Surface& S,
355 const Standard_Integer NbU,
356 const Standard_Integer NbV,
357 const Standard_Real Umin,
358 const Standard_Real Usup,
359 const Standard_Real Vmin,
360 const Standard_Real Vsup,
361 const Standard_Real TolU,
362 const Standard_Real TolV,
363 const Extrema_ExtFlag F,
364 const Extrema_ExtAlgo A)
365 : myF (P,S), myFlag(F), myAlgo(A)
367 Initialize(S, NbU, NbV, Umin, Usup, Vmin, Vsup, TolU, TolV);
372 void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S,
373 const Standard_Integer NbU,
374 const Standard_Integer NbV,
375 const Standard_Real TolU,
376 const Standard_Real TolV)
378 myumin = S.FirstUParameter();
379 myusup = S.LastUParameter();
380 myvmin = S.FirstVParameter();
381 myvsup = S.LastVParameter();
382 Initialize(S,NbU,NbV,myumin,myusup,myvmin,myvsup,TolU,TolV);
386 void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S,
387 const Standard_Integer NbU,
388 const Standard_Integer NbV,
389 const Standard_Real Umin,
390 const Standard_Real Usup,
391 const Standard_Real Vmin,
392 const Standard_Real Vsup,
393 const Standard_Real TolU,
394 const Standard_Real TolV)
396 myS = (Adaptor3d_SurfacePtr)&S;
406 if ((myusample < 2) ||
407 (myvsample < 2)) { Standard_OutOfRange::Raise(); }
411 mySphereUBTree.Nullify();
414 myInit = Standard_False;
417 inline static void fillParams(const TColStd_Array1OfReal& theKnots,
418 Standard_Integer theDegree,
419 Standard_Real theParMin,
420 Standard_Real theParMax,
421 Handle_TColStd_HArray1OfReal& theParams,
422 Standard_Integer theSample)
424 NCollection_Vector<Standard_Real> aParams;
425 Standard_Integer i = 1;
426 Standard_Real aPrevPar = theParMin;
427 aParams.Append(aPrevPar);
428 //calculation the array of parametric points depending on the knots array variation and degree of given surface
429 for ( ; i < theKnots.Length() && theKnots(i) < (theParMax - Precision::PConfusion()); i++ )
431 if ( theKnots(i+1) < theParMin + Precision::PConfusion())
434 Standard_Real aStep = (theKnots(i+1) - theKnots(i))/Max(theDegree,2);
435 Standard_Integer k =1;
436 for( ; k <= theDegree ; k++)
438 Standard_Real aPar = theKnots(i) + k * aStep;
439 if(aPar > theParMax - Precision::PConfusion())
441 if(aPar > aPrevPar + Precision::PConfusion() )
443 aParams.Append(aPar);
449 aParams.Append(theParMax);
450 Standard_Integer nbPar = aParams.Length();
451 //in case of an insufficient number of points the grid will be built later
452 if (nbPar < theSample)
454 theParams = new TColStd_HArray1OfReal(1, nbPar );
455 for( i = 0; i < nbPar; i++)
456 theParams->SetValue(i+1,aParams(i));
459 void Extrema_GenExtPS::GetGridPoints( const Adaptor3d_Surface& theSurf)
461 //creation parametric points for BSpline and Bezier surfaces
462 //with taking into account of Degree and NbKnots of BSpline or Bezier geometry
463 if(theSurf.GetType() == GeomAbs_OffsetSurface)
464 GetGridPoints(theSurf.BasisSurface()->Surface());
465 //parametric points for BSpline surfaces
466 else if( theSurf.GetType() == GeomAbs_BSplineSurface)
468 Handle(Geom_BSplineSurface) aBspl = theSurf.BSpline();
471 TColStd_Array1OfReal aUKnots(1, aBspl->NbUKnots());
472 aBspl->UKnots( aUKnots);
473 TColStd_Array1OfReal aVKnots(1, aBspl->NbVKnots());
474 aBspl->VKnots( aVKnots);
475 fillParams(aUKnots,aBspl->UDegree(),myumin, myusup, myUParams, myusample);
476 fillParams(aVKnots,aBspl->VDegree(),myvmin, myvsup, myVParams, myvsample);
479 //calculation parametric points for Bezier surfaces
480 else if(theSurf.GetType() == GeomAbs_BezierSurface)
482 Handle(Geom_BezierSurface) aBezier = theSurf.Bezier();
486 TColStd_Array1OfReal aUKnots(1,2);
487 TColStd_Array1OfReal aVKnots(1,2);
488 aBezier->Bounds(aUKnots(1), aUKnots(2), aVKnots(1), aVKnots(2));
489 fillParams(aUKnots, aBezier->UDegree(), myumin, myusup, myUParams, myusample);
490 fillParams(aVKnots, aBezier->VDegree(), myvmin, myvsup, myVParams, myvsample);
493 //creation points for surfaces based on BSpline or Bezier curves
494 else if(theSurf.GetType() == GeomAbs_SurfaceOfRevolution ||
495 theSurf.GetType() == GeomAbs_SurfaceOfExtrusion)
497 Handle(TColStd_HArray1OfReal) anArrKnots;
498 Standard_Integer aDegree = 0;
499 if(theSurf.BasisCurve()->Curve().GetType() == GeomAbs_BSplineCurve)
501 Handle(Geom_BSplineCurve) aBspl = theSurf.BasisCurve()->Curve().BSpline();
504 anArrKnots = new TColStd_HArray1OfReal(1,aBspl->NbKnots());
505 aBspl->Knots( anArrKnots->ChangeArray1() );
506 aDegree = aBspl->Degree();
511 if(theSurf.BasisCurve()->Curve().GetType() == GeomAbs_BezierCurve)
513 Handle(Geom_BezierCurve) aBez = theSurf.BasisCurve()->Curve().Bezier();
516 anArrKnots = new TColStd_HArray1OfReal(1,2);
517 anArrKnots->SetValue(1, aBez->FirstParameter());
518 anArrKnots->SetValue(2, aBez->LastParameter());
519 aDegree = aBez->Degree();
523 if(anArrKnots.IsNull())
525 if( theSurf.GetType() == GeomAbs_SurfaceOfRevolution )
526 fillParams( anArrKnots->Array1(), aDegree, myvmin, myvsup, myVParams, myvsample);
528 fillParams( anArrKnots->Array1(), aDegree, myumin, myusup, myUParams, myusample);
530 //update the number of points in sample
531 if(!myUParams.IsNull())
532 myusample = myUParams->Length();
533 if( !myVParams.IsNull())
534 myvsample = myVParams->Length();
538 void Extrema_GenExtPS::BuildGrid(const gp_Pnt &thePoint)
540 Standard_Integer NoU, NoV;
542 //if grid was already built skip its creation
544 //build parametric grid in case of a complex surface geometry (BSpline and Bezier surfaces)
547 //build grid in other cases
548 if( myUParams.IsNull() )
550 Standard_Real PasU = myusup - myumin;
551 Standard_Real U0 = PasU / myusample / 100.;
552 PasU = (PasU - U0) / (myusample - 1);
554 myUParams = new TColStd_HArray1OfReal(1,myusample );
555 Standard_Integer NoU;
556 Standard_Real U = U0;
557 for ( NoU = 1 ; NoU <= myusample; NoU++, U += PasU)
558 myUParams->SetValue(NoU, U);
561 if( myVParams.IsNull())
563 Standard_Real PasV = myvsup - myvmin;
564 Standard_Real V0 = PasV / myvsample / 100.;
565 PasV = (PasV - V0) / (myvsample - 1);
568 myVParams = new TColStd_HArray1OfReal(1,myvsample );
569 Standard_Integer NoV;
570 Standard_Real V = V0;
572 for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV)
573 myVParams->SetValue(NoV, V);
576 //If flag was changed and extrema not reinitialized Extrema would fail
577 myPoints = new Extrema_HArray2OfPOnSurfParams
578 (0, myusample + 1, 0, myvsample + 1);
579 // Calculation of distances
581 for ( NoU = 1 ; NoU <= myusample; NoU++ ) {
582 for ( NoV = 1 ; NoV <= myvsample; NoV++) {
583 gp_Pnt aP1 = myS->Value(myUParams->Value(NoU), myVParams->Value(NoV));
584 Extrema_POnSurfParams aParam
585 (myUParams->Value(NoU), myVParams->Value(NoV), aP1);
587 aParam.SetElementType(Extrema_Node);
588 aParam.SetIndices(NoU, NoV);
589 myPoints->SetValue(NoU, NoV, aParam);
593 // Fill boundary with negative square distance.
594 // It is used for computation of Maximum.
595 for (NoV = 0; NoV <= myvsample + 1; NoV++) {
596 myPoints->ChangeValue(0, NoV).SetSqrDistance(-1.);
597 myPoints->ChangeValue(myusample + 1, NoV).SetSqrDistance(-1.);
600 for (NoU = 1; NoU <= myusample; NoU++) {
601 myPoints->ChangeValue(NoU, 0).SetSqrDistance(-1.);
602 myPoints->ChangeValue(NoU, myvsample + 1).SetSqrDistance(-1.);
605 myInit = Standard_True;
608 // Compute distances to mesh.
609 // Step 1. Compute distances to nodes.
610 for ( NoU = 1 ; NoU <= myusample; NoU++ ) {
611 for ( NoV = 1 ; NoV <= myvsample; NoV++) {
612 Extrema_POnSurfParams &aParam = myPoints->ChangeValue(NoU, NoV);
614 aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value()));
618 // For search of minimum compute distances to mesh.
619 if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
621 // This is the tolerance of difference of squared values.
622 // No need to set it too small.
623 const Standard_Real aDiffTol = mytolu + mytolv;
625 // Step 2. Compute distances to edges.
626 // Assume UEdge(i, j) = { Point(i, j); Point(i + 1, j ) }
627 // Assume VEdge(i, j) = { Point(i, j); Point(i, j + 1) }
628 Handle(Extrema_HArray2OfPOnSurfParams) aUEdgePntParams =
629 new Extrema_HArray2OfPOnSurfParams(1, myusample - 1, 1, myvsample);
630 Handle(Extrema_HArray2OfPOnSurfParams) aVEdgePntParams =
631 new Extrema_HArray2OfPOnSurfParams(1, myusample, 1, myvsample - 1);
633 for ( NoU = 1 ; NoU <= myusample; NoU++ ) {
634 for ( NoV = 1 ; NoV <= myvsample; NoV++) {
635 const Extrema_POnSurfParams &aParam0 = myPoints->Value(NoU, NoV);
637 if (NoU < myusample) {
638 // Compute parameters to UEdge.
639 const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU + 1, NoV);
640 Extrema_POnSurfParams aUEdgeParam = ComputeEdgeParameters
641 (Standard_True, aParam0, aParam1, myS, thePoint, aDiffTol);
643 aUEdgePntParams->SetValue(NoU, NoV, aUEdgeParam);
646 if (NoV < myvsample) {
647 // Compute parameters to VEdge.
648 const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU, NoV + 1);
649 Extrema_POnSurfParams aVEdgeParam = ComputeEdgeParameters
650 (Standard_False, aParam0, aParam1, myS, thePoint, aDiffTol);
652 aVEdgePntParams->SetValue(NoU, NoV, aVEdgeParam);
657 // Step 3. Compute distances to faces.
658 // Assume myFacePntParams(i, j) =
659 // { Point(i, j); Point(i + 1, j); Point(i + 1, j + 1); Point(i, j + 1) }
661 // { UEdge(i, j); VEdge(i + 1, j); UEdge(i, j + 1); VEdge(i, j) }
663 new Extrema_HArray2OfPOnSurfParams(0, myusample, 0, myvsample);
664 Standard_Real aSqrDist01;
665 Standard_Real aDiffDist;
666 Standard_Boolean isOut;
668 for ( NoU = 1 ; NoU < myusample; NoU++ ) {
669 for ( NoV = 1 ; NoV < myvsample; NoV++) {
670 const Extrema_POnSurfParams &aUE0 = aUEdgePntParams->Value(NoU, NoV);
671 const Extrema_POnSurfParams &aUE1 = aUEdgePntParams->Value(NoU, NoV+1);
672 const Extrema_POnSurfParams &aVE0 = aVEdgePntParams->Value(NoU, NoV);
673 const Extrema_POnSurfParams &aVE1 = aVEdgePntParams->Value(NoU+1, NoV);
675 aSqrDist01 = aUE0.Value().SquareDistance(aUE1.Value());
676 aDiffDist = Abs(aUE0.GetSqrDistance() - aUE1.GetSqrDistance());
677 isOut = Standard_False;
679 if (aDiffDist >= aSqrDist01 - aDiffTol) {
680 // The projection is outside the face.
681 isOut = Standard_True;
683 aSqrDist01 = aVE0.Value().SquareDistance(aVE1.Value());
684 aDiffDist = Abs(aVE0.GetSqrDistance() - aVE1.GetSqrDistance());
686 if (aDiffDist >= aSqrDist01 - aDiffTol) {
687 // The projection is outside the face.
688 isOut = Standard_True;
693 // Get the closest point on an edge.
694 const Extrema_POnSurfParams &aUEMin =
695 aUE0.GetSqrDistance() < aUE1.GetSqrDistance() ? aUE0 : aUE1;
696 const Extrema_POnSurfParams &aVEMin =
697 aVE0.GetSqrDistance() < aVE1.GetSqrDistance() ? aVE0 : aVE1;
698 const Extrema_POnSurfParams &aEMin =
699 aUEMin.GetSqrDistance() < aVEMin.GetSqrDistance() ? aUEMin : aVEMin;
701 myFacePntParams->SetValue(NoU, NoV, aEMin);
703 // Find closest point inside the face.
709 // Compute U parameter.
710 aUE0.Parameter(aU[0], aV[0]);
711 aUE1.Parameter(aU[1], aV[1]);
712 aUPar = 0.5*(aU[0] + aU[1]);
713 // Compute V parameter.
714 aVE0.Parameter(aU[0], aV[0]);
715 aVE1.Parameter(aU[1], aV[1]);
716 aVPar = 0.5*(aV[0] + aV[1]);
718 Extrema_POnSurfParams aParam(aUPar, aVPar, myS->Value(aUPar, aVPar));
720 aParam.SetElementType(Extrema_Face);
721 aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value()));
722 aParam.SetIndices(NoU, NoV);
723 myFacePntParams->SetValue(NoU, NoV, aParam);
728 // Fill boundary with RealLast square distance.
729 for (NoV = 0; NoV <= myvsample; NoV++) {
730 myFacePntParams->ChangeValue(0, NoV).SetSqrDistance(RealLast());
731 myFacePntParams->ChangeValue(myusample, NoV).SetSqrDistance(RealLast());
734 for (NoU = 1; NoU < myusample; NoU++) {
735 myFacePntParams->ChangeValue(NoU, 0).SetSqrDistance(RealLast());
736 myFacePntParams->ChangeValue(NoU, myvsample).SetSqrDistance(RealLast());
742 a- Constitution of the table of distances (TbDist(0,myusample+1,0,myvsample+1)):
743 ---------------------------------------------------------------
746 // Parameterisation of the sample
748 void Extrema_GenExtPS::BuildTree()
750 // if tree already exists, assume it is already correctly filled
751 if ( ! mySphereUBTree.IsNull() )
754 Standard_Real PasU = myusup - myumin;
755 Standard_Real PasV = myvsup - myvmin;
756 Standard_Real U0 = PasU / myusample / 100.;
757 Standard_Real V0 = PasV / myvsample / 100.;
759 PasU = (PasU - U0) / (myusample - 1);
760 PasV = (PasV - V0) / (myvsample - 1);
764 //build grid of parametric points
765 myUParams = new TColStd_HArray1OfReal(1,myusample );
766 myVParams = new TColStd_HArray1OfReal(1,myvsample );
767 Standard_Integer NoU, NoV;
768 Standard_Real U = U0, V = V0;
769 for ( NoU = 1 ; NoU <= myusample; NoU++, U += PasU)
770 myUParams->SetValue(NoU, U);
771 for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV)
772 myVParams->SetValue(NoV, V);
774 // Calculation of distances
775 mySphereUBTree = new Extrema_UBTreeOfSphere;
776 Extrema_UBTreeFillerOfSphere aFiller(*mySphereUBTree);
777 Standard_Integer i = 0;
779 mySphereArray = new Bnd_HArray1OfSphere(0, myusample * myvsample);
781 for ( NoU = 1; NoU <= myusample; NoU++ ) {
782 for ( NoV = 1; NoV <= myvsample; NoV++) {
783 P1 = myS->Value(myUParams->Value(NoU), myVParams->Value(NoV));
784 Bnd_Sphere aSph(P1.XYZ(), 0/*mytolu < mytolv ? mytolu : mytolv*/, NoU, NoV);
785 aFiller.Add(i, aSph);
786 mySphereArray->SetValue( i, aSph );
793 void Extrema_GenExtPS::FindSolution(const gp_Pnt& /*P*/,
794 const Extrema_POnSurfParams &theParams)
796 math_Vector Tol(1,2);
800 math_Vector UV(1, 2);
802 theParams.Parameter(UV(1), UV(2));
804 math_Vector UVinf(1,2), UVsup(1,2);
810 math_Vector errors(1,2);
811 math_Vector root(1, 2);
813 Standard_Integer aNbMaxIter = 100;
815 gp_Pnt PStart = theParams.Value();
817 math_FunctionSetRoot S (myF,UV,Tol,UVinf,UVsup, aNbMaxIter);
819 myDone = Standard_True;
822 void Extrema_GenExtPS::SetFlag(const Extrema_ExtFlag F)
827 void Extrema_GenExtPS::SetAlgo(const Extrema_ExtAlgo A)
830 myInit = Standard_False;
835 void Extrema_GenExtPS::Perform(const gp_Pnt& P)
837 myDone = Standard_False;
840 if(myAlgo == Extrema_ExtAlgo_Grad)
843 Standard_Integer NoU,NoV;
845 if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
847 Extrema_ElementType anElemType;
850 Standard_Integer iU2;
851 Standard_Integer iV2;
852 Standard_Boolean isMin;
855 for (NoU = 1; NoU < myusample; NoU++) {
856 for (NoV = 1; NoV < myvsample; NoV++) {
857 const Extrema_POnSurfParams &aParam =
858 myFacePntParams->Value(NoU, NoV);
860 isMin = Standard_False;
861 anElemType = aParam.GetElementType();
863 if (anElemType == Extrema_Face) {
864 isMin = Standard_True;
866 // Check if it is a boundary edge or corner vertex.
867 aParam.GetIndices(iU, iV);
869 if (anElemType == Extrema_UIsoEdge) {
870 isMin = (iV == 1 || iV == myvsample);
871 } else if (anElemType == Extrema_VIsoEdge) {
872 isMin = (iU == 1 || iU == myusample);
873 } else if (anElemType == Extrema_Node) {
874 isMin = (iU == 1 || iU == myusample) &&
875 (iV == 1 || iV == myvsample);
879 // This is a middle element.
880 if (anElemType == Extrema_UIsoEdge ||
881 (anElemType == Extrema_Node && (iU == 1 || iU == myusample))) {
882 // Check the down face.
883 const Extrema_POnSurfParams &aDownParam =
884 myFacePntParams->Value(NoU, NoV - 1);
886 if (aDownParam.GetElementType() == anElemType) {
887 aDownParam.GetIndices(iU2, iV2);
888 isMin = (iU == iU2 && iV == iV2);
890 } else if (anElemType == Extrema_VIsoEdge ||
891 (anElemType == Extrema_Node && (iV == 1 || iV == myvsample))) {
892 // Check the right face.
893 const Extrema_POnSurfParams &aRightParam =
894 myFacePntParams->Value(NoU - 1, NoV);
896 if (aRightParam.GetElementType() == anElemType) {
897 aRightParam.GetIndices(iU2, iV2);
898 isMin = (iU == iU2 && iV == iV2);
900 } else if (iU == NoU && iV == NoV) {
901 // Check the lower-left node. For this purpose it is necessary
902 // to check lower-left, lower and left faces.
903 isMin = Standard_True;
905 const Extrema_POnSurfParams *anOtherParam[3] =
906 { &myFacePntParams->Value(NoU, NoV - 1), // Down
907 &myFacePntParams->Value(NoU - 1, NoV - 1), // Lower-left
908 &myFacePntParams->Value(NoU - 1, NoV) }; // Left
910 for (i = 0; i < 3 && isMin; i++) {
911 if (anOtherParam[i]->GetElementType() == Extrema_Node) {
912 anOtherParam[i]->GetIndices(iU2, iV2);
913 isMin = (iU == iU2 && iV == iV2);
915 isMin = Standard_False;
923 FindSolution(P, aParam);
929 if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX)
933 for (NoU = 1; NoU <= myusample; NoU++) {
934 for (NoV = 1; NoV <= myvsample; NoV++) {
935 Dist = myPoints->Value(NoU, NoV).GetSqrDistance();
936 if ((myPoints->Value(NoU-1,NoV-1).GetSqrDistance() <= Dist) &&
937 (myPoints->Value(NoU-1,NoV ).GetSqrDistance() <= Dist) &&
938 (myPoints->Value(NoU-1,NoV+1).GetSqrDistance() <= Dist) &&
939 (myPoints->Value(NoU ,NoV-1).GetSqrDistance() <= Dist) &&
940 (myPoints->Value(NoU ,NoV+1).GetSqrDistance() <= Dist) &&
941 (myPoints->Value(NoU+1,NoV-1).GetSqrDistance() <= Dist) &&
942 (myPoints->Value(NoU+1,NoV ).GetSqrDistance() <= Dist) &&
943 (myPoints->Value(NoU+1,NoV+1).GetSqrDistance() <= Dist)) {
945 FindSolution(P, myPoints->Value(NoU, NoV));
954 if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
956 Bnd_Sphere aSol = mySphereArray->Value(0);
957 Bnd_SphereUBTreeSelectorMin aSelector(mySphereArray, aSol);
958 //aSelector.SetMaxDist( RealLast() );
959 aSelector.DefineCheckPoint( P );
960 mySphereUBTree->Select( aSelector );
961 //TODO: check if no solution in binary tree
962 Bnd_Sphere& aSph = aSelector.Sphere();
963 Standard_Real aU = myUParams->Value(aSph.U());
964 Standard_Real aV = myVParams->Value(aSph.V());
965 Extrema_POnSurfParams aParams(aU, aV, myS->Value(aU, aV));
967 aParams.SetSqrDistance(P.SquareDistance(aParams.Value()));
968 aParams.SetIndices(aSph.U(), aSph.V());
969 FindSolution(P, aParams);
971 if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX)
973 Bnd_Sphere aSol = mySphereArray->Value(0);
974 Bnd_SphereUBTreeSelectorMax aSelector(mySphereArray, aSol);
975 //aSelector.SetMaxDist( RealLast() );
976 aSelector.DefineCheckPoint( P );
977 mySphereUBTree->Select( aSelector );
978 //TODO: check if no solution in binary tree
979 Bnd_Sphere& aSph = aSelector.Sphere();
980 Standard_Real aU = myUParams->Value(aSph.U());
981 Standard_Real aV = myVParams->Value(aSph.V());
982 Extrema_POnSurfParams aParams(aU, aV, myS->Value(aU, aV));
983 aParams.SetSqrDistance(P.SquareDistance(aParams.Value()));
984 aParams.SetIndices(aSph.U(), aSph.V());
986 FindSolution(P, aParams);
990 //=============================================================================
992 Standard_Boolean Extrema_GenExtPS::IsDone () const { return myDone; }
993 //=============================================================================
995 Standard_Integer Extrema_GenExtPS::NbExt () const
997 if (!IsDone()) { StdFail_NotDone::Raise(); }
1000 //=============================================================================
1002 Standard_Real Extrema_GenExtPS::SquareDistance (const Standard_Integer N) const
1004 if (!IsDone()) { StdFail_NotDone::Raise(); }
1005 return myF.SquareDistance(N);
1007 //=============================================================================
1009 Extrema_POnSurf Extrema_GenExtPS::Point (const Standard_Integer N) const
1011 if (!IsDone()) { StdFail_NotDone::Raise(); }
1012 return myF.Point(N);
1014 //=============================================================================