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

0032882: Modeling Data - Extrema curve/curve cannot find all solutions

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 smoskvin
parent 7090725e2b
commit d7f5072158
11 changed files with 144 additions and 381 deletions

51
src/Extrema/Extrema_GenExtCC.gxx

@ -163,7 +163,6 @@ static Standard_Real ProjPOnC(const Pnt& theP,
if (aD < aDist) if (aD < aDist)
aDist = aD; aDist = aD;
} }
aDist = sqrt(aDist);
} }
return aDist; return aDist;
} }
@ -340,10 +339,20 @@ void Extrema_GenExtCC::Perform()
aNbInter[1] = anIntervals2->Length() - 1; 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. // Lipchitz constant computation.
const Standard_Real aMaxLC = 10000.; 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 aMaxDer1 = 1.0 / C1.Resolution(1.0);
const Standard_Real aMaxDer2 = 1.0 / C2.Resolution(1.0); const Standard_Real aMaxDer2 = 1.0 / C2.Resolution(1.0);
Standard_Real aMaxDer = Max(aMaxDer1, aMaxDer2) * Sqrt(2.0); Standard_Real aMaxDer = Max(aMaxDer1, aMaxDer2) * Sqrt(2.0);
@ -383,6 +392,43 @@ void Extrema_GenExtCC::Perform()
} }
Extrema_GlobOptFuncCCC2 aFunc (C1, C2); 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); math_GlobOptMin aFinder(&aFunc, myLowBorder, myUppBorder, aLC);
aFinder.SetLipConstState(isConstLockedFlag); aFinder.SetLipConstState(isConstLockedFlag);
aFinder.SetContinuity(aContinuity == GeomAbs_C2 ? 2 : 1); aFinder.SetContinuity(aContinuity == GeomAbs_C2 ? 2 : 1);
@ -518,7 +564,6 @@ void Extrema_GenExtCC::Perform()
aVec(2) = (aCurrent.Y() + aNext.Y()) * 0.5; aVec(2) = (aCurrent.Y() + aNext.Y()) * 0.5;
aFunc.Value(aVec, aVal); aFunc.Value(aVec, aVal);
if (Abs(aVal - aF) < Precision::Confusion()) if (Abs(aVal - aF) < Precision::Confusion())
{ {
// It seems the parallel segment is found. // It seems the parallel segment is found.

15
src/Extrema/Extrema_GlobOptFuncCC.cxx

@ -43,7 +43,7 @@ static Standard_Boolean _Value(const Adaptor3d_Curve& C1,
return Standard_False; return Standard_False;
} }
F = C2.Value(v).Distance(C1.Value(u)); F = C2.Value(v).SquareDistance(C1.Value(u));
return Standard_True; return Standard_True;
} }
@ -64,7 +64,7 @@ static Standard_Boolean _Value(const Adaptor2d_Curve2d& C1,
return Standard_False; return Standard_False;
} }
F = C2.Value(v).Distance(C1.Value(u)); F = C2.Value(v).SquareDistance(C1.Value(u));
return Standard_True; return Standard_True;
} }
@ -89,13 +89,14 @@ static Standard_Boolean _Gradient(const Adaptor3d_Curve& C1,
C1.D1(X(1), C1D0, C1D1); C1.D1(X(1), C1D0, C1D1);
C2.D1(X(2), C2D0, C2D1); C2.D1(X(2), C2D0, C2D1);
G(1) = - (C2D0.X() - C1D0.X()) * C1D1.X() G(1) = - (C2D0.X() - C1D0.X()) * C1D1.X()
- (C2D0.Y() - C1D0.Y()) * C1D1.Y() - (C2D0.Y() - C1D0.Y()) * C1D1.Y()
- (C2D0.Z() - C1D0.Z()) * C1D1.Z(); - (C2D0.Z() - C1D0.Z()) * C1D1.Z();
G(2) = (C2D0.X() - C1D0.X()) * C2D1.X() G(2) = (C2D0.X() - C1D0.X()) * C2D1.X()
+ (C2D0.Y() - C1D0.Y()) * C2D1.Y() + (C2D0.Y() - C1D0.Y()) * C2D1.Y()
+ (C2D0.Z() - C1D0.Z()) * C2D1.Z(); + (C2D0.Z() - C1D0.Z()) * C2D1.Z();
G *= 2.;
return Standard_True; return Standard_True;
} }
@ -121,8 +122,11 @@ static Standard_Boolean _Gradient(const Adaptor2d_Curve2d& C1,
G(1) = - (C2D0.X() - C1D0.X()) * C1D1.X() G(1) = - (C2D0.X() - C1D0.X()) * C1D1.X()
- (C2D0.Y() - C1D0.Y()) * C1D1.Y(); - (C2D0.Y() - C1D0.Y()) * C1D1.Y();
G(2) = (C2D0.X() - C1D0.X()) * C2D1.X() G(2) = (C2D0.X() - C1D0.X()) * C2D1.X()
+ (C2D0.Y() - C1D0.Y()) * C2D1.Y(); + (C2D0.Y() - C1D0.Y()) * C2D1.Y();
G *= 2.;
return Standard_True; return Standard_True;
} }
@ -166,6 +170,7 @@ static Standard_Boolean _Hessian (const Adaptor3d_Curve& C1,
+ (C2D0.X() - C1D0.X()) * C2D2.X() + (C2D0.X() - C1D0.X()) * C2D2.X()
+ (C2D0.Y() - C1D0.Y()) * C2D2.Y() + (C2D0.Y() - C1D0.Y()) * C2D2.Y()
+ (C2D0.Z() - C1D0.Z()) * C2D2.Z(); + (C2D0.Z() - C1D0.Z()) * C2D2.Z();
H *= 2.;
return Standard_True; return Standard_True;
} }
@ -204,10 +209,11 @@ static Standard_Boolean _Hessian (const Adaptor2d_Curve2d& C1,
+ C2D1.Y() * C2D1.Y() + C2D1.Y() * C2D1.Y()
+ (C2D0.X() - C1D0.X()) * C2D2.X() + (C2D0.X() - C1D0.X()) * C2D2.X()
+ (C2D0.Y() - C1D0.Y()) * C2D2.Y(); + (C2D0.Y() - C1D0.Y()) * C2D2.Y();
H *= 2.;
return Standard_True; return Standard_True;
} }
// C0 //C0
//======================================================================= //=======================================================================
//function : Extrema_GlobOptFuncCCC0 //function : Extrema_GlobOptFuncCCC0
@ -415,6 +421,5 @@ Standard_Boolean Extrema_GlobOptFuncCCC2::Values(const math_Vector& X,Standard_R
else else
isHessianComputed = _Hessian(*myC1_2d, *myC2_2d, X, H); isHessianComputed = _Hessian(*myC1_2d, *myC2_2d, X, H);
return (Value(X, F) && Gradient(X, G) && isHessianComputed); return (Value(X, F) && Gradient(X, G) && isHessianComputed);
} }

4
src/IntTools/IntTools_FaceFace.cxx

@ -1169,7 +1169,6 @@ void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
tol2d = myTolApprox; tol2d = myTolApprox;
} }
Standard_Real aReachedTol = Precision::Confusion();
bIsDecomposited = IntTools_WLineTool:: bIsDecomposited = IntTools_WLineTool::
DecompositionOfWLine(WL, DecompositionOfWLine(WL,
myHS1, myHS1,
@ -1180,7 +1179,6 @@ void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
bAvoidLineConstructor, bAvoidLineConstructor,
myTol, myTol,
aSeqOfL, aSeqOfL,
aReachedTol,
myContext); myContext);
// //
aNbSeqOfL=aSeqOfL.Length(); aNbSeqOfL=aSeqOfL.Length();
@ -1188,7 +1186,7 @@ void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
Standard_Real aTolC = 0.; Standard_Real aTolC = 0.;
if (bIsDecomposited) { if (bIsDecomposited) {
nbiter=aNbSeqOfL; nbiter=aNbSeqOfL;
aTolC = aReachedTol; aTolC = Precision::Confusion();
} }
else { else {
nbiter=1; nbiter=1;

354
src/IntTools/IntTools_WLineTool.cxx

@ -228,203 +228,6 @@ Standard_Boolean IntTools_WLineTool::NotUseSurfacesForApprox(const TopoDS_Face&
/////////////////////// DecompositionOfWLine //////////////////////////// /////////////////////// DecompositionOfWLine ////////////////////////////
//=======================================================================
//function : CheckTangentZonesExist
//purpose : static subfunction in ComputeTangentZones
//=======================================================================
static
Standard_Boolean CheckTangentZonesExist(const Handle(GeomAdaptor_Surface)& theSurface1,
const Handle(GeomAdaptor_Surface)& 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_Surface)& theSurface1,
const Handle(GeomAdaptor_Surface)& 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 //function : IsPointOnBoundary
//purpose : static subfunction in DecompositionOfWLine //purpose : static subfunction in DecompositionOfWLine
@ -456,27 +259,6 @@ Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
return bRet; 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_Surface) 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 //function : AdjustByNeighbour
//purpose : static subfunction in DecompositionOfWLine //purpose : static subfunction in DecompositionOfWLine
@ -651,72 +433,6 @@ Standard_Boolean FindPoint(const gp_Pnt2d& theFirstPoint,
return Standard_False; 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_Surface) 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 //function : DecompositionOfWLine
//purpose : //purpose :
@ -731,7 +447,6 @@ Standard_Boolean IntTools_WLineTool::
const Standard_Boolean theAvoidLConstructor, const Standard_Boolean theAvoidLConstructor,
const Standard_Real theTol, const Standard_Real theTol,
IntPatch_SequenceOfLine& theNewLines, IntPatch_SequenceOfLine& theNewLines,
Standard_Real& theReachedTol3d,
const Handle(IntTools_Context)& aContext) const Handle(IntTools_Context)& aContext)
{ {
Standard_Boolean bRet, bAvoidLineConstructor; Standard_Boolean bRet, bAvoidLineConstructor;
@ -757,13 +472,7 @@ Standard_Boolean IntTools_WLineTool::
TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts); TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts);
TColStd_Array1OfInteger anArrayOfLineType(1, aNbPnts); TColStd_Array1OfInteger anArrayOfLineType(1, aNbPnts);
TColStd_ListOfInteger aListOfPointIndex; 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; nblines=0;
aTol=Precision::Confusion(); aTol=Precision::Confusion();
@ -834,24 +543,6 @@ Standard_Boolean IntTools_WLineTool::
bIsCurrentPointOnBoundary = Standard_True; bIsCurrentPointOnBoundary = Standard_True;
break; 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++) { }//for(j = 0; j < 2; j++) {
if(bIsCurrentPointOnBoundary){ if(bIsCurrentPointOnBoundary){
@ -930,7 +621,7 @@ Standard_Boolean IntTools_WLineTool::
Standard_Integer nbboundaries = 0; Standard_Integer nbboundaries = 0;
Standard_Boolean bIsNearBoundary = Standard_False; Standard_Boolean bIsNearBoundary = Standard_False;
Standard_Integer aZoneIndex = 0; //Standard_Integer aZoneIndex = 0;
Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1 Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1 Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
@ -980,32 +671,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; Standard_Boolean bComputeLineEnd = Standard_False;
if(nbboundaries == 2) { if(nbboundaries == 2) {
@ -1144,20 +809,7 @@ Standard_Boolean IntTools_WLineTool::
gp_Pnt2d ap1(nU1, nV1); gp_Pnt2d ap1(nU1, nV1);
gp_Pnt2d ap2; gp_Pnt2d ap2;
if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) {
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() ) {
// re-compute point near boundary if shifted on a period // re-compute point near boundary if shifted on a period
ap2 = AdjustByNeighbour( ap1, anewpoint, aGASurface ); ap2 = AdjustByNeighbour( ap1, anewpoint, aGASurface );

1
src/IntTools/IntTools_WLineTool.hxx

@ -46,7 +46,6 @@ public:
const Standard_Boolean theAvoidLConstructor, const Standard_Boolean theAvoidLConstructor,
const Standard_Real theTol, const Standard_Real theTol,
IntPatch_SequenceOfLine& theNewLines, IntPatch_SequenceOfLine& theNewLines,
Standard_Real& theReachedTol3d,
const Handle(IntTools_Context)& ); const Handle(IntTools_Context)& );
}; };

34
src/LocOpe/LocOpe_WiresOnShape.cxx

@ -1295,6 +1295,10 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge,
Standard_Real aTolV[2]; Standard_Real aTolV[2];
aTolV[0] =BRep_Tool::Tolerance(theVertices[0]); aTolV[0] =BRep_Tool::Tolerance(theVertices[0]);
aTolV[1] =BRep_Tool::Tolerance(theVertices[1]); 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); BRepAdaptor_Curve2d thePCurve(theEdge, theFace);
Bnd_Box2d theBox; Bnd_Box2d theBox;
@ -1304,6 +1308,9 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge,
Standard_Real aFpar, aLpar; Standard_Real aFpar, aLpar;
const Handle(Geom_Curve)& theCurve = BRep_Tool::Curve(theEdge, thePar[0], thePar[1]); const Handle(Geom_Curve)& theCurve = BRep_Tool::Curve(theEdge, thePar[0], thePar[1]);
GeomAdaptor_Curve theGAcurve(theCurve, 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(); Standard_Real aDistMax = Precision::Confusion() * Precision::Confusion();
TopExp_Explorer Explo(theFace, TopAbs_EDGE); TopExp_Explorer Explo(theFace, TopAbs_EDGE);
for (; Explo.More(); Explo.Next()) for (; Explo.More(); Explo.Next())
@ -1330,6 +1337,23 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge,
isOverlapped = Standard_True; isOverlapped = Standard_True;
return; 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++) for (i = 1; i <= aNbExt; i++)
{ {
Standard_Real aDist = anExtrema.SquareDistance(i); Standard_Real aDist = anExtrema.SquareDistance(i);
@ -1342,7 +1366,7 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge,
Standard_Real anIntPar = aPOnC2.Parameter(); Standard_Real anIntPar = aPOnC2.Parameter();
for (j = 0; j < 2; j++) //try to find intersection on an extremity of "theEdge" 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; break;
} }
//intersection found in the middle of the edge //intersection found in the middle of the edge
@ -1351,10 +1375,10 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge,
gp_Pnt aPoint = aCurve->Value(anIntPar); gp_Pnt aPoint = aCurve->Value(anIntPar);
gp_Pnt aPointInt = theCurve->Value(theIntPar); gp_Pnt aPointInt = theCurve->Value(theIntPar);
if (aPointInt.SquareDistance(thePnt[0]) > aTolV[0] * aTolV[0] && if (aPointInt.SquareDistance(thePnt[0]) > aTolVExt[0] &&
aPointInt.SquareDistance(thePnt[1]) > aTolV[1] * aTolV[1] && aPointInt.SquareDistance(thePnt[1]) > aTolVExt[1] &&
aPoint.SquareDistance(thePnt[0]) > aTolV[0] * aTolV[0] && aPoint.SquareDistance(thePnt[0]) > aTolVExt[0] &&
aPoint.SquareDistance(thePnt[1]) > aTolV[1] * aTolV[1]) aPoint.SquareDistance(thePnt[1]) > aTolVExt[1])
{ {
SplitPars.Append(theIntPar); SplitPars.Append(theIntPar);
if( aDist > aDistMax) if( aDist > aDistMax)

2
tests/bugs/modalg_5/bug23706_10

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

9
tests/bugs/modalg_5/bug23706_11

@ -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 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] set info [extrema r1 r2]
if { [llength $info] != 1 } { if {[regexp "ext_1" $info]} {
puts "Error : Extrema is wrong" set dist [lindex [length ext_1] end]
if { $dist > 1.0e-8 } {
puts "Error: Extrema distance is too big"
}
} else { } else {
puts "OK: Extrema is valid" puts "Error: Extrema is not found"
} }

2
tests/bugs/modalg_6/bug27665

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

17
tests/lowalgos/extcc/bug29858_03

@ -18,11 +18,12 @@ mkcurve c2 e2
set info [extrema c1 c2] set info [extrema c1 c2]
# Check result # Check result
regexp {Extrema 1 is point : +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $info full x y z if {[regexp "ext_1" $info]} {
# Point check set dist [lindex [length ext_1] end]
set good_x 0.0 if { $dist > 1.0e-10 } {
set good_y 0.070710562195021642 puts "Error: Extrema distance is too big"
set good_z -0.65305318986891325 }
checkreal "Intersection point x:" ${x} ${good_x} 0.01 0.01 } else {
checkreal "Intersection point y:" ${y} ${good_y} 0.01 0.01 puts "Error: Extrema is not found"
checkreal "Intersection point z:" ${z} ${good_z} 0.01 0.01 }

36
tests/lowalgos/extcc/bug32882 Normal 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