1 // Created on: 1995-07-18
2 // Created by: Modelistation
3 // Copyright (c) 1995-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 // Modified by skv - Thu Sep 30 15:21:07 2004 OCC593
19 #include <Adaptor3d_Curve.hxx>
20 #include <Adaptor3d_HCurve.hxx>
21 #include <Adaptor3d_HSurface.hxx>
22 #include <Adaptor3d_Surface.hxx>
23 #include <Bnd_Array1OfSphere.hxx>
24 #include <Bnd_HArray1OfSphere.hxx>
25 #include <Bnd_Sphere.hxx>
26 #include <Extrema_ExtFlag.hxx>
27 #include <Extrema_GenExtPS.hxx>
28 #include <Extrema_HUBTreeOfSphere.hxx>
29 #include <Extrema_POnSurf.hxx>
30 #include <Extrema_POnSurfParams.hxx>
31 #include <Geom_BezierCurve.hxx>
32 #include <Geom_BezierSurface.hxx>
33 #include <Geom_BSplineCurve.hxx>
34 #include <Geom_BSplineSurface.hxx>
35 #include <Geom_OffsetSurface.hxx>
36 #include <Geom_RectangularTrimmedSurface.hxx>
37 #include <GeomAbs_IsoType.hxx>
39 #include <math_FunctionSetRoot.hxx>
40 #include <math_NewtonFunctionSetRoot.hxx>
41 #include <math_Vector.hxx>
42 #include <Precision.hxx>
43 #include <Standard_OutOfRange.hxx>
44 #include <Standard_TypeMismatch.hxx>
45 #include <StdFail_NotDone.hxx>
46 #include <TColgp_Array2OfPnt.hxx>
47 #include <TColStd_Array2OfInteger.hxx>
48 #include <TColStd_Array2OfReal.hxx>
50 //IMPLEMENT_HARRAY1(Extrema_HArray1OfSphere)
51 class Bnd_SphereUBTreeSelector : public Extrema_UBTreeOfSphere::Selector
55 Bnd_SphereUBTreeSelector (const Handle(Bnd_HArray1OfSphere)& theSphereArray,
58 mySphereArray(theSphereArray),
63 void DefineCheckPoint( const gp_Pnt& theXYZ )
66 Bnd_Sphere& Sphere() const
69 virtual Standard_Boolean Reject( const Bnd_Sphere &theBnd ) const = 0;
71 virtual Standard_Boolean Accept(const Standard_Integer& theObj) = 0;
74 const Handle(Bnd_HArray1OfSphere)& mySphereArray;
77 void operator= (const Bnd_SphereUBTreeSelector&);
81 class Bnd_SphereUBTreeSelectorMin : public Bnd_SphereUBTreeSelector
84 Bnd_SphereUBTreeSelectorMin (const Handle(Bnd_HArray1OfSphere)& theSphereArray,
86 : Bnd_SphereUBTreeSelector(theSphereArray, theSol),
90 void SetMinDist( const Standard_Real theMinDist )
91 { myMinDist = theMinDist; }
93 Standard_Real MinDist() const
96 Standard_Boolean Reject( const Bnd_Sphere &theBnd ) const
98 Bnd_SphereUBTreeSelectorMin* me =
99 const_cast<Bnd_SphereUBTreeSelectorMin*>(this);
100 // myMinDist is decreased each time a nearer object is found
101 return theBnd.IsOut( myXYZ.XYZ(), me->myMinDist );
104 Standard_Boolean Accept(const Standard_Integer&);
107 Standard_Real myMinDist;
110 Standard_Boolean Bnd_SphereUBTreeSelectorMin::Accept(const Standard_Integer& theInd)
112 const Bnd_Sphere& aSph = mySphereArray->Value(theInd);
113 Standard_Real aCurDist;
115 // if ( (aCurDist = aSph.SquareDistance(myXYZ.XYZ())) < mySol.SquareDistance(myXYZ.XYZ()) )
116 if ( (aCurDist = aSph.Distance(myXYZ.XYZ())) < mySol.Distance(myXYZ.XYZ()) )
119 if ( aCurDist < myMinDist )
120 myMinDist = aCurDist;
122 return Standard_True;
125 return Standard_False;
128 class Bnd_SphereUBTreeSelectorMax : public Bnd_SphereUBTreeSelector
131 Bnd_SphereUBTreeSelectorMax (const Handle(Bnd_HArray1OfSphere)& theSphereArray,
133 : Bnd_SphereUBTreeSelector(theSphereArray, theSol),
137 void SetMaxDist( const Standard_Real theMaxDist )
138 { myMaxDist = theMaxDist; }
140 Standard_Real MaxDist() const
141 { return myMaxDist; }
143 Standard_Boolean Reject( const Bnd_Sphere &theBnd ) const
145 Bnd_SphereUBTreeSelectorMax* me =
146 const_cast<Bnd_SphereUBTreeSelectorMax*>(this);
147 // myMaxDist is decreased each time a nearer object is found
148 return theBnd.IsOut( myXYZ.XYZ(), me->myMaxDist );
151 Standard_Boolean Accept(const Standard_Integer&);
154 Standard_Real myMaxDist;
157 Standard_Boolean Bnd_SphereUBTreeSelectorMax::Accept(const Standard_Integer& theInd)
159 const Bnd_Sphere& aSph = mySphereArray->Value(theInd);
160 Standard_Real aCurDist;
162 // if ( (aCurDist = aSph.SquareDistance(myXYZ.XYZ())) > mySol.SquareDistance(myXYZ.XYZ()) )
163 if ( (aCurDist = aSph.Distance(myXYZ.XYZ())) > mySol.Distance(myXYZ.XYZ()) )
166 if ( aCurDist > myMaxDist )
167 myMaxDist = aCurDist;
169 return Standard_True;
172 return Standard_False;
175 //=============================================================================
177 /*-----------------------------------------------------------------------------
179 Find all extremum distances between point P and surface
180 S using sampling (NbU,NbV).
183 The algorithm bases on the hypothesis that sampling is precise enough,
184 if there exist N extreme distances between the point and the surface,
185 so there also exist N extrema between the point and the grid.
186 So, the algorithm consists in starting from extrema of the grid to find the
187 extrema of the surface.
188 The extrema are calculated by the algorithm math_FunctionSetRoot with the
190 - F: Extrema_FuncExtPS created from P and S,
191 - UV: math_Vector the components which of are parameters of the extremum on the
193 - Tol: Min(TolU,TolV), (Prov.:math_FunctionSetRoot does not autorize a vector)
194 - UVinf: math_Vector the components which of are lower limits of u and v,
195 - UVsup: math_Vector the components which of are upper limits of u and v.
198 a- Creation of the table of distances (TbDist(0,NbU+1,0,NbV+1)):
199 The table is expanded at will; lines 0 and NbU+1 and
200 columns 0 and NbV+1 are initialized at RealFirst() or RealLast()
201 to simplify the tests carried out at stage b
202 (there is no need to test if the point is on border of the grid).
203 b- Calculation of extrema:
204 First the minimums and then the maximums are found. These 2 procedured
205 pass in a similar way:
207 - 'borders' of table TbDist (RealLast() in case of minimums
208 and RealLast() in case of maximums),
209 - table TbSel(0,NbU+1,0,NbV+1) of selection of points for
210 calculation of local extremum (0). When a point will selected,
211 it will not be selectable, as well as the ajacent points
212 (8 at least). The corresponding addresses will be set to 1.
213 b.b- Calculation of minimums (or maximums):
214 All distances from table TbDist are parsed in a loop:
215 - search minimum (or maximum) in the grid,
216 - calculate extremum on the surface,
217 - update table TbSel.
218 -----------------------------------------------------------------------------*/
220 Extrema_GenExtPS::Extrema_GenExtPS()
222 myDone = Standard_False;
223 myInit = Standard_False;
224 myFlag = Extrema_ExtFlag_MINMAX;
225 myAlgo = Extrema_ExtAlgo_Grad;
229 Extrema_GenExtPS::Extrema_GenExtPS (const gp_Pnt& P,
230 const Adaptor3d_Surface& S,
231 const Standard_Integer NbU,
232 const Standard_Integer NbV,
233 const Standard_Real TolU,
234 const Standard_Real TolV,
235 const Extrema_ExtFlag F,
236 const Extrema_ExtAlgo A)
241 Initialize(S, NbU, NbV, TolU, TolV);
245 Extrema_GenExtPS::Extrema_GenExtPS (const gp_Pnt& P,
246 const Adaptor3d_Surface& S,
247 const Standard_Integer NbU,
248 const Standard_Integer NbV,
249 const Standard_Real Umin,
250 const Standard_Real Usup,
251 const Standard_Real Vmin,
252 const Standard_Real Vsup,
253 const Standard_Real TolU,
254 const Standard_Real TolV,
255 const Extrema_ExtFlag F,
256 const Extrema_ExtAlgo A)
261 Initialize(S, NbU, NbV, Umin, Usup, Vmin, Vsup, TolU, TolV);
266 void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S,
267 const Standard_Integer NbU,
268 const Standard_Integer NbV,
269 const Standard_Real TolU,
270 const Standard_Real TolV)
272 myumin = S.FirstUParameter();
273 myusup = S.LastUParameter();
274 myvmin = S.FirstVParameter();
275 myvsup = S.LastVParameter();
276 Initialize(S,NbU,NbV,myumin,myusup,myvmin,myvsup,TolU,TolV);
280 void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S,
281 const Standard_Integer NbU,
282 const Standard_Integer NbV,
283 const Standard_Real Umin,
284 const Standard_Real Usup,
285 const Standard_Real Vmin,
286 const Standard_Real Vsup,
287 const Standard_Real TolU,
288 const Standard_Real TolV)
290 myS = (Adaptor3d_SurfacePtr)&S;
300 if ((myusample < 2) ||
301 (myvsample < 2)) { throw Standard_OutOfRange(); }
305 mySphereUBTree.Nullify();
308 myInit = Standard_False;
311 inline static void fillParams(const TColStd_Array1OfReal& theKnots,
312 Standard_Integer theDegree,
313 Standard_Real theParMin,
314 Standard_Real theParMax,
315 Handle(TColStd_HArray1OfReal)& theParams,
316 Standard_Integer theSample)
318 NCollection_Vector<Standard_Real> aParams;
319 Standard_Integer i = 1;
320 Standard_Real aPrevPar = theParMin;
321 aParams.Append(aPrevPar);
322 //calculation the array of parametric points depending on the knots array variation and degree of given surface
323 for ( ; i < theKnots.Length() && theKnots(i) < (theParMax - Precision::PConfusion()); i++ )
325 if ( theKnots(i+1) < theParMin + Precision::PConfusion())
328 Standard_Real aStep = (theKnots(i+1) - theKnots(i))/Max(theDegree,2);
329 Standard_Integer k =1;
330 for( ; k <= theDegree ; k++)
332 Standard_Real aPar = theKnots(i) + k * aStep;
333 if(aPar > theParMax - Precision::PConfusion())
335 if(aPar > aPrevPar + Precision::PConfusion() )
337 aParams.Append(aPar);
343 aParams.Append(theParMax);
344 Standard_Integer nbPar = aParams.Length();
345 //in case of an insufficient number of points the grid will be built later
346 if (nbPar < theSample)
348 theParams = new TColStd_HArray1OfReal(1, nbPar );
349 for( i = 0; i < nbPar; i++)
350 theParams->SetValue(i+1,aParams(i));
353 void Extrema_GenExtPS::GetGridPoints( const Adaptor3d_Surface& theSurf)
355 //creation parametric points for BSpline and Bezier surfaces
356 //with taking into account of Degree and NbKnots of BSpline or Bezier geometry
357 if(theSurf.GetType() == GeomAbs_OffsetSurface)
358 GetGridPoints(theSurf.BasisSurface()->Surface());
359 //parametric points for BSpline surfaces
360 else if( theSurf.GetType() == GeomAbs_BSplineSurface)
362 Handle(Geom_BSplineSurface) aBspl = theSurf.BSpline();
365 TColStd_Array1OfReal aUKnots(1, aBspl->NbUKnots());
366 aBspl->UKnots( aUKnots);
367 TColStd_Array1OfReal aVKnots(1, aBspl->NbVKnots());
368 aBspl->VKnots( aVKnots);
369 fillParams(aUKnots,aBspl->UDegree(),myumin, myusup, myUParams, myusample);
370 fillParams(aVKnots,aBspl->VDegree(),myvmin, myvsup, myVParams, myvsample);
373 //calculation parametric points for Bezier surfaces
374 else if(theSurf.GetType() == GeomAbs_BezierSurface)
376 Handle(Geom_BezierSurface) aBezier = theSurf.Bezier();
380 TColStd_Array1OfReal aUKnots(1,2);
381 TColStd_Array1OfReal aVKnots(1,2);
382 aBezier->Bounds(aUKnots(1), aUKnots(2), aVKnots(1), aVKnots(2));
383 fillParams(aUKnots, aBezier->UDegree(), myumin, myusup, myUParams, myusample);
384 fillParams(aVKnots, aBezier->VDegree(), myvmin, myvsup, myVParams, myvsample);
387 //creation points for surfaces based on BSpline or Bezier curves
388 else if(theSurf.GetType() == GeomAbs_SurfaceOfRevolution ||
389 theSurf.GetType() == GeomAbs_SurfaceOfExtrusion)
391 Handle(TColStd_HArray1OfReal) anArrKnots;
392 Standard_Integer aDegree = 0;
393 if(theSurf.BasisCurve()->Curve().GetType() == GeomAbs_BSplineCurve)
395 Handle(Geom_BSplineCurve) aBspl = theSurf.BasisCurve()->Curve().BSpline();
398 anArrKnots = new TColStd_HArray1OfReal(1,aBspl->NbKnots());
399 aBspl->Knots( anArrKnots->ChangeArray1() );
400 aDegree = aBspl->Degree();
405 if(theSurf.BasisCurve()->Curve().GetType() == GeomAbs_BezierCurve)
407 Handle(Geom_BezierCurve) aBez = theSurf.BasisCurve()->Curve().Bezier();
410 anArrKnots = new TColStd_HArray1OfReal(1,2);
411 anArrKnots->SetValue(1, aBez->FirstParameter());
412 anArrKnots->SetValue(2, aBez->LastParameter());
413 aDegree = aBez->Degree();
417 if(anArrKnots.IsNull())
419 if( theSurf.GetType() == GeomAbs_SurfaceOfRevolution )
420 fillParams( anArrKnots->Array1(), aDegree, myvmin, myvsup, myVParams, myvsample);
422 fillParams( anArrKnots->Array1(), aDegree, myumin, myusup, myUParams, myusample);
424 //update the number of points in sample
425 if(!myUParams.IsNull())
426 myusample = myUParams->Length();
427 if( !myVParams.IsNull())
428 myvsample = myVParams->Length();
433 * This function computes the point on surface parameters on edge.
434 * if it coincides with theParam0 or theParam1, it is returned.
436 const Extrema_POnSurfParams& Extrema_GenExtPS::ComputeEdgeParameters
437 (const Standard_Boolean IsUEdge,
438 const Extrema_POnSurfParams &theParam0,
439 const Extrema_POnSurfParams &theParam1,
440 const gp_Pnt &thePoint,
441 const Standard_Real theDiffTol)
443 const Standard_Real aSqrDist01 =
444 theParam0.Value().SquareDistance(theParam1.Value());
446 if (aSqrDist01 <= theDiffTol)
448 // The points are confused. Get the first point and change its type.
453 const Standard_Real aDiffDist =
454 Abs(theParam0.GetSqrDistance() - theParam1.GetSqrDistance());
456 if (aDiffDist >= aSqrDist01 - theDiffTol)
458 // The shortest distance is one of the nodes.
459 if (theParam0.GetSqrDistance() > theParam1.GetSqrDistance())
461 // The shortest distance is the point 1.
466 // The shortest distance is the point 0.
472 // The shortest distance is inside the edge.
473 gp_XYZ aPoP(thePoint.XYZ().Subtracted(theParam0.Value().XYZ()));
474 gp_XYZ aPoP1(theParam1.Value().XYZ().Subtracted(theParam0.Value().XYZ()));
475 Standard_Real aRatio = aPoP.Dot(aPoP1)/aSqrDist01;
479 theParam0.Parameter(aU[0], aV[0]);
480 theParam1.Parameter(aU[1], aV[1]);
482 Standard_Real aUPar = aU[0];
483 Standard_Real aVPar = aV[0];
487 aUPar += aRatio*(aU[1] - aU[0]);
491 aVPar += aRatio*(aV[1] - aV[0]);
494 myGridParam.SetParameters(aUPar, aVPar, myS->Value(aUPar, aVPar));
495 Standard_Integer anIndices[2];
497 theParam0.GetIndices(anIndices[0], anIndices[1]);
498 myGridParam.SetElementType(IsUEdge ? Extrema_UIsoEdge : Extrema_VIsoEdge);
499 myGridParam.SetSqrDistance(thePoint.SquareDistance(myGridParam.Value()));
500 myGridParam.SetIndices(anIndices[0], anIndices[1]);
506 void Extrema_GenExtPS::BuildGrid(const gp_Pnt &thePoint)
508 Standard_Integer NoU, NoV;
510 //if grid was already built skip its creation
512 //build parametric grid in case of a complex surface geometry (BSpline and Bezier surfaces)
515 //build grid in other cases
516 if( myUParams.IsNull() )
518 Standard_Real PasU = myusup - myumin;
519 Standard_Real U0 = PasU / myusample / 100.;
520 PasU = (PasU - U0) / (myusample - 1);
522 myUParams = new TColStd_HArray1OfReal(1,myusample );
523 Standard_Real U = U0;
524 for ( NoU = 1 ; NoU <= myusample; NoU++, U += PasU)
525 myUParams->SetValue(NoU, U);
528 if( myVParams.IsNull())
530 Standard_Real PasV = myvsup - myvmin;
531 Standard_Real V0 = PasV / myvsample / 100.;
532 PasV = (PasV - V0) / (myvsample - 1);
535 myVParams = new TColStd_HArray1OfReal(1,myvsample );
536 Standard_Real V = V0;
538 for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV)
539 myVParams->SetValue(NoV, V);
542 //If flag was changed and extrema not reinitialized Extrema would fail
543 myPoints = new Extrema_HArray2OfPOnSurfParams
544 (0, myusample + 1, 0, myvsample + 1);
545 // Calculation of distances
547 for ( NoU = 1 ; NoU <= myusample; NoU++ ) {
548 for ( NoV = 1 ; NoV <= myvsample; NoV++) {
549 gp_Pnt aP1 = myS->Value(myUParams->Value(NoU), myVParams->Value(NoV));
550 Extrema_POnSurfParams aParam
551 (myUParams->Value(NoU), myVParams->Value(NoV), aP1);
553 aParam.SetElementType(Extrema_Node);
554 aParam.SetIndices(NoU, NoV);
555 myPoints->SetValue(NoU, NoV, aParam);
560 new Extrema_HArray2OfPOnSurfParams(0, myusample, 0, myvsample);
563 new Extrema_HArray2OfPOnSurfParams(1, myusample - 1, 1, myvsample);
565 new Extrema_HArray2OfPOnSurfParams(1, myusample, 1, myvsample - 1);
567 // Fill boundary with negative square distance.
568 // It is used for computation of Maximum.
569 for (NoV = 0; NoV <= myvsample + 1; NoV++) {
570 myPoints->ChangeValue(0, NoV).SetSqrDistance(-1.);
571 myPoints->ChangeValue(myusample + 1, NoV).SetSqrDistance(-1.);
574 for (NoU = 1; NoU <= myusample; NoU++) {
575 myPoints->ChangeValue(NoU, 0).SetSqrDistance(-1.);
576 myPoints->ChangeValue(NoU, myvsample + 1).SetSqrDistance(-1.);
579 myInit = Standard_True;
582 // Compute distances to mesh.
583 // Step 1. Compute distances to nodes.
584 for ( NoU = 1 ; NoU <= myusample; NoU++ ) {
585 for ( NoV = 1 ; NoV <= myvsample; NoV++) {
586 Extrema_POnSurfParams &aParam = myPoints->ChangeValue(NoU, NoV);
588 aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value()));
592 // For search of minimum compute distances to mesh.
593 if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
595 // This is the tolerance of difference of squared values.
596 // No need to set it too small.
597 const Standard_Real aDiffTol = mytolu + mytolv;
599 // Step 2. Compute distances to edges.
600 // Assume UEdge(i, j) = { Point(i, j); Point(i + 1, j ) }
601 // Assume VEdge(i, j) = { Point(i, j); Point(i, j + 1) }
602 for ( NoU = 1 ; NoU <= myusample; NoU++ )
604 for ( NoV = 1 ; NoV <= myvsample; NoV++)
606 const Extrema_POnSurfParams &aParam0 = myPoints->Value(NoU, NoV);
610 // Compute parameters to UEdge.
611 const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU + 1, NoV);
612 const Extrema_POnSurfParams &anEdgeParam = ComputeEdgeParameters(Standard_True, aParam0, aParam1, thePoint, aDiffTol);
614 myUEdgePntParams->SetValue(NoU, NoV, anEdgeParam);
619 // Compute parameters to VEdge.
620 const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU, NoV + 1);
621 const Extrema_POnSurfParams &anEdgeParam = ComputeEdgeParameters(Standard_False, aParam0, aParam1, thePoint, aDiffTol);
623 myVEdgePntParams->SetValue(NoU, NoV, anEdgeParam);
628 // Step 3. Compute distances to faces.
629 // Assume myFacePntParams(i, j) =
630 // { Point(i, j); Point(i + 1, j); Point(i + 1, j + 1); Point(i, j + 1) }
632 // { UEdge(i, j); VEdge(i + 1, j); UEdge(i, j + 1); VEdge(i, j) }
633 Standard_Real aSqrDist01;
634 Standard_Real aDiffDist;
635 Standard_Boolean isOut;
637 for ( NoU = 1 ; NoU < myusample; NoU++ ) {
638 for ( NoV = 1 ; NoV < myvsample; NoV++) {
639 const Extrema_POnSurfParams &aUE0 = myUEdgePntParams->Value(NoU, NoV);
640 const Extrema_POnSurfParams &aUE1 = myUEdgePntParams->Value(NoU, NoV+1);
641 const Extrema_POnSurfParams &aVE0 = myVEdgePntParams->Value(NoU, NoV);
642 const Extrema_POnSurfParams &aVE1 = myVEdgePntParams->Value(NoU+1, NoV);
644 aSqrDist01 = aUE0.Value().SquareDistance(aUE1.Value());
645 aDiffDist = Abs(aUE0.GetSqrDistance() - aUE1.GetSqrDistance());
646 isOut = Standard_False;
648 if (aDiffDist >= aSqrDist01 - aDiffTol) {
649 // The projection is outside the face.
650 isOut = Standard_True;
652 aSqrDist01 = aVE0.Value().SquareDistance(aVE1.Value());
653 aDiffDist = Abs(aVE0.GetSqrDistance() - aVE1.GetSqrDistance());
655 if (aDiffDist >= aSqrDist01 - aDiffTol) {
656 // The projection is outside the face.
657 isOut = Standard_True;
662 // Get the closest point on an edge.
663 const Extrema_POnSurfParams &aUEMin =
664 aUE0.GetSqrDistance() < aUE1.GetSqrDistance() ? aUE0 : aUE1;
665 const Extrema_POnSurfParams &aVEMin =
666 aVE0.GetSqrDistance() < aVE1.GetSqrDistance() ? aVE0 : aVE1;
667 const Extrema_POnSurfParams &aEMin =
668 aUEMin.GetSqrDistance() < aVEMin.GetSqrDistance() ? aUEMin : aVEMin;
670 myFacePntParams->SetValue(NoU, NoV, aEMin);
672 // Find closest point inside the face.
678 // Compute U parameter.
679 aUE0.Parameter(aU[0], aV[0]);
680 aUE1.Parameter(aU[1], aV[1]);
681 aUPar = 0.5*(aU[0] + aU[1]);
682 // Compute V parameter.
683 aVE0.Parameter(aU[0], aV[0]);
684 aVE1.Parameter(aU[1], aV[1]);
685 aVPar = 0.5*(aV[0] + aV[1]);
687 Extrema_POnSurfParams aParam(aUPar, aVPar, myS->Value(aUPar, aVPar));
689 aParam.SetElementType(Extrema_Face);
690 aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value()));
691 aParam.SetIndices(NoU, NoV);
692 myFacePntParams->SetValue(NoU, NoV, aParam);
697 // Fill boundary with RealLast square distance.
698 for (NoV = 0; NoV <= myvsample; NoV++) {
699 myFacePntParams->ChangeValue(0, NoV).SetSqrDistance(RealLast());
700 myFacePntParams->ChangeValue(myusample, NoV).SetSqrDistance(RealLast());
703 for (NoU = 1; NoU < myusample; NoU++) {
704 myFacePntParams->ChangeValue(NoU, 0).SetSqrDistance(RealLast());
705 myFacePntParams->ChangeValue(NoU, myvsample).SetSqrDistance(RealLast());
710 // Parameterisation of the sample
711 void Extrema_GenExtPS::BuildTree()
713 // if tree already exists, assume it is already correctly filled
714 if ( ! mySphereUBTree.IsNull() )
717 if (myS->GetType() == GeomAbs_BSplineSurface) {
718 Handle(Geom_BSplineSurface) aBspl = myS->BSpline();
719 Standard_Integer aUValue = aBspl->UDegree() * aBspl->NbUKnots();
720 Standard_Integer aVValue = aBspl->VDegree() * aBspl->NbVKnots();
721 if (aUValue > myusample)
723 if (aVValue > myvsample)
727 Standard_Real PasU = myusup - myumin;
728 Standard_Real PasV = myvsup - myvmin;
729 Standard_Real U0 = PasU / myusample / 100.;
730 Standard_Real V0 = PasV / myvsample / 100.;
732 PasU = (PasU - U0) / (myusample - 1);
733 PasV = (PasV - V0) / (myvsample - 1);
737 //build grid of parametric points
738 myUParams = new TColStd_HArray1OfReal(1,myusample );
739 myVParams = new TColStd_HArray1OfReal(1,myvsample );
740 Standard_Integer NoU, NoV;
741 Standard_Real U = U0, V = V0;
742 for ( NoU = 1 ; NoU <= myusample; NoU++, U += PasU)
743 myUParams->SetValue(NoU, U);
744 for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV)
745 myVParams->SetValue(NoV, V);
747 // Calculation of distances
748 mySphereUBTree = new Extrema_UBTreeOfSphere;
749 Extrema_UBTreeFillerOfSphere aFiller(*mySphereUBTree);
750 Standard_Integer i = 0;
752 mySphereArray = new Bnd_HArray1OfSphere(0, myusample * myvsample);
754 for ( NoU = 1; NoU <= myusample; NoU++ ) {
755 for ( NoV = 1; NoV <= myvsample; NoV++) {
756 P1 = myS->Value(myUParams->Value(NoU), myVParams->Value(NoV));
757 Bnd_Sphere aSph(P1.XYZ(), 0/*mytolu < mytolv ? mytolu : mytolv*/, NoU, NoV);
758 aFiller.Add(i, aSph);
759 mySphereArray->SetValue( i, aSph );
766 void Extrema_GenExtPS::FindSolution(const gp_Pnt& /*P*/,
767 const Extrema_POnSurfParams &theParams)
769 math_Vector Tol(1,2);
773 math_Vector UV(1, 2);
774 theParams.Parameter(UV(1), UV(2));
776 math_Vector UVinf(1,2), UVsup(1,2);
782 math_FunctionSetRoot S(myF, Tol);
783 S.Perform(myF, UV, UVinf, UVsup);
785 myDone = Standard_True;
788 void Extrema_GenExtPS::SetFlag(const Extrema_ExtFlag F)
793 void Extrema_GenExtPS::SetAlgo(const Extrema_ExtAlgo A)
796 myInit = Standard_False;
801 void Extrema_GenExtPS::Perform(const gp_Pnt& P)
803 myDone = Standard_False;
806 if(myAlgo == Extrema_ExtAlgo_Grad)
809 Standard_Integer NoU,NoV;
811 if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
813 Extrema_ElementType anElemType;
816 Standard_Integer iU2;
817 Standard_Integer iV2;
818 Standard_Boolean isMin;
821 for (NoU = 1; NoU < myusample; NoU++) {
822 for (NoV = 1; NoV < myvsample; NoV++) {
823 const Extrema_POnSurfParams &aParam =
824 myFacePntParams->Value(NoU, NoV);
826 isMin = Standard_False;
827 anElemType = aParam.GetElementType();
829 if (anElemType == Extrema_Face) {
830 isMin = Standard_True;
832 // Check if it is a boundary edge or corner vertex.
833 aParam.GetIndices(iU, iV);
835 if (anElemType == Extrema_UIsoEdge) {
836 isMin = (iV == 1 || iV == myvsample);
837 } else if (anElemType == Extrema_VIsoEdge) {
838 isMin = (iU == 1 || iU == myusample);
839 } else if (anElemType == Extrema_Node) {
840 isMin = (iU == 1 || iU == myusample) &&
841 (iV == 1 || iV == myvsample);
845 // This is a middle element.
846 if (anElemType == Extrema_UIsoEdge ||
847 (anElemType == Extrema_Node && (iU == 1 || iU == myusample))) {
848 // Check the down face.
849 const Extrema_POnSurfParams &aDownParam =
850 myFacePntParams->Value(NoU, NoV - 1);
852 if (aDownParam.GetElementType() == anElemType) {
853 aDownParam.GetIndices(iU2, iV2);
854 isMin = (iU == iU2 && iV == iV2);
856 } else if (anElemType == Extrema_VIsoEdge ||
857 (anElemType == Extrema_Node && (iV == 1 || iV == myvsample))) {
858 // Check the right face.
859 const Extrema_POnSurfParams &aRightParam =
860 myFacePntParams->Value(NoU - 1, NoV);
862 if (aRightParam.GetElementType() == anElemType) {
863 aRightParam.GetIndices(iU2, iV2);
864 isMin = (iU == iU2 && iV == iV2);
866 } else if (iU == NoU && iV == NoV) {
867 // Check the lower-left node. For this purpose it is necessary
868 // to check lower-left, lower and left faces.
869 isMin = Standard_True;
871 const Extrema_POnSurfParams *anOtherParam[3] =
872 { &myFacePntParams->Value(NoU, NoV - 1), // Down
873 &myFacePntParams->Value(NoU - 1, NoV - 1), // Lower-left
874 &myFacePntParams->Value(NoU - 1, NoV) }; // Left
876 for (i = 0; i < 3 && isMin; i++) {
877 if (anOtherParam[i]->GetElementType() == Extrema_Node) {
878 anOtherParam[i]->GetIndices(iU2, iV2);
879 isMin = (iU == iU2 && iV == iV2);
881 isMin = Standard_False;
889 FindSolution(P, aParam);
895 if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX)
899 for (NoU = 1; NoU <= myusample; NoU++)
901 for (NoV = 1; NoV <= myvsample; NoV++)
903 const Extrema_POnSurfParams &aParamMain = myPoints->Value(NoU, NoV);
904 const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU - 1, NoV - 1);
905 const Extrema_POnSurfParams &aParam2 = myPoints->Value(NoU - 1, NoV);
906 const Extrema_POnSurfParams &aParam3 = myPoints->Value(NoU - 1, NoV + 1);
907 const Extrema_POnSurfParams &aParam4 = myPoints->Value(NoU, NoV - 1);
908 const Extrema_POnSurfParams &aParam5 = myPoints->Value(NoU, NoV + 1);
909 const Extrema_POnSurfParams &aParam6 = myPoints->Value(NoU + 1, NoV - 1);
910 const Extrema_POnSurfParams &aParam7 = myPoints->Value(NoU + 1, NoV);
911 const Extrema_POnSurfParams &aParam8 = myPoints->Value(NoU + 1, NoV + 1);
913 Dist = aParamMain.GetSqrDistance();
915 if ((aParam1.GetSqrDistance() <= Dist) &&
916 (aParam2.GetSqrDistance() <= Dist) &&
917 (aParam3.GetSqrDistance() <= Dist) &&
918 (aParam4.GetSqrDistance() <= Dist) &&
919 (aParam5.GetSqrDistance() <= Dist) &&
920 (aParam6.GetSqrDistance() <= Dist) &&
921 (aParam7.GetSqrDistance() <= Dist) &&
922 (aParam8.GetSqrDistance() <= Dist))
925 FindSolution(P, myPoints->Value(NoU, NoV));
934 if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX)
936 Bnd_Sphere aSol = mySphereArray->Value(0);
937 Bnd_SphereUBTreeSelectorMin aSelector(mySphereArray, aSol);
938 //aSelector.SetMaxDist( RealLast() );
939 aSelector.DefineCheckPoint( P );
940 mySphereUBTree->Select( aSelector );
941 //TODO: check if no solution in binary tree
942 Bnd_Sphere& aSph = aSelector.Sphere();
943 Standard_Real aU = myUParams->Value(aSph.U());
944 Standard_Real aV = myVParams->Value(aSph.V());
945 Extrema_POnSurfParams aParams(aU, aV, myS->Value(aU, aV));
947 aParams.SetSqrDistance(P.SquareDistance(aParams.Value()));
948 aParams.SetIndices(aSph.U(), aSph.V());
949 FindSolution(P, aParams);
951 if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX)
953 Bnd_Sphere aSol = mySphereArray->Value(0);
954 Bnd_SphereUBTreeSelectorMax aSelector(mySphereArray, aSol);
955 //aSelector.SetMaxDist( RealLast() );
956 aSelector.DefineCheckPoint( P );
957 mySphereUBTree->Select( aSelector );
958 //TODO: check if no solution in binary tree
959 Bnd_Sphere& aSph = aSelector.Sphere();
960 Standard_Real aU = myUParams->Value(aSph.U());
961 Standard_Real aV = myVParams->Value(aSph.V());
962 Extrema_POnSurfParams aParams(aU, aV, myS->Value(aU, aV));
963 aParams.SetSqrDistance(P.SquareDistance(aParams.Value()));
964 aParams.SetIndices(aSph.U(), aSph.V());
966 FindSolution(P, aParams);
970 //=============================================================================
972 Standard_Boolean Extrema_GenExtPS::IsDone () const { return myDone; }
973 //=============================================================================
975 Standard_Integer Extrema_GenExtPS::NbExt () const
977 if (!IsDone()) { throw StdFail_NotDone(); }
980 //=============================================================================
982 Standard_Real Extrema_GenExtPS::SquareDistance (const Standard_Integer N) const
984 if (!IsDone()) { throw StdFail_NotDone(); }
985 return myF.SquareDistance(N);
987 //=============================================================================
989 const Extrema_POnSurf& Extrema_GenExtPS::Point (const Standard_Integer N) const
991 if (!IsDone()) { throw StdFail_NotDone(); }
994 //=============================================================================