X-Git-Url: http://git.dev.opencascade.org/gitweb/?p=occt.git;a=blobdiff_plain;f=src%2FExtrema%2FExtrema_GenExtPS.cxx;h=078de1cb780a5ec6bec7629b3862459baba58ed0;hp=56b211841b5d65af87ac2cdba4cde471e261940e;hb=4e18e72a22679f14d021c5fe7d5f9628107369b8;hpb=6062bf3e4a03270eebc8059d22bd9780bb45661d diff --git a/src/Extrema/Extrema_GenExtPS.cxx b/src/Extrema/Extrema_GenExtPS.cxx index 56b211841b..078de1cb78 100755 --- a/src/Extrema/Extrema_GenExtPS.cxx +++ b/src/Extrema/Extrema_GenExtPS.cxx @@ -1,7 +1,23 @@ -// File: Extrema_GenExtPS.cxx -// Created: Tue Jul 18 08:21:34 1995 -// Author: Modelistation -// +// Created on: 1995-07-18 +// Created by: Modelistation +// Copyright (c) 1995-1999 Matra Datavision +// Copyright (c) 1999-2012 OPEN CASCADE SAS +// +// The content of this file is subject to the Open CASCADE Technology Public +// License Version 6.5 (the "License"). You may not use the content of this file +// except in compliance with the License. Please obtain a copy of the License +// at http://www.opencascade.org and read it completely before using this file. +// +// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its +// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France. +// +// The Original Code and all software distributed under the License is +// distributed on an "AS IS" basis, without warranty of any kind, and the +// Initial Developer hereby disclaims all such warranties, including without +// limitation, any warranties of merchantability, fitness for a particular +// purpose or non-infringement. Please see the License for the specific terms +// and conditions governing the rights and limitations under the License. + // Modified by skv - Thu Sep 30 15:21:07 2004 OCC593 @@ -21,7 +37,16 @@ #include #include #include - +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include //IMPLEMENT_HARRAY1(Extrema_HArray1OfSphere) @@ -147,7 +172,69 @@ Standard_Boolean Bnd_SphereUBTreeSelectorMax::Accept(const Standard_Integer& the return Standard_False; } +/* + * This function computes the point on surface parameters on edge. + * if it coincides with theParam0 or theParam1, it is returned. + */ +static Extrema_POnSurfParams ComputeEdgeParameters + (const Standard_Boolean IsUEdge, + const Extrema_POnSurfParams &theParam0, + const Extrema_POnSurfParams &theParam1, + const Adaptor3d_SurfacePtr &theSurf, + const gp_Pnt &thePoint, + const Standard_Real theDiffTol) +{ + const Standard_Real aSqrDist01 = + theParam0.Value().SquareDistance(theParam1.Value()); + + if (aSqrDist01 <= theDiffTol) { + // The points are confused. Get the first point and change its type. + return theParam0; + } else { + const Standard_Real aDiffDist = + Abs(theParam0.GetSqrDistance() - theParam1.GetSqrDistance()); + + if (aDiffDist >= aSqrDist01 - theDiffTol) { + // The shortest distance is one of the nodes. + if (theParam0.GetSqrDistance() > theParam1.GetSqrDistance()) { + // The shortest distance is the point 1. + return theParam1; + } else { + // The shortest distance is the point 0. + return theParam0; + } + } else { + // The shortest distance is inside the edge. + gp_XYZ aPoP(thePoint.XYZ().Subtracted(theParam0.Value().XYZ())); + gp_XYZ aPoP1(theParam1.Value().XYZ().Subtracted(theParam0.Value().XYZ())); + Standard_Real aRatio = aPoP.Dot(aPoP1)/aSqrDist01; + Standard_Real aU[2]; + Standard_Real aV[2]; + + theParam0.Parameter(aU[0], aV[0]); + theParam1.Parameter(aU[1], aV[1]); + + Standard_Real aUPar = aU[0]; + Standard_Real aVPar = aV[0]; + if (IsUEdge) { + aUPar += aRatio*(aU[1] - aU[0]); + } else { + aVPar += aRatio*(aV[1] - aV[0]); + } + + Extrema_POnSurfParams aParam(aUPar, aVPar, theSurf->Value(aUPar, aVPar)); + Standard_Integer anIndices[2]; + + theParam0.GetIndices(anIndices[0], anIndices[1]); + aParam.SetElementType(IsUEdge ? Extrema_UIsoEdge : Extrema_VIsoEdge); + aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value())); + aParam.SetIndices(anIndices[0], anIndices[1]); + + return aParam; + } + } +} //============================================================================= @@ -194,52 +281,6 @@ Processing: - update table TbSel. -----------------------------------------------------------------------------*/ -static Standard_Boolean IsoIsDeg (const Adaptor3d_Surface& S, - const Standard_Real Param, - const GeomAbs_IsoType IT, - const Standard_Real TolMin, - const Standard_Real TolMax) -{ - Standard_Real U1=0.,U2=0.,V1=0.,V2=0.,T; - Standard_Boolean Along = Standard_True; - U1 = S.FirstUParameter(); - U2 = S.LastUParameter(); - V1 = S.FirstVParameter(); - V2 = S.LastVParameter(); - gp_Vec D1U,D1V; - gp_Pnt P; - Standard_Real Step,D1NormMax; - if (IT == GeomAbs_IsoV) - { - Step = (U2 - U1)/10; - D1NormMax=0.; - for (T=U1;T<=U2;T=T+Step) - { - S.D1(T,Param,P,D1U,D1V); - D1NormMax=Max(D1NormMax,D1U.Magnitude()); - } - - if (D1NormMax >TolMax || D1NormMax < TolMin ) - Along = Standard_False; - } - else - { - Step = (V2 - V1)/10; - D1NormMax=0.; - for (T=V1;T<=V2;T=T+Step) - { - S.D1(Param,T,P,D1U,D1V); - D1NormMax=Max(D1NormMax,D1V.Magnitude()); - } - - if (D1NormMax >TolMax || D1NormMax < TolMin ) - Along = Standard_False; - - - } - return Along; -} -//---------------------------------------------------------- Extrema_GenExtPS::Extrema_GenExtPS() { myDone = Standard_False; @@ -306,7 +347,6 @@ void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S, const Standard_Real TolU, const Standard_Real TolV) { - myInit = Standard_True; myS = (Adaptor3d_SurfacePtr)&S; myusample = NbU; myvsample = NbV; @@ -323,34 +363,334 @@ void Extrema_GenExtPS::Initialize(const Adaptor3d_Surface& S, myF.Initialize(S); mySphereUBTree.Nullify(); + myUParams.Nullify(); + myVParams.Nullify(); + myInit = Standard_False; +} - if(myAlgo == Extrema_ExtAlgo_Grad) +inline static void fillParams(const TColStd_Array1OfReal& theKnots, + Standard_Integer theDegree, + Standard_Real theParMin, + Standard_Real theParMax, + Handle_TColStd_HArray1OfReal& theParams, + Standard_Integer theSample) +{ + NCollection_Vector aParams; + Standard_Integer i = 1; + Standard_Real aPrevPar = theParMin; + aParams.Append(aPrevPar); + //calculation the array of parametric points depending on the knots array variation and degree of given surface + for ( ; i < theKnots.Length() && theKnots(i) < (theParMax - Precision::PConfusion()); i++ ) { - //If flag was changed and extrema not reinitialized Extrema would fail - mypoints = new TColgp_HArray2OfPnt(0,myusample+1,0,myvsample+1); - Standard_Real PasU = myusup - myumin; - Standard_Real PasV = myvsup - myvmin; - Standard_Real U0 = PasU / myusample / 100.; - Standard_Real V0 = PasV / myvsample / 100.; - gp_Pnt P1; - PasU = (PasU - U0) / (myusample - 1); - PasV = (PasV - V0) / (myvsample - 1); - U0 = U0/2. + myumin; - V0 = V0/2. + myvmin; + if ( theKnots(i+1) < theParMin + Precision::PConfusion()) + continue; + + Standard_Real aStep = (theKnots(i+1) - theKnots(i))/Max(theDegree,2); + Standard_Integer k =1; + for( ; k <= theDegree ; k++) + { + Standard_Real aPar = theKnots(i) + k * aStep; + if(aPar > theParMax - Precision::PConfusion()) + break; + if(aPar > aPrevPar + Precision::PConfusion() ) + { + aParams.Append(aPar); + aPrevPar = aPar; + } + } + + } + aParams.Append(theParMax); + Standard_Integer nbPar = aParams.Length(); + //in case of an insufficient number of points the grid will be built later + if (nbPar < theSample) + return; + theParams = new TColStd_HArray1OfReal(1, nbPar ); + for( i = 0; i < nbPar; i++) + theParams->SetValue(i+1,aParams(i)); +} + +void Extrema_GenExtPS::GetGridPoints( const Adaptor3d_Surface& theSurf) +{ + //creation parametric points for BSpline and Bezier surfaces + //with taking into account of Degree and NbKnots of BSpline or Bezier geometry + if(theSurf.GetType() == GeomAbs_OffsetSurface) + GetGridPoints(theSurf.BasisSurface()->Surface()); + //parametric points for BSpline surfaces + else if( theSurf.GetType() == GeomAbs_BSplineSurface) + { + Handle(Geom_BSplineSurface) aBspl = theSurf.BSpline(); + if(!aBspl.IsNull()) + { + TColStd_Array1OfReal aUKnots(1, aBspl->NbUKnots()); + aBspl->UKnots( aUKnots); + TColStd_Array1OfReal aVKnots(1, aBspl->NbVKnots()); + aBspl->VKnots( aVKnots); + fillParams(aUKnots,aBspl->UDegree(),myumin, myusup, myUParams, myusample); + fillParams(aVKnots,aBspl->VDegree(),myvmin, myvsup, myVParams, myvsample); + } + } + //calculation parametric points for Bezier surfaces + else if(theSurf.GetType() == GeomAbs_BezierSurface) + { + Handle(Geom_BezierSurface) aBezier = theSurf.Bezier(); + if(aBezier.IsNull()) + return; -// Calculation of distances + TColStd_Array1OfReal aUKnots(1,2); + TColStd_Array1OfReal aVKnots(1,2); + aBezier->Bounds(aUKnots(1), aUKnots(2), aVKnots(1), aVKnots(2)); + fillParams(aUKnots, aBezier->UDegree(), myumin, myusup, myUParams, myusample); + fillParams(aVKnots, aBezier->VDegree(), myvmin, myvsup, myVParams, myvsample); + } + //creation points for surfaces based on BSpline or Bezier curves + else if(theSurf.GetType() == GeomAbs_SurfaceOfRevolution || + theSurf.GetType() == GeomAbs_SurfaceOfExtrusion) + { + Handle(TColStd_HArray1OfReal) anArrKnots; + Standard_Integer aDegree = 0; + if(theSurf.BasisCurve()->Curve().GetType() == GeomAbs_BSplineCurve) + { + Handle(Geom_BSplineCurve) aBspl = theSurf.BasisCurve()->Curve().BSpline(); + if(!aBspl.IsNull()) + { + anArrKnots = new TColStd_HArray1OfReal(1,aBspl->NbKnots()); + aBspl->Knots( anArrKnots->ChangeArray1() ); + aDegree = aBspl->Degree(); + + } + + } + if(theSurf.BasisCurve()->Curve().GetType() == GeomAbs_BezierCurve) + { + Handle(Geom_BezierCurve) aBez = theSurf.BasisCurve()->Curve().Bezier(); + if(!aBez.IsNull()) + { + anArrKnots = new TColStd_HArray1OfReal(1,2); + anArrKnots->SetValue(1, aBez->FirstParameter()); + anArrKnots->SetValue(2, aBez->LastParameter()); + aDegree = aBez->Degree(); + + } + } + if(anArrKnots.IsNull()) + return; + if( theSurf.GetType() == GeomAbs_SurfaceOfRevolution ) + fillParams( anArrKnots->Array1(), aDegree, myvmin, myvsup, myVParams, myvsample); + else + fillParams( anArrKnots->Array1(), aDegree, myumin, myusup, myUParams, myusample); + } + //update the number of points in sample + if(!myUParams.IsNull()) + myusample = myUParams->Length(); + if( !myVParams.IsNull()) + myvsample = myVParams->Length(); + +} + +void Extrema_GenExtPS::BuildGrid(const gp_Pnt &thePoint) +{ Standard_Integer NoU, NoV; - Standard_Real U, V; - for ( NoU = 1, U = U0; NoU <= myusample; NoU++, U += PasU) { - for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) { - P1 = myS->Value(U, V); - mypoints->SetValue(NoU,NoV,P1); + + //if grid was already built skip its creation + if (!myInit) { + //build parametric grid in case of a complex surface geometry (BSpline and Bezier surfaces) + GetGridPoints(*myS); + + //build grid in other cases + if( myUParams.IsNull() ) + { + Standard_Real PasU = myusup - myumin; + Standard_Real U0 = PasU / myusample / 100.; + PasU = (PasU - U0) / (myusample - 1); + U0 = U0/2. + myumin; + myUParams = new TColStd_HArray1OfReal(1,myusample ); + Standard_Integer NoU; + Standard_Real U = U0; + for ( NoU = 1 ; NoU <= myusample; NoU++, U += PasU) + myUParams->SetValue(NoU, U); + } + + if( myVParams.IsNull()) + { + Standard_Real PasV = myvsup - myvmin; + Standard_Real V0 = PasV / myvsample / 100.; + PasV = (PasV - V0) / (myvsample - 1); + V0 = V0/2. + myvmin; + + myVParams = new TColStd_HArray1OfReal(1,myvsample ); + Standard_Integer NoV; + Standard_Real V = V0; + + for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) + myVParams->SetValue(NoV, V); + } + + //If flag was changed and extrema not reinitialized Extrema would fail + myPoints = new Extrema_HArray2OfPOnSurfParams + (0, myusample + 1, 0, myvsample + 1); + // Calculation of distances + + for ( NoU = 1 ; NoU <= myusample; NoU++ ) { + for ( NoV = 1 ; NoV <= myvsample; NoV++) { + gp_Pnt aP1 = myS->Value(myUParams->Value(NoU), myVParams->Value(NoV)); + Extrema_POnSurfParams aParam + (myUParams->Value(NoU), myVParams->Value(NoV), aP1); + + aParam.SetElementType(Extrema_Node); + aParam.SetIndices(NoU, NoV); + myPoints->SetValue(NoU, NoV, aParam); + } + } + + // Fill boundary with negative square distance. + // It is used for computation of Maximum. + for (NoV = 0; NoV <= myvsample + 1; NoV++) { + myPoints->ChangeValue(0, NoV).SetSqrDistance(-1.); + myPoints->ChangeValue(myusample + 1, NoV).SetSqrDistance(-1.); } + + for (NoU = 1; NoU <= myusample; NoU++) { + myPoints->ChangeValue(NoU, 0).SetSqrDistance(-1.); + myPoints->ChangeValue(NoU, myvsample + 1).SetSqrDistance(-1.); + } + + myInit = Standard_True; } + + // Compute distances to mesh. + // Step 1. Compute distances to nodes. + for ( NoU = 1 ; NoU <= myusample; NoU++ ) { + for ( NoV = 1 ; NoV <= myvsample; NoV++) { + Extrema_POnSurfParams &aParam = myPoints->ChangeValue(NoU, NoV); + + aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value())); + } } - //mypoints = new TColgp_HArray2OfPnt(0,myusample+1,0,myvsample+1); + // For search of minimum compute distances to mesh. + if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX) + { + // This is the tolerance of difference of squared values. + // No need to set it too small. + const Standard_Real aDiffTol = mytolu + mytolv; + + // Step 2. Compute distances to edges. + // Assume UEdge(i, j) = { Point(i, j); Point(i + 1, j ) } + // Assume VEdge(i, j) = { Point(i, j); Point(i, j + 1) } + Handle(Extrema_HArray2OfPOnSurfParams) aUEdgePntParams = + new Extrema_HArray2OfPOnSurfParams(1, myusample - 1, 1, myvsample); + Handle(Extrema_HArray2OfPOnSurfParams) aVEdgePntParams = + new Extrema_HArray2OfPOnSurfParams(1, myusample, 1, myvsample - 1); + + for ( NoU = 1 ; NoU <= myusample; NoU++ ) { + for ( NoV = 1 ; NoV <= myvsample; NoV++) { + const Extrema_POnSurfParams &aParam0 = myPoints->Value(NoU, NoV); + + if (NoU < myusample) { + // Compute parameters to UEdge. + const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU + 1, NoV); + Extrema_POnSurfParams aUEdgeParam = ComputeEdgeParameters + (Standard_True, aParam0, aParam1, myS, thePoint, aDiffTol); + + aUEdgePntParams->SetValue(NoU, NoV, aUEdgeParam); + } + + if (NoV < myvsample) { + // Compute parameters to VEdge. + const Extrema_POnSurfParams &aParam1 = myPoints->Value(NoU, NoV + 1); + Extrema_POnSurfParams aVEdgeParam = ComputeEdgeParameters + (Standard_False, aParam0, aParam1, myS, thePoint, aDiffTol); + + aVEdgePntParams->SetValue(NoU, NoV, aVEdgeParam); + } + } + } + + // Step 3. Compute distances to faces. + // Assume myFacePntParams(i, j) = + // { Point(i, j); Point(i + 1, j); Point(i + 1, j + 1); Point(i, j + 1) } + // Or + // { UEdge(i, j); VEdge(i + 1, j); UEdge(i, j + 1); VEdge(i, j) } + myFacePntParams = + new Extrema_HArray2OfPOnSurfParams(0, myusample, 0, myvsample); + Standard_Real aSqrDist01; + Standard_Real aDiffDist; + Standard_Boolean isOut; + + for ( NoU = 1 ; NoU < myusample; NoU++ ) { + for ( NoV = 1 ; NoV < myvsample; NoV++) { + const Extrema_POnSurfParams &aUE0 = aUEdgePntParams->Value(NoU, NoV); + const Extrema_POnSurfParams &aUE1 = aUEdgePntParams->Value(NoU, NoV+1); + const Extrema_POnSurfParams &aVE0 = aVEdgePntParams->Value(NoU, NoV); + const Extrema_POnSurfParams &aVE1 = aVEdgePntParams->Value(NoU+1, NoV); + + aSqrDist01 = aUE0.Value().SquareDistance(aUE1.Value()); + aDiffDist = Abs(aUE0.GetSqrDistance() - aUE1.GetSqrDistance()); + isOut = Standard_False; + + if (aDiffDist >= aSqrDist01 - aDiffTol) { + // The projection is outside the face. + isOut = Standard_True; + } else { + aSqrDist01 = aVE0.Value().SquareDistance(aVE1.Value()); + aDiffDist = Abs(aVE0.GetSqrDistance() - aVE1.GetSqrDistance()); + + if (aDiffDist >= aSqrDist01 - aDiffTol) { + // The projection is outside the face. + isOut = Standard_True; + } + } + + if (isOut) { + // Get the closest point on an edge. + const Extrema_POnSurfParams &aUEMin = + aUE0.GetSqrDistance() < aUE1.GetSqrDistance() ? aUE0 : aUE1; + const Extrema_POnSurfParams &aVEMin = + aVE0.GetSqrDistance() < aVE1.GetSqrDistance() ? aVE0 : aVE1; + const Extrema_POnSurfParams &aEMin = + aUEMin.GetSqrDistance() < aVEMin.GetSqrDistance() ? aUEMin : aVEMin; + + myFacePntParams->SetValue(NoU, NoV, aEMin); + } else { + // Find closest point inside the face. + Standard_Real aU[2]; + Standard_Real aV[2]; + Standard_Real aUPar; + Standard_Real aVPar; + + // Compute U parameter. + aUE0.Parameter(aU[0], aV[0]); + aUE1.Parameter(aU[1], aV[1]); + aUPar = 0.5*(aU[0] + aU[1]); + // Compute V parameter. + aVE0.Parameter(aU[0], aV[0]); + aVE1.Parameter(aU[1], aV[1]); + aVPar = 0.5*(aV[0] + aV[1]); + + Extrema_POnSurfParams aParam(aUPar, aVPar, myS->Value(aUPar, aVPar)); + + aParam.SetElementType(Extrema_Face); + aParam.SetSqrDistance(thePoint.SquareDistance(aParam.Value())); + aParam.SetIndices(NoU, NoV); + myFacePntParams->SetValue(NoU, NoV, aParam); + } + } + } + + // Fill boundary with RealLast square distance. + for (NoV = 0; NoV <= myvsample; NoV++) { + myFacePntParams->ChangeValue(0, NoV).SetSqrDistance(RealLast()); + myFacePntParams->ChangeValue(myusample, NoV).SetSqrDistance(RealLast()); + } + + for (NoU = 1; NoU < myusample; NoU++) { + myFacePntParams->ChangeValue(NoU, 0).SetSqrDistance(RealLast()); + myFacePntParams->ChangeValue(NoU, myvsample).SetSqrDistance(RealLast()); + } + } +} /* a- Constitution of the table of distances (TbDist(0,myusample+1,0,myvsample+1)): @@ -359,9 +699,6 @@ a- Constitution of the table of distances (TbDist(0,myusample+1,0,myvsample+1)): // Parameterisation of the sample - -} - void Extrema_GenExtPS::BuildTree() { // if tree already exists, assume it is already correctly filled @@ -378,16 +715,26 @@ void Extrema_GenExtPS::BuildTree() U0 = U0/2. + myumin; V0 = V0/2. + myvmin; + //build grid of parametric points + myUParams = new TColStd_HArray1OfReal(1,myusample ); + myVParams = new TColStd_HArray1OfReal(1,myvsample ); + Standard_Integer NoU, NoV; + Standard_Real U = U0, V = V0; + for ( NoU = 1 ; NoU <= myusample; NoU++, U += PasU) + myUParams->SetValue(NoU, U); + for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) + myVParams->SetValue(NoV, V); + // Calculation of distances mySphereUBTree = new Extrema_UBTreeOfSphere; Extrema_UBTreeFillerOfSphere aFiller(*mySphereUBTree); Standard_Integer i = 0; - Standard_Real U, V; + mySphereArray = new Bnd_HArray1OfSphere(0, myusample * myvsample); - Standard_Integer NoU, NoV; - for ( NoU = 1, U = U0; NoU <= myusample; NoU++, U += PasU) { - for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) { - P1 = myS->Value(U, V); + + for ( NoU = 1; NoU <= myusample; NoU++ ) { + for ( NoV = 1; NoV <= myvsample; NoV++) { + P1 = myS->Value(myUParams->Value(NoU), myVParams->Value(NoV)); Bnd_Sphere aSph(P1.XYZ(), 0/*mytolu < mytolv ? mytolu : mytolv*/, NoU, NoV); aFiller.Add(i, aSph); mySphereArray->SetValue( i, aSph ); @@ -397,12 +744,17 @@ void Extrema_GenExtPS::BuildTree() aFiller.Fill(); } -void Extrema_GenExtPS::FindSolution(const gp_Pnt& P, const math_Vector& UV, const Standard_Real PasU, const Standard_Real PasV, const Extrema_ExtFlag f) +void Extrema_GenExtPS::FindSolution(const gp_Pnt& /*P*/, + const Extrema_POnSurfParams &theParams) { math_Vector Tol(1,2); Tol(1) = mytolu; Tol(2) = mytolv; + math_Vector UV(1, 2); + + theParams.Parameter(UV(1), UV(2)); + math_Vector UVinf(1,2), UVsup(1,2); UVinf(1) = myumin; UVinf(2) = myvmin; @@ -411,86 +763,12 @@ void Extrema_GenExtPS::FindSolution(const gp_Pnt& P, const math_Vector& UV, cons math_Vector errors(1,2); math_Vector root(1, 2); - Standard_Real eps = 1.e-9; - Standard_Integer nbsubsample = 11; Standard_Integer aNbMaxIter = 100; - gp_Pnt PStart = myS->Value(UV(1), UV(2)); - Standard_Real DistStart = P.SquareDistance(PStart); - Standard_Real DistSol = DistStart; + gp_Pnt PStart = theParams.Value(); math_FunctionSetRoot S (myF,UV,Tol,UVinf,UVsup, aNbMaxIter); - Standard_Boolean ToResolveOnSubgrid = Standard_False; - if (f == Extrema_ExtFlag_MIN) - { - if(S.IsDone()) - { - root = S.Root(); - myF.Value(root, errors); - gp_Pnt PSol = myS->Value(root(1), root(2)); - DistSol = P.SquareDistance(PSol); - if(Abs(errors(1)) > eps || Abs(errors(2)) > eps || DistStart < DistSol) - //try to improve solution on subgrid of sample points - ToResolveOnSubgrid = Standard_True; - } - else - ToResolveOnSubgrid = Standard_True; - - if (ToResolveOnSubgrid) - { - Standard_Real u1 = Max(UV(1) - PasU, myumin), u2 = Min(UV(1) + PasU, myusup); - Standard_Real v1 = Max(UV(2) - PasV, myvmin), v2 = Min(UV(2) + PasV, myvsup); - - if(u2 - u1 < 2.*PasU) { - if(Abs(u1 - myumin) < 1.e-9) { - u2 = u1 + 2.*PasU; - u2 = Min(u2, myusup); - } - if(Abs(u2 - myusup) < 1.e-9) { - u1 = u2 - 2.*PasU; - u1 = Max(u1, myumin); - } - } - - if(v2 - v1 < 2.*PasV) { - if(Abs(v1 - myvmin) < 1.e-9) { - v2 = v1 + 2.*PasV; - v2 = Min(v2, myvsup); - } - if(Abs(v2 - myvsup) < 1.e-9) { - v1 = v2 - 2.*PasV; - v1 = Max(v1, myvmin); - } - } - - Standard_Real du = (u2 - u1)/(nbsubsample-1); - Standard_Real dv = (v2 - v1)/(nbsubsample-1); - Standard_Real u, v; - Standard_Real dist; - - Standard_Boolean NewSolution = Standard_False; - Standard_Integer Nu, Nv; - for (Nu = 1, u = u1; Nu < nbsubsample; Nu++, u += du) { - for (Nv = 1, v = v1; Nv < nbsubsample; Nv++, v += dv) { - gp_Pnt Puv = myS->Value(u, v); - dist = P.SquareDistance(Puv); - - if(dist < DistSol) { - UV(1) = u; - UV(2) = v; - NewSolution = Standard_True; - DistSol = dist; - } - } - } - - if(NewSolution) { - //try to precise - math_FunctionSetRoot S (myF,UV,Tol,UVinf,UVsup, aNbMaxIter); - } - } //end of if (ToResolveOnSubgrid) - } //end of if (f == Extrema_ExtFlag_MIN) myDone = Standard_True; } @@ -502,70 +780,101 @@ void Extrema_GenExtPS::SetFlag(const Extrema_ExtFlag F) void Extrema_GenExtPS::SetAlgo(const Extrema_ExtAlgo A) { + if(myAlgo != A) + myInit = Standard_False; myAlgo = A; + } void Extrema_GenExtPS::Perform(const gp_Pnt& P) { - //BuildTree(); myDone = Standard_False; myF.SetPoint(P); - - Standard_Real PasU = myusup - myumin; - Standard_Real PasV = myvsup - myvmin; - Standard_Real U, U0 = PasU / myusample / 100.; - Standard_Real V, V0 = PasV / myvsample / 100.; - PasU = (PasU - U0) / (myusample - 1); - PasV = (PasV - V0) / (myvsample - 1); - U0 = U0/2.+myumin; - V0 = V0/2.+myvmin; - - //math_Vector Tol(1, 2); - math_Vector UV(1,2); - + if(myAlgo == Extrema_ExtAlgo_Grad) { + BuildGrid(P); Standard_Integer NoU,NoV; - TColStd_Array2OfReal TheDist(0, myusample+1, 0, myvsample+1); - for ( NoU = 1, U = U0; NoU <= myusample; NoU++, U += PasU) { - for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) { - TheDist(NoU, NoV) = P.SquareDistance(mypoints->Value(NoU, NoV)); - } - } - - TColStd_Array2OfInteger TbSel(0,myusample+1,0,myvsample+1); - TbSel.Init(0); - - Standard_Real Dist; - if(myFlag == Extrema_ExtFlag_MIN || myFlag == Extrema_ExtFlag_MINMAX) { - for (NoV = 0; NoV <= myvsample+1; NoV++) { - TheDist(0,NoV) = RealLast(); - TheDist(myusample+1,NoV) = RealLast(); - } - for (NoU = 1; NoU <= myusample; NoU++) { - TheDist(NoU,0) = RealLast(); - TheDist(NoU,myvsample+1) = RealLast(); - } - for (NoU = 1; NoU <= myusample; NoU++) { - for (NoV = 1; NoV <= myvsample; NoV++) { - if (TbSel(NoU,NoV) == 0) { - Dist = TheDist(NoU,NoV); - if ((TheDist(NoU-1,NoV-1) >= Dist) && - (TheDist(NoU-1,NoV ) >= Dist) && - (TheDist(NoU-1,NoV+1) >= Dist) && - (TheDist(NoU ,NoV-1) >= Dist) && - (TheDist(NoU ,NoV+1) >= Dist) && - (TheDist(NoU+1,NoV-1) >= Dist) && - (TheDist(NoU+1,NoV ) >= Dist) && - (TheDist(NoU+1,NoV+1) >= Dist)) { - //Create array of UV vectors to calculate min - UV(1) = U0 + (NoU - 1) * PasU; - UV(2) = V0 + (NoV - 1) * PasV; - FindSolution(P, UV, PasU, PasV, Extrema_ExtFlag_MIN); + Extrema_ElementType anElemType; + Standard_Integer iU; + Standard_Integer iV; + Standard_Integer iU2; + Standard_Integer iV2; + Standard_Boolean isMin; + Standard_Integer i; + + for (NoU = 1; NoU < myusample; NoU++) { + for (NoV = 1; NoV < myvsample; NoV++) { + const Extrema_POnSurfParams &aParam = + myFacePntParams->Value(NoU, NoV); + + isMin = Standard_False; + anElemType = aParam.GetElementType(); + + if (anElemType == Extrema_Face) { + isMin = Standard_True; + } else { + // Check if it is a boundary edge or corner vertex. + aParam.GetIndices(iU, iV); + + if (anElemType == Extrema_UIsoEdge) { + isMin = (iV == 1 || iV == myvsample); + } else if (anElemType == Extrema_VIsoEdge) { + isMin = (iU == 1 || iU == myusample); + } else if (anElemType == Extrema_Node) { + isMin = (iU == 1 || iU == myusample) && + (iV == 1 || iV == myvsample); } + + if (!isMin) { + // This is a middle element. + if (anElemType == Extrema_UIsoEdge || + (anElemType == Extrema_Node && (iU == 1 || iU == myusample))) { + // Check the down face. + const Extrema_POnSurfParams &aDownParam = + myFacePntParams->Value(NoU, NoV - 1); + + if (aDownParam.GetElementType() == anElemType) { + aDownParam.GetIndices(iU2, iV2); + isMin = (iU == iU2 && iV == iV2); + } + } else if (anElemType == Extrema_VIsoEdge || + (anElemType == Extrema_Node && (iV == 1 || iV == myvsample))) { + // Check the right face. + const Extrema_POnSurfParams &aRightParam = + myFacePntParams->Value(NoU - 1, NoV); + + if (aRightParam.GetElementType() == anElemType) { + aRightParam.GetIndices(iU2, iV2); + isMin = (iU == iU2 && iV == iV2); + } + } else if (iU == NoU && iV == NoV) { + // Check the lower-left node. For this purpose it is necessary + // to check lower-left, lower and left faces. + isMin = Standard_True; + + const Extrema_POnSurfParams *anOtherParam[3] = + { &myFacePntParams->Value(NoU, NoV - 1), // Down + &myFacePntParams->Value(NoU - 1, NoV - 1), // Lower-left + &myFacePntParams->Value(NoU - 1, NoV) }; // Left + + for (i = 0; i < 3 && isMin; i++) { + if (anOtherParam[i]->GetElementType() == Extrema_Node) { + anOtherParam[i]->GetIndices(iU2, iV2); + isMin = (iU == iU2 && iV == iV2); + } else { + isMin = Standard_False; + } + } + } + } + } + + if (isMin) { + FindSolution(P, aParam); } } } @@ -573,31 +882,21 @@ void Extrema_GenExtPS::Perform(const gp_Pnt& P) if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX) { - for (NoV = 0; NoV <= myvsample+1; NoV++) { - TheDist(0,NoV) = RealFirst(); - TheDist(myusample+1,NoV) = RealFirst(); - } - for (NoU = 1; NoU <= myusample; NoU++) { - TheDist(NoU,0) = RealFirst(); - TheDist(NoU,myvsample+1) = RealFirst(); - } + Standard_Real Dist; + for (NoU = 1; NoU <= myusample; NoU++) { for (NoV = 1; NoV <= myvsample; NoV++) { - if (TbSel(NoU,NoV) == 0) { - Dist = TheDist(NoU,NoV); - if ((TheDist(NoU-1,NoV-1) <= Dist) && - (TheDist(NoU-1,NoV ) <= Dist) && - (TheDist(NoU-1,NoV+1) <= Dist) && - (TheDist(NoU ,NoV-1) <= Dist) && - (TheDist(NoU ,NoV+1) <= Dist) && - (TheDist(NoU+1,NoV-1) <= Dist) && - (TheDist(NoU+1,NoV ) <= Dist) && - (TheDist(NoU+1,NoV+1) <= Dist)) { - //Create array of UV vectors to calculate max - UV(1) = U0 + (NoU - 1) * PasU; - UV(2) = V0 + (NoV - 1) * PasV; - FindSolution(P, UV, PasU, PasV, Extrema_ExtFlag_MAX); - } + Dist = myPoints->Value(NoU, NoV).GetSqrDistance(); + if ((myPoints->Value(NoU-1,NoV-1).GetSqrDistance() <= Dist) && + (myPoints->Value(NoU-1,NoV ).GetSqrDistance() <= Dist) && + (myPoints->Value(NoU-1,NoV+1).GetSqrDistance() <= Dist) && + (myPoints->Value(NoU ,NoV-1).GetSqrDistance() <= Dist) && + (myPoints->Value(NoU ,NoV+1).GetSqrDistance() <= Dist) && + (myPoints->Value(NoU+1,NoV-1).GetSqrDistance() <= Dist) && + (myPoints->Value(NoU+1,NoV ).GetSqrDistance() <= Dist) && + (myPoints->Value(NoU+1,NoV+1).GetSqrDistance() <= Dist)) { + // Find maximum. + FindSolution(P, myPoints->Value(NoU, NoV)); } } } @@ -612,14 +911,16 @@ void Extrema_GenExtPS::Perform(const gp_Pnt& P) Bnd_SphereUBTreeSelectorMin aSelector(mySphereArray, aSol); //aSelector.SetMaxDist( RealLast() ); aSelector.DefineCheckPoint( P ); - Standard_Integer aNbSel = mySphereUBTree->Select( aSelector ); + mySphereUBTree->Select( aSelector ); //TODO: check if no solution in binary tree Bnd_Sphere& aSph = aSelector.Sphere(); + Standard_Real aU = myUParams->Value(aSph.U()); + Standard_Real aV = myVParams->Value(aSph.V()); + Extrema_POnSurfParams aParams(aU, aV, myS->Value(aU, aV)); - UV(1) = U0 + (aSph.U() - 1) * PasU; - UV(2) = V0 + (aSph.V() - 1) * PasV; - - FindSolution(P, UV, PasU, PasV, Extrema_ExtFlag_MIN); + aParams.SetSqrDistance(P.SquareDistance(aParams.Value())); + aParams.SetIndices(aSph.U(), aSph.V()); + FindSolution(P, aParams); } if(myFlag == Extrema_ExtFlag_MAX || myFlag == Extrema_ExtFlag_MINMAX) { @@ -627,14 +928,16 @@ void Extrema_GenExtPS::Perform(const gp_Pnt& P) Bnd_SphereUBTreeSelectorMax aSelector(mySphereArray, aSol); //aSelector.SetMaxDist( RealLast() ); aSelector.DefineCheckPoint( P ); - Standard_Integer aNbSel = mySphereUBTree->Select( aSelector ); + mySphereUBTree->Select( aSelector ); //TODO: check if no solution in binary tree Bnd_Sphere& aSph = aSelector.Sphere(); + Standard_Real aU = myUParams->Value(aSph.U()); + Standard_Real aV = myVParams->Value(aSph.V()); + Extrema_POnSurfParams aParams(aU, aV, myS->Value(aU, aV)); + aParams.SetSqrDistance(P.SquareDistance(aParams.Value())); + aParams.SetIndices(aSph.U(), aSph.V()); - UV(1) = U0 + (aSph.U() - 1) * PasU; - UV(2) = V0 + (aSph.V() - 1) * PasV; - - FindSolution(P, UV, PasU, PasV, Extrema_ExtFlag_MAX); + FindSolution(P, aParams); } } }