OpenCascade/occt
1
0
Fork 0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-04-04 18:06:22 +03:00
Code Packages Releases Activity

0032882: Modeling Data - Extrema curve/curve cannot find all solutions (OCCT 7.2 backport)

Browse Source 
Extrema/Extrema_GenExtCC.gxx - estimation of Lipchitz constant is improved
Extrema_GlobOptFuncCC.cxx - function value is changed

LocOpe/LocOpe_WiresOnShape.cxx - small correction to fix regression

lowalgos/extcc/bug32882 - new test case is added

some test were updated according new behavior of extrema algo
This commit is contained in:
ifv 2022-03-21 16:55:55 +03:00 committed by ika
parent 906412a658
commit a0e35c9556
10 changed files with 293 additions and 386 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

193
src/Extrema/Extrema_GenExtCC.gxx
View File

@ -26,6 +26,7 @@
#include <Precision.hxx>
#include <NCollection_Vector.hxx>
#include <NCollection_CellFilter.hxx>
#include <GCPnts_AbscissaPoint.hxx>
// Comparator, used in std::sort.
static Standard_Boolean comp(const gp_XY& theA,
@ -47,6 +48,57 @@ static Standard_Boolean comp(const gp_XY& theA,
return Standard_False;
}
static void ChangeIntervals(Handle(TColStd_HArray1OfReal)& theInts, const Standard_Integer theNbInts)
{
Standard_Integer aNbInts = theInts->Length() - 1;
Standard_Integer aNbAdd = theNbInts - aNbInts;
Handle(TColStd_HArray1OfReal) aNewInts = new TColStd_HArray1OfReal(1, theNbInts + 1);
Standard_Integer aNbLast = theInts->Length();
Standard_Integer i;
if (aNbInts == 1)
{
aNewInts->SetValue(1, theInts->First());
aNewInts->SetValue(theNbInts + 1, theInts->Last());
Standard_Real dt = (theInts->Last() - theInts->First()) / theNbInts;
Standard_Real t = theInts->First() + dt;
for (i = 2; i <= theNbInts; ++i, t += dt)
{
aNewInts->SetValue(i, t);
}
theInts = aNewInts;
return;
}
//
for (i = 1; i <= aNbLast; ++i)
{
aNewInts->SetValue(i, theInts->Value(i));
}
//
while (aNbAdd > 0)
{
Standard_Real anLIntMax = -1.;
Standard_Integer aMaxInd = -1;
for (i = 1; i < aNbLast; ++i)
{
Standard_Real anL = aNewInts->Value(i + 1) - aNewInts->Value(i);
if (anL > anLIntMax)
{
anLIntMax = anL;
aMaxInd = i;
}
}
Standard_Real t = (aNewInts->Value(aMaxInd + 1) + aNewInts->Value(aMaxInd)) / 2.;
for (i = aNbLast; i > aMaxInd; --i)
{
aNewInts->SetValue(i + 1, aNewInts->Value(i));
}
aNbLast++;
aNbAdd--;
aNewInts->SetValue(aMaxInd + 1, t);
}
theInts = aNewInts;
}
class Extrema_CCPointsInspector : public NCollection_CellFilter_InspectorXY
{
public:
@ -111,7 +163,6 @@ static Standard_Real ProjPOnC(const Pnt& theP,
if (aD < aDist)
aDist = aD;
}
aDist = sqrt(aDist);
}
return aDist;
}
@ -229,14 +280,79 @@ void Extrema_GenExtCC::Perform()
aNbInter[1] = C2.NbIntervals(aContinuity);
}
TColStd_Array1OfReal anIntervals1(1, aNbInter[0] + 1);
TColStd_Array1OfReal anIntervals2(1, aNbInter[1] + 1);
C1.Intervals(anIntervals1, aContinuity);
C2.Intervals(anIntervals2, aContinuity);
Standard_Real anL[2];
Standard_Integer indmax = -1, indmin = -1;
const Standard_Real mult = 20.;
if (!(Precision::IsInfinite(C1.FirstParameter()) || Precision::IsInfinite(C1.LastParameter()) ||
Precision::IsInfinite(C2.FirstParameter()) || Precision::IsInfinite(C2.LastParameter())))
{
anL[0] = GCPnts_AbscissaPoint::Length(C1);
anL[1] = GCPnts_AbscissaPoint::Length(C2);
if (anL[0] / aNbInter[0] > mult * anL[1] / aNbInter[1])
{
indmax = 0;
indmin = 1;
}
else if (anL[1] / aNbInter[1] > mult * anL[0] / aNbInter[0])
{
indmax = 1;
indmin = 0;
}
}
Standard_Integer aNbIntOpt = 0;
if (indmax >= 0)
{
aNbIntOpt = RealToInt(anL[indmax] * aNbInter[indmin] / anL[indmin] / (mult / 4.)) + 1;
if (aNbIntOpt > 100 || aNbIntOpt < aNbInter[indmax])
{
indmax = -1;
}
else
{
if (aNbIntOpt * aNbInter[indmin] > 100)
{
aNbIntOpt = 100 / aNbInter[indmin];
if (aNbIntOpt < aNbInter[indmax])
{
indmax = -1;
}
}
}
}
Handle(TColStd_HArray1OfReal) anIntervals1 = new TColStd_HArray1OfReal(1, aNbInter[0] + 1);
Handle(TColStd_HArray1OfReal) anIntervals2 = new TColStd_HArray1OfReal(1, aNbInter[1] + 1);
C1.Intervals(anIntervals1->ChangeArray1(), aContinuity);
C2.Intervals(anIntervals2->ChangeArray1(), aContinuity);
if (indmax >= 0)
{
if (indmax == 0)
{
//Change anIntervals1
ChangeIntervals(anIntervals1, aNbIntOpt);
aNbInter[0] = anIntervals1->Length() - 1;
}
else
{
//Change anIntervals2;
ChangeIntervals(anIntervals2, aNbIntOpt);
aNbInter[1] = anIntervals2->Length() - 1;
}
}
if (C1.IsClosed() && aNbInter[0] == 1)
{
ChangeIntervals(anIntervals1, 3);
aNbInter[0] = anIntervals1->Length() - 1;
}
if (C2.IsClosed() && aNbInter[1] == 1)
{
ChangeIntervals(anIntervals2, 3);
aNbInter[1] = anIntervals2->Length() - 1;
}
// Lipchitz constant computation.
const Standard_Real aMaxLC = 10000.;
Standard_Real aLC = 9.0; // Default value.
Standard_Real aLC = 100.0; // Default value.
const Standard_Real aMaxDer1 = 1.0 / C1.Resolution(1.0);
const Standard_Real aMaxDer2 = 1.0 / C2.Resolution(1.0);
Standard_Real aMaxDer = Max(aMaxDer1, aMaxDer2) * Sqrt(2.0);
@ -276,6 +392,43 @@ void Extrema_GenExtCC::Perform()
}
Extrema_GlobOptFuncCCC2 aFunc (C1, C2);
if (aLC < aMaxLC || aMaxDer > aMaxLC)
{
//Estimation of Lipschitz constant by gradient of optimization function
//using sampling in parameter space.
math_Vector aT(1, 2), aG(1, 2);
Standard_Real aF, aMaxG = 0.;
Standard_Real t1, t2, dt1, dt2;
Standard_Integer n1 = 21, n2 = 21, i1, i2;
dt1 = (C1.LastParameter() - C1.FirstParameter()) / (n1 - 1);
dt2 = (C2.LastParameter() - C2.FirstParameter()) / (n2 - 1);
for (i1 = 1, t1 = C1.FirstParameter(); i1 <= n1; ++i1, t1 += dt1)
{
aT(1) = t1;
for (i2 = 1, t2 = C2.FirstParameter(); i2 <= n2; ++i2, t2 += dt2)
{
aT(2) = t2;
aFunc.Values(aT, aF, aG);
Standard_Real aMod = aG(1)*aG(1) + aG(2)*aG(2);
aMaxG = Max(aMaxG, aMod);
}
}
aMaxG = Sqrt(aMaxG);
if (aMaxG > aMaxDer)
{
aLC = Min(aMaxG, aMaxLC);
isConstLockedFlag = Standard_True;
}
if (aMaxG > 100. * aMaxLC)
{
aLC = 100. * aMaxLC;
isConstLockedFlag = Standard_True;
}
else if (aMaxG < 0.1 * aMaxDer)
{
isConstLockedFlag = Standard_True;
}
}
math_GlobOptMin aFinder(&aFunc, myLowBorder, myUppBorder, aLC);
aFinder.SetLipConstState(isConstLockedFlag);
aFinder.SetContinuity(aContinuity == GeomAbs_C2 ? 2 : 1);
@ -286,8 +439,8 @@ void Extrema_GenExtCC::Perform()
aFinder.SetFunctionalMinimalValue(0.0); // Best distance cannot be lower than 0.0.
// Size computed to have cell index inside of int32 value.
const Standard_Real aCellSize = Max(Max(anIntervals1.Last() - anIntervals1.First(),
anIntervals2.Last() - anIntervals2.First())
const Standard_Real aCellSize = Max(Max(anIntervals1->Last() - anIntervals1->First(),
anIntervals2->Last() - anIntervals2->First())
* Precision::PConfusion() / (2.0 * Sqrt(2.0)),
Precision::PConfusion());
Extrema_CCPointsInspector anInspector(aCellSize);
@ -303,10 +456,10 @@ void Extrema_GenExtCC::Perform()
{
for(j = 1; j <= aNbInter[1]; j++)
{
aFirstBorderInterval(1) = anIntervals1(i);
aFirstBorderInterval(2) = anIntervals2(j);
aSecondBorderInterval(1) = anIntervals1(i + 1);
aSecondBorderInterval(2) = anIntervals2(j + 1);
aFirstBorderInterval(1) = anIntervals1->Value(i);
aFirstBorderInterval(2) = anIntervals2->Value(j);
aSecondBorderInterval(1) = anIntervals1->Value(i + 1);
aSecondBorderInterval(2) = anIntervals2->Value(j + 1);
aFinder.SetLocalParams(aFirstBorderInterval, aSecondBorderInterval);
aFinder.Perform(GetSingleSolutionFlag());
@ -515,6 +668,7 @@ Standard_Boolean Extrema_GenExtCC::IsDone() const
//=======================================================================
Standard_Boolean Extrema_GenExtCC::IsParallel() const
{
if (!IsDone()) throw StdFail_NotDone();
return myParallel;
}
@ -524,8 +678,7 @@ Standard_Boolean Extrema_GenExtCC::IsParallel() const
//=======================================================================
Standard_Integer Extrema_GenExtCC::NbExt() const
{
StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::NbExt()")
if (!IsDone()) throw StdFail_NotDone();
return myPoints1.Length();
}
@ -535,8 +688,10 @@ Standard_Integer Extrema_GenExtCC::NbExt() const
//=======================================================================
Standard_Real Extrema_GenExtCC::SquareDistance(const Standard_Integer N) const
{
StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::SquareDistance()")
Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::SquareDistance()")
if (N < 1 || N > NbExt())
{
throw Standard_OutOfRange();
}
return Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)).SquareDistance(Tool2::Value(*((Curve2*)myC[1]), myPoints2(N)));
}
@ -549,8 +704,10 @@ void Extrema_GenExtCC::Points(const Standard_Integer N,
POnC& P1,
POnC& P2) const
{
StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::Points()")
Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::Points()")
if (N < 1 || N > NbExt())
{
throw Standard_OutOfRange();
}
P1.SetValues(myPoints1(N), Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)));
P2.SetValues(myPoints2(N), Tool2::Value(*((Curve2*)myC[1]), myPoints2(N)));

15
src/Extrema/Extrema_GlobOptFuncCC.cxx
View File

@ -45,7 +45,7 @@ static Standard_Boolean _Value(const Adaptor3d_Curve& C1,
return Standard_False;
}
F = C2.Value(v).Distance(C1.Value(u));
F = C2.Value(v).SquareDistance(C1.Value(u));
return Standard_True;
}
@ -66,7 +66,7 @@ static Standard_Boolean _Value(const Adaptor2d_Curve2d& C1,
return Standard_False;
}
F = C2.Value(v).Distance(C1.Value(u));
F = C2.Value(v).SquareDistance(C1.Value(u));
return Standard_True;
}
@ -91,13 +91,14 @@ static Standard_Boolean _Gradient(const Adaptor3d_Curve& C1,
C1.D1(X(1), C1D0, C1D1);
C2.D1(X(2), C2D0, C2D1);
G(1) = - (C2D0.X() - C1D0.X()) * C1D1.X()
- (C2D0.Y() - C1D0.Y()) * C1D1.Y()
- (C2D0.Z() - C1D0.Z()) * C1D1.Z();
G(2) = (C2D0.X() - C1D0.X()) * C2D1.X()
+ (C2D0.Y() - C1D0.Y()) * C2D1.Y()
+ (C2D0.Z() - C1D0.Z()) * C2D1.Z();
G *= 2.;
return Standard_True;
}
@ -123,8 +124,11 @@ static Standard_Boolean _Gradient(const Adaptor2d_Curve2d& C1,
G(1) = - (C2D0.X() - C1D0.X()) * C1D1.X()
- (C2D0.Y() - C1D0.Y()) * C1D1.Y();
G(2) = (C2D0.X() - C1D0.X()) * C2D1.X()
+ (C2D0.Y() - C1D0.Y()) * C2D1.Y();
G *= 2.;
return Standard_True;
}
@ -168,6 +172,7 @@ static Standard_Boolean _Hessian (const Adaptor3d_Curve& C1,
+ (C2D0.X() - C1D0.X()) * C2D2.X()
+ (C2D0.Y() - C1D0.Y()) * C2D2.Y()
+ (C2D0.Z() - C1D0.Z()) * C2D2.Z();
H *= 2.;
return Standard_True;
}
@ -206,10 +211,11 @@ static Standard_Boolean _Hessian (const Adaptor2d_Curve2d& C1,
+ C2D1.Y() * C2D1.Y()
+ (C2D0.X() - C1D0.X()) * C2D2.X()
+ (C2D0.Y() - C1D0.Y()) * C2D2.Y();
H *= 2.;
return Standard_True;
}
// C0
//C0
//=======================================================================
//function : Extrema_GlobOptFuncCCC0
@ -405,6 +411,5 @@ Standard_Boolean Extrema_GlobOptFuncCCC2::Values(const math_Vector& X,Standard_R
else
isHessianComputed = _Hessian(*myC1_2d, *myC2_2d, X, H);
return (Value(X, F) && Gradient(X, G) && isHessianComputed);
}

4
src/IntTools/IntTools_FaceFace.cxx
View File

@ -1188,7 +1188,7 @@ void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
bAvoidLineConstructor,
myTol,
aSeqOfL,
aReachedTol,
aReachedTol, // obsolete parameter
myContext);
//
aNbSeqOfL=aSeqOfL.Length();
@ -1196,7 +1196,7 @@ void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
Standard_Real aTolC = 0.;
if (bIsDecomposited) {
nbiter=aNbSeqOfL;
aTolC = aReachedTol;
aTolC = Precision::Confusion();
}
else {
nbiter=1;

355
src/IntTools/IntTools_WLineTool.cxx
View File

@ -229,203 +229,6 @@ Standard_Boolean IntTools_WLineTool::NotUseSurfacesForApprox(const TopoDS_Face&
/////////////////////// DecompositionOfWLine ////////////////////////////
//=======================================================================
//function : CheckTangentZonesExist
//purpose : static subfunction in ComputeTangentZones
//=======================================================================
static
Standard_Boolean CheckTangentZonesExist(const Handle(GeomAdaptor_HSurface)& theSurface1,
const Handle(GeomAdaptor_HSurface)& theSurface2)
{
if ( ( theSurface1->GetType() != GeomAbs_Torus ) ||
( theSurface2->GetType() != GeomAbs_Torus ) )
return Standard_False;
gp_Torus aTor1 = theSurface1->Torus();
gp_Torus aTor2 = theSurface2->Torus();
if ( aTor1.Location().Distance( aTor2.Location() ) > Precision::Confusion() )
return Standard_False;
if ( ( fabs( aTor1.MajorRadius() - aTor2.MajorRadius() ) > Precision::Confusion() ) ||
( fabs( aTor1.MinorRadius() - aTor2.MinorRadius() ) > Precision::Confusion() ) )
return Standard_False;
if ( ( aTor1.MajorRadius() < aTor1.MinorRadius() ) ||
( aTor2.MajorRadius() < aTor2.MinorRadius() ) )
return Standard_False;
return Standard_True;
}
//=======================================================================
//function : ComputeTangentZones
//purpose : static subfunction in DecompositionOfWLine
//=======================================================================
static
Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
const Handle(GeomAdaptor_HSurface)& theSurface2,
const TopoDS_Face& theFace1,
const TopoDS_Face& theFace2,
Handle(TColgp_HArray1OfPnt2d)& theResultOnS1,
Handle(TColgp_HArray1OfPnt2d)& theResultOnS2,
Handle(TColStd_HArray1OfReal)& theResultRadius,
const Handle(IntTools_Context)& aContext)
{
Standard_Integer aResult = 0;
if ( !CheckTangentZonesExist( theSurface1, theSurface2 ) )
return aResult;
TColgp_SequenceOfPnt2d aSeqResultS1, aSeqResultS2;
TColStd_SequenceOfReal aSeqResultRad;
gp_Torus aTor1 = theSurface1->Torus();
gp_Torus aTor2 = theSurface2->Torus();
gp_Ax2 anax1( aTor1.Location(), aTor1.Axis().Direction() );
gp_Ax2 anax2( aTor2.Location(), aTor2.Axis().Direction() );
Standard_Integer j = 0;
for ( j = 0; j < 2; j++ ) {
Standard_Real aCoef = ( j == 0 ) ? -1 : 1;
Standard_Real aRadius1 = fabs(aTor1.MajorRadius() + aCoef * aTor1.MinorRadius());
Standard_Real aRadius2 = fabs(aTor2.MajorRadius() + aCoef * aTor2.MinorRadius());
gp_Circ aCircle1( anax1, aRadius1 );
gp_Circ aCircle2( anax2, aRadius2 );
// roughly compute radius of tangent zone for perpendicular case
Standard_Real aCriteria = Precision::Confusion() * 0.5;
Standard_Real aT1 = aCriteria;
Standard_Real aT2 = aCriteria;
if ( j == 0 ) {
// internal tangency
Standard_Real aR = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
//aT1 = aCriteria * aCriteria + aR * aR - ( aR - aCriteria ) * ( aR - aCriteria );
aT1 = 2. * aR * aCriteria;
aT2 = aT1;
}
else {
// external tangency
Standard_Real aRb = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
Standard_Real aRm = ( aRadius1 < aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
Standard_Real aDelta = aRb - aCriteria;
aDelta *= aDelta;
aDelta -= aRm * aRm;
aDelta /= 2. * (aRb - aRm);
aDelta -= 0.5 * (aRb - aRm);
aT1 = 2. * aRm * (aRm - aDelta);
aT2 = aT1;
}
aCriteria = ( aT1 > aT2) ? aT1 : aT2;
if ( aCriteria > 0 )
aCriteria = sqrt( aCriteria );
if ( aCriteria > 0.5 * aTor1.MinorRadius() ) {
// too big zone -> drop to minimum
aCriteria = Precision::Confusion();
}
GeomAdaptor_Curve aC1( new Geom_Circle(aCircle1) );
GeomAdaptor_Curve aC2( new Geom_Circle(aCircle2) );
Extrema_ExtCC anExtrema(aC1, aC2, 0, 2. * M_PI, 0, 2. * M_PI,
Precision::PConfusion(), Precision::PConfusion());
if ( anExtrema.IsDone() ) {
Standard_Integer i = 0;
for ( i = 1; i <= anExtrema.NbExt(); i++ ) {
if ( anExtrema.SquareDistance(i) > aCriteria * aCriteria )
continue;
Extrema_POnCurv P1, P2;
anExtrema.Points( i, P1, P2 );
Standard_Boolean bFoundResult = Standard_True;
gp_Pnt2d pr1, pr2;
Standard_Integer surfit = 0;
for ( surfit = 0; surfit < 2; surfit++ ) {
GeomAPI_ProjectPointOnSurf& aProjector =
(surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
gp_Pnt aP3d = (surfit == 0) ? P1.Value() : P2.Value();
aProjector.Perform(aP3d);
if(!aProjector.IsDone())
bFoundResult = Standard_False;
else {
if(aProjector.LowerDistance() > aCriteria) {
bFoundResult = Standard_False;
}
else {
Standard_Real foundU = 0, foundV = 0;
aProjector.LowerDistanceParameters(foundU, foundV);
if ( surfit == 0 )
pr1 = gp_Pnt2d( foundU, foundV );
else
pr2 = gp_Pnt2d( foundU, foundV );
}
}
}
if ( bFoundResult ) {
aSeqResultS1.Append( pr1 );
aSeqResultS2.Append( pr2 );
aSeqResultRad.Append( aCriteria );
// torus is u and v periodic
const Standard_Real twoPI = M_PI + M_PI;
Standard_Real arr1tmp[2] = {pr1.X(), pr1.Y()};
Standard_Real arr2tmp[2] = {pr2.X(), pr2.Y()};
// iteration on period bounds
for ( Standard_Integer k1 = 0; k1 < 2; k1++ ) {
Standard_Real aBound = ( k1 == 0 ) ? 0 : twoPI;
Standard_Real aShift = ( k1 == 0 ) ? twoPI : -twoPI;
// iteration on surfaces
for ( Standard_Integer k2 = 0; k2 < 2; k2++ ) {
Standard_Real* arr1 = ( k2 == 0 ) ? arr1tmp : arr2tmp;
Standard_Real* arr2 = ( k2 != 0 ) ? arr1tmp : arr2tmp;
TColgp_SequenceOfPnt2d& aSeqS1 = ( k2 == 0 ) ? aSeqResultS1 : aSeqResultS2;
TColgp_SequenceOfPnt2d& aSeqS2 = ( k2 != 0 ) ? aSeqResultS1 : aSeqResultS2;
if (fabs(arr1[0] - aBound) < Precision::PConfusion()) {
aSeqS1.Append( gp_Pnt2d( arr1[0] + aShift, arr1[1] ) );
aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
aSeqResultRad.Append( aCriteria );
}
if (fabs(arr1[1] - aBound) < Precision::PConfusion()) {
aSeqS1.Append( gp_Pnt2d( arr1[0], arr1[1] + aShift) );
aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
aSeqResultRad.Append( aCriteria );
}
}
} //
}
}
}
}
aResult = aSeqResultRad.Length();
if ( aResult > 0 ) {
theResultOnS1 = new TColgp_HArray1OfPnt2d( 1, aResult );
theResultOnS2 = new TColgp_HArray1OfPnt2d( 1, aResult );
theResultRadius = new TColStd_HArray1OfReal( 1, aResult );
for ( Standard_Integer i = 1 ; i <= aResult; i++ ) {
theResultOnS1->SetValue( i, aSeqResultS1.Value(i) );
theResultOnS2->SetValue( i, aSeqResultS2.Value(i) );
theResultRadius->SetValue( i, aSeqResultRad.Value(i) );
}
}
return aResult;
}
//=======================================================================
//function : IsPointOnBoundary
//purpose : static subfunction in DecompositionOfWLine
@ -457,27 +260,6 @@ Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
return bRet;
}
//=======================================================================
//function : IsInsideTanZone
//purpose : Check if point is inside a radial tangent zone.
// static subfunction in DecompositionOfWLine and FindPoint
//=======================================================================
static
Standard_Boolean IsInsideTanZone(const gp_Pnt2d& thePoint,
const gp_Pnt2d& theTanZoneCenter,
const Standard_Real theZoneRadius,
Handle(GeomAdaptor_HSurface) theGASurface)
{
Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
Standard_Real aRadiusSQR = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
aRadiusSQR *= aRadiusSQR;
if ( thePoint.SquareDistance( theTanZoneCenter ) <= aRadiusSQR )
return Standard_True;
return Standard_False;
}
//=======================================================================
//function : AdjustByNeighbour
//purpose : static subfunction in DecompositionOfWLine
@ -652,72 +434,6 @@ Standard_Boolean FindPoint(const gp_Pnt2d& theFirstPoint,
return Standard_False;
}
//=======================================================================
//function : FindPoint
//purpose : Find point on the boundary of radial tangent zone
// static subfunction in DecompositionOfWLine
//=======================================================================
static
Standard_Boolean FindPoint(const gp_Pnt2d& theFirstPoint,
const gp_Pnt2d& theLastPoint,
const Standard_Real theUmin,
const Standard_Real theUmax,
const Standard_Real theVmin,
const Standard_Real theVmax,
const gp_Pnt2d& theTanZoneCenter,
const Standard_Real theZoneRadius,
Handle(GeomAdaptor_HSurface) theGASurface,
gp_Pnt2d& theNewPoint) {
theNewPoint = theLastPoint;
if ( !IsInsideTanZone( theLastPoint, theTanZoneCenter, theZoneRadius, theGASurface) )
return Standard_False;
Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
Standard_Real aRadius = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
gp_Ax22d anAxis( theTanZoneCenter, gp_Dir2d(1, 0), gp_Dir2d(0, 1) );
gp_Circ2d aCircle( anAxis, aRadius );
//
gp_Vec2d aDir( theLastPoint.XY() - theFirstPoint.XY() );
Standard_Real aLength = aDir.Magnitude();
if ( aLength <= gp::Resolution() )
return Standard_False;
gp_Lin2d aLine( theFirstPoint, aDir );
//
Handle(Geom2d_Line) aCLine = new Geom2d_Line( aLine );
Handle(Geom2d_TrimmedCurve) aC1 = new Geom2d_TrimmedCurve( aCLine, 0, aLength );
Handle(Geom2d_Circle) aC2 = new Geom2d_Circle( aCircle );
Standard_Real aTol = aRadius * 0.001;
aTol = ( aTol < Precision::PConfusion() ) ? Precision::PConfusion() : aTol;
Geom2dAPI_InterCurveCurve anIntersector;
anIntersector.Init( aC1, aC2, aTol );
if ( anIntersector.NbPoints() == 0 )
return Standard_False;
Standard_Boolean aFound = Standard_False;
Standard_Real aMinDist = aLength * aLength;
Standard_Integer i = 0;
for ( i = 1; i <= anIntersector.NbPoints(); i++ ) {
gp_Pnt2d aPInt = anIntersector.Point( i );
if ( aPInt.SquareDistance( theFirstPoint ) < aMinDist ) {
if ( ( aPInt.X() >= theUmin ) && ( aPInt.X() <= theUmax ) &&
( aPInt.Y() >= theVmin ) && ( aPInt.Y() <= theVmax ) ) {
theNewPoint = aPInt;
aFound = Standard_True;
}
}
}
return aFound;
}
//=======================================================================
//function : DecompositionOfWLine
//purpose :
@ -732,7 +448,7 @@ Standard_Boolean IntTools_WLineTool::
const Standard_Boolean theAvoidLConstructor,
const Standard_Real theTol,
IntPatch_SequenceOfLine& theNewLines,
Standard_Real& theReachedTol3d,
Standard_Real& /*theReachedTol3d*/,
const Handle(IntTools_Context)& aContext)
{
Standard_Boolean bRet, bAvoidLineConstructor;
@ -758,13 +474,7 @@ Standard_Boolean IntTools_WLineTool::
TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts);
TColStd_Array1OfInteger anArrayOfLineType(1, aNbPnts);
TColStd_ListOfInteger aListOfPointIndex;
Handle(TColgp_HArray1OfPnt2d) aTanZoneS1;
Handle(TColgp_HArray1OfPnt2d) aTanZoneS2;
Handle(TColStd_HArray1OfReal) aTanZoneRadius;
Standard_Integer aNbZone = ComputeTangentZones( theSurface1, theSurface2, theFace1, theFace2,
aTanZoneS1, aTanZoneS2, aTanZoneRadius, aContext);
//
nblines=0;
aTol=Precision::Confusion();
@ -835,24 +545,6 @@ Standard_Boolean IntTools_WLineTool::
bIsCurrentPointOnBoundary = Standard_True;
break;
}
else {
// check if a point belong to a tangent zone. Begin
Standard_Integer zIt = 0;
for ( zIt = 1; zIt <= aNbZone; zIt++ ) {
gp_Pnt2d aPZone = (i == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
if ( IsInsideTanZone(gp_Pnt2d( U, V ), aPZone, aZoneRadius, aGASurface ) ) {
// set boundary flag to split the curve by a tangent zone
bIsPointOnBoundary = Standard_True;
bIsCurrentPointOnBoundary = Standard_True;
if ( theReachedTol3d < aZoneRadius ) {
theReachedTol3d = aZoneRadius;
}
break;
}
}
}
}//for(j = 0; j < 2; j++) {
if(bIsCurrentPointOnBoundary){
@ -931,7 +623,7 @@ Standard_Boolean IntTools_WLineTool::
Standard_Integer nbboundaries = 0;
Standard_Boolean bIsNearBoundary = Standard_False;
Standard_Integer aZoneIndex = 0;
//Standard_Integer aZoneIndex = 0;
Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
@ -981,32 +673,6 @@ Standard_Boolean IntTools_WLineTool::
}
}
// check if a point belong to a tangent zone. Begin
for ( Standard_Integer zIt = 1; zIt <= aNbZone; zIt++ ) {
gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
Standard_Real nU1, nV1;
if(surfit == 0)
aNeighbourPoint.ParametersOnS1(nU1, nV1);
else
aNeighbourPoint.ParametersOnS2(nU1, nV1);
gp_Pnt2d ap1(nU1, nV1);
gp_Pnt2d ap2 = AdjustByNeighbour( ap1, gp_Pnt2d( U, V ), aGASurface );
if ( IsInsideTanZone( ap2, aPZone, aZoneRadius, aGASurface ) ) {
aZoneIndex = zIt;
bIsNearBoundary = Standard_True;
if ( theReachedTol3d < aZoneRadius ) {
theReachedTol3d = aZoneRadius;
}
}
}
// check if a point belong to a tangent zone. End
Standard_Boolean bComputeLineEnd = Standard_False;
if(nbboundaries == 2) {
@ -1145,20 +811,7 @@ Standard_Boolean IntTools_WLineTool::
gp_Pnt2d ap1(nU1, nV1);
gp_Pnt2d ap2;
if ( aZoneIndex ) {
// exclude point from a tangent zone
anewpoint = AdjustByNeighbour( ap1, anewpoint, aGASurface );
gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(aZoneIndex) : aTanZoneS2->Value(aZoneIndex);
Standard_Real aZoneRadius = aTanZoneRadius->Value(aZoneIndex);
if ( FindPoint(ap1, anewpoint, umin, umax, vmin, vmax,
aPZone, aZoneRadius, aGASurface, ap2) ) {
anewpoint = ap2;
found = Standard_True;
}
}
else if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) {
if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) {
// re-compute point near boundary if shifted on a period
ap2 = AdjustByNeighbour( ap1, anewpoint, aGASurface );

34
src/LocOpe/LocOpe_WiresOnShape.cxx
View File

@ -1300,6 +1300,10 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge,
Standard_Real aTolV[2];
aTolV[0] =BRep_Tool::Tolerance(theVertices[0]);
aTolV[1] =BRep_Tool::Tolerance(theVertices[1]);
Standard_Real ext = 16.; // = 4 * 4 - to avoid creating microedges, area around vertices is increased
// up to 4 vertex tolerance. Such approach is usual for other topological
// algorithms, for example, Boolean Operations.
Standard_Real aTolVExt[2] = { ext * aTolV[0] * aTolV[0], ext * aTolV[1] * aTolV[1] };
BRepAdaptor_Curve2d thePCurve(theEdge, theFace);
Bnd_Box2d theBox;
@ -1309,6 +1313,9 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge,
Standard_Real aFpar, aLpar;
const Handle(Geom_Curve)& theCurve = BRep_Tool::Curve(theEdge, thePar[0], thePar[1]);
GeomAdaptor_Curve theGAcurve(theCurve, thePar[0], thePar[1]);
Standard_Real aTolV2d[2] = { theGAcurve.Resolution(aTolV[0]), theGAcurve.Resolution(aTolV[1]) };
aTolV2d[0] = Max(aTolV2d[0], Precision::PConfusion());
aTolV2d[1] = Max(aTolV2d[1], Precision::PConfusion());
Standard_Real aDistMax = Precision::Confusion() * Precision::Confusion();
TopExp_Explorer Explo(theFace, TopAbs_EDGE);
for (; Explo.More(); Explo.Next())
@ -1335,6 +1342,23 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge,
isOverlapped = Standard_True;
return;
}
// Check extremity distances
Standard_Real dists[4];
gp_Pnt aP11, aP12, aP21, aP22;
anExtrema.TrimmedSquareDistances(dists[0], dists[1], dists[2], dists[3],
aP11, aP12, aP21, aP22);
for (i = 0; i < 4; ++i)
{
if (i < 2)
j = 0;
else
j = 1;
if (dists[i] < aTolVExt[j] / ext)
{
return;
}
}
for (i = 1; i <= aNbExt; i++)
{
Standard_Real aDist = anExtrema.SquareDistance(i);
@ -1347,7 +1371,7 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge,
Standard_Real anIntPar = aPOnC2.Parameter();
for (j = 0; j < 2; j++) //try to find intersection on an extremity of "theEdge"
{
if (Abs(theIntPar - thePar[j]) <= Precision::PConfusion())
if (Abs(theIntPar - thePar[j]) <= aTolV2d[j])
break;
}
//intersection found in the middle of the edge
@ -1356,10 +1380,10 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge,
gp_Pnt aPoint = aCurve->Value(anIntPar);
gp_Pnt aPointInt = theCurve->Value(theIntPar);
if (aPointInt.SquareDistance(thePnt[0]) > aTolV[0] * aTolV[0] &&
aPointInt.SquareDistance(thePnt[1]) > aTolV[1] * aTolV[1] &&
aPoint.SquareDistance(thePnt[0]) > aTolV[0] * aTolV[0] &&
aPoint.SquareDistance(thePnt[1]) > aTolV[1] * aTolV[1])
if (aPointInt.SquareDistance(thePnt[0]) > aTolVExt[0] &&
aPointInt.SquareDistance(thePnt[1]) > aTolVExt[1] &&
aPoint.SquareDistance(thePnt[0]) > aTolVExt[0] &&
aPoint.SquareDistance(thePnt[1]) > aTolVExt[1])
{
SplitPars.Append(theIntPar);
if( aDist > aDistMax)

2
tests/bugs/modalg_5/bug23706_10
View File

@ -13,7 +13,7 @@ set info [extrema r3 r4]
if {[regexp "ext_1" $info]} {
set dist [lindex [length ext_1] end]
if { $dist > 4.0e-13 } {
if { $dist > 5.0e-11 } {
puts "Error: Extrema distance is too big"
}
} else {

9
tests/bugs/modalg_5/bug23706_11
View File

@ -10,8 +10,11 @@ bsplinecurve r1 2 5 1 3 2 1 3 1 4 1 5 3 2 5 3 1 3 7 3 1 4 8 3 1 4 8 3 1 5 9 3 1
bsplinecurve r2 2 5 2 3 2.5 1 3 1 3.5 1 4 3 -1 2 3 1 1 11 3 1 3 9 3 1 3 9 3 1 3 9 3 1 5 7 3 1 7 4 3 1
set info [extrema r1 r2]
if { [llength $info] != 1 } {
puts "Error : Extrema is wrong"
if {[regexp "ext_1" $info]} {
set dist [lindex [length ext_1] end]
if { $dist > 1.0e-8 } {
puts "Error: Extrema distance is too big"
}
} else {
puts "OK: Extrema is valid"
puts "Error: Extrema is not found"
}

2
tests/bugs/modalg_6/bug27665
View File

@ -5,7 +5,7 @@ puts ""
#################################################
# BrepExrtrema_DistShapeShape bad performance on OCCT 6.7.0
#################################################
cpulimit 100
cpulimit 500
restore [locate_data_file bug27665_wircmpd.brep] w
explode w

29
tests/lowalgos/extcc/bug29858_03 Normal file
View File

@ -0,0 +1,29 @@
puts "========"
puts "OCC29858"
puts "========"
puts ""
#################################################
# Regression in GeomAPI_ExtremaCurveCurve
#################################################
# Read input
restore [locate_data_file bug29858_03_e1.brep] e1
restore [locate_data_file bug29858_03_e2.brep] e2
# Extract geometry from topology
mkcurve c1 e1
mkcurve c2 e2
# Run extrema
set info [extrema c1 c2]
# Check result
if {[regexp "ext_1" $info]} {
set dist [lindex [length ext_1] end]
if { $dist > 1.0e-10 } {
puts "Error: Extrema distance is too big"
}
} else {
puts "Error: Extrema is not found"
}

36
tests/lowalgos/extcc/bug32882 Normal file
View File

@ -0,0 +1,36 @@
puts "========"
puts "OCC32882: Modeling Data - Extrema curve/curve cannot find all solutions"
puts "========"
puts ""
# Read input
restore [locate_data_file bug32882.brep] cc
explode cc
# Extract geometry from topology
mkcurve c1 cc_1
mkcurve c2 cc_2
mkcurve c3 cc_3
# Run extrema c1/c3
set info [extrema c1 c3]
# Check number of solution
if { [llength $info] != 3 } {
puts "Error: Invalid extrema number in extrema c1-c3 output"
}
# Check result
checklength ext_1 -l 2.929642751e-14 -eps .01
checklength ext_2 -l 3.480934286e-14 -eps .01
checklength ext_3 -l 3.177643716e-14 -eps .01
# Run extrema c3/c2
set info [extrema c3 c2]
# Check number of solutions
if { [llength $info] != 2 } {
puts "Error: Invalid extrema number in extrema c3-c2 output"
}
checklength ext_1 -l 5.684341886e-14 -eps .01
checklength ext_2 -l 2.929642751e-14 -eps .01