From 493c359ffdad3f204dc092a18581b1cb4d206794 Mon Sep 17 00:00:00 2001 From: nbv Date: Tue, 16 May 2017 15:14:50 +0300 Subject: [PATCH] Patch including (particular) fixes for the issues #23178, #26586, #27766, #26852 and #24023. --- src/BSplCLib/BSplCLib_2.cxx | 5 +- src/Bnd/Bnd.cdl | 3 +- src/Bnd/Bnd_Range.cxx | 37 + src/Bnd/Bnd_Range.hxx | 124 ++ src/Bnd/FILES | 2 + src/Geom2dConvert/Geom2dConvert.cxx | 49 +- src/IntPatch/FILES | 2 + src/IntPatch/IntPatch.cdl | 1 + src/IntPatch/IntPatch_ImpImpIntersection.cdl | 5 +- .../IntPatch_ImpImpIntersection_1.gxx | 13 +- .../IntPatch_ImpImpIntersection_2.gxx | 65 +- .../IntPatch_ImpImpIntersection_4.gxx | 1699 +++++++++-------- src/IntPatch/IntPatch_Intersection.cdl | 18 - src/IntPatch/IntPatch_Intersection.cxx | 508 +---- src/IntPatch/IntPatch_WLineTool.cxx | 1688 ++++++++++++++++ src/IntPatch/IntPatch_WLineTool.hxx | 104 + tests/blend/simple/Q6 | 2 +- tests/bugs/modalg_5/bug25742_2 | 6 +- tests/bugs/modalg_6/bug23178 | 43 + tests/bugs/modalg_6/bug27766 | 30 + 20 files changed, 3024 insertions(+), 1380 deletions(-) create mode 100644 src/Bnd/Bnd_Range.cxx create mode 100644 src/Bnd/Bnd_Range.hxx create mode 100644 src/Bnd/FILES create mode 100644 src/IntPatch/IntPatch_WLineTool.cxx create mode 100644 src/IntPatch/IntPatch_WLineTool.hxx create mode 100644 tests/bugs/modalg_6/bug23178 create mode 100644 tests/bugs/modalg_6/bug27766 diff --git a/src/BSplCLib/BSplCLib_2.cxx b/src/BSplCLib/BSplCLib_2.cxx index 204755cf1d..8ecd0036ce 100644 --- a/src/BSplCLib/BSplCLib_2.cxx +++ b/src/BSplCLib/BSplCLib_2.cxx @@ -42,8 +42,7 @@ struct BSplCLib_DataContainer BSplCLib_DataContainer(Standard_Integer Degree) { (void)Degree; // avoid compiler warning - Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() || - BSplCLib::MaxDegree() > 25, + Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree(), "BSplCLib: bspline degree is greater than maximum supported"); } @@ -324,7 +323,7 @@ BSplCLib::BuildBSpMatrix(const TColStd_Array1OfReal& Parameters, ReturnCode = 0, FirstNonZeroBsplineIndex, BandWidth, - MaxOrder = 21, + MaxOrder = BSplCLib::MaxDegree() + 1, Order ; math_Matrix BSplineBasis(1, MaxOrder, diff --git a/src/Bnd/Bnd.cdl b/src/Bnd/Bnd.cdl index b7653aca11..075ebc861d 100644 --- a/src/Bnd/Bnd.cdl +++ b/src/Bnd/Bnd.cdl @@ -105,5 +105,6 @@ is class Box; class B3f instantiates B3x from Bnd (ShortReal from Standard); -- 3D box with single-precision coordinates - + + imported Range; end Bnd; diff --git a/src/Bnd/Bnd_Range.cxx b/src/Bnd/Bnd_Range.cxx new file mode 100644 index 0000000000..d79b211dea --- /dev/null +++ b/src/Bnd/Bnd_Range.cxx @@ -0,0 +1,37 @@ +// Created on: 2016-06-07 +// Created by: Nikolai BUKHALOV +// Copyright (c) 2016 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + + +#include + +//======================================================================= +//function : Common +//purpose : +//======================================================================= +void Bnd_Range::Common(const Bnd_Range& theOther) +{ + if(theOther.IsVoid()) + { + SetVoid(); + } + + if(IsVoid()) + { + return; + } + + myFirst = Max(myFirst, theOther.myFirst); + myLast = Min(myLast, theOther.myLast); +} diff --git a/src/Bnd/Bnd_Range.hxx b/src/Bnd/Bnd_Range.hxx new file mode 100644 index 0000000000..439008f1f8 --- /dev/null +++ b/src/Bnd/Bnd_Range.hxx @@ -0,0 +1,124 @@ +// Created on: 2016-06-07 +// Created by: Nikolai BUKHALOV +// Copyright (c) 2016 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _Bnd_Range_HeaderFile +#define _Bnd_Range_HeaderFile + +#include +#include + +//! This class describes a range in 1D space restricted +//! by two real values. +//! A range can be void indicating there is no point included in the range. + +class Bnd_Range +{ +public: + + //! Default constructor. Creates VOID range. + Bnd_Range() : myFirst(0.0), myLast(-1.0) + { + }; + + //! Constructor. Never creates VOID range. + Bnd_Range(const Standard_Real theMin, const Standard_Real theMax) : + myFirst(theMin), myLast(theMax) + { + if(myLast < myFirst) + Standard_ConstructionError::Raise("Last < First"); + }; + + //! Replaces with common-part of and theOther + Standard_EXPORT void Common(const Bnd_Range& theOther); + + //! Extends to include theParameter + void Add(const Standard_Real theParameter) + { + if(IsVoid()) + { + myFirst = myLast = theParameter; + return; + } + + myFirst = Min(myFirst, theParameter); + myLast = Max(myLast, theParameter); + } + + //! Obtain MIN boundary of . + //! If is VOID the method returns false. + Standard_Boolean GetMin(Standard_Real& thePar) const + { + if(IsVoid()) + { + return Standard_False; + } + + thePar = myFirst; + return Standard_True; + } + + //! Obtain MAX boundary of . + //! If is VOID the method returns false. + Standard_Boolean GetMAX(Standard_Real& thePar) const + { + if(IsVoid()) + { + return Standard_False; + } + + thePar = myLast; + return Standard_True; + } + + //! Returns range value (MAX-MIN). Returns negative value for VOID range. + Standard_Real Delta() const + { + return (myLast - myFirst); + } + + //! Is initialized. + Standard_Boolean IsVoid() const + { + return (myLast < myFirst); + } + + //! Initializes by default parameters. Makes VOID. + void SetVoid() + { + myLast = -1.0; + myFirst = 0.0; + } + + //! Extends this to the given value (in both side) + void Enlarge(const Standard_Real theDelta) + { + if (IsVoid()) + { + return; + } + + myFirst -= theDelta; + myLast += theDelta; + } + +private: + //! Start of range + Standard_Real myFirst; + + //! End of range + Standard_Real myLast; +}; + +#endif \ No newline at end of file diff --git a/src/Bnd/FILES b/src/Bnd/FILES new file mode 100644 index 0000000000..6284d7369b --- /dev/null +++ b/src/Bnd/FILES @@ -0,0 +1,2 @@ +Bnd_Range.cxx +Bnd_Range.hxx diff --git a/src/Geom2dConvert/Geom2dConvert.cxx b/src/Geom2dConvert/Geom2dConvert.cxx index 911c171937..285a70bac6 100644 --- a/src/Geom2dConvert/Geom2dConvert.cxx +++ b/src/Geom2dConvert/Geom2dConvert.cxx @@ -1053,6 +1053,7 @@ void Geom2dConvert::ConcatG1(TColGeom2d_Array1OfBSplineCurve& ArrayOf for (j=1;j<=nb_curve-1;j++){ //boucle secondaire a l'interieur de chaque groupe Curve1=ArrayOfCurves(j); if ( (j==(nb_curve-1)) &&(Need2DegRepara(ArrayOfCurves))){ + const Standard_Integer aNewCurveDegree = Min(2 * Curve1->Degree(), Geom2d_BSplineCurve::MaxDegree()); Curve2->D1(Curve2->LastParameter(),Pint,Vec1); Curve1->D1(Curve1->FirstParameter(),Pint,Vec2); lambda=Vec2.Magnitude()/Vec1.Magnitude(); @@ -1102,7 +1103,7 @@ void Geom2dConvert::ConcatG1(TColGeom2d_Array1OfBSplineCurve& ArrayOf Curve1FlatKnots, Curve1Poles, FlatKnots, - 2*Curve1->Degree(), + aNewCurveDegree, NewPoles, Status ); @@ -1112,7 +1113,7 @@ void Geom2dConvert::ConcatG1(TColGeom2d_Array1OfBSplineCurve& ArrayOf Curve1FlatKnots, Curve1Weights, FlatKnots, - 2*Curve1->Degree(), + aNewCurveDegree, NewWeights, Status ); @@ -1138,7 +1139,7 @@ void Geom2dConvert::ConcatG1(TColGeom2d_Array1OfBSplineCurve& ArrayOf for (ii=1;ii<=NewPoles.Length();ii++) for (jj=1;jj<=2;jj++) NewPoles(ii).SetCoord(jj,NewPoles(ii).Coord(jj)/NewWeights(ii)); - Curve1= new Geom2d_BSplineCurve(NewPoles,NewWeights,KnotC1,KnotC1Mults,2*Curve1->Degree()); + Curve1 = new Geom2d_BSplineCurve(NewPoles, NewWeights, KnotC1, KnotC1Mults, aNewCurveDegree); } Geom2dConvert_CompCurveToBSplineCurve C(Handle(Geom2d_BSplineCurve)::DownCast(Curve2)); fusion=C.Add(Curve1, @@ -1304,6 +1305,8 @@ void Geom2dConvert::ConcatC1(TColGeom2d_Array1OfBSplineCurve& ArrayOf else Curve1=ArrayOfCurves(j); + const Standard_Integer aNewCurveDegree = Min(2 * Curve1->Degree(), Geom2d_BSplineCurve::MaxDegree()); + if (j==0) //initialisation en debut de groupe Curve2=Curve1; else{ @@ -1343,7 +1346,7 @@ void Geom2dConvert::ConcatC1(TColGeom2d_Array1OfBSplineCurve& ArrayOf TColStd_Array1OfReal FlatKnots(1,Curve1FlatKnots.Length()+(Curve1->Degree()*Curve1->NbKnots())); BSplCLib::KnotSequence(KnotC1,KnotC1Mults,FlatKnots); - TColgp_Array1OfPnt2d NewPoles(1,FlatKnots.Length()-(2*Curve1->Degree()+1)); + TColgp_Array1OfPnt2d NewPoles(1, FlatKnots.Length() - (aNewCurveDegree + 1)); Standard_Integer Status; TColStd_Array1OfReal Curve1Weights(1,Curve1->NbPoles()); Curve1->Weights(Curve1Weights); @@ -1358,18 +1361,18 @@ void Geom2dConvert::ConcatC1(TColGeom2d_Array1OfBSplineCurve& ArrayOf Curve1FlatKnots, Curve1Poles, FlatKnots, - 2*Curve1->Degree(), + aNewCurveDegree, NewPoles, Status ); - TColStd_Array1OfReal NewWeights(1,FlatKnots.Length()-(2*Curve1->Degree()+1)); + TColStd_Array1OfReal NewWeights(1, FlatKnots.Length() - (aNewCurveDegree + 1)); // BSplCLib::FunctionReparameterise(reparameterise_evaluator, BSplCLib::FunctionReparameterise(ev, Curve1->Degree(), Curve1FlatKnots, Curve1Weights, FlatKnots, - 2*Curve1->Degree(), + aNewCurveDegree, NewWeights, Status ); @@ -1377,7 +1380,7 @@ void Geom2dConvert::ConcatC1(TColGeom2d_Array1OfBSplineCurve& ArrayOf for (jj=1;jj<=2;jj++) NewPoles(ii).SetCoord(jj,NewPoles(ii).Coord(jj)/NewWeights(ii)); } - Curve1= new Geom2d_BSplineCurve(NewPoles,NewWeights,KnotC1,KnotC1Mults,2*Curve1->Degree()); + Curve1 = new Geom2d_BSplineCurve(NewPoles, NewWeights, KnotC1, KnotC1Mults, aNewCurveDegree); } Geom2dConvert_CompCurveToBSplineCurve C(Handle(Geom2d_BSplineCurve)::DownCast(Curve2)); fusion=C.Add(Curve1, @@ -1448,7 +1451,7 @@ void Geom2dConvert::C0BSplineToC1BSplineCurve(Handle(Geom2d_BSplineCurve)& BS, Standard_Integer i,j,nbcurveC1=1; Standard_Real U1,U2; Standard_Boolean closed_flag = Standard_False ; - gp_Pnt2d point; + gp_Pnt2d point1, point2; gp_Vec2d V1,V2; Standard_Boolean fusion; @@ -1481,15 +1484,20 @@ void Geom2dConvert::C0BSplineToC1BSplineCurve(Handle(Geom2d_BSplineCurve)& BS, BSbis->Segment(U1,U2); ArrayOfCurves(i)=BSbis; } + + const Standard_Real anAngularToler = 1.0e-7; Handle(TColStd_HArray1OfInteger) ArrayOfIndices; Handle(TColGeom2d_HArray1OfBSplineCurve) ArrayOfConcatenated; - BS->D1(BS->FirstParameter(),point,V1); //a verifier - BS->D1(BS->LastParameter(),point,V2); + BS->D1(BS->FirstParameter(),point1,V1); //a verifier + BS->D1(BS->LastParameter(),point2,V2); - if ((BS->IsClosed())&&(V1.IsParallel(V2,Precision::Confusion()))) - closed_flag = Standard_True ; - + if ((point1.SquareDistance(point2) < gp::Resolution()) && + (V1.IsParallel(V2, anAngularToler))) + { + closed_flag = Standard_True; + } + Geom2dConvert::ConcatC1(ArrayOfCurves, ArrayOfToler, ArrayOfIndices, @@ -1539,7 +1547,7 @@ void Geom2dConvert::C0BSplineToArrayOfC1BSplineCurve(const Handle(Geom2d_BSpline Standard_Integer i,j,nbcurveC1=1; Standard_Real U1,U2; Standard_Boolean closed_flag = Standard_False ; - gp_Pnt2d point; + gp_Pnt2d point1, point2; gp_Vec2d V1,V2; // Standard_Boolean fusion; @@ -1573,11 +1581,14 @@ void Geom2dConvert::C0BSplineToArrayOfC1BSplineCurve(const Handle(Geom2d_BSpline Handle(TColStd_HArray1OfInteger) ArrayOfIndices; - BS->D1(BS->FirstParameter(),point,V1); - BS->D1(BS->LastParameter(),point,V2); + BS->D1(BS->FirstParameter(),point1,V1); + BS->D1(BS->LastParameter(),point2,V2); - if ((BS->IsClosed())&&(V1.IsParallel(V2,AngularTolerance))) - closed_flag = Standard_True ; + if (((point1.SquareDistance(point2) < gp::Resolution())) && + (V1.IsParallel(V2, AngularTolerance))) + { + closed_flag = Standard_True; + } Geom2dConvert::ConcatC1(ArrayOfCurves, ArrayOfToler, diff --git a/src/IntPatch/FILES b/src/IntPatch/FILES index 446e204fee..662f8f3191 100755 --- a/src/IntPatch/FILES +++ b/src/IntPatch/FILES @@ -5,3 +5,5 @@ IntPatch_ImpImpIntersection_3.gxx IntPatch_ImpImpIntersection_4.gxx IntPatch_ImpImpIntersection_5.gxx IntPatch_ImpImpIntersection_6.gxx +IntPatch_WLineTool.cxx +IntPatch_WLineTool.hxx \ No newline at end of file diff --git a/src/IntPatch/IntPatch.cdl b/src/IntPatch/IntPatch.cdl index f4a06df072..c3c557d5ae 100644 --- a/src/IntPatch/IntPatch.cdl +++ b/src/IntPatch/IntPatch.cdl @@ -152,4 +152,5 @@ is HCurve2dTool from IntPatch, CSFunction from IntPatch); + imported WLineTool; end IntPatch; diff --git a/src/IntPatch/IntPatch_ImpImpIntersection.cdl b/src/IntPatch/IntPatch_ImpImpIntersection.cdl index 15421c8e86..e6f9ea36d4 100644 --- a/src/IntPatch/IntPatch_ImpImpIntersection.cdl +++ b/src/IntPatch/IntPatch_ImpImpIntersection.cdl @@ -63,7 +63,6 @@ is S1: HSurface from Adaptor3d; D1: TopolTool from Adaptor3d; S2: HSurface from Adaptor3d; D2: TopolTool from Adaptor3d; TolArc,TolTang: Real from Standard; - isTheTrimmed: Boolean from Standard = Standard_False; theIsReqToKeepRLine: Boolean from Standard = Standard_False) raises ConstructionError from Standard @@ -103,8 +102,8 @@ is TangentFaces(me) ---Purpose: Returns True if the two patches are considered as - -- entierly tangent, i-e every restriction arc of one - -- patch is inside the geometric base of the otehr patch. + -- entirely tangent, i.e. every restriction arc of one + -- patch is inside the geometric base of the other patch. returns Boolean from Standard ---C++: inline diff --git a/src/IntPatch/IntPatch_ImpImpIntersection_1.gxx b/src/IntPatch/IntPatch_ImpImpIntersection_1.gxx index 18146a1684..e2f049e236 100644 --- a/src/IntPatch/IntPatch_ImpImpIntersection_1.gxx +++ b/src/IntPatch/IntPatch_ImpImpIntersection_1.gxx @@ -70,16 +70,7 @@ static void ProcessBounds(const Handle(IntPatch_ALine)&, const Standard_Real); -static Standard_Boolean IntCyCy(const IntSurf_Quadric&, - const IntSurf_Quadric&, - const Standard_Real, - Standard_Boolean&, - Standard_Boolean&, - Standard_Boolean&, - IntPatch_SequenceOfLine&, - IntPatch_SequenceOfPoint&); - -static Standard_Boolean IntCyCyTrim(const IntSurf_Quadric& theQuad1, +static Standard_Boolean IntCyCy(const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, const Standard_Real theTol3D, const Standard_Real theTol2D, @@ -87,6 +78,8 @@ static Standard_Boolean IntCyCyTrim(const IntSurf_Quadric& theQuad1, const Bnd_Box2d& theUVSurf2, const Standard_Boolean isTheReverse, Standard_Boolean& isTheEmpty, + Standard_Boolean& isTheSameSurface, + Standard_Boolean& isTheMultiplePoint, IntPatch_SequenceOfLine& theSlin, IntPatch_SequenceOfPoint& theSPnt); diff --git a/src/IntPatch/IntPatch_ImpImpIntersection_2.gxx b/src/IntPatch/IntPatch_ImpImpIntersection_2.gxx index f94f8096d2..145b17d3b2 100644 --- a/src/IntPatch/IntPatch_ImpImpIntersection_2.gxx +++ b/src/IntPatch/IntPatch_ImpImpIntersection_2.gxx @@ -14,6 +14,8 @@ // Alternatively, this file may be used under the terms of Open CASCADE // commercial license or contractual agreement. +#include + static Standard_Integer SetQuad(const Handle(Adaptor3d_HSurface)& theS, GeomAbs_SurfaceType& theTS, @@ -40,7 +42,7 @@ IntPatch_ImpImpIntersection::IntPatch_ImpImpIntersection const Standard_Real TolTang, const Standard_Boolean theIsReqToKeepRLine) { - Perform(S1,D1,S2,D2,TolArc,TolTang, Standard_False, theIsReqToKeepRLine); + Perform(S1,D1,S2,D2,TolArc,TolTang, theIsReqToKeepRLine); } //======================================================================= //function : Perform @@ -52,13 +54,14 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_TopolTool)& D2, const Standard_Real TolArc, const Standard_Real TolTang, - const Standard_Boolean isTheTrimmed, - const Standard_Boolean theIsReqToKeepRLine) { + const Standard_Boolean theIsReqToKeepRLine) +{ done = Standard_False; - Standard_Boolean isTrimmed = isTheTrimmed; spnt.Clear(); slin.Clear(); + Standard_Boolean isPostProcessingRequired = Standard_True; + empt = Standard_True; tgte = Standard_False; oppo = Standard_False; @@ -149,14 +152,6 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)& S1, case 22: { // Cylinder/Cylinder Standard_Boolean isDONE = Standard_False; - - if(!isTrimmed) - { - isDONE = IntCyCy(quad1, quad2, TolTang, empt, - SameSurf, multpoint, slin, spnt); - } - else - { Bnd_Box2d aBox1, aBox2; const Standard_Real aU1f = S1->FirstUParameter(); @@ -177,36 +172,52 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)& S1, aBox2.Add(gp_Pnt2d(aU2f, S2->FirstVParameter())); aBox2.Add(gp_Pnt2d(aU2l, S2->LastVParameter())); - const Standard_Real a2DTol = Min( S1->UResolution(TolTang), - S2->UResolution(TolTang)); - + // Resolution is too big if the cylinder radius is + // too small. Therefore, we shall bounde its value above. + // Here, we use simple constant. + const Standard_Real a2DTol = Min(1.0e-4, Min( S1->UResolution(TolTang), + S2->UResolution(TolTang))); + + //The bigger range the bigger nunber of points in Walking-line (WLine) + //we will be able to add and, consequently, we will obtain more + //precise intersection line. + //Every point of WLine is determined as function from U1-parameter, + //where U1 is U-parameter on 1st quadric. + //Therefore, we should use quadric with bigger range as 1st parameter + //in IntCyCy() function. + //On the other hand, there is no point in reversing in case of + //analytical intersection (when result is line, ellipse, point...). + //This result is independent of the arguments order. Standard_Boolean isReversed = ((aU2l - aU2f) < (aU1l - aU1f)); if(isReversed) { - isDONE = IntCyCyTrim(quad2, quad1, TolTang, a2DTol, aBox2, aBox1, - isReversed, empt, slin, spnt); + isDONE = IntCyCy(quad2, quad1, TolTang, a2DTol, aBox2, aBox1, + Standard_True, empt, SameSurf, multpoint, slin, spnt); } else { - isDONE = IntCyCyTrim(quad1, quad2, TolTang, a2DTol, aBox1, aBox2, - isReversed, empt, slin, spnt); + isDONE = IntCyCy(quad1, quad2, TolTang, a2DTol, aBox1, aBox2, + Standard_False, empt, SameSurf, multpoint, slin, spnt); } if(!isDONE) { - isDONE = IntCyCy(quad1, quad2, TolTang, empt, - SameSurf, multpoint, slin, spnt); - isTrimmed = Standard_False; + return; } - } - if (!isDONE) + bEmpty = empt; + if(!slin.IsEmpty()) { - return; + const Handle(IntPatch_WLine)& aWLine = + Handle(IntPatch_WLine)::DownCast(slin.Value(1)); + + if(!aWLine.IsNull()) + {//No geometric solution + isPostProcessingRequired = Standard_False; + } } - bEmpty = empt; break; } // @@ -299,7 +310,7 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)& S1, } // - if(!isTrimmed) + if(isPostProcessingRequired) { if (!SameSurf) { AFunc.SetQuadric(quad2); diff --git a/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx b/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx index fda766fd82..34fe6c5841 100644 --- a/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx +++ b/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx @@ -15,6 +15,7 @@ // commercial license or contractual agreement. #include +#include #include #include #include @@ -25,7 +26,10 @@ //If Abs(a) <= aNulValue then it is considered that a = 0. static const Standard_Real aNulValue = 1.0e-11; -struct stCoeffsValue; +static void ShortCosForm( const Standard_Real theCosFactor, + const Standard_Real theSinFactor, + Standard_Real& theCoeff, + Standard_Real& theAngle); // static Standard_Boolean ExploreCurve(const gp_Cylinder& aCy, @@ -33,10 +37,6 @@ static IntAna_Curve& aC, const Standard_Real aTol, IntAna_ListOfCurve& aLC); -static - Standard_Boolean IsToReverse(const gp_Cylinder& Cy1, - const gp_Cylinder& Cy2, - const Standard_Real Tol); static Standard_Boolean InscribePoint(const Standard_Real theUfTarget, const Standard_Real theUlTarget, @@ -45,16 +45,399 @@ static Standard_Boolean InscribePoint(const Standard_Real theUfTarget, const Standard_Real thePeriod, const Standard_Boolean theFlForce); + +class ComputationMethods +{ +public: + //Stores equations coefficients + struct stCoeffsValue + { + stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&); + + math_Vector mVecA1; + math_Vector mVecA2; + math_Vector mVecB1; + math_Vector mVecB2; + math_Vector mVecC1; + math_Vector mVecC2; + math_Vector mVecD; + + Standard_Real mK21; //sinU2 + Standard_Real mK11; //sinU1 + Standard_Real mL21; //cosU2 + Standard_Real mL11; //cosU1 + Standard_Real mM1; //Free member + + Standard_Real mK22; //sinU2 + Standard_Real mK12; //sinU1 + Standard_Real mL22; //cosU2 + Standard_Real mL12; //cosU1 + Standard_Real mM2; //Free member + + Standard_Real mK1; + Standard_Real mL1; + Standard_Real mK2; + Standard_Real mL2; + + Standard_Real mFIV1; + Standard_Real mPSIV1; + Standard_Real mFIV2; + Standard_Real mPSIV2; + + Standard_Real mB; + Standard_Real mC; + Standard_Real mFI1; + Standard_Real mFI2; + }; + + + //! Determines, if U2(U1) function is increasing. + static Standard_Boolean CylCylMonotonicity(const Standard_Real theU1par, + const Standard_Integer theWLIndex, + const stCoeffsValue& theCoeffs, + const Standard_Real thePeriod, + Standard_Boolean& theIsIncreasing); + + //! Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0, + //! esimates the tolerance of U2-computing (estimation result is + //! assigned to *theDelta value). + static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par, + const Standard_Integer theWLIndex, + const stCoeffsValue& theCoeffs, + Standard_Real& theU2, + Standard_Real* const theDelta = 0); + + static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1, + const Standard_Real theU2, + const stCoeffsValue& theCoeffs, + Standard_Real& theV1, + Standard_Real& theV2); + + static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par, + const Standard_Integer theWLIndex, + const stCoeffsValue& theCoeffs, + Standard_Real& theU2, + Standard_Real& theV1, + Standard_Real& theV2); + +}; + +ComputationMethods::stCoeffsValue::stCoeffsValue(const gp_Cylinder& theCyl1, + const gp_Cylinder& theCyl2): + mVecA1(-theCyl1.Radius()*theCyl1.XAxis().Direction().XYZ()), + mVecA2(theCyl2.Radius()*theCyl2.XAxis().Direction().XYZ()), + mVecB1(-theCyl1.Radius()*theCyl1.YAxis().Direction().XYZ()), + mVecB2(theCyl2.Radius()*theCyl2.YAxis().Direction().XYZ()), + mVecC1(theCyl1.Axis().Direction().XYZ()), + mVecC2(theCyl2.Axis().Direction().XYZ().Reversed()), + mVecD(theCyl2.Location().XYZ() - theCyl1.Location().XYZ()) +{ + enum CoupleOfEquation + { + COENONE = 0, + COE12 = 1, + COE23 = 2, + COE13 = 3 + }aFoundCouple = COENONE; + + + Standard_Real aDetV1V2 = 0.0; + + const Standard_Real aDelta1 = mVecC1(1)*mVecC2(2)-mVecC1(2)*mVecC2(1); //1-2 + const Standard_Real aDelta2 = mVecC1(2)*mVecC2(3)-mVecC1(3)*mVecC2(2); //2-3 + const Standard_Real aDelta3 = mVecC1(1)*mVecC2(3)-mVecC1(3)*mVecC2(1); //1-3 + const Standard_Real anAbsD1 = Abs(aDelta1); //1-2 + const Standard_Real anAbsD2 = Abs(aDelta2); //2-3 + const Standard_Real anAbsD3 = Abs(aDelta3); //1-3 + + if(anAbsD1 >= anAbsD2) + { + if(anAbsD3 > anAbsD1) + { + aFoundCouple = COE13; + aDetV1V2 = aDelta3; + } + else + { + aFoundCouple = COE12; + aDetV1V2 = aDelta1; + } + } + else + { + if(anAbsD3 > anAbsD2) + { + aFoundCouple = COE13; + aDetV1V2 = aDelta3; + } + else + { + aFoundCouple = COE23; + aDetV1V2 = aDelta2; + } + } + + // In point of fact, every determinant (aDelta1, aDelta2 and aDelta3) is + // cross-product between directions (i.e. sine of angle). + // If sine is too small then sine is (approx.) equal to angle itself. + // Therefore, in this case we should compare sine with angular tolerance. + // This constant is used for check if axes are parallel (see constructor + // AxeOperator::AxeOperator(...) in IntAna_QuadQuadGeo.cxx file). + if(Abs(aDetV1V2) < Precision::Angular()) + { + Standard_Failure::Raise("Error. Exception in divide by zerro (IntCyCyTrim)!!!!"); + } + + switch(aFoundCouple) + { + case COE12: + break; + case COE23: + { + math_Vector aVTemp(mVecA1); + mVecA1(1) = aVTemp(2); + mVecA1(2) = aVTemp(3); + mVecA1(3) = aVTemp(1); + + aVTemp = mVecA2; + mVecA2(1) = aVTemp(2); + mVecA2(2) = aVTemp(3); + mVecA2(3) = aVTemp(1); + + aVTemp = mVecB1; + mVecB1(1) = aVTemp(2); + mVecB1(2) = aVTemp(3); + mVecB1(3) = aVTemp(1); + + aVTemp = mVecB2; + mVecB2(1) = aVTemp(2); + mVecB2(2) = aVTemp(3); + mVecB2(3) = aVTemp(1); + + aVTemp = mVecC1; + mVecC1(1) = aVTemp(2); + mVecC1(2) = aVTemp(3); + mVecC1(3) = aVTemp(1); + + aVTemp = mVecC2; + mVecC2(1) = aVTemp(2); + mVecC2(2) = aVTemp(3); + mVecC2(3) = aVTemp(1); + + aVTemp = mVecD; + mVecD(1) = aVTemp(2); + mVecD(2) = aVTemp(3); + mVecD(3) = aVTemp(1); + + break; + } + case COE13: + { + math_Vector aVTemp = mVecA1; + mVecA1(2) = aVTemp(3); + mVecA1(3) = aVTemp(2); + + aVTemp = mVecA2; + mVecA2(2) = aVTemp(3); + mVecA2(3) = aVTemp(2); + + aVTemp = mVecB1; + mVecB1(2) = aVTemp(3); + mVecB1(3) = aVTemp(2); + + aVTemp = mVecB2; + mVecB2(2) = aVTemp(3); + mVecB2(3) = aVTemp(2); + + aVTemp = mVecC1; + mVecC1(2) = aVTemp(3); + mVecC1(3) = aVTemp(2); + + aVTemp = mVecC2; + mVecC2(2) = aVTemp(3); + mVecC2(3) = aVTemp(2); + + aVTemp = mVecD; + mVecD(2) = aVTemp(3); + mVecD(3) = aVTemp(2); + + break; + } + default: + break; + } + + //------- For V1 (begin) + //sinU2 + mK21 = (mVecC2(2)*mVecB2(1)-mVecC2(1)*mVecB2(2))/aDetV1V2; + //sinU1 + mK11 = (mVecC2(2)*mVecB1(1)-mVecC2(1)*mVecB1(2))/aDetV1V2; + //cosU2 + mL21 = (mVecC2(2)*mVecA2(1)-mVecC2(1)*mVecA2(2))/aDetV1V2; + //cosU1 + mL11 = (mVecC2(2)*mVecA1(1)-mVecC2(1)*mVecA1(2))/aDetV1V2; + //Free member + mM1 = (mVecC2(2)*mVecD(1)-mVecC2(1)*mVecD(2))/aDetV1V2; + //------- For V1 (end) + + //------- For V2 (begin) + //sinU2 + mK22 = (mVecC1(1)*mVecB2(2)-mVecC1(2)*mVecB2(1))/aDetV1V2; + //sinU1 + mK12 = (mVecC1(1)*mVecB1(2)-mVecC1(2)*mVecB1(1))/aDetV1V2; + //cosU2 + mL22 = (mVecC1(1)*mVecA2(2)-mVecC1(2)*mVecA2(1))/aDetV1V2; + //cosU1 + mL12 = (mVecC1(1)*mVecA1(2)-mVecC1(2)*mVecA1(1))/aDetV1V2; + //Free member + mM2 = (mVecC1(1)*mVecD(2)-mVecC1(2)*mVecD(1))/aDetV1V2; + //------- For V1 (end) + + ShortCosForm(mL11, mK11, mK1, mFIV1); + ShortCosForm(mL21, mK21, mL1, mPSIV1); + ShortCosForm(mL12, mK12, mK2, mFIV2); + ShortCosForm(mL22, mK22, mL2, mPSIV2); + + const Standard_Real aA1=mVecC1(3)*mK21+mVecC2(3)*mK22-mVecB2(3), //sinU2 + aA2=mVecC1(3)*mL21+mVecC2(3)*mL22-mVecA2(3), //cosU2 + aB1=mVecB1(3)-mVecC1(3)*mK11-mVecC2(3)*mK12, //sinU1 + aB2=mVecA1(3)-mVecC1(3)*mL11-mVecC2(3)*mL12; //cosU1 + + mC =mVecD(3) - mVecC1(3)*mM1 -mVecC2(3)*mM2; //Free + + Standard_Real aA = 0.0; + + ShortCosForm(aB2,aB1,mB,mFI1); + ShortCosForm(aA2,aA1,aA,mFI2); + + mB /= aA; + mC /= aA; +} + +class WorkWithBoundaries +{ +public: + enum SearchBoundType + { + SearchNONE = 0, + SearchV1 = 1, + SearchV2 = 2 + }; + + struct StPInfo + { + StPInfo() + { + mySurfID = 0; + myU1 = RealLast(); + myV1 = RealLast(); + myU2 = RealLast(); + myV2 = RealLast(); + } + + //Equal to 0 for 1st surface non-zerro for 2nd one. + Standard_Integer mySurfID; + + Standard_Real myU1; + Standard_Real myV1; + Standard_Real myU2; + Standard_Real myV2; + + bool operator>(const StPInfo& theOther) const + { + return myU1 > theOther.myU1; + } + + bool operator<(const StPInfo& theOther) const + { + return myU1 < theOther.myU1; + } + + bool operator==(const StPInfo& theOther) const + { + return myU1 == theOther.myU1; + } + }; + + WorkWithBoundaries(const IntSurf_Quadric& theQuad1, + const IntSurf_Quadric& theQuad2, + const ComputationMethods::stCoeffsValue& theCoeffs, + const Bnd_Box2d& theUVSurf1, + const Bnd_Box2d& theUVSurf2, + const Standard_Integer theNbWLines, + const Standard_Real thePeriod, + const Standard_Real theTol3D, + const Standard_Real theTol2D, + const Standard_Boolean isTheReverse) : + myQuad1(theQuad1), myQuad2(theQuad2), myCoeffs(theCoeffs), + myUVSurf1(theUVSurf1), myUVSurf2(theUVSurf2), myNbWLines(theNbWLines), + myPeriod(thePeriod), myTol3D(theTol3D), myTol2D(theTol2D), + myIsReverse(isTheReverse) + { + }; + + void AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL, + const Standard_Real theU1, + const Standard_Real theU2, + const Standard_Real theV1, + const Standard_Real theV1Prev, + const Standard_Real theV2, + const Standard_Real theV2Prev, + const Standard_Integer theWLIndex, + const Standard_Boolean theFlForce, + Standard_Boolean& isTheFound1, + Standard_Boolean& isTheFound2) const; + + Standard_Boolean BoundariesComputing(Standard_Real theU1f[], + Standard_Real theU1l[]) const; + + void BoundaryEstimation(const gp_Cylinder& theCy1, + const gp_Cylinder& theCy2, + Bnd_Range& theOutBoxS1, + Bnd_Range& theOutBoxS2) const; + +protected: + + Standard_Boolean + SearchOnVBounds(const SearchBoundType theSBType, + const Standard_Real theVzad, + const Standard_Real theVInit, + const Standard_Real theInitU2, + const Standard_Real theInitMainVar, + Standard_Real& theMainVariableValue) const; + + const WorkWithBoundaries& operator=(const WorkWithBoundaries&); + +private: + friend class ComputationMethods; + + const IntSurf_Quadric& myQuad1; + const IntSurf_Quadric& myQuad2; + const ComputationMethods::stCoeffsValue& myCoeffs; + const Bnd_Box2d& myUVSurf1; + const Bnd_Box2d& myUVSurf2; + const Standard_Integer myNbWLines; + const Standard_Real myPeriod; + const Standard_Real myTol3D; + const Standard_Real myTol2D; + const Standard_Boolean myIsReverse; +}; + +static + Standard_Boolean ExploreCurve(const gp_Cylinder& aCy, + const gp_Cone& aCo, + IntAna_Curve& aC, + const Standard_Real aTol, + IntAna_ListOfCurve& aLC); + static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, const Handle(IntSurf_LineOn2S)& theLine, - const stCoeffsValue& theCoeffs, + const ComputationMethods::stCoeffsValue& theCoeffs, const Standard_Integer theWLIndex, const Standard_Integer theMinNbPoints, const Standard_Integer theStartPointOnLine, const Standard_Integer theEndPointOnLine, - const Standard_Real theU2f, - const Standard_Real theU2l, const Standard_Real theTol2D, const Standard_Real thePeriodOfSurf2, const Standard_Boolean isTheReverse); @@ -74,6 +457,29 @@ static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX) } } +//======================================================================= +//function : ExtremaLineLine +//purpose : Computes extrema between the given lines. Returns parameters +// on correspond curve. +//======================================================================= +static inline void ExtremaLineLine(const gp_Ax1& theC1, + const gp_Ax1& theC2, + const Standard_Real theCosA, + const Standard_Real theSqSinA, + Standard_Real& thePar1, + Standard_Real& thePar2) +{ + const gp_Dir &aD1 = theC1.Direction(), + &aD2 = theC2.Direction(); + + const gp_XYZ aL1L2 = theC2.Location().XYZ() - theC1.Location().XYZ(); + const Standard_Real aD1L = aD1.XYZ().Dot(aL1L2), + aD2L = aD2.XYZ().Dot(aL1L2); + + thePar1 = (aD1L - theCosA * aD2L) / theSqSinA; + thePar2 = (theCosA * aD1L - aD2L) / theSqSinA; +} + //======================================================================= //function : VBoundaryPrecise //purpose : By default, we shall consider, that V1 and V2 will increase @@ -165,7 +571,7 @@ static Standard_Boolean StepComputing(const math_Matrix& theMatr, Standard_Real& theV1Found, Standard_Real& theV2Found*/) { -#ifdef OCCT_DEBUG +#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG bool flShow = false; if(flShow) @@ -426,19 +832,21 @@ void ProcessBounds(const Handle(IntPatch_ALine)& alig, //-- ligne coura alig->SetLastPoint(alig->NbVertex()); } } + //======================================================================= -//function : IntCyCy +//function : CyCyAnalyticalIntersect //purpose : //======================================================================= -Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, +Standard_Boolean CyCyAnalyticalIntersect( const IntSurf_Quadric& Quad1, const IntSurf_Quadric& Quad2, + const IntAna_QuadQuadGeo& theInter, const Standard_Real Tol, + const Standard_Boolean isTheReverse, Standard_Boolean& Empty, Standard_Boolean& Same, Standard_Boolean& Multpoint, IntPatch_SequenceOfLine& slin, IntPatch_SequenceOfPoint& spnt) - { IntPatch_Point ptsol; @@ -453,19 +861,12 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, gp_Cylinder Cy1(Quad1.Cylinder()); gp_Cylinder Cy2(Quad2.Cylinder()); - IntAna_QuadQuadGeo inter(Cy1,Cy2,Tol); + typint = theInter.TypeInter(); + Standard_Integer NbSol = theInter.NbSolutions(); + Empty = Standard_False; + Same = Standard_False; - if (!inter.IsDone()) - { - return Standard_False; - } - - typint = inter.TypeInter(); - Standard_Integer NbSol = inter.NbSolutions(); - Empty = Standard_False; - Same = Standard_False; - - switch (typint) + switch (typint) { case IntAna_Empty: { @@ -481,11 +882,22 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, case IntAna_Point: { - gp_Pnt psol(inter.Point(1)); + gp_Pnt psol(theInter.Point(1)); + ptsol.SetValue(psol,Tol,Standard_True); + Standard_Real U1,V1,U2,V2; + + if(isTheReverse) + { + Quad2.Parameters(psol, U1, V1); + Quad1.Parameters(psol, U2, V2); + } + else + { Quad1.Parameters(psol,U1,V1); Quad2.Parameters(psol,U2,V2); - ptsol.SetValue(psol,Tol,Standard_True); + } + ptsol.SetParameters(U1,V1,U2,V2); spnt.Append(ptsol); } @@ -496,10 +908,14 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, gp_Pnt ptref; if (NbSol == 1) { // Cylinders are tangent to each other by line - linsol = inter.Line(1); + linsol = theInter.Line(1); ptref = linsol.Location(); + + //Radius-vectors gp_Dir crb1(gp_Vec(ptref,Cy1.Location())); gp_Dir crb2(gp_Vec(ptref,Cy2.Location())); + + //outer normal lines gp_Vec norm1(Quad1.Normale(ptref)); gp_Vec norm2(Quad2.Normale(ptref)); IntSurf_Situation situcyl1; @@ -507,6 +923,11 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, if (crb1.Dot(crb2) < 0.) { // centre de courbures "opposes" + //ATTENTION!!! + // Normal and Radius-vector of the 1st(!) cylinder + // is used for judging what the situation of the 2nd(!) + // cylinder is. + if (norm1.Dot(crb1) > 0.) { situcyl2 = IntSurf_Inside; @@ -576,9 +997,11 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, { for (i=1; i <= NbSol; i++) { - linsol = inter.Line(i); + linsol = theInter.Line(i); ptref = linsol.Location(); gp_Vec lsd = linsol.Direction(); + + //Theoretically, qwe = +/- 1.0. Standard_Real qwe = lsd.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)); if (qwe >0.00000001) { @@ -608,7 +1031,7 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, gp_Pnt ptref; IntPatch_Point pmult1, pmult2; - elipsol = inter.Ellipse(1); + elipsol = theInter.Ellipse(1); gp_Pnt pttang1(ElCLib::Value(0.5*M_PI, elipsol)); gp_Pnt pttang2(ElCLib::Value(1.5*M_PI, elipsol)); @@ -620,18 +1043,39 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, pmult2.SetMultiple(Standard_True); Standard_Real oU1,oV1,oU2,oV2; + + if(isTheReverse) + { + Quad2.Parameters(pttang1,oU1,oV1); + Quad1.Parameters(pttang1,oU2,oV2); + } + else + { Quad1.Parameters(pttang1,oU1,oV1); Quad2.Parameters(pttang1,oU2,oV2); + } + pmult1.SetParameters(oU1,oV1,oU2,oV2); + if(isTheReverse) + { + Quad2.Parameters(pttang2,oU1,oV1); + Quad1.Parameters(pttang2,oU2,oV2); + } + else + { Quad1.Parameters(pttang2,oU1,oV1); Quad2.Parameters(pttang2,oU2,oV2); + } + pmult2.SetParameters(oU1,oV1,oU2,oV2); // on traite la premiere ellipse //-- Calcul de la Transition de la ligne ElCLib::D1(0.,elipsol,ptref,Tgt); + + //Theoretically, qwe = +/- |Tgt|. Standard_Real qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)); if (qwe>0.00000001) { @@ -660,8 +1104,18 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, aIP.SetValue(aP,Tol,Standard_False); aIP.SetMultiple(Standard_False); // + + if(isTheReverse) + { + Quad2.Parameters(aP, aU1, aV1); + Quad1.Parameters(aP, aU2, aV2); + } + else + { Quad1.Parameters(aP, aU1, aV1); Quad2.Parameters(aP, aU2, aV2); + } + aIP.SetParameters(aU1, aV1, aU2, aV2); // aIP.SetParameter(0.); @@ -685,7 +1139,7 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, //-- Transitions calculee au point 0 OK // // on traite la deuxieme ellipse - elipsol = inter.Ellipse(2); + elipsol = theInter.Ellipse(2); Standard_Real param1 = ElCLib::Parameter(elipsol,pttang1); Standard_Real param2 = ElCLib::Parameter(elipsol,pttang2); @@ -704,6 +1158,8 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, //-- Calcul des transitions de ligne pour la premiere ligne ElCLib::D1(parampourtransition,elipsol,ptref,Tgt); + + //Theoretically, qwe = +/- |Tgt|. qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)); if(qwe> 0.00000001) { @@ -731,439 +1187,112 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, aIP.SetValue(aP,Tol,Standard_False); aIP.SetMultiple(Standard_False); // - Quad1.Parameters(aP, aU1, aV1); - Quad2.Parameters(aP, aU2, aV2); - aIP.SetParameters(aU1, aV1, aU2, aV2); - // - aIP.SetParameter(0.); - glig->AddVertex(aIP); - glig->SetFirstPoint(1); - // - aIP.SetParameter(2.*M_PI); - glig->AddVertex(aIP); - glig->SetLastPoint(2); - } - // - glig->AddVertex(pmult1); - glig->AddVertex(pmult2); - // - slin.Append(glig); - } - break; - - case IntAna_NoGeometricSolution: - { - Standard_Boolean bReverse; - Standard_Real U1,V1,U2,V2; - gp_Pnt psol; - // - bReverse=IsToReverse(Cy1, Cy2, Tol); - if (bReverse) - { - Cy2=Quad1.Cylinder(); - Cy1=Quad2.Cylinder(); - } - // - IntAna_IntQuadQuad anaint(Cy1,Cy2,Tol); - if (!anaint.IsDone()) - { - return Standard_False; - } - - if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) - { - Empty = Standard_True; - } - else - { - NbSol = anaint.NbPnt(); - for (i = 1; i <= NbSol; i++) - { - psol = anaint.Point(i); - Quad1.Parameters(psol,U1,V1); - Quad2.Parameters(psol,U2,V2); - ptsol.SetValue(psol,Tol,Standard_True); - ptsol.SetParameters(U1,V1,U2,V2); - spnt.Append(ptsol); - } - - gp_Pnt ptvalid, ptf, ptl; - gp_Vec tgvalid; - - Standard_Real first,last,para; - IntAna_Curve curvsol; - Standard_Boolean tgfound; - Standard_Integer kount; - - NbSol = anaint.NbCurve(); - for (i = 1; i <= NbSol; i++) - { - curvsol = anaint.Curve(i); - curvsol.Domain(first,last); - ptf = curvsol.Value(first); - ptl = curvsol.Value(last); - - para = last; - kount = 1; - tgfound = Standard_False; - - while (!tgfound) - { - para = (1.123*first + para)/2.123; - tgfound = curvsol.D1u(para,ptvalid,tgvalid); - if (!tgfound) - { - kount ++; - tgfound = kount > 5; - } - } - - Handle(IntPatch_ALine) alig; - if (kount <= 5) - { - Standard_Real qwe = tgvalid.DotCross( Quad2.Normale(ptvalid), - Quad1.Normale(ptvalid)); - if(qwe>0.00000001) - { - trans1 = IntSurf_Out; - trans2 = IntSurf_In; - } - else if(qwe<-0.00000001) - { - trans1 = IntSurf_In; - trans2 = IntSurf_Out; - } - else - { - trans1=trans2=IntSurf_Undecided; - } - alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2); - } - else - { - alig = new IntPatch_ALine(curvsol,Standard_False); - //-- cout << "Transition indeterminee" << endl; - } - - Standard_Boolean TempFalse1 = Standard_False; - Standard_Boolean TempFalse2 = Standard_False; - - ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1,ptf,first, - TempFalse2,ptl,last,Multpoint,Tol); - slin.Append(alig); - } - } - } - break; - - default: - return Standard_False; - } - - return Standard_True; -} - -//======================================================================= -//function : ShortCosForm -//purpose : Represents theCosFactor*cosA+theSinFactor*sinA as -// theCoeff*cos(A-theAngle) if it is possibly (all angles are -// in radians). -//======================================================================= -static void ShortCosForm( const Standard_Real theCosFactor, - const Standard_Real theSinFactor, - Standard_Real& theCoeff, - Standard_Real& theAngle) -{ - theCoeff = sqrt(theCosFactor*theCosFactor+theSinFactor*theSinFactor); - theAngle = 0.0; - if(IsEqual(theCoeff, 0.0)) - { - theAngle = 0.0; - return; - } - - theAngle = acos(Abs(theCosFactor/theCoeff)); - - if(theSinFactor > 0.0) - { - if(IsEqual(theCosFactor, 0.0)) - { - theAngle = M_PI/2.0; - } - else if(theCosFactor < 0.0) - { - theAngle = M_PI-theAngle; - } - } - else if(IsEqual(theSinFactor, 0.0)) - { - if(theCosFactor < 0.0) - { - theAngle = M_PI; - } - } - if(theSinFactor < 0.0) - { - if(theCosFactor > 0.0) - { - theAngle = 2.0*M_PI-theAngle; - } - else if(IsEqual(theCosFactor, 0.0)) - { - theAngle = 3.0*M_PI/2.0; - } - else if(theCosFactor < 0.0) - { - theAngle = M_PI+theAngle; - } - } -} - -enum SearchBoundType -{ - SearchNONE = 0, - SearchV1 = 1, - SearchV2 = 2 -}; - -//Stores equations coefficients -struct stCoeffsValue -{ - stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&); - - math_Vector mVecA1; - math_Vector mVecA2; - math_Vector mVecB1; - math_Vector mVecB2; - math_Vector mVecC1; - math_Vector mVecC2; - math_Vector mVecD; - - Standard_Real mK21; //sinU2 - Standard_Real mK11; //sinU1 - Standard_Real mL21; //cosU2 - Standard_Real mL11; //cosU1 - Standard_Real mM1; //Free member - - Standard_Real mK22; //sinU2 - Standard_Real mK12; //sinU1 - Standard_Real mL22; //cosU2 - Standard_Real mL12; //cosU1 - Standard_Real mM2; //Free member - - Standard_Real mK1; - Standard_Real mL1; - Standard_Real mK2; - Standard_Real mL2; - - Standard_Real mFIV1; - Standard_Real mPSIV1; - Standard_Real mFIV2; - Standard_Real mPSIV2; - - Standard_Real mB; - Standard_Real mC; - Standard_Real mFI1; - Standard_Real mFI2; -}; - -stCoeffsValue::stCoeffsValue( const gp_Cylinder& theCyl1, - const gp_Cylinder& theCyl2): - mVecA1(-theCyl1.Radius()*theCyl1.XAxis().Direction().XYZ()), - mVecA2(theCyl2.Radius()*theCyl2.XAxis().Direction().XYZ()), - mVecB1(-theCyl1.Radius()*theCyl1.YAxis().Direction().XYZ()), - mVecB2(theCyl2.Radius()*theCyl2.YAxis().Direction().XYZ()), - mVecC1(theCyl1.Axis().Direction().XYZ()), - mVecC2(theCyl2.Axis().Direction().XYZ().Reversed()), - mVecD(theCyl2.Location().XYZ() - theCyl1.Location().XYZ()) -{ - enum CoupleOfEquation - { - COENONE = 0, - COE12 = 1, - COE23 = 2, - COE13 = 3 - }aFoundCouple = COENONE; - - - Standard_Real aDetV1V2 = 0.0; - - const Standard_Real aDelta1 = mVecC1(1)*mVecC2(2)-mVecC1(2)*mVecC2(1); //1-2 - const Standard_Real aDelta2 = mVecC1(2)*mVecC2(3)-mVecC1(3)*mVecC2(2); //2-3 - const Standard_Real aDelta3 = mVecC1(1)*mVecC2(3)-mVecC1(3)*mVecC2(1); //1-3 - const Standard_Real anAbsD1 = Abs(aDelta1); //1-2 - const Standard_Real anAbsD2 = Abs(aDelta2); //2-3 - const Standard_Real anAbsD3 = Abs(aDelta3); //1-3 - - if(anAbsD1 >= anAbsD2) - { - if(anAbsD3 > anAbsD1) - { - aFoundCouple = COE13; - aDetV1V2 = aDelta3; - } - else - { - aFoundCouple = COE12; - aDetV1V2 = aDelta1; - } - } - else - { - if(anAbsD3 > anAbsD2) - { - aFoundCouple = COE13; - aDetV1V2 = aDelta3; - } - else - { - aFoundCouple = COE23; - aDetV1V2 = aDelta2; - } - } - - // In point of fact, every determinant (aDelta1, aDelta2 and aDelta3) is - // cross-product between directions (i.e. sine of angle). - // If sine is too small then sine is (approx.) equal to angle itself. - // Therefore, in this case we should compare sine with angular tolerance. - // This constant is used for check if axes are parallel (see constructor - // AxeOperator::AxeOperator(...) in IntAna_QuadQuadGeo.cxx file). - if(Abs(aDetV1V2) < Precision::Angular()) - { - Standard_Failure::Raise("Error. Exception in divide by zerro (IntCyCyTrim)!!!!"); - } - - switch(aFoundCouple) - { - case COE12: - break; - case COE23: - { - math_Vector aVTemp(mVecA1); - mVecA1(1) = aVTemp(2); - mVecA1(2) = aVTemp(3); - mVecA1(3) = aVTemp(1); - - aVTemp = mVecA2; - mVecA2(1) = aVTemp(2); - mVecA2(2) = aVTemp(3); - mVecA2(3) = aVTemp(1); - - aVTemp = mVecB1; - mVecB1(1) = aVTemp(2); - mVecB1(2) = aVTemp(3); - mVecB1(3) = aVTemp(1); - - aVTemp = mVecB2; - mVecB2(1) = aVTemp(2); - mVecB2(2) = aVTemp(3); - mVecB2(3) = aVTemp(1); - - aVTemp = mVecC1; - mVecC1(1) = aVTemp(2); - mVecC1(2) = aVTemp(3); - mVecC1(3) = aVTemp(1); - - aVTemp = mVecC2; - mVecC2(1) = aVTemp(2); - mVecC2(2) = aVTemp(3); - mVecC2(3) = aVTemp(1); - - aVTemp = mVecD; - mVecD(1) = aVTemp(2); - mVecD(2) = aVTemp(3); - mVecD(3) = aVTemp(1); - - break; - } - case COE13: - { - math_Vector aVTemp = mVecA1; - mVecA1(2) = aVTemp(3); - mVecA1(3) = aVTemp(2); - - aVTemp = mVecA2; - mVecA2(2) = aVTemp(3); - mVecA2(3) = aVTemp(2); - - aVTemp = mVecB1; - mVecB1(2) = aVTemp(3); - mVecB1(3) = aVTemp(2); - - aVTemp = mVecB2; - mVecB2(2) = aVTemp(3); - mVecB2(3) = aVTemp(2); - - aVTemp = mVecC1; - mVecC1(2) = aVTemp(3); - mVecC1(3) = aVTemp(2); - - aVTemp = mVecC2; - mVecC2(2) = aVTemp(3); - mVecC2(3) = aVTemp(2); - aVTemp = mVecD; - mVecD(2) = aVTemp(3); - mVecD(3) = aVTemp(2); + if(isTheReverse) + { + Quad2.Parameters(aP, aU1, aV1); + Quad1.Parameters(aP, aU2, aV2); + } + else + { + Quad1.Parameters(aP, aU1, aV1); + Quad2.Parameters(aP, aU2, aV2); + } - break; + aIP.SetParameters(aU1, aV1, aU2, aV2); + // + aIP.SetParameter(0.); + glig->AddVertex(aIP); + glig->SetFirstPoint(1); + // + aIP.SetParameter(2.*M_PI); + glig->AddVertex(aIP); + glig->SetLastPoint(2); + } + // + glig->AddVertex(pmult1); + glig->AddVertex(pmult2); + // + slin.Append(glig); } - default: break; - } - - //------- For V1 (begin) - //sinU2 - mK21 = (mVecC2(2)*mVecB2(1)-mVecC2(1)*mVecB2(2))/aDetV1V2; - //sinU1 - mK11 = (mVecC2(2)*mVecB1(1)-mVecC2(1)*mVecB1(2))/aDetV1V2; - //cosU2 - mL21 = (mVecC2(2)*mVecA2(1)-mVecC2(1)*mVecA2(2))/aDetV1V2; - //cosU1 - mL11 = (mVecC2(2)*mVecA1(1)-mVecC2(1)*mVecA1(2))/aDetV1V2; - //Free member - mM1 = (mVecC2(2)*mVecD(1)-mVecC2(1)*mVecD(2))/aDetV1V2; - //------- For V1 (end) - - //------- For V2 (begin) - //sinU2 - mK22 = (mVecC1(1)*mVecB2(2)-mVecC1(2)*mVecB2(1))/aDetV1V2; - //sinU1 - mK12 = (mVecC1(1)*mVecB1(2)-mVecC1(2)*mVecB1(1))/aDetV1V2; - //cosU2 - mL22 = (mVecC1(1)*mVecA2(2)-mVecC1(2)*mVecA2(1))/aDetV1V2; - //cosU1 - mL12 = (mVecC1(1)*mVecA1(2)-mVecC1(2)*mVecA1(1))/aDetV1V2; - //Free member - mM2 = (mVecC1(1)*mVecD(2)-mVecC1(2)*mVecD(1))/aDetV1V2; - //------- For V1 (end) - - ShortCosForm(mL11, mK11, mK1, mFIV1); - ShortCosForm(mL21, mK21, mL1, mPSIV1); - ShortCosForm(mL12, mK12, mK2, mFIV2); - ShortCosForm(mL22, mK22, mL2, mPSIV2); - const Standard_Real aA1=mVecC1(3)*mK21+mVecC2(3)*mK22-mVecB2(3), //sinU2 - aA2=mVecC1(3)*mL21+mVecC2(3)*mL22-mVecA2(3), //cosU2 - aB1=mVecB1(3)-mVecC1(3)*mK11-mVecC2(3)*mK12, //sinU1 - aB2=mVecA1(3)-mVecC1(3)*mL11-mVecC2(3)*mL12; //cosU1 + case IntAna_Parabola: + case IntAna_Hyperbola: + Standard_Failure::Raise("IntCyCy(): Wrong intersection type!"); + + case IntAna_Circle: + // Circle is useful when we will work with trimmed surfaces + // (two cylinders can be tangent by their basises, e.g. circle) + case IntAna_NoGeometricSolution: + default: + return Standard_False; + } - mC =mVecD(3) - mVecC1(3)*mM1 -mVecC2(3)*mM2; //Free + return Standard_True; +} - Standard_Real aA = 0.0; +//======================================================================= +//function : ShortCosForm +//purpose : Represents theCosFactor*cosA+theSinFactor*sinA as +// theCoeff*cos(A-theAngle) if it is possibly (all angles are +// in radians). +//======================================================================= +static void ShortCosForm( const Standard_Real theCosFactor, + const Standard_Real theSinFactor, + Standard_Real& theCoeff, + Standard_Real& theAngle) +{ + theCoeff = sqrt(theCosFactor*theCosFactor+theSinFactor*theSinFactor); + theAngle = 0.0; + if(IsEqual(theCoeff, 0.0)) + { + theAngle = 0.0; + return; + } - ShortCosForm(aB2,aB1,mB,mFI1); - ShortCosForm(aA2,aA1,aA,mFI2); + theAngle = acos(Abs(theCosFactor/theCoeff)); - mB /= aA; - mC /= aA; + if(theSinFactor > 0.0) + { + if(IsEqual(theCosFactor, 0.0)) + { + theAngle = M_PI/2.0; + } + else if(theCosFactor < 0.0) + { + theAngle = M_PI-theAngle; + } + } + else if(IsEqual(theSinFactor, 0.0)) + { + if(theCosFactor < 0.0) + { + theAngle = M_PI; + } + } + if(theSinFactor < 0.0) + { + if(theCosFactor > 0.0) + { + theAngle = 2.0*M_PI-theAngle; + } + else if(IsEqual(theCosFactor, 0.0)) + { + theAngle = 3.0*M_PI/2.0; + } + else if(theCosFactor < 0.0) + { + theAngle = M_PI+theAngle; + } + } } //======================================================================= //function : CylCylMonotonicity //purpose : Determines, if U2(U1) function is increasing. //======================================================================= -static Standard_Boolean CylCylMonotonicity( const Standard_Real theU1par, +Standard_Boolean ComputationMethods::CylCylMonotonicity(const Standard_Real theU1par, const Standard_Integer theWLIndex, const stCoeffsValue& theCoeffs, const Standard_Real thePeriod, @@ -1235,11 +1364,11 @@ static Standard_Boolean CylCylMonotonicity( const Standard_Real theU1par, // esimates the tolerance of U2-computing (estimation result is // assigned to *theDelta value). //======================================================================= -static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par, +Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par, const Standard_Integer theWLIndex, const stCoeffsValue& theCoeffs, Standard_Real& theU2, - Standard_Real* const theDelta = 0) + Standard_Real* const theDelta) { //This formula is got from some experience and can be changed. const Standard_Real aTol0 = Min(10.0*Epsilon(1.0)*theCoeffs.mB, aNulValue); @@ -1308,7 +1437,7 @@ static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par, return Standard_True; } -static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1, +Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1, const Standard_Real theU2, const stCoeffsValue& theCoeffs, Standard_Real& theV1, @@ -1328,7 +1457,7 @@ static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1, } -static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par, +Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par, const Standard_Integer theWLIndex, const stCoeffsValue& theCoeffs, Standard_Real& theU2, @@ -1348,13 +1477,13 @@ static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par, //function : SearchOnVBounds //purpose : //======================================================================= -static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType, - const stCoeffsValue& theCoeffs, +Standard_Boolean WorkWithBoundaries:: + SearchOnVBounds(const SearchBoundType theSBType, const Standard_Real theVzad, const Standard_Real theVInit, const Standard_Real theInitU2, const Standard_Real theInitMainVar, - Standard_Real& theMainVariableValue) + Standard_Real& theMainVariableValue) const { const Standard_Integer aNbDim = 3; const Standard_Real aMaxError = 4.0*M_PI; // two periods @@ -1381,40 +1510,40 @@ static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType, aSinU2 = sin(aU2Prev), aCosU2 = cos(aU2Prev); - math_Vector aVecFreeMem = (theCoeffs.mVecA2 * aU2Prev + - theCoeffs.mVecB2) * aSinU2 - - (theCoeffs.mVecB2 * aU2Prev - - theCoeffs.mVecA2) * aCosU2 + - (theCoeffs.mVecA1 * aMainVarPrev + - theCoeffs.mVecB1) * aSinU1 - - (theCoeffs.mVecB1 * aMainVarPrev - - theCoeffs.mVecA1) * aCosU1 + - theCoeffs.mVecD; + math_Vector aVecFreeMem = (myCoeffs.mVecA2 * aU2Prev + + myCoeffs.mVecB2) * aSinU2 - + (myCoeffs.mVecB2 * aU2Prev - + myCoeffs.mVecA2) * aCosU2 + + (myCoeffs.mVecA1 * aMainVarPrev + + myCoeffs.mVecB1) * aSinU1 - + (myCoeffs.mVecB1 * aMainVarPrev - + myCoeffs.mVecA1) * aCosU1 + + myCoeffs.mVecD; math_Vector aMSum(1, 3); switch(theSBType) { case SearchV1: - aMatr.SetCol(3, theCoeffs.mVecC2); - aMSum = theCoeffs.mVecC1 * theVzad; + aMatr.SetCol(3, myCoeffs.mVecC2); + aMSum = myCoeffs.mVecC1 * theVzad; aVecFreeMem -= aMSum; - aMSum += theCoeffs.mVecC2*anOtherVar; + aMSum += myCoeffs.mVecC2*anOtherVar; break; case SearchV2: - aMatr.SetCol(3, theCoeffs.mVecC1); - aMSum = theCoeffs.mVecC2 * theVzad; + aMatr.SetCol(3, myCoeffs.mVecC1); + aMSum = myCoeffs.mVecC2 * theVzad; aVecFreeMem -= aMSum; - aMSum += theCoeffs.mVecC1*anOtherVar; + aMSum += myCoeffs.mVecC1*anOtherVar; break; default: return Standard_False; } - aMatr.SetCol(1, theCoeffs.mVecA1 * aSinU1 - theCoeffs.mVecB1 * aCosU1); - aMatr.SetCol(2, theCoeffs.mVecA2 * aSinU2 - theCoeffs.mVecB2 * aCosU2); + aMatr.SetCol(1, myCoeffs.mVecA1 * aSinU1 - myCoeffs.mVecB1 * aCosU1); + aMatr.SetCol(2, myCoeffs.mVecA2 * aSinU2 - myCoeffs.mVecB2 * aCosU2); Standard_Real aDetMainSyst = aMatr.Determinant(); @@ -1456,15 +1585,15 @@ static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType, if(anError > anErrorPrev) {//Method diverges. Keep the best result - const Standard_Real aSinU1 = sin(aMainVarPrev), - aCosU1 = cos(aMainVarPrev), - aSinU2 = sin(aU2Prev), - aCosU2 = cos(aU2Prev); - aMSum -= (theCoeffs.mVecA1*aCosU1 + - theCoeffs.mVecB1*aSinU1 + - theCoeffs.mVecA2*aCosU2 + - theCoeffs.mVecB2*aSinU2 + - theCoeffs.mVecD); + const Standard_Real aSinU1Last = sin(aMainVarPrev), + aCosU1Last = cos(aMainVarPrev), + aSinU2Last = sin(aU2Prev), + aCosU2Last = cos(aU2Prev); + aMSum -= (myCoeffs.mVecA1*aCosU1Last + + myCoeffs.mVecB1*aSinU1Last + + myCoeffs.mVecA2*aCosU2Last + + myCoeffs.mVecB2*aSinU2Last + + myCoeffs.mVecD); const Standard_Real aSQNorm = aMSum.Norm2(); return (aSQNorm < aTol2); } @@ -1531,28 +1660,11 @@ Standard_Boolean InscribePoint( const Standard_Real theUfTarget, return Standard_True; } - if(IsEqual(thePeriod, 0.0)) - {//it cannot be inscribed - return Standard_False; - } - - //Make theUGiven nearer to theUfTarget (in order to avoid - //excess iterations) - theUGiven += thePeriod*Floor((theUfTarget-theUGiven)/thePeriod); - Standard_Real aDU = theUGiven - theUfTarget; + const Standard_Real aUf = theUfTarget - theTol2D; + const Standard_Real aUl = aUf + thePeriod; - if(aDU > 0.0) - aDU = -thePeriod; - else - aDU = +thePeriod; + theUGiven = ElCLib::InPeriod(theUGiven, aUf, aUl); - while(((theUGiven - theUfTarget)*aDU < 0.0) && - !((theUfTarget - theUGiven <= theTol2D) && - (theUGiven - theUlTarget <= theTol2D))) - { - theUGiven += aDU; - } - return ((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D)); } @@ -1645,7 +1757,7 @@ static Standard_Boolean ExcludeNearElements(Standard_Real theArr[], //======================================================================= static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, - const stCoeffsValue& theCoeffs, + const ComputationMethods::stCoeffsValue& theCoeffs, const Standard_Boolean isTheReverse, const Standard_Boolean isThePrecise, const gp_Pnt2d& thePntOnSurf1, @@ -1713,9 +1825,22 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1, aPlast.ParametersOnS1(aUl, aVl); if(!theFlBefore && (aU1par <= aUl)) - {//Parameter value must be increased if theFlBefore == FALSE. + { + //Parameter value must be increased if theFlBefore == FALSE. + + aU1par += thePeriodOfSurf1; + + //The condition is as same as in + //InscribePoint(...) function + if((theUfSurf1 - aU1par > theTol2D) || + (aU1par - theUlSurf1 > theTol2D)) + { + //New aU1par is out of target interval. + //Go back to old value. + return Standard_False; } + } //theTol2D is minimal step along parameter changed. //Therefore, if we apply this minimal step two @@ -1760,8 +1885,7 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1, { SeekAdditionalPoints( theQuad1, theQuad2, theLine, theCoeffs, theWLIndex, aDeltaStep, aNbPnts-2, - aNbPnts-1, theUfSurf2, theUlSurf2, - theTol2D, thePeriodOfSurf1, isTheReverse); + aNbPnts-1, theTol2D, thePeriodOfSurf1, isTheReverse); } } @@ -1772,15 +1896,7 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1, //function : AddBoundaryPoint //purpose : //======================================================================= -static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1, - const IntSurf_Quadric& theQuad2, - const Handle(IntPatch_WLine)& theWL, - const stCoeffsValue& theCoeffs, - const Bnd_Box2d& theUVSurf1, - const Bnd_Box2d& theUVSurf2, - const Standard_Real theTol3D, - const Standard_Real theTol2D, - const Standard_Real thePeriod, +void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL, const Standard_Real theU1, const Standard_Real theU2, const Standard_Real theV1, @@ -1788,10 +1904,9 @@ static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1, const Standard_Real theV2, const Standard_Real theV2Prev, const Standard_Integer theWLIndex, - const Standard_Boolean isTheReverse, const Standard_Boolean theFlForce, Standard_Boolean& isTheFound1, - Standard_Boolean& isTheFound2) + Standard_Boolean& isTheFound2) const { Standard_Real aUSurf1f = 0.0, //const aUSurf1l = 0.0, @@ -1802,176 +1917,97 @@ static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1, aVSurf2f = 0.0, aVSurf2l = 0.0; - theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l); - theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l); + myUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l); + myUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l); - SearchBoundType aTS1 = SearchNONE, aTS2 = SearchNONE; - Standard_Real aV1zad = aVSurf1f, aV2zad = aVSurf2f; + const Standard_Integer aSize = 4; + const Standard_Real anArrVzad[aSize] = {aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l}; - Standard_Real anUpar1 = theU1, anUpar2 = theU1; - Standard_Real aVf = theV1, aVl = theV1Prev; + StPInfo aUVPoint[aSize]; - if( (Abs(aVf-aVSurf1f) <= theTol2D) || - ((aVf-aVSurf1f)*(aVl-aVSurf1f) <= 0.0)) - { - aTS1 = SearchV1; - aV1zad = aVSurf1f; - isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1f, theV2, theU2, theU1, anUpar1); - } - else if((Abs(aVf-aVSurf1l) <= theTol2D) || - ((aVf-aVSurf1l)*(aVl-aVSurf1l) <= 0.0)) + for(Standard_Integer anIDSurf = 0; anIDSurf < 4; anIDSurf+=2) { - aTS1 = SearchV1; - aV1zad = aVSurf1l; - isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1l, theV2, theU2, theU1, anUpar1); - } + const Standard_Real aVf = (anIDSurf == 0) ? theV1 : theV2, + aVl = (anIDSurf == 0) ? theV1Prev : theV2Prev; - aVf = theV2; - aVl = theV2Prev; + const SearchBoundType aTS = (anIDSurf == 0) ? SearchV1 : SearchV2; - if((Abs(aVf-aVSurf2f) <= theTol2D) || - ((aVf-aVSurf2f)*(aVl-aVSurf2f) <= 0.0)) - { - aTS2 = SearchV2; - aV2zad = aVSurf2f; - isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2f, theV1, theU2, theU1, anUpar2); - } - else if((Abs(aVf-aVSurf2l) <= theTol2D) || - ((aVf-aVSurf2l)*(aVl-aVSurf2l) <= 0.0)) + for(Standard_Integer anIDBound = 0; anIDBound < 2; anIDBound++) { - aTS2 = SearchV2; - aV2zad = aVSurf2l; - isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2l, theV1, theU2, theU1, anUpar2); - } - if(!isTheFound1 && !isTheFound2) - return Standard_True; + const Standard_Integer anIndex = anIDSurf+anIDBound; - //We increase U1 parameter. Therefore, added point must have U1 parameter less than - //or equal to current (conditions "(anUpar1 < theU1)" and "(anUpar2 < theU1)"). + aUVPoint[anIndex].mySurfID = anIDSurf; - if(anUpar1 < anUpar2) + if((Abs(aVf-anArrVzad[anIndex]) > myTol2D) && + ((aVf-anArrVzad[anIndex])*(aVl-anArrVzad[anIndex]) > 0.0)) { - if(isTheFound1 && (anUpar1 < theU1)) - { - Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2; - if(!CylCylComputeParameters(anUpar1, theWLIndex, theCoeffs, aU2, aV1, aV2)) - { - isTheFound1 = Standard_False; + continue; } - if(aTS1 == SearchV1) - aV1 = aV1zad; - - if(aTS1 == SearchV2) - aV2 = aV2zad; - - if(isTheFound1 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False, - gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2), - aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l, - aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod, - theWL->Curve(), theWLIndex, theTol3D, - theTol2D, theFlForce)) + const Standard_Boolean aRes = SearchOnVBounds(aTS, anArrVzad[anIndex], + (anIDSurf == 0) ? theV2 : theV1, + theU2, theU1, + aUVPoint[anIndex].myU1); + + if(!aRes || aUVPoint[anIndex].myU1 >= theU1) { - isTheFound1 = Standard_False; + aUVPoint[anIndex].myU1 = RealLast(); + continue; } - } else { - isTheFound1 = Standard_False; - } + Standard_Real &aU1 = aUVPoint[anIndex].myU1, + &aU2 = aUVPoint[anIndex].myU2, + &aV1 = aUVPoint[anIndex].myV1, + &aV2 = aUVPoint[anIndex].myV2; + aU2 = theU2; + aV1 = theV1; + aV2 = theV2; - if(isTheFound2 && (anUpar2 < theU1)) + if(!ComputationMethods:: + CylCylComputeParameters(aU1, theWLIndex, myCoeffs, aU2, aV1, aV2)) { - Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2; - if(!CylCylComputeParameters(anUpar2, theWLIndex, theCoeffs, aU2, aV1, aV2)) - { - isTheFound2 = Standard_False; + aU1 = RealLast(); + continue; } - if(aTS2 == SearchV1) - aV1 = aV1zad; - - if(aTS2 == SearchV2) - aV2 = aV2zad; - - if(isTheFound2 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False, - gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2), - aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l, - aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod, - theWL->Curve(), theWLIndex, theTol3D, - theTol2D, theFlForce)) - { - isTheFound2 = Standard_False; + if(aTS == SearchV1) + aV1 = anArrVzad[anIndex]; + else //if(aTS[anIndex] == SearchV2) + aV2 = anArrVzad[anIndex]; } } - else - { - isTheFound2 = Standard_False; } - } - else - { - if(isTheFound2 && (anUpar2 < theU1)) - { - Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2; - if(!CylCylComputeParameters(anUpar2, theWLIndex, theCoeffs, aU2, aV1, aV2)) - { - isTheFound2 = Standard_False; - } - if(aTS2 == SearchV1) - aV1 = aV1zad; + std::sort(aUVPoint, aUVPoint+aSize); - if(aTS2 == SearchV2) - aV2 = aV2zad; + isTheFound1 = isTheFound2 = Standard_False; - if(isTheFound2 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False, - gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2), - aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l, - aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod, - theWL->Curve(), theWLIndex, theTol3D, - theTol2D, theFlForce)) + for(Standard_Integer i = 0; i < aSize; i++) { - isTheFound2 = Standard_False; - } - } - else - { - isTheFound2 = Standard_False; - } + if(aUVPoint[i].myU1 == RealLast()) + break; - if(isTheFound1 && (anUpar1 < theU1)) + if(!AddPointIntoWL(myQuad1, myQuad2, myCoeffs, myIsReverse, Standard_False, + gp_Pnt2d(aUVPoint[i].myU1, aUVPoint[i].myV1), + gp_Pnt2d(aUVPoint[i].myU2, aUVPoint[i].myV2), + aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l, + aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, myPeriod, + theWL->Curve(), theWLIndex, myTol3D, myTol2D, theFlForce)) { - Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2; - if(!CylCylComputeParameters(anUpar1, theWLIndex, theCoeffs, aU2, aV1, aV2)) - { - isTheFound1 = Standard_False; + continue; } - if(aTS1 == SearchV1) - aV1 = aV1zad; - - if(aTS1 == SearchV2) - aV2 = aV2zad; - - if(isTheFound1 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False, - gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2), - aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l, - aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod, - theWL->Curve(), theWLIndex, theTol3D, - theTol2D, theFlForce)) + if(aUVPoint[i].mySurfID == 0) { - isTheFound1 = Standard_False; + isTheFound1 = Standard_True; } - } else { - isTheFound1 = Standard_False; + isTheFound2 = Standard_True; } } - - return Standard_True; } //======================================================================= @@ -1981,13 +2017,11 @@ static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1, static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, const Handle(IntSurf_LineOn2S)& theLine, - const stCoeffsValue& theCoeffs, + const ComputationMethods::stCoeffsValue& theCoeffs, const Standard_Integer theWLIndex, const Standard_Integer theMinNbPoints, const Standard_Integer theStartPointOnLine, const Standard_Integer theEndPointOnLine, - const Standard_Real theU2f, - const Standard_Real theU2l, const Standard_Real theTol2D, const Standard_Real thePeriodOfSurf2, const Standard_Boolean isTheReverse) @@ -2062,16 +2096,23 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1, U1prec = 0.5*(U1f+U1l); - if(!CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec)) + if(!ComputationMethods:: + CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec)) + { continue; + } - InscribePoint(theU2f, theU2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False); + MinMax(U2f, U2l); + if(!InscribePoint(U2f, U2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False)) + { + continue; + } const gp_Pnt aP1(theQuad1.Value(U1prec, V1prec)), aP2(theQuad2.Value(U2prec, V2prec)); const gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ())); -#ifdef OCCT_DEBUG - //cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << endl; +#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG + std::cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << std::endl; #endif IntSurf_PntOn2S anIP; @@ -2102,27 +2143,24 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1, //purpose : Computes true domain of future intersection curve (allows // avoiding excess iterations) //======================================================================= -//======================================================================= -static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, - const Standard_Real thePeriod, - Standard_Real theU1f[], - Standard_Real theU1l[]) +Standard_Boolean WorkWithBoundaries::BoundariesComputing(Standard_Real theU1f[], + Standard_Real theU1l[]) const { - if(theCoeffs.mB > 0.0) + if(myCoeffs.mB > 0.0) { - if(theCoeffs.mB + Abs(theCoeffs.mC) < -1.0) + if(myCoeffs.mB + Abs(myCoeffs.mC) < -1.0) {//There is NOT intersection return Standard_False; } - else if(theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0) + else if(myCoeffs.mB + Abs(myCoeffs.mC) <= 1.0) {//U=[0;2*PI]+aFI1 - theU1f[0] = theCoeffs.mFI1; - theU1l[0] = thePeriod + theCoeffs.mFI1; + theU1f[0] = myCoeffs.mFI1; + theU1l[0] = myPeriod + myCoeffs.mFI1; } - else if((1 + theCoeffs.mC <= theCoeffs.mB) && - (theCoeffs.mB <= 1 - theCoeffs.mC)) + else if((1 + myCoeffs.mC <= myCoeffs.mB) && + (myCoeffs.mB <= 1 - myCoeffs.mC)) { - Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB; + Standard_Real anArg = -(myCoeffs.mC + 1) / myCoeffs.mB; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) @@ -2130,15 +2168,15 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, const Standard_Real aDAngle = acos(anArg); //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1) - theU1f[0] = theCoeffs.mFI1; - theU1l[0] = aDAngle + theCoeffs.mFI1; - theU1f[1] = thePeriod - aDAngle + theCoeffs.mFI1; - theU1l[1] = thePeriod + theCoeffs.mFI1; + theU1f[0] = myCoeffs.mFI1; + theU1l[0] = aDAngle + myCoeffs.mFI1; + theU1f[1] = myPeriod - aDAngle + myCoeffs.mFI1; + theU1l[1] = myPeriod + myCoeffs.mFI1; } - else if((1 - theCoeffs.mC <= theCoeffs.mB) && - (theCoeffs.mB <= 1 + theCoeffs.mC)) + else if((1 - myCoeffs.mC <= myCoeffs.mB) && + (myCoeffs.mB <= 1 + myCoeffs.mC)) { - Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB; + Standard_Real anArg = (1 - myCoeffs.mC) / myCoeffs.mB; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) @@ -2147,13 +2185,13 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, const Standard_Real aDAngle = acos(anArg); //U=[aDAngle;2*PI-aDAngle]+aFI1 - theU1f[0] = aDAngle + theCoeffs.mFI1; - theU1l[0] = thePeriod - aDAngle + theCoeffs.mFI1; + theU1f[0] = aDAngle + myCoeffs.mFI1; + theU1l[0] = myPeriod - aDAngle + myCoeffs.mFI1; } - else if(theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0) + else if(myCoeffs.mB - Abs(myCoeffs.mC) >= 1.0) { - Standard_Real anArg1 = (1 - theCoeffs.mC) / theCoeffs.mB, - anArg2 = -(theCoeffs.mC + 1) / theCoeffs.mB; + Standard_Real anArg1 = (1 - myCoeffs.mC) / myCoeffs.mB, + anArg2 = -(myCoeffs.mC + 1) / myCoeffs.mB; if(anArg1 > 1.0) anArg1 = 1.0; if(anArg1 < -1.0) @@ -2168,31 +2206,31 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, //(U=[aDAngle1;aDAngle2]+aFI1) || //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1) - theU1f[0] = aDAngle1 + theCoeffs.mFI1; - theU1l[0] = aDAngle2 + theCoeffs.mFI1; - theU1f[1] = thePeriod - aDAngle2 + theCoeffs.mFI1; - theU1l[1] = thePeriod - aDAngle1 + theCoeffs.mFI1; + theU1f[0] = aDAngle1 + myCoeffs.mFI1; + theU1l[0] = aDAngle2 + myCoeffs.mFI1; + theU1f[1] = myPeriod - aDAngle2 + myCoeffs.mFI1; + theU1l[1] = myPeriod - aDAngle1 + myCoeffs.mFI1; } else { return Standard_False; } } - else if(theCoeffs.mB < 0.0) + else if(myCoeffs.mB < 0.0) { - if(theCoeffs.mB + Abs(theCoeffs.mC) > 1.0) + if(myCoeffs.mB + Abs(myCoeffs.mC) > 1.0) {//There is NOT intersection return Standard_False; } - else if(-theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0) + else if(-myCoeffs.mB + Abs(myCoeffs.mC) <= 1.0) {//U=[0;2*PI]+aFI1 - theU1f[0] = theCoeffs.mFI1; - theU1l[0] = thePeriod + theCoeffs.mFI1; + theU1f[0] = myCoeffs.mFI1; + theU1l[0] = myPeriod + myCoeffs.mFI1; } - else if((-theCoeffs.mC - 1 <= theCoeffs.mB) && - ( theCoeffs.mB <= theCoeffs.mC - 1)) + else if((-myCoeffs.mC - 1 <= myCoeffs.mB) && + ( myCoeffs.mB <= myCoeffs.mC - 1)) { - Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB; + Standard_Real anArg = (1 - myCoeffs.mC) / myCoeffs.mB; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) @@ -2201,15 +2239,15 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, const Standard_Real aDAngle = acos(anArg); //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1) - theU1f[0] = theCoeffs.mFI1; - theU1l[0] = aDAngle + theCoeffs.mFI1; - theU1f[1] = thePeriod - aDAngle + theCoeffs.mFI1; - theU1l[1] = thePeriod + theCoeffs.mFI1; + theU1f[0] = myCoeffs.mFI1; + theU1l[0] = aDAngle + myCoeffs.mFI1; + theU1f[1] = myPeriod - aDAngle + myCoeffs.mFI1; + theU1l[1] = myPeriod + myCoeffs.mFI1; } - else if((theCoeffs.mC - 1 <= theCoeffs.mB) && - (theCoeffs.mB <= -theCoeffs.mB - 1)) + else if((myCoeffs.mC - 1 <= myCoeffs.mB) && + (myCoeffs.mB <= -myCoeffs.mB - 1)) { - Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB; + Standard_Real anArg = -(myCoeffs.mC + 1) / myCoeffs.mB; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) @@ -2218,13 +2256,13 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, const Standard_Real aDAngle = acos(anArg); //U=[aDAngle;2*PI-aDAngle]+aFI1 - theU1f[0] = aDAngle + theCoeffs.mFI1; - theU1l[0] = thePeriod - aDAngle + theCoeffs.mFI1; + theU1f[0] = aDAngle + myCoeffs.mFI1; + theU1l[0] = myPeriod - aDAngle + myCoeffs.mFI1; } - else if(-theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0) + else if(-myCoeffs.mB - Abs(myCoeffs.mC) >= 1.0) { - Standard_Real anArg1 = -(theCoeffs.mC + 1) / theCoeffs.mB, - anArg2 = (1 - theCoeffs.mC) / theCoeffs.mB; + Standard_Real anArg1 = -(myCoeffs.mC + 1) / myCoeffs.mB, + anArg2 = (1 - myCoeffs.mC) / myCoeffs.mB; if(anArg1 > 1.0) anArg1 = 1.0; if(anArg1 < -1.0) @@ -2239,10 +2277,10 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, //(U=[aDAngle1;aDAngle2]+aFI1) || //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1) - theU1f[0] = aDAngle1 + theCoeffs.mFI1; - theU1l[0] = aDAngle2 + theCoeffs.mFI1; - theU1f[1] = thePeriod - aDAngle2 + theCoeffs.mFI1; - theU1l[1] = thePeriod - aDAngle1 + theCoeffs.mFI1; + theU1f[0] = aDAngle1 + myCoeffs.mFI1; + theU1l[0] = aDAngle2 + myCoeffs.mFI1; + theU1f[1] = myPeriod - aDAngle2 + myCoeffs.mFI1; + theU1l[1] = myPeriod - aDAngle1 + myCoeffs.mFI1; } else { @@ -2262,7 +2300,7 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, //purpose : theNbCritPointsMax contains true number of critical points. // It must be initialized correctly before function calling //======================================================================= -static void CriticalPointsComputing(const stCoeffsValue& theCoeffs, +static void CriticalPointsComputing(const ComputationMethods::stCoeffsValue& theCoeffs, const Standard_Real theUSurf1f, const Standard_Real theUSurf1l, const Standard_Real theUSurf2f, @@ -2385,10 +2423,72 @@ static void CriticalPointsComputing(const stCoeffsValue& theCoeffs, } //======================================================================= -//function : IntCyCyTrim +//function : BoundaryEstimation +//purpose : Rough estimation of the parameter range. +//======================================================================= +void WorkWithBoundaries::BoundaryEstimation(const gp_Cylinder& theCy1, + const gp_Cylinder& theCy2, + Bnd_Range& theOutBoxS1, + Bnd_Range& theOutBoxS2) const +{ + const gp_Dir &aD1 = theCy1.Axis().Direction(), + &aD2 = theCy2.Axis().Direction(); + const Standard_Real aR1 = theCy1.Radius(), + aR2 = theCy2.Radius(); + + //Let consider a parallelogram. Its edges are parallel to aD1 and aD2. + //Its altitudes are equal to 2*aR1 and 2*aR2 (diameters of the cylinders). + //In fact, this parallelogram is a projection of the cylinders to the plane + //created by the intersected axes aD1 and aD2 (if the axes are skewed then + //one axis can be translated by parallel shifting till intersection). + + const Standard_Real aCosA = aD1.Dot(aD2); + const Standard_Real aSqSinA = 1.0 - aCosA * aCosA; + + //If sine is small then it can be compared with angle. + if (aSqSinA < Precision::Angular()*Precision::Angular()) + return; + + //Half of delta V. Delta V is a distance between + //projections of two opposite parallelogram vertices + //(joined by the maximal diagonal) to the cylinder axis. + const Standard_Real aSinA = sqrt(aSqSinA); + const Standard_Real aHDV1 = (aR1 * aCosA + aR2)/aSinA, + aHDV2 = (aR2 * aCosA + aR1)/aSinA; + +#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG + //The code in this block is created for test only.It is stupidly to create + //OCCT-test for the method, which will be changed possibly never. + std::cout << "Reference values: aHDV1 = " << aHDV1 << "; aHDV2 = " << aHDV2 << std::endl; +#endif + + //V-parameters of intersection point of the axes (in case of skewed axes, + //see comment above). + Standard_Real aV01 = 0.0, aV02 = 0.0; + ExtremaLineLine(theCy1.Axis(), theCy2.Axis(), aCosA, aSqSinA, aV01, aV02); + + theOutBoxS1.Add(aV01 - aHDV1); + theOutBoxS1.Add(aV01 + aHDV1); + + theOutBoxS2.Add(aV02 - aHDV2); + theOutBoxS2.Add(aV02 + aHDV2); + + theOutBoxS1.Enlarge(Precision::Confusion()); + theOutBoxS2.Enlarge(Precision::Confusion()); + + Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0; + myUVSurf1.Get(aU1, aV1, aU2, aV2); + theOutBoxS1.Common(Bnd_Range(aV1, aV2)); + + myUVSurf2.Get(aU1, aV1, aU2, aV2); + theOutBoxS2.Common(Bnd_Range(aV1, aV2)); +} + +//======================================================================= +//function : IntCyCy //purpose : //======================================================================= -Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, +Standard_Boolean IntCyCy( const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, const Standard_Real theTol3D, const Standard_Real theTol2D, @@ -2396,23 +2496,19 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, const Bnd_Box2d& theUVSurf2, const Standard_Boolean isTheReverse, Standard_Boolean& isTheEmpty, + Standard_Boolean& isTheSameSurface, + Standard_Boolean& isTheMultiplePoint, IntPatch_SequenceOfLine& theSlin, IntPatch_SequenceOfPoint& theSPnt) { - Standard_Real aUSurf1f = 0.0, //const - aUSurf1l = 0.0, - aVSurf1f = 0.0, - aVSurf1l = 0.0; - Standard_Real aUSurf2f = 0.0, //const - aUSurf2l = 0.0, - aVSurf2f = 0.0, - aVSurf2l = 0.0; - - theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l); - theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l); + isTheEmpty = Standard_True; + isTheSameSurface = Standard_False; + isTheMultiplePoint = Standard_False; + theSlin.Clear(); + theSPnt.Clear(); const gp_Cylinder& aCyl1 = theQuad1.Cylinder(), - aCyl2 = theQuad2.Cylinder(); + &aCyl2 = theQuad2.Cylinder(); IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D); @@ -2421,17 +2517,26 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, return Standard_False; } - IntAna_ResultType aTypInt = anInter.TypeInter(); + if(anInter.TypeInter() != IntAna_NoGeometricSolution) + { + return CyCyAnalyticalIntersect( theQuad1, theQuad2, anInter, + theTol3D, isTheReverse, isTheEmpty, + isTheSameSurface, isTheMultiplePoint, + theSlin, theSPnt); + } - if(aTypInt != IntAna_NoGeometricSolution) - { //It is not necessary (because result is an analytic curve) or - //it is impossible to make Walking-line. + Standard_Real aUSurf1f = 0.0, //const + aUSurf1l = 0.0, + aVSurf1f = 0.0, + aVSurf1l = 0.0; + Standard_Real aUSurf2f = 0.0, //const + aUSurf2l = 0.0, + aVSurf2f = 0.0, + aVSurf2l = 0.0; - return Standard_False; - } + theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l); + theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l); - theSlin.Clear(); - theSPnt.Clear(); const Standard_Integer aNbMaxPoints = 2000; const Standard_Integer aNbMinPoints = 200; const Standard_Integer aNbPoints = Min(Max(aNbMinPoints, @@ -2444,14 +2549,23 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, const Standard_Integer aNbWLines = 2; - const stCoeffsValue anEquationCoeffs(aCyl1, aCyl2); + const ComputationMethods::stCoeffsValue anEquationCoeffs(aCyl1, aCyl2); + + const WorkWithBoundaries aBoundWork(theQuad1, theQuad2, anEquationCoeffs, + theUVSurf1, theUVSurf2, aNbWLines, + aPeriod, theTol3D, theTol2D, isTheReverse); + + Bnd_Range aRangeS1, aRangeS2; + aBoundWork.BoundaryEstimation(aCyl1, aCyl2, aRangeS1, aRangeS2); + if (aRangeS1.IsVoid() || aRangeS2.IsVoid()) + return Standard_True; //Boundaries const Standard_Integer aNbOfBoundaries = 2; Standard_Real aU1f[aNbOfBoundaries] = {-Precision::Infinite(), -Precision::Infinite()}; Standard_Real aU1l[aNbOfBoundaries] = {Precision::Infinite(), Precision::Infinite()}; - if(!BoundariesComputing(anEquationCoeffs, aPeriod, aU1f, aU1l)) + if(!aBoundWork.BoundariesComputing(aU1f, aU1l)) return Standard_True; for(Standard_Integer i = 0; i < aNbOfBoundaries; i++) @@ -2545,7 +2659,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, anUexpect[i] = anUf; } - Standard_Real aCriticalDelta[aNbCritPointsMax]; + Standard_Real aCriticalDelta[aNbCritPointsMax] = {0}; for(Standard_Integer aCritPID = 0; aCritPID < aNbCritPoints; aCritPID++) { //We are not intersted in elements of aCriticalDelta array //if their index is greater than or equal to aNbCritPoints @@ -2564,10 +2678,10 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, { anU1 = anU1crit[i]; - for(Standard_Integer i = 0; i < aNbWLines; i++) + for(Standard_Integer j = 0; j < aNbWLines; j++) { - aWLFindStatus[i] = WLFStatus_Broken; - anUexpect[i] = anU1; + aWLFindStatus[j] = WLFStatus_Broken; + anUexpect[j] = anU1; } break; @@ -2620,7 +2734,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, //or on seam-edge strictly (if it is possible). Standard_Real aTol = theTol2D; - CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i], &aTol); + ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs, + aU2[i], &aTol); InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False); aTol = Max(aTol, theTol2D); @@ -2636,7 +2751,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, } else { - CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]); + ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]); InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False); } @@ -2651,7 +2766,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, //correct value. Standard_Boolean isIncreasing = Standard_True; - CylCylMonotonicity(anU1+aStepMin, i, anEquationCoeffs, aPeriod, isIncreasing); + ComputationMethods::CylCylMonotonicity(anU1+aStepMin, i, anEquationCoeffs, + aPeriod, isIncreasing); //If U2(U1) is increasing and U2 is considered to be equal aUSurf2l //then after the next step (when U1 will be increased) U2 will be @@ -2690,7 +2806,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, } } - CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs, aV1[i], aV2[i]); + ComputationMethods::CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs, + aV1[i], aV2[i]); if(isFirst) { @@ -2763,7 +2880,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False; Standard_Boolean isForce = Standard_False; - if((aWLFindStatus[i] == WLFStatus_Absent)) + if (aWLFindStatus[i] == WLFStatus_Absent) { if(((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1-aUSurf1l) < theTol2D)) { @@ -2771,11 +2888,9 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, } } - AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs, - theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod, - anU1, aU2[i], aV1[i], aV1Prev[i], - aV2[i], aV2Prev[i], i, isTheReverse, - isForce, isFound1, isFound2); + aBoundWork.AddBoundaryPoint(aWLine[i], anU1, aU2[i], aV1[i], aV1Prev[i], + aV2[i], aV2Prev[i], i, isForce, + isFound1, isFound2); const Standard_Boolean isPrevVBound = !isVIntersect && ((Abs(aV1Prev[i] - aVSurf1f) <= theTol2D) || @@ -2783,7 +2898,6 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, (Abs(aV2Prev[i] - aVSurf2f) <= theTol2D) || (Abs(aV2Prev[i] - aVSurf2l) <= theTol2D)); - aV1Prev[i] = aV1[i]; aV2Prev[i] = aV2[i]; @@ -2798,7 +2912,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, aWLFindStatus[i] = WLFStatus_Exist; } - if(( aWLFindStatus[i] != WLFStatus_Broken) || (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl)) + if( (aWLFindStatus[i] != WLFStatus_Broken) || + (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl)) { if(aWLine[i]->NbPnts() > 0) { @@ -2872,10 +2987,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False; - AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs, - theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod, - anU1, aU2[i], aV1[i], aV1Prev[i], - aV2[i], aV2Prev[i], i, isTheReverse, + aBoundWork.AddBoundaryPoint(aWLine[i], anU1, aU2[i], aV1[i], aV1Prev[i], + aV2[i], aV2Prev[i], i, Standard_False, isFound1, isFound2); if(isFound1 || isFound2) @@ -2952,7 +3065,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, continue; } - CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs, aU2[i], aV1[i], aV2[i]); + ComputationMethods::CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs, + aU2[i], aV1[i], aV2[i]); AddPointIntoWL( theQuad1, theQuad2, anEquationCoeffs, isTheReverse, Standard_True, gp_Pnt2d(anUmaxAdded, aV1[i]), @@ -2969,8 +3083,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, //Step computing { - const Standard_Real aDeltaV1 = (aVSurf1l - aVSurf1f)/IntToReal(aNbPoints), - aDeltaV2 = (aVSurf2l - aVSurf2f)/IntToReal(aNbPoints); + const Standard_Real aDeltaV1 = aRangeS1.Delta()/IntToReal(aNbPoints), + aDeltaV2 = aRangeS2.Delta()/IntToReal(aNbPoints); math_Matrix aMatr(1, 3, 1, 5); @@ -3078,8 +3192,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, isAddedIntoWL[i] = Standard_True; SeekAdditionalPoints( theQuad1, theQuad2, aWLine[i]->Curve(), anEquationCoeffs, i, aNbPoints, 1, - aWLine[i]->NbPnts(), aUSurf2f, aUSurf2l, - theTol2D, aPeriod, isTheReverse); + aWLine[i]->NbPnts(), theTol2D, aPeriod, + isTheReverse); aWLine[i]->ComputeVertexParameters(theTol3D); theSlin.Append(aWLine[i]); @@ -3090,7 +3204,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, isAddedIntoWL[i] = Standard_False; } -#ifdef OCCT_DEBUG +#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG //aWLine[i]->Dump(); #endif } @@ -3104,8 +3218,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, { for(Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); aNbLin++) { - const Handle(IntPatch_WLine)& aWLine1 = - Handle(IntPatch_WLine)::DownCast(theSlin.Value(aNbLin)); + Handle(IntPatch_WLine) aWLine1 (Handle(IntPatch_WLine):: + DownCast(theSlin.Value(aNbLin))); const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1); const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aWLine1->NbPnts()); @@ -3153,13 +3267,13 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, for (Standard_Integer i = 0; i < aNbWLines; i++) { Standard_Real anU2t = 0.0; - if(!CylCylComputeParameters(anUl, i, anEquationCoeffs, anU2t)) + if(!ComputationMethods::CylCylComputeParameters(anUl, i, anEquationCoeffs, anU2t)) continue; - const Standard_Real aDU = Abs(anU2t - aCurU2); - if(aDU < aDelta) + const Standard_Real aDU2 = Abs(anU2t - aCurU2); + if(aDU2 < aDelta) { - aDelta = aDU; + aDelta = aDU2; anIndex = i; } } @@ -3169,7 +3283,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, { Standard_Real anU2 = 0.0, anV1 = 0.0, anV2 = 0.0; Standard_Boolean isDone = - CylCylComputeParameters(anUf, anIndex, anEquationCoeffs, + ComputationMethods::CylCylComputeParameters(anUf, anIndex, anEquationCoeffs, anU2, anV1, anV2); anUf += aDU; @@ -3178,8 +3292,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, continue; } - if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse, Standard_True, - gp_Pnt2d(anUf, anV1), gp_Pnt2d(anU2, anV2), + if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse, + Standard_True, gp_Pnt2d(anUf, anV1), gp_Pnt2d(anU2, anV2), aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod, aWLine->Curve(), anIndex, theTol3D, @@ -3193,8 +3307,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, { SeekAdditionalPoints( theQuad1, theQuad2, aWLine->Curve(), anEquationCoeffs, anIndex, aNbMinPoints, - 1, aWLine->NbPnts(), aUSurf2f, aUSurf2l, - theTol2D, aPeriod, isTheReverse); + 1, aWLine->NbPnts(), theTol2D, aPeriod, + isTheReverse); aWLine->ComputeVertexParameters(theTol3D); theSlin.Append(aWLine); @@ -3615,7 +3729,7 @@ Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1, ptvalid = curvsol.Value(para); alig = new IntPatch_ALine(curvsol,Standard_False); kept = Standard_True; - //-- cout << "Transition indeterminee" << endl; + //-- std::cout << "Transition indeterminee" << std::endl; } if (kept) { Standard_Boolean Nfirstp = !firstp; @@ -3697,46 +3811,3 @@ Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy, // return bFind; } -//======================================================================= -//function : IsToReverse -//purpose : -//======================================================================= -Standard_Boolean IsToReverse(const gp_Cylinder& Cy1, - const gp_Cylinder& Cy2, - const Standard_Real Tol) -{ - Standard_Boolean bRet; - Standard_Real aR1,aR2, dR, aSc1, aSc2; - // - bRet=Standard_False; - // - aR1=Cy1.Radius(); - aR2=Cy2.Radius(); - dR=aR1-aR2; - if (dR<0.) { - dR=-dR; - } - if (dR>Tol) { - bRet=aR1>aR2; - } - else { - gp_Dir aDZ(0.,0.,1.); - // - const gp_Dir& aDir1=Cy1.Axis().Direction(); - aSc1=aDir1*aDZ; - if (aSc1<0.) { - aSc1=-aSc1; - } - // - const gp_Dir& aDir2=Cy2.Axis().Direction(); - aSc2=aDir2*aDZ; - if (aSc2<0.) { - aSc2=-aSc2; - } - // - if (aSc2 #include +#include + #include #define DEBUG 0 static const Standard_Integer aNbPointsInALine = 200; -//======================================================================= -//function : MinMax -//purpose : Replaces theParMIN = MIN(theParMIN, theParMAX), -// theParMAX = MAX(theParMIN, theParMAX). -//======================================================================= -static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX) -{ - if(theParMIN > theParMAX) - { - const Standard_Real aTmp = theParMAX; - theParMAX = theParMIN; - theParMIN = aTmp; - } -} - -//======================================================================= -//function : IsSeam -//purpose : Returns: -// 0 - if interval [theU1, theU2] does not intersect the "seam-edge" -// or if "seam-edge" do not exist; -// 1 - if interval (theU1, theU2) intersect the "seam-edge". -// 2 - if theU1 or/and theU2 lie ON the "seam-edge" -// -//ATTENTION!!! -// If (theU1 == theU2) then this function will return only both 0 or 2. -//======================================================================= -static Standard_Integer IsSeam( const Standard_Real theU1, - const Standard_Real theU2, - const Standard_Real thePeriod) -{ - if(IsEqual(thePeriod, 0.0)) - return 0; - - //If interval [theU1, theU2] intersect seam-edge then there exists an integer - //number N such as - // (theU1 <= T*N <= theU2) <=> (theU1/T <= N <= theU2/T), - //where T is the period. - //I.e. the inerval [theU1/T, theU2/T] must contain at least one - //integer number. In this case, Floor(theU1/T) and Floor(theU2/T) - //return different values or theU1/T is strictly integer number. - //Examples: - // 1. theU1/T==2.8, theU2/T==3.5 => Floor(theU1/T) == 2, Floor(theU2/T) == 3. - // 2. theU1/T==2.0, theU2/T==2.6 => Floor(theU1/T) == Floor(theU2/T) == 2. - - const Standard_Real aVal1 = theU1/thePeriod, - aVal2 = theU2/thePeriod; - const Standard_Integer aPar1 = static_cast(Floor(aVal1)); - const Standard_Integer aPar2 = static_cast(Floor(aVal2)); - if(aPar1 != aPar2) - {//Interval (theU1, theU2] intersects seam-edge - if(IsEqual(aVal2, static_cast(aPar2))) - {//aVal2 is an integer number => theU2 lies ON the "seam-edge" - return 2; - } - - return 1; - } - - //Here, aPar1 == aPar2. - - if(IsEqual(aVal1, static_cast(aPar1))) - {//aVal1 is an integer number => theU1 lies ON the "seam-edge" - return 2; - } - - //If aVal2 is a true integer number then always (aPar1 != aPar2). - - return 0; -} - -//======================================================================= -//function : IsSeamOrBound -//purpose : Returns TRUE if segment [thePtf, thePtl] intersects "seam-edge" -// (if it exist) or surface boundaries and both thePtf and thePtl do -// not match "seam-edge" or boundaries. -// Point thePtmid lies in this segment. If thePtmid match -// "seam-edge" or boundaries strictly (without any tolerance) then -// the function will return TRUE. -// See comments in function body for detail information. -//======================================================================= -static Standard_Boolean IsSeamOrBound(const IntSurf_PntOn2S& thePtf, - const IntSurf_PntOn2S& thePtl, - const IntSurf_PntOn2S& thePtmid, - const Standard_Real theU1Period, - const Standard_Real theU2Period, - const Standard_Real theV1Period, - const Standard_Real theV2Period, - const Standard_Real theUfSurf1, - const Standard_Real theUlSurf1, - const Standard_Real theVfSurf1, - const Standard_Real theVlSurf1, - const Standard_Real theUfSurf2, - const Standard_Real theUlSurf2, - const Standard_Real theVfSurf2, - const Standard_Real theVlSurf2) -{ - Standard_Real aU11 = 0.0, aU12 = 0.0, aV11 = 0.0, aV12 = 0.0; - Standard_Real aU21 = 0.0, aU22 = 0.0, aV21 = 0.0, aV22 = 0.0; - thePtf.Parameters(aU11, aV11, aU12, aV12); - thePtl.Parameters(aU21, aV21, aU22, aV22); - - MinMax(aU11, aU21); - MinMax(aV11, aV21); - MinMax(aU12, aU22); - MinMax(aV12, aV22); - - if((aU11 - theUfSurf1)*(aU21 - theUfSurf1) < 0.0) - {//Interval [aU11, aU21] intersects theUfSurf1 - return Standard_True; - } - - if((aU11 - theUlSurf1)*(aU21 - theUlSurf1) < 0.0) - {//Interval [aU11, aU21] intersects theUlSurf1 - return Standard_True; - } - - if((aV11 - theVfSurf1)*(aV21 - theVfSurf1) < 0.0) - {//Interval [aV11, aV21] intersects theVfSurf1 - return Standard_True; - } - - if((aV11 - theVlSurf1)*(aV21 - theVlSurf1) < 0.0) - {//Interval [aV11, aV21] intersects theVlSurf1 - return Standard_True; - } - - if((aU12 - theUfSurf2)*(aU22 - theUfSurf2) < 0.0) - {//Interval [aU12, aU22] intersects theUfSurf2 - return Standard_True; - } - - if((aU12 - theUlSurf2)*(aU22 - theUlSurf2) < 0.0) - {//Interval [aU12, aU22] intersects theUlSurf2 - return Standard_True; - } - - if((aV12 - theVfSurf2)*(aV22 - theVfSurf2) < 0.0) - {//Interval [aV12, aV22] intersects theVfSurf2 - return Standard_True; - } - - if((aV12 - theVlSurf2)*(aV22 - theVlSurf2) < 0.0) - {//Interval [aV12, aV22] intersects theVlSurf2 - return Standard_True; - } - - if(IsSeam(aU11, aU21, theU1Period)) - return Standard_True; - - if(IsSeam(aV11, aV21, theV1Period)) - return Standard_True; - - if(IsSeam(aU12, aU22, theU2Period)) - return Standard_True; - - if(IsSeam(aV12, aV22, theV2Period)) - return Standard_True; - - /* - The segment [thePtf, thePtl] does not intersect the boundaries and - the seam-edge of the surfaces. - Nevertheless, following situation is possible: - - seam or - bound - | - thePtf * | - | - * thePtmid - thePtl * | - | - - This case must be processed, too. - */ - - Standard_Real aU1 = 0.0, aU2 = 0.0, aV1 = 0.0, aV2 = 0.0; - thePtmid.Parameters(aU1, aV1, aU2, aV2); - - if(IsEqual(aU1, theUfSurf1) || IsEqual(aU1, theUlSurf1)) - return Standard_True; - - if(IsEqual(aU2, theUfSurf2) || IsEqual(aU2, theUlSurf2)) - return Standard_True; - - if(IsEqual(aV1, theVfSurf1) || IsEqual(aV1, theVlSurf1)) - return Standard_True; - - if(IsEqual(aV2, theVfSurf2) || IsEqual(aV2, theVlSurf2)) - return Standard_True; - - if(IsSeam(aU1, aU1, theU1Period)) - return Standard_True; - - if(IsSeam(aU2, aU2, theU2Period)) - return Standard_True; - - if(IsSeam(aV1, aV1, theV1Period)) - return Standard_True; - - if(IsSeam(aV2, aV2, theV2Period)) - return Standard_True; - - return Standard_False; -} - -//======================================================================= -//function : JoinWLines -//purpose : joins all WLines from theSlin to one if it is possible and -// records the result into theSlin again. -// Lines will be kept to be splitted if: -// a) they are separated (has no common points); -// b) resulted line (after joining) go through -// seam-edges or surface boundaries. -// -// In addition, if points in theSPnt lies at least in one of -// the line in theSlin, this point will be deleted. -//======================================================================= -static void JoinWLines(IntPatch_SequenceOfLine& theSlin, - IntPatch_SequenceOfPoint& theSPnt, - const Standard_Real theTol3D, - const Standard_Real theU1Period, - const Standard_Real theU2Period, - const Standard_Real theV1Period, - const Standard_Real theV2Period, - const Standard_Real theUfSurf1, - const Standard_Real theUlSurf1, - const Standard_Real theVfSurf1, - const Standard_Real theVlSurf1, - const Standard_Real theUfSurf2, - const Standard_Real theUlSurf2, - const Standard_Real theVfSurf2, - const Standard_Real theVlSurf2) -{ - if(theSlin.Length() == 0) - return; - - for(Standard_Integer aNumOfLine1 = 1; aNumOfLine1 <= theSlin.Length(); aNumOfLine1++) - { - const Handle(IntPatch_WLine)& aWLine1 = Handle(IntPatch_WLine)::DownCast(theSlin.Value(aNumOfLine1)); - - if(aWLine1.IsNull()) - {//We must have failed to join not-point-lines - return; - } - - const Standard_Integer aNbPntsWL1 = aWLine1->NbPnts(); - const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1); - const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aNbPntsWL1); - - for(Standard_Integer aNPt = 1; aNPt <= theSPnt.Length(); aNPt++) - { - const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNPt).PntOn2S(); - - if( aPntCur.IsSame(aPntFWL1, Precision::Confusion()) || - aPntCur.IsSame(aPntLWL1, Precision::Confusion())) - { - theSPnt.Remove(aNPt); - aNPt--; - } - } - - Standard_Boolean hasBeenRemoved = Standard_False; - for(Standard_Integer aNumOfLine2 = aNumOfLine1 + 1; aNumOfLine2 <= theSlin.Length(); aNumOfLine2++) - { - const Handle(IntPatch_WLine)& aWLine2 = Handle(IntPatch_WLine)::DownCast(theSlin.Value(aNumOfLine2)); - - const Standard_Integer aNbPntsWL2 = aWLine2->NbPnts(); - - const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1); - const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aNbPntsWL1); - - const IntSurf_PntOn2S& aPntFWL2 = aWLine2->Point(1); - const IntSurf_PntOn2S& aPntLWL2 = aWLine2->Point(aNbPntsWL2); - - if(aPntFWL1.IsSame(aPntFWL2, Precision::Confusion())) - { - const IntSurf_PntOn2S& aPt1 = aWLine1->Point(2); - const IntSurf_PntOn2S& aPt2 = aWLine2->Point(2); - if(!IsSeamOrBound(aPt1, aPt2, aPntFWL1, theU1Period, theU2Period, - theV1Period, theV2Period, theUfSurf1, theUlSurf1, - theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2, - theVfSurf2, theVlSurf2)) - { - aWLine1->ClearVertexes(); - for(Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++) - { - const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt); - aWLine1->Curve()->InsertBefore(1, aPt); - } - - aWLine1->ComputeVertexParameters(theTol3D); - - theSlin.Remove(aNumOfLine2); - aNumOfLine2--; - hasBeenRemoved = Standard_True; - - continue; - } - } - - if(aPntFWL1.IsSame(aPntLWL2, Precision::Confusion())) - { - const IntSurf_PntOn2S& aPt1 = aWLine1->Point(2); - const IntSurf_PntOn2S& aPt2 = aWLine2->Point(aNbPntsWL2-1); - if(!IsSeamOrBound(aPt1, aPt2, aPntFWL1, theU1Period, theU2Period, - theV1Period, theV2Period, theUfSurf1, theUlSurf1, - theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2, - theVfSurf2, theVlSurf2)) - { - aWLine1->ClearVertexes(); - for(Standard_Integer aNPt = aNbPntsWL2; aNPt >= 1; aNPt--) - { - const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt); - aWLine1->Curve()->InsertBefore(1, aPt); - } - - aWLine1->ComputeVertexParameters(theTol3D); - - theSlin.Remove(aNumOfLine2); - aNumOfLine2--; - hasBeenRemoved = Standard_True; - - continue; - } - } - - if(aPntLWL1.IsSame(aPntFWL2, Precision::Confusion())) - { - const IntSurf_PntOn2S& aPt1 = aWLine1->Point(aNbPntsWL1-1); - const IntSurf_PntOn2S& aPt2 = aWLine2->Point(2); - if(!IsSeamOrBound(aPt1, aPt2, aPntLWL1, theU1Period, theU2Period, - theV1Period, theV2Period, theUfSurf1, theUlSurf1, - theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2, - theVfSurf2, theVlSurf2)) - { - aWLine1->ClearVertexes(); - for(Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++) - { - const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt); - aWLine1->Curve()->Add(aPt); - } - - aWLine1->ComputeVertexParameters(theTol3D); - - theSlin.Remove(aNumOfLine2); - aNumOfLine2--; - hasBeenRemoved = Standard_True; - - continue; - } - } - - if(aPntLWL1.IsSame(aPntLWL2, Precision::Confusion())) - { - const IntSurf_PntOn2S& aPt1 = aWLine1->Point(aNbPntsWL1-1); - const IntSurf_PntOn2S& aPt2 = aWLine2->Point(aNbPntsWL2-1); - if(!IsSeamOrBound(aPt1, aPt2, aPntLWL1, theU1Period, theU2Period, - theV1Period, theV2Period, theUfSurf1, theUlSurf1, - theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2, - theVfSurf2, theVlSurf2)) - { - aWLine1->ClearVertexes(); - for(Standard_Integer aNPt = aNbPntsWL2; aNPt >= 1; aNPt--) - { - const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt); - aWLine1->Curve()->Add(aPt); - } - - aWLine1->ComputeVertexParameters(theTol3D); - - theSlin.Remove(aNumOfLine2); - aNumOfLine2--; - hasBeenRemoved = Standard_True; - - continue; - } - } - } - - if(hasBeenRemoved) - aNumOfLine1--; - } -} - //====================================================================== // function: SequenceOfLine //====================================================================== @@ -1293,18 +911,9 @@ void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)& theS1, ListOfPnts.Clear(); if(isGeomInt) { - if(theD1->DomainIsInfinite() || theD2->DomainIsInfinite()) - { - GeomGeomPerfom( theS1, theD1, theS2, theD2, TolArc, - TolTang, ListOfPnts, RestrictLine, - typs1, typs2, theIsReqToKeepRLine); - } - else - { - GeomGeomPerfomTrimSurf( theS1, theD1, theS2, theD2, - TolArc, TolTang, ListOfPnts, RestrictLine, - typs1, typs2, theIsReqToKeepRLine); - } + GeomGeomPerfom( theS1, theD1, theS2, theD2, TolArc, + TolTang, ListOfPnts, RestrictLine, + typs1, typs2, theIsReqToKeepRLine); } else { @@ -1537,16 +1146,8 @@ void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)& theS1, } else if(ts1 == 1) { - if(theD1->DomainIsInfinite() || theD2->DomainIsInfinite()) - { - GeomGeomPerfom(theS1, theD1, theS2, theD2, TolArc, - TolTang, ListOfPnts, RestrictLine, typs1, typs2); - } - else - { - GeomGeomPerfomTrimSurf(theS1, theD1, theS2, theD2, - TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2); - } + GeomGeomPerfom(theS1, theD1, theS2, theD2, TolArc, + TolTang, ListOfPnts, RestrictLine, typs1, typs2); } } @@ -1754,13 +1355,32 @@ void IntPatch_Intersection::GeomGeomPerfom(const Handle(Adaptor3d_HSurface)& the slin.Append(wlin); } else - slin.Append(interii.Line(i)); + { + slin.Append(line); + } } for (Standard_Integer i = 1; i <= interii.NbPnts(); i++) { spnt.Append(interii.Point(i)); } + + if ((typs1 == GeomAbs_Cylinder) && (typs2 == GeomAbs_Cylinder)) + { + IntPatch_WLineTool::JoinWLines(slin, spnt, TolTang, + theS1->IsUPeriodic() ? theS1->UPeriod() : 0.0, + theS2->IsUPeriodic() ? theS2->UPeriod() : 0.0, + theS1->IsVPeriodic() ? theS1->VPeriod() : 0.0, + theS2->IsVPeriodic() ? theS2->VPeriod() : 0.0, + theS1->FirstUParameter(), + theS1->LastUParameter(), + theS1->FirstVParameter(), + theS1->LastVParameter(), + theS2->FirstUParameter(), + theS2->LastUParameter(), + theS2->FirstVParameter(), + theS2->LastVParameter()); + } } } else @@ -1852,80 +1472,6 @@ void IntPatch_Intersection:: } } -//======================================================================= -//function : GeomGeomPerfomTrimSurf -//purpose : This function returns ready walking-line (which is not need -// in convertation) as an intersection line between two -// trimmed surfaces. -//======================================================================= -void IntPatch_Intersection:: - GeomGeomPerfomTrimSurf( const Handle(Adaptor3d_HSurface)& theS1, - const Handle(Adaptor3d_TopolTool)& theD1, - const Handle(Adaptor3d_HSurface)& theS2, - const Handle(Adaptor3d_TopolTool)& theD2, - const Standard_Real theTolArc, - const Standard_Real theTolTang, - IntSurf_ListOfPntOn2S& theListOfPnts, - const Standard_Boolean RestrictLine, - const GeomAbs_SurfaceType theTyps1, - const GeomAbs_SurfaceType theTyps2, - const Standard_Boolean theIsReqToKeepRLine) -{ - IntSurf_Quadric Quad1,Quad2; - - if((theTyps1 == GeomAbs_Cylinder) && (theTyps2 == GeomAbs_Cylinder)) - { - IntPatch_ImpImpIntersection anInt; - anInt.Perform(theS1, theD1, theS2, theD2, myTolArc, - myTolTang, Standard_True, theIsReqToKeepRLine); - - done = anInt.IsDone(); - - if(done) - { - empt = anInt.IsEmpty(); - if (!empt) - { - tgte = anInt.TangentFaces(); - if (tgte) - oppo = anInt.OppositeFaces(); - - const Standard_Integer aNbLin = anInt.NbLines(); - const Standard_Integer aNbPts = anInt.NbPnts(); - - for(Standard_Integer aLID = 1; aLID <= aNbLin; aLID++) - { - const Handle(IntPatch_Line)& aLine = anInt.Line(aLID); - slin.Append(aLine); - } - - for(Standard_Integer aPID = 1; aPID <= aNbPts; aPID++) - { - const IntPatch_Point& aPoint = anInt.Point(aPID); - spnt.Append(aPoint); - } - - JoinWLines( slin, spnt, theTolTang, - theS1->IsUPeriodic()? theS1->UPeriod() : 0.0, - theS2->IsUPeriodic()? theS2->UPeriod() : 0.0, - theS1->IsVPeriodic()? theS1->VPeriod() : 0.0, - theS2->IsVPeriodic()? theS2->VPeriod() : 0.0, - theS1->FirstUParameter(), theS1->LastUParameter(), - theS1->FirstVParameter(), theS1->LastVParameter(), - theS2->FirstUParameter(), theS2->LastUParameter(), - theS2->FirstVParameter(), theS2->LastVParameter()); - } - } - } - else - { - GeomGeomPerfom(theS1, theD1, theS2, theD2, - theTolArc, theTolTang, theListOfPnts, - RestrictLine, theTyps1, theTyps2, theIsReqToKeepRLine); - } -} - - void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_TopolTool)& D1, const Handle(Adaptor3d_HSurface)& S2, diff --git a/src/IntPatch/IntPatch_WLineTool.cxx b/src/IntPatch/IntPatch_WLineTool.cxx new file mode 100644 index 0000000000..25bf5aaf59 --- /dev/null +++ b/src/IntPatch/IntPatch_WLineTool.cxx @@ -0,0 +1,1688 @@ +// Copyright (c) 1999-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include +#include +#include +#include + +// It is pure empirical value. +const Standard_Real IntPatch_WLineTool::myMaxConcatAngle = M_PI/6; + +//Bit-mask is used for information about +//the operation made in +//IntPatch_WLineTool::ExtendTwoWlinesToEachOther() method. +enum +{ + IntPatchWT_EnAll = 0x00, + IntPatchWT_DisLastLast = 0x01, + IntPatchWT_DisLastFirst = 0x02, + IntPatchWT_DisFirstLast = 0x04, + IntPatchWT_DisFirstFirst = 0x08 +}; + +//======================================================================= +//function : MinMax +//purpose : Replaces theParMIN = MIN(theParMIN, theParMAX), +// theParMAX = MAX(theParMIN, theParMAX). +// +// Static subfunction in IsSeamOrBound. +//======================================================================= +static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX) +{ + if(theParMIN > theParMAX) + { + const Standard_Real aTmp = theParMAX; + theParMAX = theParMIN; + theParMIN = aTmp; + } +} + +//========================================================================= +// function : FillPointsHash +// purpose : Fill points hash by input data. +// Static subfunction in ComputePurgedWLine. +//========================================================================= +static void FillPointsHash(const Handle(IntPatch_WLine) &theWLine, + NCollection_Array1 &thePointsHash) +{ + // 1 - Delete point. + // 0 - Store point. + // -1 - Vertex point (not delete). + Standard_Integer i, v; + + for(i = 1; i <= theWLine->NbPnts(); i++) + thePointsHash.SetValue(i, 0); + + for(v = 1; v <= theWLine->NbVertex(); v++) + { + IntPatch_Point aVertex = theWLine->Vertex(v); + Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine(); + thePointsHash.SetValue(avertexindex, -1); + } +} + +//========================================================================= +// function : MakeNewWLine +// purpose : Makes new walking line according to the points hash +// Static subfunction in ComputePurgedWLine and DeleteOuter. +//========================================================================= +static Handle(IntPatch_WLine) MakeNewWLine(const Handle(IntPatch_WLine) &theWLine, + const NCollection_Array1 &thePointsHash) +{ + Standard_Integer i; + + Handle(IntSurf_LineOn2S) aPurgedLineOn2S = new IntSurf_LineOn2S(); + Handle(IntPatch_WLine) aLocalWLine = new IntPatch_WLine(aPurgedLineOn2S, Standard_False); + Standard_Integer anOldLineIdx = 1, aVertexIdx = 1; + for(i = 1; i <= thePointsHash.Upper(); i++) + { + if (thePointsHash(i) == 0) + { + // Store this point. + aPurgedLineOn2S->Add(theWLine->Point(i)); + anOldLineIdx++; + } + else if (thePointsHash(i) == -1) + { + // Add vertex. + IntPatch_Point aVertex = theWLine->Vertex(aVertexIdx++); + aVertex.SetParameter(anOldLineIdx++); + aLocalWLine->AddVertex(aVertex); + aPurgedLineOn2S->Add(theWLine->Point(i)); + } + } + + return aLocalWLine; +} + +//========================================================================= +// function : MovePoint +// purpose : Move point into surface param space. No interpolation used +// because walking algorithm should care for closeness to the param space. +// Static subfunction in ComputePurgedWLine. +//========================================================================= +static void MovePoint(const Handle(Adaptor3d_HSurface) &theS1, + Standard_Real &U1, Standard_Real &V1) +{ + if (U1 < theS1->FirstUParameter()) + U1 = theS1->FirstUParameter(); + + if (U1 > theS1->LastUParameter()) + U1 = theS1->LastUParameter(); + + if (V1 < theS1->FirstVParameter()) + V1 = theS1->FirstVParameter(); + + if (V1 > theS1->LastVParameter()) + V1 = theS1->LastVParameter(); +} + +//========================================================================= +// function : DeleteOuterPoints +// purpose : Check and delete out of bounds points on walking line. +// Static subfunction in ComputePurgedWLine. +//========================================================================= +static Handle(IntPatch_WLine) + DeleteOuterPoints(const Handle(IntPatch_WLine) &theWLine, + const Handle(Adaptor3d_HSurface) &theS1, + const Handle(Adaptor3d_HSurface) &theS2, + const Handle(Adaptor3d_TopolTool) &theDom1, + const Handle(Adaptor3d_TopolTool) &theDom2) +{ + Standard_Integer i; + + NCollection_Array1 aDelOuterPointsHash(1, theWLine->NbPnts()); + FillPointsHash(theWLine, aDelOuterPointsHash); + + if (theS1->IsUPeriodic() || theS1->IsVPeriodic() || + theS2->IsUPeriodic() || theS2->IsVPeriodic() ) + return theWLine; + + gp_Pnt2d aPntOnF1, aPntOnF2; + Standard_Real aX1, aY1, aX2, aY2; + + // Iterate over points in walking line and delete which are out of bounds. + // Forward. + Standard_Boolean isAllDeleted = Standard_True; + Standard_Boolean aChangedFirst = Standard_False; + Standard_Integer aFirstGeomIdx = 1; + for(i = 1; i <= theWLine->NbPnts(); i++) + { + theWLine->Point(i).Parameters(aX1, aY1, aX2, aY2); + aPntOnF1.SetCoord(aX1, aY1); + aPntOnF2.SetCoord(aX2, aY2); + + TopAbs_State aState1 = theDom1->Classify(aPntOnF1, Precision::Confusion()); + TopAbs_State aState2 = theDom2->Classify(aPntOnF2, Precision::Confusion()); + + if (aState1 == TopAbs_OUT || + aState2 == TopAbs_OUT ) + { + aDelOuterPointsHash(i) = 1; + aChangedFirst = Standard_True; + } + else + { + isAllDeleted = Standard_False; + + aFirstGeomIdx = Max (i - 1, 1); + if (aDelOuterPointsHash(i) == -1) + aFirstGeomIdx = i; // Use data what lies in (i) point / vertex. + + aDelOuterPointsHash(i) = -1; + break; + } + } + + if (isAllDeleted) + { + // ALL points are out of bounds: + // case boolean bcut_complex F5 and similar. + return theWLine; + } + + // Backward. + Standard_Boolean aChangedLast = Standard_False; + Standard_Integer aLastGeomIdx = theWLine->NbPnts(); + for(i = theWLine->NbPnts(); i >= 1; i--) + { + theWLine->Point(i).Parameters(aX1, aY1, aX2, aY2); + aPntOnF1.SetCoord(aX1, aY1); + aPntOnF2.SetCoord(aX2, aY2); + + TopAbs_State aState1 = theDom1->Classify(aPntOnF1, Precision::Confusion()); + TopAbs_State aState2 = theDom2->Classify(aPntOnF2, Precision::Confusion()); + + if (aState1 == TopAbs_OUT || + aState2 == TopAbs_OUT ) + { + aDelOuterPointsHash(i) = 1; + aChangedLast = Standard_True; // Move vertex to first good point + } + else + { + aLastGeomIdx = Min (i + 1, theWLine->NbPnts()); + if (aDelOuterPointsHash(i) == -1) + aLastGeomIdx = i; // Use data what lies in (i) point / vertex. + + aDelOuterPointsHash(i) = -1; + break; + } + } + + if (!aChangedFirst && !aChangedLast) + { + // Nothing is done, return input. + return theWLine; + } + + // Build new line and modify geometry of necessary vertexes. + Handle(IntPatch_WLine) aLocalWLine = MakeNewWLine(theWLine, aDelOuterPointsHash); + + if (aChangedFirst) + { + // Vertex geometry. + IntPatch_Point aVertex = aLocalWLine->Vertex(1); + aVertex.SetValue(theWLine->Point(aFirstGeomIdx).Value()); + Standard_Real aU1, aU2, aV1, aV2; + theWLine->Point(aFirstGeomIdx).Parameters(aU1, aV1, aU2, aV2); + MovePoint(theS1, aU1, aV1); + MovePoint(theS2, aU2, aV2); + aVertex.SetParameters(aU1, aV1, aU2, aV2); + aLocalWLine->Replace(1, aVertex); + // Change point in walking line. + aLocalWLine->SetPoint(1, aVertex); + } + + if (aChangedLast) + { + // Vertex geometry. + IntPatch_Point aVertex = aLocalWLine->Vertex(aLocalWLine->NbVertex()); + aVertex.SetValue(theWLine->Point(aLastGeomIdx).Value()); + Standard_Real aU1, aU2, aV1, aV2; + theWLine->Point(aLastGeomIdx).Parameters(aU1, aV1, aU2, aV2); + MovePoint(theS1, aU1, aV1); + MovePoint(theS2, aU2, aV2); + aVertex.SetParameters(aU1, aV1, aU2, aV2); + aLocalWLine->Replace(aLocalWLine->NbVertex(), aVertex); + // Change point in walking line. + aLocalWLine->SetPoint(aLocalWLine->NbPnts(), aVertex); + } + + + return aLocalWLine; +} + +//========================================================================= +// function : IsInsideIn2d +// purpose : Check if aNextPnt lies inside of tube build on aBasePnt and aBaseVec. +// In 2d space. Static subfunction in DeleteByTube. +//========================================================================= +static Standard_Boolean IsInsideIn2d(const gp_Pnt2d& aBasePnt, + const gp_Vec2d& aBaseVec, + const gp_Pnt2d& aNextPnt, + const Standard_Real aSquareMaxDist) +{ + gp_Vec2d aVec2d(aBasePnt, aNextPnt); + + //d*d = (basevec^(nextpnt-basepnt))**2 / basevec**2 + Standard_Real aCross = aVec2d.Crossed(aBaseVec); + Standard_Real aSquareDist = aCross * aCross + / aBaseVec.SquareMagnitude(); + + return (aSquareDist <= aSquareMaxDist); +} + +//========================================================================= +// function : IsInsideIn3d +// purpose : Check if aNextPnt lies inside of tube build on aBasePnt and aBaseVec. +// In 3d space. Static subfunction in DeleteByTube. +//========================================================================= +static Standard_Boolean IsInsideIn3d(const gp_Pnt& aBasePnt, + const gp_Vec& aBaseVec, + const gp_Pnt& aNextPnt, + const Standard_Real aSquareMaxDist) +{ + gp_Vec aVec(aBasePnt, aNextPnt); + + //d*d = (basevec^(nextpnt-basepnt))**2 / basevec**2 + Standard_Real aSquareDist = aVec.CrossSquareMagnitude(aBaseVec) + / aBaseVec.SquareMagnitude(); + + return (aSquareDist <= aSquareMaxDist); +} + +static const Standard_Integer aMinNbBadDistr = 15; +static const Standard_Integer aNbSingleBezier = 30; + +//========================================================================= +// function : DeleteByTube +// purpose : Check and delete points using tube criteria. +// Static subfunction in ComputePurgedWLine. +//========================================================================= +static Handle(IntPatch_WLine) + DeleteByTube(const Handle(IntPatch_WLine) &theWLine, + const Handle(Adaptor3d_HSurface) &theS1, + const Handle(Adaptor3d_HSurface) &theS2) +{ + // III: Check points for tube criteria: + // Workaround to handle case of small amount points after purge. + // Test "boolean boptuc_complex B5" and similar. + Standard_Integer aNbPnt = 0 , i; + + if (theWLine->NbPnts() <= 2) + return theWLine; + + NCollection_Array1 aNewPointsHash(1, theWLine->NbPnts()); + FillPointsHash(theWLine, aNewPointsHash); + + // Inital computations. + Standard_Real UonS1[3], VonS1[3], UonS2[3], VonS2[3]; + theWLine->Point(1).ParametersOnS1(UonS1[0], VonS1[0]); + theWLine->Point(2).ParametersOnS1(UonS1[1], VonS1[1]); + theWLine->Point(1).ParametersOnS2(UonS2[0], VonS2[0]); + theWLine->Point(2).ParametersOnS2(UonS2[1], VonS2[1]); + + gp_Pnt2d aBase2dPnt1(UonS1[0], VonS1[0]); + gp_Pnt2d aBase2dPnt2(UonS2[0], VonS2[0]); + gp_Vec2d aBase2dVec1(UonS1[1] - UonS1[0], VonS1[1] - VonS1[0]); + gp_Vec2d aBase2dVec2(UonS2[1] - UonS2[0], VonS2[1] - VonS2[0]); + gp_Pnt aBase3dPnt = theWLine->Point(1).Value(); + gp_Vec aBase3dVec(theWLine->Point(1).Value(), theWLine->Point(2).Value()); + + // Choose base tolerance and scale it to pipe algorithm. + const Standard_Real aBaseTolerance = Precision::Approximation(); + Standard_Real aResS1Tol = Min(theS1->UResolution(aBaseTolerance), + theS1->VResolution(aBaseTolerance)); + Standard_Real aResS2Tol = Min(theS2->UResolution(aBaseTolerance), + theS2->VResolution(aBaseTolerance)); + Standard_Real aTol1 = aResS1Tol * aResS1Tol; + Standard_Real aTol2 = aResS2Tol * aResS2Tol; + Standard_Real aTol3d = aBaseTolerance * aBaseTolerance; + + const Standard_Real aLimitCoeff = 0.99 * 0.99; + for(i = 3; i <= theWLine->NbPnts(); i++) + { + Standard_Boolean isDeleteState = Standard_False; + + theWLine->Point(i).ParametersOnS1(UonS1[2], VonS1[2]); + theWLine->Point(i).ParametersOnS2(UonS2[2], VonS2[2]); + gp_Pnt2d aPnt2dOnS1(UonS1[2], VonS1[2]); + gp_Pnt2d aPnt2dOnS2(UonS2[2], VonS2[2]); + const gp_Pnt& aPnt3d = theWLine->Point(i).Value(); + + if (aNewPointsHash(i - 1) != - 1 && + IsInsideIn2d(aBase2dPnt1, aBase2dVec1, aPnt2dOnS1, aTol1) && + IsInsideIn2d(aBase2dPnt2, aBase2dVec2, aPnt2dOnS2, aTol2) && + IsInsideIn3d(aBase3dPnt, aBase3dVec, aPnt3d, aTol3d) ) + { + // Handle possible uneven parametrization on one of 2d subspaces. + // Delete point only when expected lengths are close to each other (aLimitCoeff). + // Example: + // c2d1 - line + // c3d - line + // c2d2 - geometrically line, but have uneven parametrization -> c2d2 is bspline. + gp_XY aPntOnS1[2]= { gp_XY(UonS1[1] - UonS1[0], VonS1[1] - VonS1[0]) + , gp_XY(UonS1[2] - UonS1[1], VonS1[2] - VonS1[1])}; + gp_XY aPntOnS2[2]= { gp_XY(UonS2[1] - UonS2[0], VonS2[1] - VonS2[0]) + , gp_XY(UonS2[2] - UonS2[1], VonS2[2] - VonS2[1])}; + + Standard_Real aStepOnS1 = aPntOnS1[0].SquareModulus() / aPntOnS1[1].SquareModulus(); + Standard_Real aStepOnS2 = aPntOnS2[0].SquareModulus() / aPntOnS2[1].SquareModulus(); + + // Check very rare case when wline fluctuates nearly one point and some of them may be equal. + // Middle point will be deleted when such situation occurs. + // bugs moddata_2 bug469. + if (Min(aStepOnS1, aStepOnS2) >= aLimitCoeff * Max(aStepOnS1, aStepOnS2)) + { + // Set hash flag to "Delete" state. + isDeleteState = Standard_True; + aNewPointsHash.SetValue(i - 1, 1); + + // Change middle point. + UonS1[1] = UonS1[2]; + UonS2[1] = UonS2[2]; + VonS1[1] = VonS1[2]; + VonS2[1] = VonS2[2]; + } + } + + if (!isDeleteState) + { + // Compute new pipe parameters. + UonS1[0] = UonS1[1]; + VonS1[0] = VonS1[1]; + UonS2[0] = UonS2[1]; + VonS2[0] = VonS2[1]; + + UonS1[1] = UonS1[2]; + VonS1[1] = VonS1[2]; + UonS2[1] = UonS2[2]; + VonS2[1] = VonS2[2]; + + aBase2dPnt1.SetCoord(UonS1[0], VonS1[0]); + aBase2dPnt2.SetCoord(UonS2[0], VonS2[0]); + aBase2dVec1.SetCoord(UonS1[1] - UonS1[0], VonS1[1] - VonS1[0]); + aBase2dVec2.SetCoord(UonS2[1] - UonS2[0], VonS2[1] - VonS2[0]); + aBase3dPnt = theWLine->Point(i - 1).Value(); + aBase3dVec = gp_Vec(theWLine->Point(i - 1).Value(), theWLine->Point(i).Value()); + + aNbPnt++; + } + } + + // Workaround to handle case of small amount of points after purge. + // Test "boolean boptuc_complex B5" and similar. + // This is possible since there are at least two points. + if (aNewPointsHash(1) == -1 && + aNewPointsHash(2) == -1 && + aNbPnt <= 3) + { + // Delete first. + aNewPointsHash(1) = 1; + } + if (aNewPointsHash(theWLine->NbPnts() - 1) == -1 && + aNewPointsHash(theWLine->NbPnts() ) == -1 && + aNbPnt <= 3) + { + // Delete last. + aNewPointsHash(theWLine->NbPnts()) = 1; + } + + // Purgre when too small amount of points left. + if (aNbPnt <= 2) + { + for(i = aNewPointsHash.Lower(); i <= aNewPointsHash.Upper(); i++) + { + if (aNewPointsHash(i) != -1) + { + aNewPointsHash(i) = 1; + } + } + } + + // Handle possible bad distribution of points, + // which are will converted into one single bezier curve (less than 30 points). + // Make distribution more even: + // max step will be nearly to 0.1 of param distance. + if (aNbPnt + 2 > aMinNbBadDistr && + aNbPnt + 2 < aNbSingleBezier ) + { + for(Standard_Integer anIdx = 1; anIdx <= 8; anIdx++) + { + Standard_Integer aHashIdx = + Standard_Integer(anIdx * theWLine->NbPnts() / 9); + + //Vertex must be stored as VERTEX (HASH = -1) + if (aNewPointsHash(aHashIdx) != -1) + aNewPointsHash(aHashIdx) = 0; + } + } + + return MakeNewWLine(theWLine, aNewPointsHash); +} + +//======================================================================= +//function : IsOnPeriod +//purpose : Checks, if [theU1, theU2] intersects the period-value +// (k*thePeriod, where k is an integer number (k = 0, +/-1, +/-2, ...). +// +// Returns: +// 0 - if interval [theU1, theU2] does not intersect the "period-value" +// or if thePeriod == 0.0; +// 1 - if interval (theU1, theU2) intersect the "period-value". +// 2 - if theU1 or/and theU2 lie ON the "period-value" +// +//ATTENTION!!! +// If (theU1 == theU2) then this function will return only both 0 or 2. +//======================================================================= +static Standard_Integer IsOnPeriod(const Standard_Real theU1, + const Standard_Real theU2, + const Standard_Real thePeriod) +{ + if(thePeriod < RealSmall()) + return 0; + + //If interval [theU1, theU2] intersect seam-edge then there exists an integer + //number N such as + // (theU1 <= T*N <= theU2) <=> (theU1/T <= N <= theU2/T), + //where T is the period. + //I.e. the inerval [theU1/T, theU2/T] must contain at least one + //integer number. In this case, Floor(theU1/T) and Floor(theU2/T) + //return different values or theU1/T is strictly integer number. + //Examples: + // 1. theU1/T==2.8, theU2/T==3.5 => Floor(theU1/T) == 2, Floor(theU2/T) == 3. + // 2. theU1/T==2.0, theU2/T==2.6 => Floor(theU1/T) == Floor(theU2/T) == 2. + + const Standard_Real aVal1 = theU1/thePeriod, + aVal2 = theU2/thePeriod; + const Standard_Integer aPar1 = static_cast(Floor(aVal1)); + const Standard_Integer aPar2 = static_cast(Floor(aVal2)); + if(aPar1 != aPar2) + {//Interval (theU1, theU2] intersects seam-edge + if(IsEqual(aVal2, static_cast(aPar2))) + {//aVal2 is an integer number => theU2 lies ON the "seam-edge" + return 2; + } + + return 1; + } + + //Here, aPar1 == aPar2. + + if(IsEqual(aVal1, static_cast(aPar1))) + {//aVal1 is an integer number => theU1 lies ON the "seam-edge" + return 2; + } + + //If aVal2 is a true integer number then always (aPar1 != aPar2). + + return 0; +} + +//======================================================================= +//function : IsSeamOrBound +//purpose : Returns TRUE if segment [thePtf, thePtl] intersects "seam-edge" +// (if it exist) or surface boundaries and both thePtf and thePtl do +// not match "seam-edge" or boundaries. +// Point thePtmid lies in this segment (in both 3D and 2D-space). +// If thePtmid match "seam-edge" or boundaries strictly +// (without any tolerance) then the function will return TRUE. +// See comments in function body for detail information. +//======================================================================= +static Standard_Boolean IsSeamOrBound(const IntSurf_PntOn2S& thePtf, + const IntSurf_PntOn2S& thePtl, + const IntSurf_PntOn2S& thePtmid, + const Standard_Real theU1Period, + const Standard_Real theU2Period, + const Standard_Real theV1Period, + const Standard_Real theV2Period, + const Standard_Real theUfSurf1, + const Standard_Real theUlSurf1, + const Standard_Real theVfSurf1, + const Standard_Real theVlSurf1, + const Standard_Real theUfSurf2, + const Standard_Real theUlSurf2, + const Standard_Real theVfSurf2, + const Standard_Real theVlSurf2) +{ + Standard_Real aU11 = 0.0, aU12 = 0.0, aV11 = 0.0, aV12 = 0.0; + Standard_Real aU21 = 0.0, aU22 = 0.0, aV21 = 0.0, aV22 = 0.0; + thePtf.Parameters(aU11, aV11, aU12, aV12); + thePtl.Parameters(aU21, aV21, aU22, aV22); + + MinMax(aU11, aU21); + MinMax(aV11, aV21); + MinMax(aU12, aU22); + MinMax(aV12, aV22); + + if((aU11 - theUfSurf1)*(aU21 - theUfSurf1) < 0.0) + {//Interval [aU11, aU21] intersects theUfSurf1 + return Standard_True; + } + + if((aU11 - theUlSurf1)*(aU21 - theUlSurf1) < 0.0) + {//Interval [aU11, aU21] intersects theUlSurf1 + return Standard_True; + } + + if((aV11 - theVfSurf1)*(aV21 - theVfSurf1) < 0.0) + {//Interval [aV11, aV21] intersects theVfSurf1 + return Standard_True; + } + + if((aV11 - theVlSurf1)*(aV21 - theVlSurf1) < 0.0) + {//Interval [aV11, aV21] intersects theVlSurf1 + return Standard_True; + } + + if((aU12 - theUfSurf2)*(aU22 - theUfSurf2) < 0.0) + {//Interval [aU12, aU22] intersects theUfSurf2 + return Standard_True; + } + + if((aU12 - theUlSurf2)*(aU22 - theUlSurf2) < 0.0) + {//Interval [aU12, aU22] intersects theUlSurf2 + return Standard_True; + } + + if((aV12 - theVfSurf2)*(aV22 - theVfSurf2) < 0.0) + {//Interval [aV12, aV22] intersects theVfSurf2 + return Standard_True; + } + + if((aV12 - theVlSurf2)*(aV22 - theVlSurf2) < 0.0) + {//Interval [aV12, aV22] intersects theVlSurf2 + return Standard_True; + } + + if(IsOnPeriod(aU11, aU21, theU1Period)) + return Standard_True; + + if(IsOnPeriod(aV11, aV21, theV1Period)) + return Standard_True; + + if(IsOnPeriod(aU12, aU22, theU2Period)) + return Standard_True; + + if(IsOnPeriod(aV12, aV22, theV2Period)) + return Standard_True; + + /* + The segment [thePtf, thePtl] does not intersect the boundaries and + the seam-edge of the surfaces. + Nevertheless, following situation is possible: + + seam or + bound + | + thePtf * | + | + * thePtmid + thePtl * | + | + + This case must be processed, too. + */ + + Standard_Real aU1 = 0.0, aU2 = 0.0, aV1 = 0.0, aV2 = 0.0; + thePtmid.Parameters(aU1, aV1, aU2, aV2); + + if(IsEqual(aU1, theUfSurf1) || IsEqual(aU1, theUlSurf1)) + return Standard_True; + + if(IsEqual(aU2, theUfSurf2) || IsEqual(aU2, theUlSurf2)) + return Standard_True; + + if(IsEqual(aV1, theVfSurf1) || IsEqual(aV1, theVlSurf1)) + return Standard_True; + + if(IsEqual(aV2, theVfSurf2) || IsEqual(aV2, theVlSurf2)) + return Standard_True; + + if(IsOnPeriod(aU1, aU1, theU1Period)) + return Standard_True; + + if(IsOnPeriod(aU2, aU2, theU2Period)) + return Standard_True; + + if(IsOnPeriod(aV1, aV1, theV1Period)) + return Standard_True; + + if(IsOnPeriod(aV2, aV2, theV2Period)) + return Standard_True; + + return Standard_False; +} + +#if 0 +//======================================================================= +//function : AbjustPeriodicToPrevPoint +//purpose : Returns theCurrentParam in order to the distance betwen +// theRefParam and theCurrentParam is less than 0.5*thePeriod. +//======================================================================= +static void AbjustPeriodicToPrevPoint(const Standard_Real theRefParam, + const Standard_Real thePeriod, + Standard_Real& theCurrentParam) +{ + if(thePeriod == 0.0) + return; + + Standard_Real aDeltaPar = 2.0*(theRefParam - theCurrentParam); + const Standard_Real anIncr = Sign(thePeriod, aDeltaPar); + while(Abs(aDeltaPar) > thePeriod) + { + theCurrentParam += anIncr; + aDeltaPar = 2.0*(theRefParam-theCurrentParam); + } +} + +//======================================================================= +//function : IsIntersectionPoint +//purpose : Returns True if thePmid is intersection point +// between theS1 and theS2 with given tolerance. +// In this case, parameters of thePmid on every quadric +// will be recomputed and returned. +//======================================================================= +static Standard_Boolean IsIntersectionPoint(const gp_Pnt& thePmid, + const IntSurf_Quadric& theS1, + const IntSurf_Quadric& theS2, + const IntSurf_PntOn2S& theRefPt, + const Standard_Real theTol, + const Standard_Real theU1Period, + const Standard_Real theU2Period, + const Standard_Real theV1Period, + const Standard_Real theV2Period, + Standard_Real &theU1, + Standard_Real &theV1, + Standard_Real &theU2, + Standard_Real &theV2) +{ + Standard_Real aU1Ref = 0.0, aV1Ref = 0.0, aU2Ref = 0.0, aV2Ref = 0.0; + theRefPt.Parameters(aU1Ref, aV1Ref, aU2Ref, aV2Ref); + theS1.Parameters(thePmid, theU1, theV1); + theS2.Parameters(thePmid, theU2, theV2); + + AbjustPeriodicToPrevPoint(aU1Ref, theU1Period, theU1); + AbjustPeriodicToPrevPoint(aV1Ref, theV1Period, theV1); + AbjustPeriodicToPrevPoint(aU2Ref, theU2Period, theU2); + AbjustPeriodicToPrevPoint(aV2Ref, theV2Period, theV2); + + const gp_Pnt aP1(theS1.Value(theU1, theV1)); + const gp_Pnt aP2(theS2.Value(theU2, theV2)); + + return (aP1.SquareDistance(aP2) <= theTol*theTol); +} + +//======================================================================= +//function : ExtendFirst +//purpose : Adds thePOn2S to the begin of theWline +//======================================================================= +static void ExtendFirst(const Handle(IntPatch_WLine)& theWline, + const IntSurf_PntOn2S& theAddedPt) +{ + Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0; + theAddedPt.Parameters(aU1, aV1, aU2, aV2); + + if(theAddedPt.IsSame(theWline->Point(1), Precision::Confusion())) + { + theWline->Curve()->Value(1, theAddedPt); + for(Standard_Integer i = 1; i <= theWline->NbVertex(); i++) + { + IntPatch_Point &aVert = theWline->ChangeVertex(i); + if(aVert.ParameterOnLine() != 1) + break; + + aVert.SetParameters(aU1, aV1, aU2, aV2); + aVert.SetValue(theAddedPt.Value()); + } + + return; + } + + theWline->Curve()->InsertBefore(1, theAddedPt); + + for(Standard_Integer i = 1; i <= theWline->NbVertex(); i++) + { + IntPatch_Point &aVert = theWline->ChangeVertex(i); + + if(aVert.ParameterOnLine() == 1) + { + aVert.SetParameters(aU1, aV1, aU2, aV2); + aVert.SetValue(theAddedPt.Value()); + } + else + { + aVert.SetParameter(aVert.ParameterOnLine()+1); + } + } +} + +//======================================================================= +//function : ExtendLast +//purpose : Adds thePOn2S to the end of theWline +//======================================================================= +static void ExtendLast(const Handle(IntPatch_WLine)& theWline, + const IntSurf_PntOn2S& theAddedPt) +{ + Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0; + theAddedPt.Parameters(aU1, aV1, aU2, aV2); + + const Standard_Integer aNbPnts = theWline->NbPnts(); + if(theAddedPt.IsSame(theWline->Point(aNbPnts), Precision::Confusion())) + { + theWline->Curve()->Value(aNbPnts, theAddedPt); + } + else + { + theWline->Curve()->Add(theAddedPt); + } + + for(Standard_Integer i = theWline->NbVertex(); i >= 1; i--) + { + IntPatch_Point &aVert = theWline->ChangeVertex(i); + if(aVert.ParameterOnLine() != aNbPnts) + break; + + aVert.SetParameters(aU1, aV1, aU2, aV2); + aVert.SetValue(theAddedPt.Value()); + aVert.SetParameter(theWline->NbPnts()); + } +} + +//========================================================================= +// function: IsOutOfDomain +// purpose : Checks, if 2D-representation of thePOn2S is in surfaces domain, +// defined by bounding-boxes theBoxS1 and theBoxS2 +//========================================================================= +static Standard_Boolean IsOutOfDomain(const Bnd_Box2d& theBoxS1, + const Bnd_Box2d& theBoxS2, + const IntSurf_PntOn2S &thePOn2S, + const Standard_Real theU1Period, + const Standard_Real theU2Period, + const Standard_Real theV1Period, + const Standard_Real theV2Period) +{ + Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0; + Standard_Real aU1min = 0.0, aU1max = 0.0, aV1min = 0.0, aV1max = 0.0; + Standard_Real aU2min = 0.0, aU2max = 0.0, aV2min = 0.0, aV2max = 0.0; + + thePOn2S.Parameters(aU1, aV1, aU2, aV2); + + theBoxS1.Get(aU1min, aV1min, aU1max, aV1max); + theBoxS2.Get(aU2min, aV2min, aU2max, aV2max); + + aU1 = ElCLib::InPeriod(aU1, aU1min, aU1min + theU1Period); + aV1 = ElCLib::InPeriod(aV1, aV1min, aV1min + theV1Period); + aU2 = ElCLib::InPeriod(aU2, aU2min, aU2min + theU2Period); + aV2 = ElCLib::InPeriod(aV2, aV2min, aV2min + theV2Period); + + return (theBoxS1.IsOut(gp_Pnt2d(aU1, aV1)) || + theBoxS2.IsOut(gp_Pnt2d(aU2, aV2))); +} + +//======================================================================= +//function : CheckArgumentsToExtend +//purpose : Check if extending is possible +// (see IntPatch_WLineTool::ExtendTwoWlinesToEachOther) +//======================================================================= +Standard_Boolean CheckArgumentsToExtend(const IntSurf_Quadric& theS1, + const IntSurf_Quadric& theS2, + const IntSurf_PntOn2S& thePtWL1, + const IntSurf_PntOn2S& thePtWL2, + IntSurf_PntOn2S& theNewPoint, + const gp_Vec& theVec1, + const gp_Vec& theVec2, + const gp_Vec& theVec3, + const Bnd_Box2d& theBoxS1, + const Bnd_Box2d& theBoxS2, + const Standard_Real theToler3D, + const Standard_Real theU1Period, + const Standard_Real theU2Period, + const Standard_Real theV1Period, + const Standard_Real theV2Period) +{ + const Standard_Real aSqToler = theToler3D*theToler3D; + + if(theVec3.SquareMagnitude() <= aSqToler) + { + return Standard_False; + } + + if((theVec1.Angle(theVec2) > IntPatch_WLineTool::myMaxConcatAngle) || + (theVec1.Angle(theVec3) > IntPatch_WLineTool::myMaxConcatAngle) || + (theVec2.Angle(theVec3) > IntPatch_WLineTool::myMaxConcatAngle)) + { + return Standard_False; + } + + const gp_Pnt aPmid(0.5*(thePtWL1.Value().XYZ()+thePtWL2.Value().XYZ())); + + Standard_Real aU1=0.0, aV1=0.0, aU2=0.0, aV2=0.0; + + theBoxS1.Get(aU1, aV1, aU2, aV2); + const Standard_Real aU1f = aU1, aV1f = aV1; + theBoxS2.Get(aU1, aV1, aU2, aV2); + const Standard_Real aU2f = aU1, aV2f = aV1; + + if(!IsIntersectionPoint(aPmid, theS1, theS2, thePtWL1, theToler3D, + theU1Period, theU2Period, theV1Period, theV2Period, + aU1, aV1, aU2, aV2)) + { + return Standard_False; + } + + theNewPoint.SetValue(aPmid, aU1, aV1, aU2, aV2); + + if(IsOutOfDomain(theBoxS1, theBoxS2, theNewPoint, + theU1Period, theU2Period, + theV1Period, theV2Period)) + { + return Standard_False; + } + + Standard_Real aU11 = 0.0, aV11 = 0.0, aU21 = 0.0, aV21 = 0.0, + aU12 = 0.0, aV12 = 0.0, aU22 = 0.0, aV22 = 0.0; + + thePtWL1.Parameters(aU11, aV11, aU21, aV21); + thePtWL2.Parameters(aU12, aV12, aU22, aV22); + + if(IsOnPeriod(aU11 - aU1f, aU12 - aU1f, theU1Period) || + IsOnPeriod(aV11 - aV1f, aV12 - aV1f, theV1Period) || + IsOnPeriod(aU21 - aU2f, aU22 - aU2f, theU2Period) || + IsOnPeriod(aV21 - aV2f, aV22 - aV2f, theV2Period)) + { + return Standard_False; + } + + return Standard_True; +} +#endif +//======================================================================= +//function : CheckArgumentsToJoin +//purpose : Check if joining is possible +// (see IntPatch_WLineTool::JoinWLines) +//======================================================================= +Standard_Boolean CheckArgumentsToJoin(const gp_Vec& theVec1, + const gp_Vec& theVec2) +{ + // [0, PI] - range + const Standard_Real anAngle = theVec1.Angle(theVec2); + + return (anAngle < IntPatch_WLineTool::myMaxConcatAngle); +} + +#if 0 +//======================================================================= +//function : ExtendTwoWLFirstFirst +//purpose : Performs extending theWLine1 and theWLine2 through their +// respecting start point. +//======================================================================= +static void ExtendTwoWLFirstFirst(const IntSurf_Quadric& theS1, + const IntSurf_Quadric& theS2, + const Handle(IntPatch_WLine)& theWLine1, + const Handle(IntPatch_WLine)& theWLine2, + const IntSurf_PntOn2S& thePtWL1, + const IntSurf_PntOn2S& thePtWL2, + const gp_Vec& theVec1, + const gp_Vec& theVec2, + const gp_Vec& theVec3, + const Bnd_Box2d& theBoxS1, + const Bnd_Box2d& theBoxS2, + const Standard_Real theToler3D, + const Standard_Real theU1Period, + const Standard_Real theU2Period, + const Standard_Real theV1Period, + const Standard_Real theV2Period, + unsigned int &theCheckResult, + Standard_Boolean &theHasBeenJoined) +{ + IntSurf_PntOn2S aPOn2S; + if(!CheckArgumentsToExtend(theS1, theS2, thePtWL1, thePtWL2, aPOn2S, + theVec1, theVec2, theVec3, + theBoxS1, theBoxS2, theToler3D, + theU1Period, theU2Period, theV1Period, theV2Period)) + { + return; + } + + theCheckResult |= (IntPatchWT_DisFirstLast | IntPatchWT_DisLastFirst); + + ExtendFirst(theWLine1, aPOn2S); + ExtendFirst(theWLine2, aPOn2S); + + if(theHasBeenJoined) + return; + + Standard_Real aPrm = theWLine1->Vertex(1).ParameterOnLine(); + while(theWLine1->Vertex(1).ParameterOnLine() == aPrm) + theWLine1->RemoveVertex(1); + + aPrm = theWLine2->Vertex(1).ParameterOnLine(); + while(theWLine2->Vertex(1).ParameterOnLine() == aPrm) + theWLine2->RemoveVertex(1); + + const Standard_Integer aNbPts = theWLine2->NbPnts(); + for(Standard_Integer aNPt = 2; aNPt <= aNbPts; aNPt++) + { + const IntSurf_PntOn2S& aPt = theWLine2->Point(aNPt); + theWLine1->Curve()->InsertBefore(1, aPt); + } + + for(Standard_Integer aNVtx = 1; aNVtx <= theWLine1->NbVertex(); aNVtx++) + { + IntPatch_Point &aVert = theWLine1->ChangeVertex(aNVtx); + const Standard_Real aCurParam = aVert.ParameterOnLine(); + aVert.SetParameter(aNbPts+aCurParam-1); + } + + for(Standard_Integer aNVtx = 1; aNVtx <= theWLine2->NbVertex(); aNVtx++) + { + IntPatch_Point &aVert = theWLine2->ChangeVertex(aNVtx); + const Standard_Real aCurParam = aVert.ParameterOnLine(); + aVert.SetParameter(aNbPts-aCurParam+1); + theWLine1->AddVertex(aVert, Standard_True); + } + + theHasBeenJoined = Standard_True; +} + +//======================================================================= +//function : ExtendTwoWLFirstLast +//purpose : Performs extending theWLine1 through its start point and theWLine2 +// through its end point. +//======================================================================= +static void ExtendTwoWLFirstLast(const IntSurf_Quadric& theS1, + const IntSurf_Quadric& theS2, + const Handle(IntPatch_WLine)& theWLine1, + const Handle(IntPatch_WLine)& theWLine2, + const IntSurf_PntOn2S& thePtWL1, + const IntSurf_PntOn2S& thePtWL2, + const gp_Vec& theVec1, + const gp_Vec& theVec2, + const gp_Vec& theVec3, + const Bnd_Box2d& theBoxS1, + const Bnd_Box2d& theBoxS2, + const Standard_Real theToler3D, + const Standard_Real theU1Period, + const Standard_Real theU2Period, + const Standard_Real theV1Period, + const Standard_Real theV2Period, + unsigned int &theCheckResult, + Standard_Boolean &theHasBeenJoined) +{ + IntSurf_PntOn2S aPOn2S; + if(!CheckArgumentsToExtend(theS1, theS2, thePtWL1, thePtWL2, aPOn2S, + theVec1, theVec2, theVec3, + theBoxS1, theBoxS2, theToler3D, + theU1Period, theU2Period, theV1Period, theV2Period)) + { + return; + } + + theCheckResult |= IntPatchWT_DisLastLast; + + ExtendFirst(theWLine1, aPOn2S); + ExtendLast (theWLine2, aPOn2S); + + if(theHasBeenJoined) + return; + + Standard_Real aPrm = theWLine1->Vertex(1).ParameterOnLine(); + while(theWLine1->Vertex(1).ParameterOnLine() == aPrm) + theWLine1->RemoveVertex(1); + + aPrm = theWLine2->Vertex(theWLine2->NbVertex()).ParameterOnLine(); + while(theWLine2->Vertex(theWLine2->NbVertex()).ParameterOnLine() == aPrm) + theWLine2->RemoveVertex(theWLine2->NbVertex()); + + const Standard_Integer aNbPts = theWLine2->NbPnts(); + for(Standard_Integer aNPt = aNbPts - 1; aNPt >= 1; aNPt--) + { + const IntSurf_PntOn2S& aPt = theWLine2->Point(aNPt); + theWLine1->Curve()->InsertBefore(1, aPt); + } + + for(Standard_Integer aNVtx = 1; aNVtx <= theWLine1->NbVertex(); aNVtx++) + { + IntPatch_Point &aVert = theWLine1->ChangeVertex(aNVtx); + const Standard_Real aCurParam = aVert.ParameterOnLine(); + aVert.SetParameter(aNbPts+aCurParam-1); + } + + for(Standard_Integer aNVtx = theWLine2->NbVertex(); aNVtx >= 1; aNVtx--) + { + const IntPatch_Point &aVert = theWLine2->Vertex(aNVtx); + theWLine1->AddVertex(aVert, Standard_True); + } + + theHasBeenJoined = Standard_True; +} + +//======================================================================= +//function : ExtendTwoWLLastFirst +//purpose : Performs extending theWLine1 through its end point and theWLine2 +// through its start point. +//======================================================================= +static void ExtendTwoWLLastFirst(const IntSurf_Quadric& theS1, + const IntSurf_Quadric& theS2, + const Handle(IntPatch_WLine)& theWLine1, + const Handle(IntPatch_WLine)& theWLine2, + const IntSurf_PntOn2S& thePtWL1, + const IntSurf_PntOn2S& thePtWL2, + const gp_Vec& theVec1, + const gp_Vec& theVec2, + const gp_Vec& theVec3, + const Bnd_Box2d& theBoxS1, + const Bnd_Box2d& theBoxS2, + const Standard_Real theToler3D, + const Standard_Real theU1Period, + const Standard_Real theU2Period, + const Standard_Real theV1Period, + const Standard_Real theV2Period, + unsigned int &theCheckResult, + Standard_Boolean &theHasBeenJoined) +{ + IntSurf_PntOn2S aPOn2S; + if(!CheckArgumentsToExtend(theS1, theS2, thePtWL1, thePtWL2, aPOn2S, + theVec1, theVec2, theVec3, + theBoxS1, theBoxS2, theToler3D, + theU1Period, theU2Period, theV1Period, theV2Period)) + { + return; + } + + theCheckResult |= IntPatchWT_DisLastLast; + + ExtendLast (theWLine1, aPOn2S); + ExtendFirst(theWLine2, aPOn2S); + + if(theHasBeenJoined) + { + return; + } + + Standard_Real aPrm = theWLine1->Vertex(theWLine1->NbVertex()).ParameterOnLine(); + while(theWLine1->Vertex(theWLine1->NbVertex()).ParameterOnLine() == aPrm) + theWLine1->RemoveVertex(theWLine1->NbVertex()); + + aPrm = theWLine2->Vertex(1).ParameterOnLine(); + while(theWLine2->Vertex(1).ParameterOnLine() == aPrm) + theWLine2->RemoveVertex(1); + + const Standard_Integer aNbPts = theWLine1->NbPnts(); + for(Standard_Integer aNPt = 2; aNPt <= theWLine2->NbPnts(); aNPt++) + { + const IntSurf_PntOn2S& aPt = theWLine2->Point(aNPt); + theWLine1->Curve()->Add(aPt); + } + + for(Standard_Integer aNVtx = 1; aNVtx <= theWLine2->NbVertex(); aNVtx++) + { + IntPatch_Point &aVert = theWLine2->ChangeVertex(aNVtx); + const Standard_Real aCurParam = aVert.ParameterOnLine(); + aVert.SetParameter(aNbPts+aCurParam-1); + theWLine1->AddVertex(aVert, Standard_False); + } + + theHasBeenJoined = Standard_True; +} + +//======================================================================= +//function : ExtendTwoWLLastLast +//purpose : +//======================================================================= +static void ExtendTwoWLLastLast(const IntSurf_Quadric& theS1, + const IntSurf_Quadric& theS2, + const Handle(IntPatch_WLine)& theWLine1, + const Handle(IntPatch_WLine)& theWLine2, + const IntSurf_PntOn2S& thePtWL1, + const IntSurf_PntOn2S& thePtWL2, + const gp_Vec& theVec1, + const gp_Vec& theVec2, + const gp_Vec& theVec3, + const Bnd_Box2d& theBoxS1, + const Bnd_Box2d& theBoxS2, + const Standard_Real theToler3D, + const Standard_Real theU1Period, + const Standard_Real theU2Period, + const Standard_Real theV1Period, + const Standard_Real theV2Period, + unsigned int &/*theCheckResult*/, + Standard_Boolean &theHasBeenJoined) +{ + IntSurf_PntOn2S aPOn2S; + if(!CheckArgumentsToExtend(theS1, theS2, thePtWL1, thePtWL2, aPOn2S, + theVec1, theVec2, theVec3, + theBoxS1, theBoxS2, theToler3D, + theU1Period, theU2Period, theV1Period, theV2Period)) + { + return; + } + + //theCheckResult |= IntPatchWT_DisLastLast; + + ExtendLast(theWLine1, aPOn2S); + ExtendLast(theWLine2, aPOn2S); + + if(theHasBeenJoined) + return; + + Standard_Real aPrm = theWLine1->Vertex(theWLine1->NbVertex()).ParameterOnLine(); + while(theWLine1->Vertex(theWLine1->NbVertex()).ParameterOnLine() == aPrm) + theWLine1->RemoveVertex(theWLine1->NbVertex()); + + aPrm = theWLine2->Vertex(theWLine2->NbVertex()).ParameterOnLine(); + while(theWLine2->Vertex(theWLine2->NbVertex()).ParameterOnLine() == aPrm) + theWLine2->RemoveVertex(theWLine2->NbVertex()); + + const Standard_Integer aNbPts = theWLine1->NbPnts() + theWLine2->NbPnts(); + for(Standard_Integer aNPt = theWLine2->NbPnts()-1; aNPt >= 1; aNPt--) + { + const IntSurf_PntOn2S& aPt = theWLine2->Point(aNPt); + theWLine1->Curve()->Add(aPt); + } + + for(Standard_Integer aNVtx = theWLine2->NbVertex(); aNVtx >= 1; aNVtx--) + { + IntPatch_Point &aVert = theWLine2->ChangeVertex(aNVtx); + const Standard_Real aCurParam = aVert.ParameterOnLine(); + aVert.SetParameter(aNbPts - aCurParam); + theWLine1->AddVertex(aVert, Standard_False); + } + + theHasBeenJoined = Standard_True; +} + +#endif +//========================================================================= +// function : ComputePurgedWLine +// purpose : +//========================================================================= +Handle(IntPatch_WLine) IntPatch_WLineTool:: + ComputePurgedWLine(const Handle(IntPatch_WLine) &theWLine, + const Handle(Adaptor3d_HSurface) &theS1, + const Handle(Adaptor3d_HSurface) &theS2, + const Handle(Adaptor3d_TopolTool) &theDom1, + const Handle(Adaptor3d_TopolTool) &theDom2, + const Standard_Boolean theRestrictLine) +{ + Standard_Integer i, k, v, nb, nbvtx; + Handle(IntPatch_WLine) aResult; + nbvtx = theWLine->NbVertex(); + nb = theWLine->NbPnts(); + if (nb==2) + { + const IntSurf_PntOn2S& p1 = theWLine->Point(1); + const IntSurf_PntOn2S& p2 = theWLine->Point(2); + if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) + return aResult; + } + + Handle(IntPatch_WLine) aLocalWLine; + Handle(IntPatch_WLine) aTmpWLine = theWLine; + Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S(); + aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False); + for(i = 1; i <= nb; i++) + aLineOn2S->Add(theWLine->Point(i)); + + for(v = 1; v <= nbvtx; v++) + aLocalWLine->AddVertex(theWLine->Vertex(v)); + + // I: Delete equal points + for(i = 1; i <= aLineOn2S->NbPoints(); i++) + { + Standard_Integer aStartIndex = i + 1; + Standard_Integer anEndIndex = i + 5; + nb = aLineOn2S->NbPoints(); + anEndIndex = (anEndIndex > nb) ? nb : anEndIndex; + + if((aStartIndex > nb) || (anEndIndex <= 1)) + continue; + + k = aStartIndex; + + while(k <= anEndIndex) + { + if(i != k) + { + IntSurf_PntOn2S p1 = aLineOn2S->Value(i); + IntSurf_PntOn2S p2 = aLineOn2S->Value(k); + + Standard_Real UV[8]; + p1.Parameters(UV[0], UV[1], UV[2], UV[3]); + p2.Parameters(UV[4], UV[5], UV[6], UV[7]); + + Standard_Real aMax = Abs(UV[0]); + for(Standard_Integer anIdx = 1; anIdx < 8; anIdx++) + { + if (aMax < Abs(UV[anIdx])) + aMax = Abs(UV[anIdx]); + } + + if(p1.Value().IsEqual(p2.Value(), gp::Resolution()) || + Abs(UV[0] - UV[4]) + Abs(UV[1] - UV[5]) < 1.0e-16 * aMax || + Abs(UV[2] - UV[6]) + Abs(UV[3] - UV[7]) < 1.0e-16 * aMax ) + { + aTmpWLine = aLocalWLine; + aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False); + + for(v = 1; v <= aTmpWLine->NbVertex(); v++) + { + IntPatch_Point aVertex = aTmpWLine->Vertex(v); + Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine(); + + if(avertexindex >= k) + { + aVertex.SetParameter(aVertex.ParameterOnLine() - 1.); + } + aLocalWLine->AddVertex(aVertex); + } + aLineOn2S->RemovePoint(k); + anEndIndex--; + continue; + } + } + k++; + } + } + + if (aLineOn2S->NbPoints() <= 2) + { + if (aLineOn2S->NbPoints() == 2) + return aLocalWLine; + else + return aResult; + } + + // Avoid purge in case of C0 continuity: + // Intersection approximator may produce invalid curve after purge, example: + // bugs modalg_5 bug24731. + // Do not run purger when base number of points is too small. + if (theS1->UContinuity() == GeomAbs_C0 || + theS1->VContinuity() == GeomAbs_C0 || + theS2->UContinuity() == GeomAbs_C0 || + theS2->VContinuity() == GeomAbs_C0 || + nb < aNbSingleBezier) + { + return aLocalWLine; + } + + if (theRestrictLine) + { + // II: Delete out of borders points. + aLocalWLine = DeleteOuterPoints(aLocalWLine, theS1, theS2, theDom1, theDom2); + } + + // III: Delete points by tube criteria. + Handle(IntPatch_WLine) aLocalWLineTube = + DeleteByTube(aLocalWLine, theS1, theS2); + + if(aLocalWLineTube->NbPnts() > 1) + { + aResult = aLocalWLineTube; + } + return aResult; +} + + +//======================================================================= +//function : JoinWLines +//purpose : +//======================================================================= +void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin, + IntPatch_SequenceOfPoint& theSPnt, + const Standard_Real theTol3D, + const Standard_Real theU1Period, + const Standard_Real theU2Period, + const Standard_Real theV1Period, + const Standard_Real theV2Period, + const Standard_Real theUfSurf1, + const Standard_Real theUlSurf1, + const Standard_Real theVfSurf1, + const Standard_Real theVlSurf1, + const Standard_Real theUfSurf2, + const Standard_Real theUlSurf2, + const Standard_Real theVfSurf2, + const Standard_Real theVlSurf2) +{ + if(theSlin.Length() == 0) + return; + + for(Standard_Integer aNumOfLine1 = 1; aNumOfLine1 <= theSlin.Length(); aNumOfLine1++) + { + Handle(IntPatch_WLine) aWLine1 (Handle(IntPatch_WLine)::DownCast(theSlin.Value(aNumOfLine1))); + + if(aWLine1.IsNull()) + {//We must have failed to join not-point-lines + continue; + } + + const Standard_Integer aNbPntsWL1 = aWLine1->NbPnts(); + const IntSurf_PntOn2S& aPntFW1 = aWLine1->Point(1); + const IntSurf_PntOn2S& aPntLW1 = aWLine1->Point(aNbPntsWL1); + + for(Standard_Integer aNPt = 1; aNPt <= theSPnt.Length(); aNPt++) + { + const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNPt).PntOn2S(); + + if( aPntCur.IsSame(aPntFW1, Precision::Confusion()) || + aPntCur.IsSame(aPntLW1, Precision::Confusion())) + { + theSPnt.Remove(aNPt); + aNPt--; + } + } + + Standard_Boolean hasBeenRemoved = Standard_False; + for(Standard_Integer aNumOfLine2 = aNumOfLine1 + 1; aNumOfLine2 <= theSlin.Length(); aNumOfLine2++) + { + Handle(IntPatch_WLine) aWLine2 (Handle(IntPatch_WLine)::DownCast(theSlin.Value(aNumOfLine2))); + + if(aWLine2.IsNull()) + continue; + + const Standard_Integer aNbPntsWL2 = aWLine2->NbPnts(); + + const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1); + const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aNbPntsWL1); + + const IntSurf_PntOn2S& aPntFWL2 = aWLine2->Point(1); + const IntSurf_PntOn2S& aPntLWL2 = aWLine2->Point(aNbPntsWL2); + + if(aPntFWL1.IsSame(aPntFWL2, Precision::Confusion())) + { + const IntSurf_PntOn2S& aPt1 = aWLine1->Point(2); + const IntSurf_PntOn2S& aPt2 = aWLine2->Point(2); + Standard_Boolean aCond = + CheckArgumentsToJoin(gp_Vec(aPntFWL1.Value(), aPt1.Value()), + gp_Vec(aPt2.Value(), aPntFWL2.Value())); + + aCond = aCond && !IsSeamOrBound(aPt1, aPt2, aPntFWL1, + theU1Period, theU2Period, + theV1Period, theV2Period, + theUfSurf1, theUlSurf1, + theVfSurf1, theVlSurf1, + theUfSurf2, theUlSurf2, + theVfSurf2, theVlSurf2); + + if(aCond) + { + aWLine1->ClearVertexes(); + for(Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++) + { + const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt); + aWLine1->Curve()->InsertBefore(1, aPt); + } + + aWLine1->ComputeVertexParameters(theTol3D); + + theSlin.Remove(aNumOfLine2); + aNumOfLine2--; + hasBeenRemoved = Standard_True; + + continue; + } + } + + if(aPntFWL1.IsSame(aPntLWL2, Precision::Confusion())) + { + const IntSurf_PntOn2S& aPt1 = aWLine1->Point(2); + const IntSurf_PntOn2S& aPt2 = aWLine2->Point(aNbPntsWL2-1); + + Standard_Boolean aCond = + CheckArgumentsToJoin(gp_Vec(aPntFWL1.Value(), aPt1.Value()), + gp_Vec(aPt2.Value(), aPntLWL2.Value())); + + aCond = aCond && !IsSeamOrBound(aPt1, aPt2, aPntFWL1, + theU1Period, theU2Period, + theV1Period, theV2Period, + theUfSurf1, theUlSurf1, + theVfSurf1, theVlSurf1, + theUfSurf2, theUlSurf2, + theVfSurf2, theVlSurf2); + + if(aCond) + { + aWLine1->ClearVertexes(); + for(Standard_Integer aNPt = aNbPntsWL2; aNPt >= 1; aNPt--) + { + const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt); + aWLine1->Curve()->InsertBefore(1, aPt); + } + + aWLine1->ComputeVertexParameters(theTol3D); + + theSlin.Remove(aNumOfLine2); + aNumOfLine2--; + hasBeenRemoved = Standard_True; + + continue; + } + } + + if(aPntLWL1.IsSame(aPntFWL2, Precision::Confusion())) + { + const IntSurf_PntOn2S& aPt1 = aWLine1->Point(aNbPntsWL1-1); + const IntSurf_PntOn2S& aPt2 = aWLine2->Point(2); + + Standard_Boolean aCond = + CheckArgumentsToJoin(gp_Vec(aPt1.Value(), aPntLWL1.Value()), + gp_Vec(aPntFWL2.Value(), aPt2.Value())); + + aCond = aCond && !IsSeamOrBound(aPt1, aPt2, aPntLWL1, + theU1Period, theU2Period, + theV1Period, theV2Period, + theUfSurf1, theUlSurf1, + theVfSurf1, theVlSurf1, + theUfSurf2, theUlSurf2, + theVfSurf2, theVlSurf2); + + if(aCond) + { + aWLine1->ClearVertexes(); + for(Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++) + { + const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt); + aWLine1->Curve()->Add(aPt); + } + + aWLine1->ComputeVertexParameters(theTol3D); + + theSlin.Remove(aNumOfLine2); + aNumOfLine2--; + hasBeenRemoved = Standard_True; + + continue; + } + } + + if(aPntLWL1.IsSame(aPntLWL2, Precision::Confusion())) + { + const IntSurf_PntOn2S& aPt1 = aWLine1->Point(aNbPntsWL1-1); + const IntSurf_PntOn2S& aPt2 = aWLine2->Point(aNbPntsWL2-1); + + Standard_Boolean aCond = + CheckArgumentsToJoin(gp_Vec(aPt1.Value(), aPntLWL1.Value()), + gp_Vec(aPntLWL2.Value(), aPt2.Value())); + + aCond = aCond && !IsSeamOrBound(aPt1, aPt2, aPntLWL1, + theU1Period, theU2Period, + theV1Period, theV2Period, + theUfSurf1, theUlSurf1, + theVfSurf1, theVlSurf1, + theUfSurf2, theUlSurf2, + theVfSurf2, theVlSurf2); + + if(aCond) + { + aWLine1->ClearVertexes(); + for(Standard_Integer aNPt = aNbPntsWL2; aNPt >= 1; aNPt--) + { + const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt); + aWLine1->Curve()->Add(aPt); + } + + aWLine1->ComputeVertexParameters(theTol3D); + + theSlin.Remove(aNumOfLine2); + aNumOfLine2--; + hasBeenRemoved = Standard_True; + + continue; + } + } + } + + if(hasBeenRemoved) + aNumOfLine1--; + } +} + +//======================================================================= +//function : ExtendTwoWlinesToEachOther +//purpose : Performs extending theWLine1 and theWLine2 through their +// respecting end point. +//======================================================================= +void IntPatch_WLineTool::ExtendTwoWlinesToEachOther(IntPatch_SequenceOfLine& /*theSlin*/, + const IntSurf_Quadric& /*theS1*/, + const IntSurf_Quadric& /*theS2*/, + const Standard_Real /*theToler3D*/, + const Standard_Real /*theU1Period*/, + const Standard_Real /*theU2Period*/, + const Standard_Real /*theV1Period*/, + const Standard_Real /*theV2Period*/, + const Bnd_Box2d& /*theBoxS1*/, + const Bnd_Box2d& /*theBoxS2*/) +{ +#if 0 + if(theSlin.Length() < 2) + return; + + gp_Vec aVec1, aVec2, aVec3; + + for(Standard_Integer aNumOfLine1 = 1; aNumOfLine1 <= theSlin.Length(); aNumOfLine1++) + { + Handle(IntPatch_WLine) aWLine1 (Handle(IntPatch_WLine):: + DownCast(theSlin.Value(aNumOfLine1))); + + if(aWLine1.IsNull()) + {//We must have failed to join not-point-lines + continue; + } + + const Standard_Integer aNbPntsWL1 = aWLine1->NbPnts(); + + if(aWLine1->Vertex(1).ParameterOnLine() != 1) + continue; + + if(aWLine1->Vertex(aWLine1->NbVertex()).ParameterOnLine() != aWLine1->NbPnts()) + continue; + + for(Standard_Integer aNumOfLine2 = aNumOfLine1 + 1; + aNumOfLine2 <= theSlin.Length(); aNumOfLine2++) + { + Handle(IntPatch_WLine) aWLine2 (Handle(IntPatch_WLine):: + DownCast(theSlin.Value(aNumOfLine2))); + + if(aWLine2.IsNull()) + continue; + + if(aWLine2->Vertex(1).ParameterOnLine() != 1) + continue; + + if(aWLine2->Vertex(aWLine2->NbVertex()).ParameterOnLine() != aWLine2->NbPnts()) + continue; + + //Enable/Disable of some ckeck. Bit-mask is used for it. + //E.g. if 1st point of aWLine1 matches with + //1st point of aWLine2 then we do not need in check + //1st point of aWLine1 and last point of aWLine2 etc. + unsigned int aCheckResult = IntPatchWT_EnAll; + + Standard_Boolean hasBeenJoined = Standard_False; + + const Standard_Integer aNbPntsWL2 = aWLine2->NbPnts(); + + const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1); + const IntSurf_PntOn2S& aPntFp1WL1 = aWLine1->Point(2); + + const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aNbPntsWL1); + const IntSurf_PntOn2S& aPntLm1WL1 = aWLine1->Point(aNbPntsWL1-1); + + const IntSurf_PntOn2S& aPntFWL2 = aWLine2->Point(1); + const IntSurf_PntOn2S& aPntFp1WL2 = aWLine2->Point(2); + + const IntSurf_PntOn2S& aPntLWL2 = aWLine2->Point(aNbPntsWL2); + const IntSurf_PntOn2S& aPntLm1WL2 = aWLine2->Point(aNbPntsWL2-1); + + //if(!(aCheckResult & IntPatchWT_DisFirstFirst)) + {// First/First + aVec1.SetXYZ(aPntFp1WL1.Value().XYZ() - aPntFWL1.Value().XYZ()); + aVec2.SetXYZ(aPntFWL2.Value().XYZ() - aPntFp1WL2.Value().XYZ()); + aVec3.SetXYZ(aPntFWL1.Value().XYZ() - aPntFWL2.Value().XYZ()); + + ExtendTwoWLFirstFirst(theS1, theS2, aWLine1, aWLine2, aPntFWL1, aPntFWL2, + aVec1, aVec2, aVec3, theBoxS1, theBoxS2, theToler3D, + theU1Period, theU2Period, theV1Period, theV2Period, + aCheckResult, hasBeenJoined); + } + + if(!(aCheckResult & IntPatchWT_DisFirstLast)) + {// First/Last + aVec1.SetXYZ(aPntFp1WL1.Value().XYZ() - aPntFWL1.Value().XYZ()); + aVec2.SetXYZ(aPntLWL2.Value().XYZ() - aPntLm1WL2.Value().XYZ()); + aVec3.SetXYZ(aPntFWL1.Value().XYZ() - aPntLWL2.Value().XYZ()); + + ExtendTwoWLFirstLast(theS1, theS2, aWLine1, aWLine2, aPntFWL1, aPntLWL2, + aVec1, aVec2, aVec3, theBoxS1, theBoxS2, theToler3D, + theU1Period, theU2Period, theV1Period, theV2Period, + aCheckResult, hasBeenJoined); + } + + if(!(aCheckResult & IntPatchWT_DisLastFirst)) + {// Last/First + aVec1.SetXYZ(aPntLWL1.Value().XYZ() - aPntLm1WL1.Value().XYZ()); + aVec2.SetXYZ(aPntFp1WL2.Value().XYZ() - aPntFWL2.Value().XYZ()); + aVec3.SetXYZ(aPntFWL2.Value().XYZ() - aPntLWL1.Value().XYZ()); + + ExtendTwoWLLastFirst(theS1, theS2, aWLine1, aWLine2, aPntLWL1, aPntFWL2, + aVec1, aVec2, aVec3, theBoxS1, theBoxS2, theToler3D, + theU1Period, theU2Period, theV1Period, theV2Period, + aCheckResult, hasBeenJoined); + } + + if(!(aCheckResult & IntPatchWT_DisLastLast)) + {// Last/Last + aVec1.SetXYZ(aPntLWL1.Value().XYZ() - aPntLm1WL1.Value().XYZ()); + aVec2.SetXYZ(aPntLm1WL2.Value().XYZ() - aPntLWL2.Value().XYZ()); + aVec3.SetXYZ(aPntLWL2.Value().XYZ() - aPntLWL1.Value().XYZ()); + + ExtendTwoWLLastLast(theS1, theS2, aWLine1, aWLine2, aPntLWL1, aPntLWL2, + aVec1, aVec2, aVec3, theBoxS1, theBoxS2, theToler3D, + theU1Period, theU2Period, theV1Period, theV2Period, + aCheckResult, hasBeenJoined); + } + + if(hasBeenJoined) + { + theSlin.Remove(aNumOfLine2); + aNumOfLine2--; + } + } + } +#endif +} diff --git a/src/IntPatch/IntPatch_WLineTool.hxx b/src/IntPatch/IntPatch_WLineTool.hxx new file mode 100644 index 0000000000..02e4584a10 --- /dev/null +++ b/src/IntPatch/IntPatch_WLineTool.hxx @@ -0,0 +1,104 @@ +// Copyright (c) 1999-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _IntPatch_WLineTool_HeaderFile +#define _IntPatch_WLineTool_HeaderFile + +#include +#include +#include +#include +#include +class TopoDS_Face; +class GeomAdaptor_HSurface; +class GeomInt_LineConstructor; +class IntTools_Context; +class Adaptor3d_TopolTool; +class Adaptor3d_HSurface; + +//! IntPatch_WLineTool provides set of static methods related to walking lines. +class IntPatch_WLineTool +{ +public: + + DEFINE_STANDARD_ALLOC + + //! I + //! Removes equal points (leave one of equal points) from theWLine + //! and recompute vertex parameters. + //! + //! II + //! Removes point out of borders in case of non periodic surfaces. + //! This step is done only if theRestrictLine is true. + //! + //! III + //! Removes exceed points using tube criteria: + //! delete 7D point if it lies near to expected lines in 2d and 3d. + //! Each task (2d, 2d, 3d) have its own tolerance and checked separately. + //! + //! Returns new WLine or null WLine if the number + //! of the points is less than 2. + Standard_EXPORT static + Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine) &theWLine, + const Handle(Adaptor3d_HSurface) &theS1, + const Handle(Adaptor3d_HSurface) &theS2, + const Handle(Adaptor3d_TopolTool) &theDom1, + const Handle(Adaptor3d_TopolTool) &theDom2, + const Standard_Boolean theRestrictLine); + +//! Joins all WLines from theSlin to one if it is possible and records +//! the result into theSlin again. Lines will be kept to be splitted if: +//! a) they are separated (has no common points); +//! b) resulted line (after joining) go through seam-edges or surface boundaries. +//! +//! In addition, if points in theSPnt lies at least in one of the line in theSlin, +//! this point will be deleted. + Standard_EXPORT static void JoinWLines(IntPatch_SequenceOfLine& theSlin, + IntPatch_SequenceOfPoint& theSPnt, + const Standard_Real theTol3D, + const Standard_Real theU1Period, + const Standard_Real theU2Period, + const Standard_Real theV1Period, + const Standard_Real theV2Period, + const Standard_Real theUfSurf1, + const Standard_Real theUlSurf1, + const Standard_Real theVfSurf1, + const Standard_Real theVlSurf1, + const Standard_Real theUfSurf2, + const Standard_Real theUlSurf2, + const Standard_Real theVfSurf2, + const Standard_Real theVlSurf2); + + +//! Extends every line from theSlin (if it is possible) to be started/finished +//! in strictly determined point (in the place of joint of two lines). +//! As result, some gaps between two lines will vanish. +//! The Walking lines are supposed (algorithm will do nothing for not-Walking line) +//! to be computed as a result of intersection of two quadrics. +//! The quadrics definition is accepted in input parameters. + Standard_EXPORT static void ExtendTwoWlinesToEachOther(IntPatch_SequenceOfLine& theSlin, + const IntSurf_Quadric& theS1, + const IntSurf_Quadric& theS2, + const Standard_Real theToler3D, + const Standard_Real theU1Period, + const Standard_Real theU2Period, + const Standard_Real theV1Period, + const Standard_Real theV2Period, + const Bnd_Box2d& theBoxS1, + const Bnd_Box2d& theBoxS2); + +//! Max angle to concatenate two WLines to avoid result with C0-continuity + static const Standard_Real myMaxConcatAngle; +}; + +#endif \ No newline at end of file diff --git a/tests/blend/simple/Q6 b/tests/blend/simple/Q6 index 03f132cd90..088b277e75 100644 --- a/tests/blend/simple/Q6 +++ b/tests/blend/simple/Q6 @@ -3,4 +3,4 @@ tscale s 0 0 0 SCALE explode s E blend result s SCALE*2 s_5 -set square 1.6539e+08 +set square 1.65391e+008 diff --git a/tests/bugs/modalg_5/bug25742_2 b/tests/bugs/modalg_5/bug25742_2 index b780d877d5..489ee4d63a 100755 --- a/tests/bugs/modalg_5/bug25742_2 +++ b/tests/bugs/modalg_5/bug25742_2 @@ -56,10 +56,10 @@ if { $z > ${max_time} } { mksurface s1 b1_4 mksurface s2 b2_1 -dchrono h2 stop -set q2 [dchrono h2 show] +dchrono h2 reset +dchrono h2 start -set CurveNumb [intersect i s1 s2] +set CurveNumb [intersect resi s1 s2] dchrono h2 stop set q2 [dchrono h2 show] diff --git a/tests/bugs/modalg_6/bug23178 b/tests/bugs/modalg_6/bug23178 new file mode 100644 index 0000000000..a16ee989dc --- /dev/null +++ b/tests/bugs/modalg_6/bug23178 @@ -0,0 +1,43 @@ +puts "================" +puts "OCC23178" +puts "================" +puts "" +####################################################################### +# Intersection of cylinders fails to produce results +####################################################################### + +set GoodNbCurv 6 + +foreach c [directory result*] { + unset $c +} + +cylinder s1 4.999999999813, 1.35836143386791, 0.20975890778078 -1, 0, 0 0, 0, 1 10 +cylinder s2 9.999999999813, 9.1401960568287, 0.11837300583107 0, -0.75870713028692, 0.651431877061437 0, -0.651431877061437, -0.75870713028692 5 + +intersect result s1 s2 + +foreach c [directory result*] { + bounds $c U1 U2 + + if {[dval U2-U1] < 1.0e-9} { + puts "Error: Wrong curve's range!" + } + + xdistcs $c s1 U1 U2 10 2.0e-7 + xdistcs $c s2 U1 U2 10 2.0e-7 +} + +set NbCurv [llength [directory result*]] + +if { $NbCurv == $GoodNbCurv } { + puts "OK: Number of curves is good!" +} else { + puts "Error: $GoodNbCurv is expected but $NbCurv is found!" +} + +smallview +don result* +fit +disp s1 s2 +checkview -screenshot -2d -path ${imagedir}/${test_image}.png diff --git a/tests/bugs/modalg_6/bug27766 b/tests/bugs/modalg_6/bug27766 new file mode 100644 index 0000000000..c039cd8efc --- /dev/null +++ b/tests/bugs/modalg_6/bug27766 @@ -0,0 +1,30 @@ +puts "========" +puts "OCC27766" +puts "========" +puts "" +################################################# +# Incorrect section curves between attached cylinders +################################################# + +set ExpTol 1.0e-7 +set GoodNbCurv 3 + +restore [locate_data_file bug27761_c1.brep] c1 +restore [locate_data_file bug27761_c2.brep] c2 + +set log [bopcurves c1 c2 -2d] + +regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} ${log} full Toler NbCurv + +if {${NbCurv} != ${GoodNbCurv}} { + puts "Error: Number of curves is bad!" +} + +checkreal TolReached $Toler $ExpTol 0.0 0.1 + +smallview +don c_* +fit +disp c1 c2 + +checkview -screenshot -2d -path ${imagedir}/${test_image}.png -- 2.39.5