// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.
+#include <algorithm>
+#include <Bnd_Range.hxx>
#include <IntAna_ListOfCurve.hxx>
-#include <IntAna_ListIteratorOfListOfCurve.hxx>
-#include <IntPatch_WLine.hxx>
-
#include <math_Matrix.hxx>
+#include <NCollection_IncAllocator.hxx>
+#include <Standard_DivideByZero.hxx>
//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,
- const gp_Cone& aCo,
- 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 ExploreCurve(const gp_Cone& theCo,
+ IntAna_Curve& aC,
+ const Standard_Real aTol,
+ IntAna_ListOfCurve& aLC);
static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
const Standard_Real theUlTarget,
const Standard_Real thePeriod,
const Standard_Boolean theFlForce);
+
+class ComputationMethods
+{
+ //Every cylinder can be represented by the following equation in parametric form:
+ // S(U,V) = L + R*cos(U)*Xd+R*sin(U)*Yd+V*Zd,
+ //where location L, directions Xd, Yd and Zd have type gp_XYZ.
+
+ //Intersection points between two cylinders can be found from the following system:
+ // S1(U1, V1) = S2(U2, V2)
+ //or
+ // {X01 + R1*cos(U1)*Xx1 + R1*sin(U1)*Yx1 + V1*Zx1 = X02 + R2*cos(U2)*Xx2 + R2*sin(U2)*Yx2 + V2*Zx2
+ // {Y01 + R1*cos(U1)*Xy1 + R1*sin(U1)*Yy1 + V1*Zy1 = Y02 + R2*cos(U2)*Xy2 + R2*sin(U2)*Yy2 + V2*Zy2 (1)
+ // {Z01 + R1*cos(U1)*Xz1 + R1*sin(U1)*Yz1 + V1*Zz1 = Z02 + R2*cos(U2)*Xz2 + R2*sin(U2)*Yz2 + V2*Zz2
+
+ //The formula (1) can be rewritten as follows
+ // {C11*V1+C21*V2=A11*cos(U1)+B11*sin(U1)+A21*cos(U2)+B21*sin(U2)+D1
+ // {C12*V1+C22*V2=A12*cos(U1)+B12*sin(U1)+A22*cos(U2)+B22*sin(U2)+D2 (2)
+ // {C13*V1+C23*V2=A13*cos(U1)+B13*sin(U1)+A23*cos(U2)+B23*sin(U2)+D3
+
+ //Hereafter we consider that in system
+ // {C11*V1+C21*V2=A11*cos(U1)+B11*sin(U1)+A21*cos(U2)+B21*sin(U2)+D1 (3)
+ // {C12*V1+C22*V2=A12*cos(U1)+B12*sin(U1)+A22*cos(U2)+B22*sin(U2)+D2
+ //variables V1 and V2 can be found unambiguously, i.e. determinant
+ // |C11 C21|
+ // | | != 0
+ // |C12 C22|
+ //
+ //In this case, variables V1 and V2 can be found as follows:
+ // {V1 = K11*sin(U1)+K21*sin(U2)+L11*cos(U1)+L21*cos(U2)+M1 = K1*cos(U1-FIV1)+L1*cos(U2-PSIV1)+M1 (4)
+ // {V2 = K12*sin(U1)+K22*sin(U2)+L12*cos(U1)+L22*cos(U2)+M2 = K2*cos(U2-FIV2)+L2*cos(U2-PSIV2)+M2
+
+ //Having substituted result of (4) to the 3rd equation of (2), we will obtain equation
+ // cos(U2-FI2) = B*cos(U1-FI1)+C. (5)
+
+ //I.e. when U1 is taken different given values (from domain), correspond U2 value can be computed
+ //from equation (5). After that, V1 and V2 can be computed from the system (4) (see
+ //CylCylComputeParameters(...) methods).
+
+ //It is important to remark that equation (5) (in general) has two solutions: U2=FI2 +/- f(U1).
+ //Therefore, we are getting here two intersection lines.
+
+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())
+ {
+ throw Standard_Failure("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-zero 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)
+ {
+ };
+
+ // Returns parameters of system solved while finding
+ // intersection line
+ const ComputationMethods::stCoeffsValue &SICoeffs() const
+ {
+ return myCoeffs;
+ }
+
+ // Returns quadric correspond to the index theIdx.
+ const IntSurf_Quadric& GetQSurface(const Standard_Integer theIdx) const
+ {
+ if (theIdx <= 1)
+ return myQuad1;
+
+ return myQuad2;
+ }
+
+ // Returns TRUE in case of reverting surfaces
+ Standard_Boolean IsReversed() const
+ {
+ return myIsReverse;
+ }
+
+ // Returns 2D-tolerance
+ Standard_Real Get2dTolerance() const
+ {
+ return myTol2D;
+ }
+
+ // Returns 3D-tolerance
+ Standard_Real Get3dTolerance() const
+ {
+ return myTol3D;
+ }
+
+ // Returns UV-bounds of 1st surface
+ const Bnd_Box2d& UVS1() const
+ {
+ return myUVSurf1;
+ }
+
+ // Returns UV-bounds of 2nd surface
+ const Bnd_Box2d& UVS2() const
+ {
+ return myUVSurf2;
+ }
+
+ void AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL,
+ const Standard_Real theU1,
+ const Standard_Real theU1Min,
+ 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;
+
+ static Standard_Boolean BoundariesComputing(const ComputationMethods::stCoeffsValue &theCoeffs,
+ const Standard_Real thePeriod,
+ Bnd_Range theURange[]);
+
+ void BoundaryEstimation(const gp_Cylinder& theCy1,
+ const gp_Cylinder& theCy2,
+ Bnd_Range& theOutBoxS1,
+ Bnd_Range& theOutBoxS2) const;
+
+protected:
+
+ //Solves equation (2) (see declaration of ComputationMethods class) in case,
+ //when V1 or V2 (is set by theSBType argument) is known (corresponds to the boundary
+ //and equal to theVzad) but U1 is unknown. Computation is made by numeric methods and
+ //requires initial values (theVInit, theInitU2 and theInitMainVar).
+ 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 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);
}
}
+//=======================================================================
+//function : ExtremaLineLine
+//purpose : Computes extrema between the given lines. Returns parameters
+// on correspond curve (see correspond method for Extrema_ExtElC class).
+//=======================================================================
+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
-// if U1 increases. But if it is not, new V1set and/or V2set
+//purpose : By default, we shall consider, that V1 and V2 will be increased
+// if U1 is increased. But if it is not, new V1set and/or V2set
// must be computed as [V1current - DeltaV1] (analogically
// for V2). This function processes this case.
//=======================================================================
static void VBoundaryPrecise( const math_Matrix& theMatr,
const Standard_Real theV1AfterDecrByDelta,
const Standard_Real theV2AfterDecrByDelta,
- const Standard_Real theV1f,
- const Standard_Real theV2f,
Standard_Real& theV1Set,
Standard_Real& theV2Set)
{
aSyst.SetCol(2, theMatr.Col(2));
aSyst.SetCol(3, theMatr.Col(4));
+ //We have the system (see comment to StepComputing(...) function)
+ // {a11*dV1 + a12*dV2 + a14*dU2 = -a13*dU1
+ // {a21*dV1 + a22*dV2 + a24*dU2 = -a23*dU1
+ // {a31*dV1 + a32*dV2 + a34*dU2 = -a33*dU1
+
const Standard_Real aDet = aSyst.Determinant();
aSyst.SetCol(1, theMatr.Col(3));
const Standard_Real aDet2 = aSyst.Determinant();
- if(aDet*aDet1 > 0.0)
+ //Now,
+ // dV1 = -dU1*aDet1/aDet
+ // dV2 = -dU1*aDet2/aDet
+
+ //If U1 is increased then dU1 > 0.
+ //If (aDet1/aDet > 0) then dV1 < 0 and
+ //V1 will be decreased after increasing U1.
+
+ //We have analogical situation with V2-parameter.
+
+ if(aDet*aDet1 > 0.0)
{
- theV1Set = Max(theV1AfterDecrByDelta, theV1f);
+ theV1Set = theV1AfterDecrByDelta;
}
if(aDet*aDet2 > 0.0)
{
- theV2Set = Max(theV2AfterDecrByDelta, theV2f);
+ theV2Set = theV2AfterDecrByDelta;
}
}
const Standard_Real theV2Cur,
const Standard_Real theDeltaV1,
const Standard_Real theDeltaV2,
- const Standard_Real theV1f,
- const Standard_Real theV1l,
- const Standard_Real theV2f,
- const Standard_Real theV2l,
Standard_Real& theDeltaU1Found/*,
Standard_Real& theDeltaU2Found,
Standard_Real& theV1Found,
Standard_Real& theV2Found*/)
{
+#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
+ bool flShow = false;
+
+ if(flShow)
+ {
+ printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
+ theMatr(1,1), theMatr(1,2), theMatr(1,3), theMatr(1,4), theMatr(1,5));
+ printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
+ theMatr(2,1), theMatr(2,2), theMatr(2,3), theMatr(2,4), theMatr(2,5));
+ printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
+ theMatr(3,1), theMatr(3,2), theMatr(3,3), theMatr(3,4), theMatr(3,5));
+ }
+#endif
+
Standard_Boolean isSuccess = Standard_False;
theDeltaU1Found/* = theDeltaU2Found*/ = RealLast();
//theV1Found = theV1set;
//By default, increasing V1(U1) and V2(U1) functions is
//considered
- Standard_Real aV1Set = Min(theV1Cur + theDeltaV1, theV1l),
- aV2Set = Min(theV2Cur + theDeltaV2, theV2l);
+ Standard_Real aV1Set = theV1Cur + theDeltaV1,
+ aV2Set = theV2Cur + theDeltaV2;
- //However, what is indead?
- VBoundaryPrecise( theMatr, theV1Cur - theDeltaV1, theV2Cur - theDeltaV2,
- theV1f, theV2f, aV1Set, aV2Set);
+ //However, what is indeed?
+ VBoundaryPrecise( theMatr, theV1Cur - theDeltaV1,
+ theV2Cur - theDeltaV2, aV1Set, aV2Set);
aSyst.SetCol(2, theMatr.Col(3));
aSyst.SetCol(3, theMatr.Col(4));
if (!procf) {
d=ptf.Distance(ptsol.Value());
if (d <= Tol) {
+ ptsol.SetTolerance(Tol);
if (!ptsol.IsMultiple()) {
//-- le point ptsol (de aligold) est declare multiple sur aligold
Multpoint = Standard_True;
}
if (!procl) {
if (ptl.Distance(ptsol.Value()) <= Tol) {
+ ptsol.SetTolerance(Tol);
if (!ptsol.IsMultiple()) {
Multpoint = Standard_True;
ptsol.SetMultiple(Standard_True);
}
}
}
+
+ ptsol.SetTolerance(Tol);
if (!procf && !procl) {
Quad1.Parameters(ptf,U1,V1);
Quad2.Parameters(ptf,U2,V2);
alig->SetLastPoint(alig->NbVertex());
}
}
+
//=======================================================================
-//function : IntCyCy
-//purpose :
+//function : CyCyAnalyticalIntersect
+//purpose : Checks if intersection curve is analytical (line, circle, ellipse)
+// and returns these curves.
//=======================================================================
-Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
- const IntSurf_Quadric& Quad2,
- const Standard_Real Tol,
- Standard_Boolean& Empty,
- Standard_Boolean& Same,
- Standard_Boolean& Multpoint,
- IntPatch_SequenceOfLine& slin,
- IntPatch_SequenceOfPoint& spnt)
-
+Standard_Boolean CyCyAnalyticalIntersect( const IntSurf_Quadric& Quad1,
+ const IntSurf_Quadric& Quad2,
+ const IntAna_QuadQuadGeo& theInter,
+ const Standard_Real Tol,
+ Standard_Boolean& Empty,
+ Standard_Boolean& Same,
+ Standard_Boolean& Multpoint,
+ IntPatch_SequenceOfLine& slin,
+ IntPatch_SequenceOfPoint& spnt)
{
IntPatch_Point ptsol;
-
+
Standard_Integer i;
IntSurf_TypeTrans trans1,trans2;
gp_Cylinder Cy1(Quad1.Cylinder());
gp_Cylinder Cy2(Quad2.Cylinder());
- IntAna_QuadQuadGeo inter(Cy1,Cy2,Tol);
-
- if (!inter.IsDone())
- {
- return Standard_False;
- }
-
- typint = inter.TypeInter();
- Standard_Integer NbSol = inter.NbSolutions();
+ typint = theInter.TypeInter();
+ Standard_Integer NbSol = theInter.NbSolutions();
Empty = Standard_False;
Same = Standard_False;
switch (typint)
- {
+ {
case IntAna_Empty:
{
Empty = Standard_True;
- }
+ }
break;
case IntAna_Same:
- {
+ {
Same = Standard_True;
- }
+ }
break;
-
+
case IntAna_Point:
{
- gp_Pnt psol(inter.Point(1));
- Standard_Real U1,V1,U2,V2;
- Quad1.Parameters(psol,U1,V1);
- Quad2.Parameters(psol,U2,V2);
+ gp_Pnt psol(theInter.Point(1));
ptsol.SetValue(psol,Tol,Standard_True);
+
+ Standard_Real U1,V1,U2,V2;
+ Quad1.Parameters(psol, U1, V1);
+ Quad2.Parameters(psol, U2, V2);
+
ptsol.SetParameters(U1,V1,U2,V2);
spnt.Append(ptsol);
}
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;
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;
}
}
else
- {
+ {
if (Cy1.Radius() < Cy2.Radius())
{
if (norm1.Dot(crb1) > 0.)
{
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)
{
gp_Vec Tgt;
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));
-
+
Multpoint = Standard_True;
pmult1.SetValue(pttang1,Tol,Standard_True);
pmult2.SetValue(pttang2,Tol,Standard_True);
pmult1.SetMultiple(Standard_True);
pmult2.SetMultiple(Standard_True);
-
+
Standard_Real oU1,oV1,oU2,oV2;
- Quad1.Parameters(pttang1,oU1,oV1);
- Quad2.Parameters(pttang1,oU2,oV2);
- pmult1.SetParameters(oU1,oV1,oU2,oV2);
+ Quad1.Parameters(pttang1, oU1, oV1);
+ Quad2.Parameters(pttang1, oU2, oV2);
+ pmult1.SetParameters(oU1,oV1,oU2,oV2);
Quad1.Parameters(pttang2,oU1,oV1);
Quad2.Parameters(pttang2,oU2,oV2);
+
pmult2.SetParameters(oU1,oV1,oU2,oV2);
-
- // on traite la premiere ellipse
+ // 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)
{
trans2 = IntSurf_Out;
}
else
- {
+ {
trans1=trans2=IntSurf_Undecided;
}
//
Quad1.Parameters(aP, aU1, aV1);
Quad2.Parameters(aP, aU2, aV2);
+
aIP.SetParameters(aU1, aV1, aU2, aV2);
//
aIP.SetParameter(0.);
//-- 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);
//-- 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)
{
Standard_Real aU1, aV1, aU2, aV2;
IntPatch_Point aIP;
gp_Pnt aP (ElCLib::Value(0., elipsol));
- //
- 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;
- }
- }
-
- if(Abs(aDetV1V2) < aNulValue)
- {
- 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);
+ //
+ aIP.SetValue(aP,Tol,Standard_False);
+ aIP.SetMultiple(Standard_False);
+ //
- aVTemp = mVecD;
- mVecD(2) = aVTemp(3);
- mVecD(3) = aVTemp(2);
+ 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);
+ case IntAna_Parabola:
+ case IntAna_Hyperbola:
+ throw Standard_Failure("IntCyCy(): Wrong intersection type!");
- 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_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,
- const Standard_Integer theWLIndex,
- const stCoeffsValue& theCoeffs,
- const Standard_Real thePeriod,
- Standard_Boolean& theIsIncreasing)
+Standard_Boolean ComputationMethods::CylCylMonotonicity(const Standard_Real theU1par,
+ const Standard_Integer theWLIndex,
+ const stCoeffsValue& theCoeffs,
+ const Standard_Real thePeriod,
+ Standard_Boolean& theIsIncreasing)
{
// U2(U1) = FI2 + (+/-)acos(B*cos(U1 - FI1) + C);
//Accordingly,
isPlus = Standard_False;
break;
default:
- //Standard_Failure::Raise("Error. Range Error!!!!");
+ //throw Standard_Failure("Error. Range Error!!!!");
return Standard_False;
}
return Standard_True;
}
-static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
+//=======================================================================
+//function : CylCylComputeParameters
+//purpose : Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0,
+// estimates the tolerance of U2-computing (estimation result is
+// assigned to *theDelta value).
+//=======================================================================
+Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
const Standard_Integer theWLIndex,
const stCoeffsValue& theCoeffs,
- Standard_Real& theU2)
+ Standard_Real& theU2,
+ Standard_Real* const theDelta)
{
- Standard_Real aSign = 0.0;
+ //This formula is got from some experience and can be changed.
+ const Standard_Real aTol0 = Min(10.0*Epsilon(1.0)*theCoeffs.mB, aNulValue);
+ const Standard_Real aTol = 1.0 - aTol0;
- switch(theWLIndex)
- {
- case 0:
- aSign = 1.0;
- break;
- case 1:
- aSign = -1.0;
- break;
- default:
- //Standard_Failure::Raise("Error. Range Error!!!!");
+ if(theWLIndex < 0 || theWLIndex > 1)
return Standard_False;
- }
- Standard_Real anArg = theCoeffs.mB *
- cos(theU1par - theCoeffs.mFI1) +
- theCoeffs.mC;
+ const Standard_Real aSign = theWLIndex ? -1.0 : 1.0;
+
+ Standard_Real anArg = cos(theU1par - theCoeffs.mFI1);
+ anArg = theCoeffs.mB*anArg + theCoeffs.mC;
+
+ if(anArg >= aTol)
+ {
+ if(theDelta)
+ *theDelta = 0.0;
- if(aNulValue > 1.0 - anArg)
anArg = 1.0;
+ }
+ else if(anArg <= -aTol)
+ {
+ if(theDelta)
+ *theDelta = 0.0;
- if(anArg + 1.0 < aNulValue)
anArg = -1.0;
+ }
+ else if(theDelta)
+ {
+ //There is a case, when
+ // const double aPar = 0.99999999999721167;
+ // const double aFI2 = 2.3593296083566181e-006;
+
+ //Then
+ // aPar - cos(aFI2) == -5.10703e-015 ==> cos(aFI2) == aPar.
+ //Theoreticaly, in this case
+ // aFI2 == +/- acos(aPar).
+ //However,
+ // acos(aPar) - aFI2 == 2.16362e-009.
+ //Error is quite big.
+
+ //This error should be estimated. Let use following way, which is based
+ //on expanding to Taylor series.
+
+ // acos(p)-acos(p+x) = x/sqrt(1-p*p).
+
+ //If p == (1-d) (when p > 0) or p == (-1+d) (when p < 0) then
+ // acos(p)-acos(p+x) = x/sqrt(d*(2-d)).
+
+ //Here always aTol0 <= d <= 1. Max(x) is considered (!) to be equal to aTol0.
+ //In this case
+ // 8*aTol0 <= acos(p)-acos(p+x) <= sqrt(2/(2-aTol0)-1),
+ // because 0 < aTol0 < 1.
+ //E.g. when aTol0 = 1.0e-11,
+ // 8.0e-11 <= acos(p)-acos(p+x) < 2.24e-6.
+
+ const Standard_Real aDelta = Min(1.0-anArg, 1.0+anArg);
+ Standard_DivideByZero_Raise_if((aDelta*aDelta < RealSmall()) || (aDelta >= 2.0),
+ "IntPatch_ImpImpIntersection_4.gxx, CylCylComputeParameters()");
+ *theDelta = aTol0/sqrt(aDelta*(2.0-aDelta));
+ }
- theU2 = theCoeffs.mFI2 + aSign*acos(anArg);
+ theU2 = acos(anArg);
+ theU2 = theCoeffs.mFI2 + aSign*theU2;
return Standard_True;
}
-static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1,
+//=======================================================================
+//function : CylCylComputeParameters
+//purpose : Computes V1 and V2 (V-parameters of the 1st and 2nd cylinder respectively).
+//=======================================================================
+Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1,
const Standard_Real theU2,
const stCoeffsValue& theCoeffs,
Standard_Real& theV1,
return Standard_True;
}
-
-static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
+//=======================================================================
+//function : CylCylComputeParameters
+//purpose : Computes U2 (U-parameter of the 2nd cylinder),
+// V1 and V2 (V-parameters of the 1st and 2nd cylinder respectively).
+//=======================================================================
+Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
const Standard_Integer theWLIndex,
const stCoeffsValue& theCoeffs,
Standard_Real& theU2,
//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
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();
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);
}
const Standard_Real thePeriod,
const Standard_Boolean theFlForce)
{
+ if(Precision::IsInfinite(theUGiven))
+ {
+ return Standard_False;
+ }
+
if((theUfTarget - theUGiven <= theTol2D) &&
(theUGiven - theUlTarget <= theTol2D))
{//it has already inscribed
return Standard_True;
}
- if(IsEqual(thePeriod, 0.0))
- {//it cannot be inscribed
- return Standard_False;
- }
-
- Standard_Real aDU = theUGiven - theUfTarget;
-
- if(aDU > 0.0)
- aDU = -thePeriod;
- else
- aDU = +thePeriod;
+ const Standard_Real aUf = theUfTarget - theTol2D;
+ const Standard_Real aUl = aUf + thePeriod;
- while(((theUGiven - theUfTarget)*aDU < 0.0) && !((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D)))
- {
- theUGiven += aDU;
- }
+ theUGiven = ElCLib::InPeriod(theUGiven, aUf, aUl);
- return ((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D));
+ return ((theUfTarget - theUGiven <= theTol2D) &&
+ (theUGiven - theUlTarget <= theTol2D));
}
//=======================================================================
//function : InscribeInterval
-//purpose : Adjusts theUfGiven and after that fits theUlGiven to result
+//purpose : Shifts theRange in order to make at least one of its
+// boundary in the range [theUfTarget, theUlTarget]
//=======================================================================
static Standard_Boolean InscribeInterval(const Standard_Real theUfTarget,
- const Standard_Real theUlTarget,
- Standard_Real& theUfGiven,
- Standard_Real& theUlGiven,
- const Standard_Real theTol2D,
- const Standard_Real thePeriod)
+ const Standard_Real theUlTarget,
+ Bnd_Range &theRange,
+ const Standard_Real theTol2D,
+ const Standard_Real thePeriod)
{
- Standard_Real anUpar = theUfGiven;
- const Standard_Real aDelta = theUlGiven - theUfGiven;
- if(InscribePoint(theUfTarget, theUlTarget, anUpar,
- theTol2D, thePeriod, Standard_False))
+ Standard_Real anUpar = 0.0;
+ if (!theRange.GetMin(anUpar))
+ {
+ return Standard_False;
+ }
+
+ const Standard_Real aDelta = theRange.Delta();
+ if(InscribePoint(theUfTarget, theUlTarget, anUpar,
+ theTol2D, thePeriod, (Abs(theUlTarget-anUpar) < theTol2D)))
{
- theUfGiven = anUpar;
- theUlGiven = theUfGiven + aDelta;
+ theRange.SetVoid();
+ theRange.Add(anUpar);
+ theRange.Add(anUpar + aDelta);
}
else
{
- anUpar = theUlGiven;
+ if (!theRange.GetMax (anUpar))
+ {
+ return Standard_False;
+ }
+
if(InscribePoint(theUfTarget, theUlTarget, anUpar,
- theTol2D, thePeriod, Standard_False))
+ theTol2D, thePeriod, (Abs(theUfTarget-anUpar) < theTol2D)))
{
- theUlGiven = anUpar;
- theUfGiven = theUlGiven - aDelta;
+ theRange.SetVoid();
+ theRange.Add(anUpar);
+ theRange.Add(anUpar - aDelta);
}
else
{
}
//=======================================================================
-//function : InscribeAndSortArray
-//purpose : Sort from Min to Max value
+//function : ExcludeNearElements
+//purpose : Checks if theArr contains two almost equal elements.
+// If it is true then one of equal elements will be excluded
+// (made infinite).
+// Returns TRUE if at least one element of theArr has been changed.
+//ATTENTION!!!
+// 1. Every not infinite element of theArr is considered to be
+// in [0, T] interval (where T is the period);
+// 2. theArr must be sorted in ascending order.
//=======================================================================
-static void InscribeAndSortArray( Standard_Real theArr[],
- const Standard_Integer theNOfMembers,
- const Standard_Real theUf,
- const Standard_Real theUl,
- const Standard_Real theTol2D,
- const Standard_Real thePeriod)
+static Standard_Boolean ExcludeNearElements(Standard_Real theArr[],
+ const Standard_Integer theNOfMembers,
+ const Standard_Real theUSurf1f,
+ const Standard_Real theUSurf1l,
+ const Standard_Real theTol)
{
- for(Standard_Integer i = 0; i < theNOfMembers; i++)
+ Standard_Boolean aRetVal = Standard_False;
+ for(Standard_Integer i = 1; i < theNOfMembers; i++)
{
- if(Precision::IsInfinite(theArr[i]))
- {
- if(theArr[i] < 0.0)
- theArr[i] = -theArr[i];
+ Standard_Real &anA = theArr[i],
+ &anB = theArr[i-1];
- continue;
- }
+ //Here, anA >= anB
- InscribePoint(theUf, theUl, theArr[i], theTol2D, thePeriod, Standard_False);
+ if(Precision::IsInfinite(anA))
+ break;
- for(Standard_Integer j = i - 1; j >= 0; j--)
+ if((anA-anB) < theTol)
{
+ if((anB != 0.0) && (anB != theUSurf1f) && (anB != theUSurf1l))
+ anA = (anA + anB)/2.0;
+ else
+ anA = anB;
- if(theArr[j+1] < theArr[j])
- {
- Standard_Real anUtemp = theArr[j+1];
- theArr[j+1] = theArr[j];
- theArr[j] = anUtemp;
- }
+ //Make this element infinite an forget it
+ //(we will not use it in next iterations).
+ anB = Precision::Infinite();
+ aRetVal = Standard_True;
}
}
-}
+ return aRetVal;
+}
//=======================================================================
//function : AddPointIntoWL
//purpose : Surf1 is a surface, whose U-par is variable.
+// If theFlBefore == TRUE then we enable the U1-parameter
+// of the added point to be less than U1-parameter of
+// previously added point (in general U1-parameter is
+// always increased; therefore, this situation is abnormal).
+// If theOnlyCheck==TRUE then no point will be added to theLine.
//=======================================================================
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,
const Standard_Integer theWLIndex,
const Standard_Real theTol3D,
const Standard_Real theTol2D,
- const Standard_Boolean theFlForce,
- const Standard_Boolean theFlBefore = Standard_False)
+ const Standard_Boolean theFlBefore = Standard_False,
+ const Standard_Boolean theOnlyCheck = Standard_False)
{
+ //Check if the point is in the domain or can be inscribed in the domain after adjusting.
+
gp_Pnt aPt1(theQuad1.Value(thePntOnSurf1.X(), thePntOnSurf1.Y())),
aPt2(theQuad2.Value(thePntOnSurf2.X(), thePntOnSurf2.Y()));
Standard_Real aU1par = thePntOnSurf1.X();
+
+ // aU1par always increases. Therefore, we must reduce its
+ // value in order to continue creation of WLine.
if(!InscribePoint(theUfSurf1, theUlSurf1, aU1par, theTol2D,
- thePeriodOfSurf1, theFlForce))
+ thePeriodOfSurf1, aU1par > 0.5*(theUfSurf1+theUlSurf1)))
return Standard_False;
+ if ((theLine->NbPoints() > 0) &&
+ ((theUlSurf1 - theUfSurf1) >= thePeriodOfSurf1) &&
+ (((aU1par + thePeriodOfSurf1 - theUlSurf1) <= theTol2D) ||
+ ((aU1par - thePeriodOfSurf1 - theUfSurf1) >= theTol2D)))
+ {
+ // aU1par can be adjusted to both theUlSurf1 and theUfSurf1
+ // with equal possibilities. This code fragment allows choosing
+ // correct parameter from these two variants.
+
+ Standard_Real aU1 = 0.0, aV1 = 0.0;
+ if (isTheReverse)
+ {
+ theLine->Value(theLine->NbPoints()).ParametersOnS2(aU1, aV1);
+ }
+ else
+ {
+ theLine->Value(theLine->NbPoints()).ParametersOnS1(aU1, aV1);
+ }
+
+ const Standard_Real aDelta = aU1 - aU1par;
+ if (2.0*Abs(aDelta) > thePeriodOfSurf1)
+ {
+ aU1par += Sign(thePeriodOfSurf1, aDelta);
+ }
+ }
+
Standard_Real aU2par = thePntOnSurf2.X();
if(!InscribePoint(theUfSurf2, theUlSurf2, aU2par, theTol2D,
thePeriodOfSurf1, Standard_False))
if((aV2par - theVlSurf2 > theTol2D) || (theVfSurf2 - aV2par > theTol2D))
return Standard_False;
+ //Get intersection point and add it in the WL
IntSurf_PntOn2S aPnt;
if(isTheReverse)
aPlast.ParametersOnS1(aUl, aVl);
if(!theFlBefore && (aU1par <= aUl))
- {//Parameter value must be increased if theFlBefore == FALSE.
- return Standard_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;
+ }
}
+ if (theOnlyCheck)
+ return Standard_True;
+
//theTol2D is minimal step along parameter changed.
//Therefore, if we apply this minimal step two
//neighbour points will be always "same". Consequently,
}
}
+ if (theOnlyCheck)
+ return Standard_True;
+
theLine->Add(aPnt);
if(!isThePrecise)
if((1 < aDeltaStep) && (aDeltaStep < 2000))
{
+ //Add new points in case of non-uniform distribution of existing points
SeekAdditionalPoints( theQuad1, theQuad2, theLine,
theCoeffs, theWLIndex, aDeltaStep, aNbPnts-2,
- aNbPnts-1, theUfSurf2, theUlSurf2,
- theTol2D, thePeriodOfSurf1, isTheReverse);
+ aNbPnts-1, theTol2D, thePeriodOfSurf1, isTheReverse);
}
}
//=======================================================================
//function : AddBoundaryPoint
-//purpose :
+//purpose : Find intersection point on V-boundary.
//=======================================================================
-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 theU1Min,
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 isTheReverse,
const Standard_Boolean theFlForce,
Standard_Boolean& isTheFound1,
- Standard_Boolean& isTheFound2)
+ Standard_Boolean& isTheFound2) const
{
Standard_Real aUSurf1f = 0.0, //const
aUSurf1l = 0.0,
aVSurf2f = 0.0,
aVSurf2l = 0.0;
- theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
- theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
-
- SearchBoundType aTS1 = SearchNONE, aTS2 = SearchNONE;
- Standard_Real aV1zad = aVSurf1f, aV2zad = aVSurf2f;
-
- Standard_Real anUpar1 = theU1, anUpar2 = theU1;
- Standard_Real aVf = theV1, aVl = theV1Prev;
-
- 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))
- {
- aTS1 = SearchV1;
- aV1zad = aVSurf1l;
- isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1l, theV2, theU2, theU1, anUpar1);
- }
+ myUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
+ myUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
+
+ const Standard_Integer aSize = 4;
+ const Standard_Real anArrVzad[aSize] = {aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l};
- aVf = theV2;
- aVl = theV2Prev;
+ StPInfo aUVPoint[aSize];
- if((Abs(aVf-aVSurf2f) <= theTol2D) ||
- ((aVf-aVSurf2f)*(aVl-aVSurf2f) <= 0.0))
+ for(Standard_Integer anIDSurf = 0; anIDSurf < 4; anIDSurf+=2)
{
- 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))
- {
- aTS2 = SearchV2;
- aV2zad = aVSurf2l;
- isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2l, theV1, theU2, theU1, anUpar2);
- }
+ const Standard_Real aVf = (anIDSurf == 0) ? theV1 : theV2,
+ aVl = (anIDSurf == 0) ? theV1Prev : theV2Prev;
- if(!isTheFound1 && !isTheFound2)
- return Standard_True;
+ const SearchBoundType aTS = (anIDSurf == 0) ? SearchV1 : SearchV2;
- if(anUpar1 < anUpar2)
- {
- if(isTheFound1)
+ for(Standard_Integer anIDBound = 0; anIDBound < 2; anIDBound++)
{
- Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
- if(!CylCylComputeParameters(anUpar1, theWLIndex, theCoeffs, aU2, aV1, aV2))
- {
- isTheFound1 = Standard_False;
- }
-
- 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))
- {
- isTheFound1 = Standard_False;
- }
- }
+ const Standard_Integer anIndex = anIDSurf+anIDBound;
- if(isTheFound2)
- {
- Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
- if(!CylCylComputeParameters(anUpar2, theWLIndex, theCoeffs, aU2, aV1, aV2))
+ aUVPoint[anIndex].mySurfID = anIDSurf;
+
+ if((Abs(aVf-anArrVzad[anIndex]) > myTol2D) &&
+ ((aVf-anArrVzad[anIndex])*(aVl-anArrVzad[anIndex]) > 0.0))
{
- isTheFound2 = Standard_False;
+ continue;
}
- if(aTS2 == SearchV1)
- aV1 = aV1zad;
+ //Segment [aVf, aVl] intersects at least one V-boundary (first or last)
+ // (in general, case is possible, when aVf > aVl).
- if(aTS2 == SearchV2)
- aV2 = aV2zad;
+ // Precise intersection point
+ const Standard_Boolean aRes = SearchOnVBounds(aTS, anArrVzad[anIndex],
+ (anIDSurf == 0) ? theV2 : theV1,
+ theU2, theU1,
+ aUVPoint[anIndex].myU1);
- 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))
+ // aUVPoint[anIndex].myU1 is considered to be nearer to theU1 than
+ // to theU1+/-Period
+ if (!aRes || (aUVPoint[anIndex].myU1 >= theU1) ||
+ (aUVPoint[anIndex].myU1 < theU1Min))
{
- isTheFound2 = Standard_False;
+ //Intersection point is not found or out of the domain
+ aUVPoint[anIndex].myU1 = RealLast();
+ continue;
}
- }
- }
- else
- {
- if(isTheFound2)
- {
- Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
- if(!CylCylComputeParameters(anUpar2, theWLIndex, theCoeffs, aU2, aV1, aV2))
+ else
{
- isTheFound2 = Standard_False;
- }
-
- if(aTS2 == SearchV1)
- aV1 = aV1zad;
-
- if(aTS2 == SearchV2)
- aV2 = aV2zad;
+ //intersection point is found
+
+ 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(!ComputationMethods::
+ CylCylComputeParameters(aU1, theWLIndex, myCoeffs, aU2, aV1, aV2))
+ {
+ // Found point is wrong
+ aU1 = RealLast();
+ continue;
+ }
- 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;
+ //Point on true V-boundary.
+ if(aTS == SearchV1)
+ aV1 = anArrVzad[anIndex];
+ else //if(aTS[anIndex] == SearchV2)
+ aV2 = anArrVzad[anIndex];
}
}
+ }
- if(isTheFound1)
- {
- Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
- if(!CylCylComputeParameters(anUpar1, theWLIndex, theCoeffs, aU2, aV1, aV2))
- {
- isTheFound1 = Standard_False;
- }
+ // Sort with acceding U1-parameter.
+ std::sort(aUVPoint, aUVPoint+aSize);
+
+ isTheFound1 = isTheFound2 = Standard_False;
- if(aTS1 == SearchV1)
- aV1 = aV1zad;
+ //Adding found points on boundary in the WLine.
+ for(Standard_Integer i = 0; i < aSize; i++)
+ {
+ if(aUVPoint[i].myU1 == RealLast())
+ break;
- if(aTS1 == SearchV2)
- aV2 = aV2zad;
+ 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))
+ {
+ continue;
+ }
- 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))
- {
- isTheFound1 = Standard_False;
- }
+ if(aUVPoint[i].mySurfID == 0)
+ {
+ isTheFound1 = Standard_True;
+ }
+ else
+ {
+ isTheFound2 = Standard_True;
}
}
-
- return Standard_True;
}
//=======================================================================
//function : SeekAdditionalPoints
-//purpose :
+//purpose : Inserts additional intersection points between neighbor points.
+// This action is bone with several iterations. During every iteration,
+// new point is inserted in middle of every interval.
+// The process will be finished if NbPoints >= theMinNbPoints.
//=======================================================================
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)
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;
}
//=======================================================================
-//function : CriticalPointsComputing
-//purpose :
+//function : BoundariesComputing
+//purpose : Computes true domain of future intersection curve (allows
+// avoiding excess iterations)
//=======================================================================
-static void CriticalPointsComputing(const stCoeffsValue& theCoeffs,
- const Standard_Real theUSurf1f,
- const Standard_Real theUSurf1l,
- const Standard_Real theUSurf2f,
- const Standard_Real theUSurf2l,
- const Standard_Real thePeriod,
- const Standard_Real theTol2D,
- const Standard_Integer theNbCritPointsMax,
- Standard_Real theU1crit[])
+Standard_Boolean WorkWithBoundaries::
+ BoundariesComputing(const ComputationMethods::stCoeffsValue &theCoeffs,
+ const Standard_Real thePeriod,
+ Bnd_Range theURange[])
{
- theU1crit[0] = 0.0;
- theU1crit[1] = thePeriod;
- theU1crit[2] = theUSurf1f;
- theU1crit[3] = theUSurf1l;
+ //All comments to this method is related to the comment
+ //to ComputationMethods class
+
+ //So, we have the equation
+ // cos(U2-FI2)=B*cos(U1-FI1)+C
+ //Evidently,
+ // -1 <= B*cos(U1-FI1)+C <= 1
- const Standard_Real aCOS = cos(theCoeffs.mFI2);
- const Standard_Real aBSB = Abs(theCoeffs.mB);
- if((theCoeffs.mC - aBSB <= aCOS) && (aCOS <= theCoeffs.mC + aBSB))
+ if (theCoeffs.mB > 0.0)
{
- Standard_Real anArg = (aCOS - theCoeffs.mC) / theCoeffs.mB;
- if(anArg > 1.0)
- anArg = 1.0;
- if(anArg < -1.0)
- anArg = -1.0;
-
- theU1crit[4] = -acos(anArg) + theCoeffs.mFI1;
- theU1crit[5] = acos(anArg) + theCoeffs.mFI1;
- }
-
- Standard_Real aSf = cos(theUSurf2f - theCoeffs.mFI2);
- Standard_Real aSl = cos(theUSurf2l - theCoeffs.mFI2);
- MinMax(aSf, aSl);
+ // -(1+C)/B <= cos(U1-FI1) <= (1-C)/B
- theU1crit[6] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? -acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : -Precision::Infinite();
- theU1crit[7] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? -acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : Precision::Infinite();
- theU1crit[8] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : -Precision::Infinite();
- theU1crit[9] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : Precision::Infinite();
-
- //preparative treatment of array
- InscribeAndSortArray(theU1crit, theNbCritPointsMax, 0.0, thePeriod, theTol2D, thePeriod);
- for(Standard_Integer i = 1; i < theNbCritPointsMax; i++)
- {
- Standard_Real &a = theU1crit[i],
- &b = theU1crit[i-1];
- const Standard_Real aRemain = fmod(Abs(a - b), thePeriod); // >= 0, because Abs(a - b) >= 0
- if((Abs(a - b) < theTol2D) || (aRemain < theTol2D) || (Abs(aRemain - thePeriod) < theTol2D))
+ if (theCoeffs.mB + Abs(theCoeffs.mC) < -1.0)
{
- a = (a + b)/2.0;
- b = Precision::Infinite();
- }
- }
-}
-
-//=======================================================================
-//function : IntCyCyTrim
-//purpose :
-//=======================================================================
-Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
- const IntSurf_Quadric& theQuad2,
- const Standard_Real theTol3D,
- const Standard_Real theTol2D,
- const Bnd_Box2d& theUVSurf1,
- const Bnd_Box2d& theUVSurf2,
- const Standard_Boolean isTheReverse,
- Standard_Boolean& isTheEmpty,
- 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);
-
- const gp_Cylinder& aCyl1 = theQuad1.Cylinder(),
- aCyl2 = theQuad2.Cylinder();
-
- IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
-
- if (!anInter.IsDone())
- {
- return Standard_False;
- }
-
- IntAna_ResultType aTypInt = anInter.TypeInter();
-
- if(aTypInt != IntAna_NoGeometricSolution)
- { //It is not necessary (because result is an analytic curve) or
- //it is impossible to make Walking-line.
-
- return Standard_False;
- }
-
- theSlin.Clear();
- theSPnt.Clear();
- const Standard_Integer aNbMaxPoints = 2000;
- const Standard_Integer aNbMinPoints = 200;
- const Standard_Integer aNbPoints = Min(Max(aNbMinPoints,
- RealToInt(20.0*aCyl1.Radius())), aNbMaxPoints);
- const Standard_Real aPeriod = 2.0*M_PI;
- const Standard_Real aStepMin = theTol2D,
- aStepMax = (aUSurf1l-aUSurf1f)/IntToReal(aNbPoints);
- const Standard_Integer aNbWLines = 2;
-
- const stCoeffsValue anEquationCoeffs(aCyl1, aCyl2);
-
- //Boundaries
- const Standard_Integer aNbOfBoundaries = 2;
- Standard_Real aU1f[aNbOfBoundaries] = {-Precision::Infinite(), -Precision::Infinite()};
- Standard_Real aU1l[aNbOfBoundaries] = {Precision::Infinite(), Precision::Infinite()};
-
- if(anEquationCoeffs.mB > 0.0)
- {
- if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) < -1.0)
- {//There is NOT intersection
- return Standard_True;
+ //(1-C)/B < -1 or -(1+C)/B > 1 ==> No solution
+
+ return Standard_False;
}
- else if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
- {//U=[0;2*PI]+aFI1
- aU1f[0] = anEquationCoeffs.mFI1;
- aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
+ else if (theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
+ {
+ //(1-C)/B >= 1 and -(1+C)/B <= -1 ==> U=[0;2*PI]+aFI1
+ theURange[0].Add(theCoeffs.mFI1);
+ theURange[0].Add(thePeriod + theCoeffs.mFI1);
}
- else if((1 + anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
- (anEquationCoeffs.mB <= 1 - anEquationCoeffs.mC))
+ else if ((1 + theCoeffs.mC <= theCoeffs.mB) &&
+ (theCoeffs.mB <= 1 - theCoeffs.mC))
{
- Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
+ //(1-C)/B >= 1 and -(1+C)/B >= -1 ==>
+ //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1),
+ //where aDAngle = acos(-(myCoeffs.mC + 1) / myCoeffs.mB)
+
+ Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
if(anArg > 1.0)
anArg = 1.0;
if(anArg < -1.0)
anArg = -1.0;
const Standard_Real aDAngle = acos(anArg);
- //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
- aU1f[0] = anEquationCoeffs.mFI1;
- aU1l[0] = aDAngle + anEquationCoeffs.mFI1;
- aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
- aU1l[1] = aPeriod + anEquationCoeffs.mFI1;
+ theURange[0].Add(theCoeffs.mFI1);
+ theURange[0].Add(aDAngle + theCoeffs.mFI1);
+ theURange[1].Add(thePeriod - aDAngle + theCoeffs.mFI1);
+ theURange[1].Add(thePeriod + theCoeffs.mFI1);
}
- else if((1 - anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
- (anEquationCoeffs.mB <= 1 + anEquationCoeffs.mC))
+ else if ((1 - theCoeffs.mC <= theCoeffs.mB) &&
+ (theCoeffs.mB <= 1 + theCoeffs.mC))
{
- Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
+ //(1-C)/B <= 1 and -(1+C)/B <= -1 ==> U=[aDAngle;2*PI-aDAngle]+aFI1
+ //where aDAngle = acos((1 - myCoeffs.mC) / myCoeffs.mB)
+
+ Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
if(anArg > 1.0)
anArg = 1.0;
if(anArg < -1.0)
anArg = -1.0;
const Standard_Real aDAngle = acos(anArg);
- //U=[aDAngle;2*PI-aDAngle]+aFI1
-
- aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
- aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
+ theURange[0].Add(aDAngle + theCoeffs.mFI1);
+ theURange[0].Add(thePeriod - aDAngle + theCoeffs.mFI1);
}
- else if(anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
+ else if (theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
{
- Standard_Real anArg1 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB,
- anArg2 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
+ //(1-C)/B <= 1 and -(1+C)/B >= -1 ==>
+ //(U=[aDAngle1;aDAngle2]+aFI1) ||
+ //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
+ //where aDAngle1 = acos((1 - myCoeffs.mC) / myCoeffs.mB),
+ // aDAngle2 = acos(-(myCoeffs.mC + 1) / myCoeffs.mB).
+
+ Standard_Real anArg1 = (1 - theCoeffs.mC) / theCoeffs.mB,
+ anArg2 = -(theCoeffs.mC + 1) / theCoeffs.mB;
if(anArg1 > 1.0)
anArg1 = 1.0;
if(anArg1 < -1.0)
const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
//(U=[aDAngle1;aDAngle2]+aFI1) ||
//(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
-
- aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1;
- aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1;
- aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1;
- aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1;
+ theURange[0].Add(aDAngle1 + theCoeffs.mFI1);
+ theURange[0].Add(aDAngle2 + theCoeffs.mFI1);
+ theURange[1].Add(thePeriod - aDAngle2 + theCoeffs.mFI1);
+ theURange[1].Add(thePeriod - aDAngle1 + theCoeffs.mFI1);
}
else
{
- Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
+ return Standard_False;
}
}
- else if(anEquationCoeffs.mB < 0.0)
+ else if (theCoeffs.mB < 0.0)
{
- if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) > 1.0)
- {//There is NOT intersection
- return Standard_True;
+ // (1-C)/B <= cos(U1-FI1) <= -(1+C)/B
+
+ if (theCoeffs.mB + Abs(theCoeffs.mC) > 1.0)
+ {
+ // -(1+C)/B < -1 or (1-C)/B > 1 ==> No solutions
+ return Standard_False;
}
- else if(-anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
- {//U=[0;2*PI]+aFI1
- aU1f[0] = anEquationCoeffs.mFI1;
- aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
+ else if (-theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
+ {
+ // -(1+C)/B >= 1 and (1-C)/B <= -1 ==> U=[0;2*PI]+aFI1
+ theURange[0].Add(theCoeffs.mFI1);
+ theURange[0].Add(thePeriod + theCoeffs.mFI1);
}
- else if((-anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
- ( anEquationCoeffs.mB <= anEquationCoeffs.mC - 1))
+ else if ((-theCoeffs.mC - 1 <= theCoeffs.mB) &&
+ (theCoeffs.mB <= theCoeffs.mC - 1))
{
- Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
+ // -(1+C)/B >= 1 and (1-C)/B >= -1 ==>
+ //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1),
+ //where aDAngle = acos((1 - myCoeffs.mC) / myCoeffs.mB)
+
+ Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
if(anArg > 1.0)
anArg = 1.0;
if(anArg < -1.0)
anArg = -1.0;
const Standard_Real aDAngle = acos(anArg);
- //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
-
- aU1f[0] = anEquationCoeffs.mFI1;
- aU1l[0] = aDAngle + anEquationCoeffs.mFI1;
- aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
- aU1l[1] = aPeriod + anEquationCoeffs.mFI1;
+ theURange[0].Add(theCoeffs.mFI1);
+ theURange[0].Add(aDAngle + theCoeffs.mFI1);
+ theURange[1].Add(thePeriod - aDAngle + theCoeffs.mFI1);
+ theURange[1].Add(thePeriod + theCoeffs.mFI1);
}
- else if((anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
- (anEquationCoeffs.mB <= -anEquationCoeffs.mB - 1))
+ else if ((theCoeffs.mC - 1 <= theCoeffs.mB) &&
+ (theCoeffs.mB <= -theCoeffs.mB - 1))
{
- Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
+ // -(1+C)/B <= 1 and (1-C)/B <= -1 ==> U=[aDAngle;2*PI-aDAngle]+aFI1,
+ //where aDAngle = acos(-(myCoeffs.mC + 1) / myCoeffs.mB).
+
+ Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
if(anArg > 1.0)
anArg = 1.0;
if(anArg < -1.0)
anArg = -1.0;
const Standard_Real aDAngle = acos(anArg);
- //U=[aDAngle;2*PI-aDAngle]+aFI1
-
- aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
- aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
+ theURange[0].Add(aDAngle + theCoeffs.mFI1);
+ theURange[0].Add(thePeriod - aDAngle + theCoeffs.mFI1);
}
- else if(-anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
+ else if (-theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
{
- Standard_Real anArg1 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB,
- anArg2 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
+ // -(1+C)/B <= 1 and (1-C)/B >= -1 ==>
+ //(U=[aDAngle1;aDAngle2]+aFI1) || (U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1),
+ //where aDAngle1 = acos(-(myCoeffs.mC + 1) / myCoeffs.mB),
+ // aDAngle2 = acos((1 - myCoeffs.mC) / myCoeffs.mB)
+
+ Standard_Real anArg1 = -(theCoeffs.mC + 1) / theCoeffs.mB,
+ anArg2 = (1 - theCoeffs.mC) / theCoeffs.mB;
if(anArg1 > 1.0)
anArg1 = 1.0;
if(anArg1 < -1.0)
anArg2 = -1.0;
const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
- //(U=[aDAngle1;aDAngle2]+aFI1) ||
- //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
-
- aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1;
- aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1;
- aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1;
- aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1;
+ theURange[0].Add(aDAngle1 + theCoeffs.mFI1);
+ theURange[0].Add(aDAngle2 + theCoeffs.mFI1);
+ theURange[1].Add(thePeriod - aDAngle2 + theCoeffs.mFI1);
+ theURange[1].Add(thePeriod - aDAngle1 + theCoeffs.mFI1);
}
else
{
- Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
+ return Standard_False;
}
}
else
{
- Standard_Failure::Raise("Error. Exception. Unhandled case (B-parameter computation)!");
+ return Standard_False;
+ }
+
+ return Standard_True;
+}
+
+//=======================================================================
+//function : CriticalPointsComputing
+//purpose : theNbCritPointsMax contains true number of critical points.
+// It must be initialized correctly before function calling
+//=======================================================================
+static void CriticalPointsComputing(const ComputationMethods::stCoeffsValue& theCoeffs,
+ const Standard_Real theUSurf1f,
+ const Standard_Real theUSurf1l,
+ const Standard_Real theUSurf2f,
+ const Standard_Real theUSurf2l,
+ const Standard_Real thePeriod,
+ const Standard_Real theTol2D,
+ Standard_Integer& theNbCritPointsMax,
+ Standard_Real theU1crit[])
+{
+ //[0...1] - in these points parameter U1 goes through
+ // the seam-edge of the first cylinder.
+ //[2...3] - First and last U1 parameter.
+ //[4...5] - in these points parameter U2 goes through
+ // the seam-edge of the second cylinder.
+ //[6...9] - in these points an intersection line goes through
+ // U-boundaries of the second surface.
+ //[10...11] - Boundary of monotonicity interval of U2(U1) function
+ // (see CylCylMonotonicity() function)
+
+ theU1crit[0] = 0.0;
+ theU1crit[1] = thePeriod;
+ theU1crit[2] = theUSurf1f;
+ theU1crit[3] = theUSurf1l;
+
+ const Standard_Real aCOS = cos(theCoeffs.mFI2);
+ const Standard_Real aBSB = Abs(theCoeffs.mB);
+ if((theCoeffs.mC - aBSB <= aCOS) && (aCOS <= theCoeffs.mC + aBSB))
+ {
+ Standard_Real anArg = (aCOS - theCoeffs.mC) / theCoeffs.mB;
+ if(anArg > 1.0)
+ anArg = 1.0;
+ if(anArg < -1.0)
+ anArg = -1.0;
+
+ theU1crit[4] = -acos(anArg) + theCoeffs.mFI1;
+ theU1crit[5] = acos(anArg) + theCoeffs.mFI1;
}
- for(Standard_Integer i = 0; i < aNbOfBoundaries; i++)
+ Standard_Real aSf = cos(theUSurf2f - theCoeffs.mFI2);
+ Standard_Real aSl = cos(theUSurf2l - theCoeffs.mFI2);
+ MinMax(aSf, aSl);
+
+ //In accorance with pure mathematic, theU1crit[6] and [8]
+ //must be -Precision::Infinite() instead of used +Precision::Infinite()
+ theU1crit[6] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
+ -acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
+ Precision::Infinite();
+ theU1crit[7] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
+ -acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
+ Precision::Infinite();
+ theU1crit[8] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
+ acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
+ Precision::Infinite();
+ theU1crit[9] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
+ acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
+ Precision::Infinite();
+
+ theU1crit[10] = theCoeffs.mFI1;
+ theU1crit[11] = M_PI+theCoeffs.mFI1;
+
+ //preparative treatment of array. This array must have faled to contain negative
+ //infinity number
+
+ for(Standard_Integer i = 0; i < theNbCritPointsMax; i++)
{
- if(Precision::IsInfinite(aU1f[i]) && Precision::IsInfinite(aU1l[i]))
+ if(Precision::IsInfinite(theU1crit[i]))
+ {
continue;
+ }
+
+ theU1crit[i] = fmod(theU1crit[i], thePeriod);
+ if(theU1crit[i] < 0.0)
+ theU1crit[i] += thePeriod;
+ }
+
+ //Here all not infinite elements of theU1crit are in [0, thePeriod) range
- InscribeInterval(aUSurf1f, aUSurf1l, aU1f[i], aU1l[i], theTol2D, aPeriod);
+ do
+ {
+ std::sort(theU1crit, theU1crit + theNbCritPointsMax);
}
+ while(ExcludeNearElements(theU1crit, theNbCritPointsMax,
+ theUSurf1f, theUSurf1l, theTol2D));
- if( !Precision::IsInfinite(aU1f[0]) && !Precision::IsInfinite(aU1f[1]) &&
- !Precision::IsInfinite(aU1l[0]) && !Precision::IsInfinite(aU1l[1]))
+ //Here all not infinite elements in theU1crit are different and sorted.
+
+ while(theNbCritPointsMax > 0)
{
- if( ((aU1f[1] <= aU1l[0]) || (aU1l[1] <= aU1l[0])) &&
- ((aU1f[0] <= aU1l[1]) || (aU1l[0] <= aU1l[1])))
- {//Join all intervals to one
- aU1f[0] = Min(aU1f[0], aU1f[1]);
- aU1l[0] = Max(aU1l[0], aU1l[1]);
+ Standard_Real &anB = theU1crit[theNbCritPointsMax-1];
+ if(Precision::IsInfinite(anB))
+ {
+ theNbCritPointsMax--;
+ continue;
+ }
+
+ //1st not infinte element is found
+
+ if(theNbCritPointsMax == 1)
+ break;
+
+ //Here theNbCritPointsMax > 1
- aU1f[1] = -Precision::Infinite();
- aU1l[1] = Precision::Infinite();
+ Standard_Real &anA = theU1crit[0];
+
+ //Compare 1st and last significant elements of theU1crit
+ //They may still differs by period.
+
+ if (Abs(anB - anA - thePeriod) < theTol2D)
+ {//E.g. anA == 2.0e-17, anB == (thePeriod-1.0e-18)
+ anA = (anA + anB - thePeriod)/2.0;
+ anB = Precision::Infinite();
+ theNbCritPointsMax--;
}
+
+ //Out of "while(theNbCritPointsMax > 0)" cycle.
+ break;
}
- //Critical points
- //[0...1] - in these points parameter U1 goes through
- // the seam-edge of the first cylinder.
- //[2...3] - First and last U1 parameter.
- //[4...5] - in these points parameter U2 goes through
- // the seam-edge of the second cylinder.
- //[6...9] - in these points an intersection line goes through
- // U-boundaries of the second surface.
- const Standard_Integer aNbCritPointsMax = 10;
- Standard_Real anU1crit[aNbCritPointsMax] = {Precision::Infinite(),
- Precision::Infinite(),
- Precision::Infinite(),
- Precision::Infinite(),
- Precision::Infinite(),
- Precision::Infinite(),
- Precision::Infinite(),
- Precision::Infinite(),
- Precision::Infinite(),
- Precision::Infinite()};
+ //Attention! Here theU1crit may be unsorted.
+}
- CriticalPointsComputing(anEquationCoeffs, aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
- aPeriod, theTol2D, aNbCritPointsMax, anU1crit);
+//=======================================================================
+//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 = aD1.XYZ().CrossSquareMagnitude(aD2.XYZ());
+
+ //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 anAbsCosA = Abs(aCosA);
+ const Standard_Real aHDV1 = (aR1 * anAbsCosA + aR2) / aSinA,
+ aHDV2 = (aR2 * anAbsCosA + 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 : CyCyNoGeometric
+//purpose :
+//=======================================================================
+static IntPatch_ImpImpIntersection::IntStatus
+ CyCyNoGeometric(const gp_Cylinder &theCyl1,
+ const gp_Cylinder &theCyl2,
+ const WorkWithBoundaries &theBW,
+ Bnd_Range theRange[],
+ const Standard_Integer theNbOfRanges /*=2*/,
+ Standard_Boolean& isTheEmpty,
+ IntPatch_SequenceOfLine& theSlin,
+ IntPatch_SequenceOfPoint& theSPnt)
+{
+ Standard_Real aUSurf1f = 0.0, aUSurf1l = 0.0,
+ aUSurf2f = 0.0, aUSurf2l = 0.0,
+ aVSurf1f = 0.0, aVSurf1l = 0.0,
+ aVSurf2f = 0.0, aVSurf2l = 0.0;
+
+ theBW.UVS1().Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
+ theBW.UVS2().Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
+
+ Bnd_Range aRangeS1, aRangeS2;
+ theBW.BoundaryEstimation(theCyl1, theCyl2, aRangeS1, aRangeS2);
+ if (aRangeS1.IsVoid() || aRangeS2.IsVoid())
+ return IntPatch_ImpImpIntersection::IntStatus_OK;
+
+ {
+ //Quotation of the message from issue #26894 (author MSV):
+ //"We should return fail status from intersector if the result should be an
+ //infinite curve of non-analytical type... I propose to define the limit for the
+ //extent as the radius divided by 1e+2 and multiplied by 1e+7.
+ //Thus, taking into account the number of valuable digits (15), we provide reliable
+ //computations with an error not exceeding R/100."
+ const Standard_Real aF = 1.0e+5;
+ const Standard_Real aMaxV1Range = aF*theCyl1.Radius(), aMaxV2Range = aF*theCyl2.Radius();
+ if ((aRangeS1.Delta() > aMaxV1Range) || (aRangeS2.Delta() > aMaxV2Range))
+ return IntPatch_ImpImpIntersection::IntStatus_InfiniteSectionCurve;
+ }
+
+ const ComputationMethods::stCoeffsValue &anEquationCoeffs = theBW.SICoeffs();
+ const IntSurf_Quadric& aQuad1 = theBW.GetQSurface(1);
+ const IntSurf_Quadric& aQuad2 = theBW.GetQSurface(2);
+ const Standard_Boolean isReversed = theBW.IsReversed();
+ const Standard_Real aTol2D = theBW.Get2dTolerance();
+ const Standard_Real aTol3D = theBW.Get3dTolerance();
+ const Standard_Real aPeriod = 2.0*M_PI;
+ const Standard_Integer aNbMaxPoints = 2000;
+ const Standard_Integer aNbMinPoints = 200;
+ const Standard_Integer aNbPoints = Min(Max(aNbMinPoints,
+ RealToInt(20.0*theCyl1.Radius())), aNbMaxPoints);
+ const Standard_Real aStepMin = aTol2D,
+ aStepMax = (aUSurf1l - aUSurf1f > M_PI / 100.0) ?
+ (aUSurf1l - aUSurf1f) / IntToReal(aNbPoints) :
+ aUSurf1l - aUSurf1f;
+
+ //The main idea of the algorithm is to change U1-parameter
+ //(U-parameter of theCyl1) from aU1f to aU1l with some step
+ //(step is adaptive) and to obtain set of intersection points.
+
+ for (Standard_Integer i = 0; i < theNbOfRanges; i++)
+ {
+ if (theRange[i].IsVoid())
+ continue;
+
+ InscribeInterval(aUSurf1f, aUSurf1l, theRange[i], aTol2D, aPeriod);
+ }
+
+ if (theRange[0].Union(theRange[1]))
+ {
+ // Works only if (theNbOfRanges == 2).
+ theRange[1].SetVoid();
+ }
+
+ //Critical points are the value of U1-parameter in the points
+ //where WL must be decomposed.
+
+ //When U1 goes through critical points its value is set up to this
+ //parameter forcefully and the intersection point is added in the line.
+ //After that, the WL is broken (next U1 value will be correspond to the new WL).
+
+ //See CriticalPointsComputing(...) function to get detail information about this array.
+ const Standard_Integer aNbCritPointsMax = 12;
+ Standard_Real anU1crit[aNbCritPointsMax] = { Precision::Infinite(),
+ Precision::Infinite(),
+ Precision::Infinite(),
+ Precision::Infinite(),
+ Precision::Infinite(),
+ Precision::Infinite(),
+ Precision::Infinite(),
+ Precision::Infinite(),
+ Precision::Infinite(),
+ Precision::Infinite(),
+ Precision::Infinite(),
+ Precision::Infinite() };
+
+ //This list of critical points is not full because it does not contain any points
+ //which intersection line goes through V-bounds of cylinders in.
+ //They are computed by numerical methods on - line (during algorithm working).
+ //The moment is caught, when intersection line goes through V-bounds of any cylinder.
+
+ Standard_Integer aNbCritPoints = aNbCritPointsMax;
+ CriticalPointsComputing(anEquationCoeffs, aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
+ aPeriod, aTol2D, aNbCritPoints, anU1crit);
//Getting Walking-line
enum WLFStatus
{
+ // No points have been added in WL
WLFStatus_Absent = 0,
- WLFStatus_Exist = 1,
+ // WL contains at least one point
+ WLFStatus_Exist = 1,
+ // WL has been finished in some critical point
+ // We should start new line
WLFStatus_Broken = 2
};
- for(Standard_Integer aCurInterval = 0; aCurInterval < aNbOfBoundaries; aCurInterval++)
+ const Standard_Integer aNbWLines = 2;
+ for (Standard_Integer aCurInterval = 0; aCurInterval < theNbOfRanges; aCurInterval++)
{
- if(Precision::IsInfinite(aU1f[aCurInterval]) && Precision::IsInfinite(aU1l[aCurInterval]))
- continue;
-
+ //Process every continuous region
Standard_Boolean isAddedIntoWL[aNbWLines];
- for(Standard_Integer i = 0; i < aNbWLines; i++)
+ for (Standard_Integer i = 0; i < aNbWLines; i++)
isAddedIntoWL[i] = Standard_False;
- Standard_Real anUf = aU1f[aCurInterval], anUl = aU1l[aCurInterval];
- const Standard_Boolean isDeltaPeriod = IsEqual(anUl-anUf, aPeriod);
+ Standard_Real anUf = 1.0, anUl = 0.0;
+ if (!theRange[aCurInterval].GetBounds(anUf, anUl))
+ continue;
+
+ const Standard_Boolean isDeltaPeriod = IsEqual(anUl - anUf, aPeriod);
//Inscribe and sort critical points
- InscribeAndSortArray(anU1crit, aNbCritPointsMax, anUf, anUl, theTol2D, aPeriod);
+ for (Standard_Integer i = 0; i < aNbCritPoints; i++)
+ {
+ InscribePoint(anUf, anUl, anU1crit[i], 0.0, aPeriod, Standard_False);
+ }
- while(anUf < anUl)
+ std::sort(anU1crit, anU1crit + aNbCritPoints);
+
+ while (anUf < anUl)
{
+ //Change value of U-parameter on the 1st surface from anUf to anUl
+ //(anUf will be modified in the cycle body).
+ //Step is computed adaptively (see comments below).
+
Standard_Real aU2[aNbWLines], aV1[aNbWLines], aV2[aNbWLines];
WLFStatus aWLFindStatus[aNbWLines];
Standard_Real aV1Prev[aNbWLines], aV2Prev[aNbWLines];
Handle(IntSurf_LineOn2S) aL2S[aNbWLines];
Handle(IntPatch_WLine) aWLine[aNbWLines];
- for(Standard_Integer i = 0; i < aNbWLines; i++)
+ for (Standard_Integer i = 0; i < aNbWLines; i++)
{
aL2S[i] = new IntSurf_LineOn2S();
aWLine[i] = new IntPatch_WLine(aL2S[i], Standard_False);
+ aWLine[i]->SetCreatingWayInfo(IntPatch_WLine::IntPatch_WLImpImp);
aWLFindStatus[i] = WLFStatus_Absent;
isAddingWLEnabled[i] = Standard_True;
aU2[i] = aV1[i] = aV2[i] = 0.0;
aV1Prev[i] = aV2Prev[i] = 0.0;
anUexpect[i] = anUf;
}
-
- Standard_Real anU1 = anUf;
- Standard_Real aCriticalDelta[aNbCritPointsMax];
- for(Standard_Integer i = 0; i < aNbCritPointsMax; i++)
- aCriticalDelta[i] = anU1 - anU1crit[i];
+ Standard_Real aCriticalDelta[aNbCritPointsMax] = { 0 };
+ for (Standard_Integer aCritPID = 0; aCritPID < aNbCritPoints; aCritPID++)
+ { //We are not interested in elements of aCriticalDelta array
+ //if their index is greater than or equal to aNbCritPoints
+
+ aCriticalDelta[aCritPID] = anUf - anU1crit[aCritPID];
+ }
+ Standard_Real anU1 = anUf, aMinCriticalParam = anUf;
Standard_Boolean isFirst = Standard_True;
- while(anU1 <= anUl)
+ while (anU1 <= anUl)
{
- for(Standard_Integer i = 0; i < aNbCritPointsMax; i++)
+ //Change value of U-parameter on the 1st surface from anUf to anUl
+ //(anUf will be modified in the cycle body). However, this cycle
+ //can be broken if WL goes though some critical point.
+ //Step is computed adaptively (see comments below).
+
+ for (Standard_Integer i = 0; i < aNbCritPoints; i++)
{
- if((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0)
+ if ((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0)
{
+ //WL has gone through i-th critical point
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;
}
}
- if(IsEqual(anU1, anUl))
+ if (IsEqual(anU1, anUl))
{
- for(Standard_Integer i = 0; i < aNbWLines; i++)
+ for (Standard_Integer i = 0; i < aNbWLines; i++)
{
aWLFindStatus[i] = WLFStatus_Broken;
anUexpect[i] = anU1;
- if(isDeltaPeriod)
+ if (isDeltaPeriod)
{
- //if isAddedIntoWL* == TRUE WLine contains only one point
+ //if isAddedIntoWL[i] == TRUE WLine contains only one point
//(which was end point of previous WLine). If we will
//add point found on the current step WLine will contain only
//two points. At that both these points will be equal to the
}
else
{
- isAddingWLEnabled[i] = ((theTol2D >= (anUexpect[i] - anU1)) ||
+ isAddingWLEnabled[i] = ((aTol2D >= (anUexpect[i] - anU1)) ||
(aWLFindStatus[i] == WLFStatus_Absent));
}
}//for(Standard_Integer i = 0; i < aNbWLines; i++)
}
else
{
- for(Standard_Integer i = 0; i < aNbWLines; i++)
+ for (Standard_Integer i = 0; i < aNbWLines; i++)
{
- isAddingWLEnabled[i] = ((theTol2D >= (anUexpect[i] - anU1)) ||
+ isAddingWLEnabled[i] = ((aTol2D >= (anUexpect[i] - anU1)) ||
(aWLFindStatus[i] == WLFStatus_Absent));
}//for(Standard_Integer i = 0; i < aNbWLines; i++)
}
- for(Standard_Integer i = 0; i < aNbWLines; i++)
+ for (Standard_Integer i = 0; i < aNbWLines; i++)
{
- const Standard_Integer aNbPntsWL = aWLine[i].IsNull() ? 0 :
- aWLine[i]->Curve()->NbPoints();
+ const Standard_Integer aNbPntsWL = aWLine[i].IsNull() ? 0 :
+ aWLine[i]->Curve()->NbPoints();
- CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]);
+ if ((aWLFindStatus[i] == WLFStatus_Broken) ||
+ (aWLFindStatus[i] == WLFStatus_Absent))
+ {//Begin and end of WLine must be on boundary point
+ //or on seam-edge strictly (if it is possible).
+
+ Standard_Real aTol = aTol2D;
+ ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs,
+ aU2[i], &aTol);
+ InscribePoint(aUSurf2f, aUSurf2l, aU2[i], aTol2D, aPeriod, Standard_False);
+
+ aTol = Max(aTol, aTol2D);
- InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False);
+ if (Abs(aU2[i]) <= aTol)
+ aU2[i] = 0.0;
+ else if (Abs(aU2[i] - aPeriod) <= aTol)
+ aU2[i] = aPeriod;
+ else if (Abs(aU2[i] - aUSurf2f) <= aTol)
+ aU2[i] = aUSurf2f;
+ else if (Abs(aU2[i] - aUSurf2l) <= aTol)
+ aU2[i] = aUSurf2l;
+ }
+ else
+ {
+ ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]);
+ InscribePoint(aUSurf2f, aUSurf2l, aU2[i], aTol2D, aPeriod, Standard_False);
+ }
- if(aNbPntsWL == 0)
+ if (aNbPntsWL == 0)
{//the line has not contained any points yet
- if(((aUSurf2f + aPeriod - aUSurf2l) <= 2.0*theTol2D) &&
- ((Abs(aU2[i] - aUSurf2f) < theTol2D) ||
- (Abs(aU2[i]-aUSurf2l) < theTol2D)))
+ if (((aUSurf2f + aPeriod - aUSurf2l) <= 2.0*aTol2D) &&
+ ((Abs(aU2[i] - aUSurf2f) < aTol2D) ||
+ (Abs(aU2[i] - aUSurf2l) < aTol2D)))
{
//In this case aU2[i] can have two values: current aU2[i] or
//aU2[i]+aPeriod (aU2[i]-aPeriod). It is necessary to choose
//correct value.
Standard_Boolean isIncreasing = Standard_True;
- CylCylMonotonicity(anU1, 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
//increased too. And we will go out of surface boundary.
//Therefore, If U2(U1) is increasing then U2 must be equal aUSurf2f.
//Analogically, if U2(U1) is decreasing.
- if(isIncreasing)
+ if (isIncreasing)
{
aU2[i] = aUSurf2f;
}
}
else
{//aNbPntsWL > 0
- if(((aUSurf2l - aUSurf2f) >= aPeriod) &&
- ((Abs(aU2[i] - aUSurf2f) < theTol2D) ||
- (Abs(aU2[i]-aUSurf2l) < theTol2D)))
+ if (((aUSurf2l - aUSurf2f) >= aPeriod) &&
+ ((Abs(aU2[i] - aUSurf2f) < aTol2D) ||
+ (Abs(aU2[i] - aUSurf2l) < aTol2D)))
{//end of the line
Standard_Real aU2prev = 0.0, aV2prev = 0.0;
- if(isTheReverse)
+ if (isReversed)
aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS1(aU2prev, aV2prev);
else
aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS2(aU2prev, aV2prev);
- if(2.0*Abs(aU2prev - aU2[i]) > aPeriod)
+ if (2.0*Abs(aU2prev - aU2[i]) > aPeriod)
{
- if(aU2prev > aU2[i])
+ if (aU2prev > aU2[i])
aU2[i] += aPeriod;
else
aU2[i] -= aPeriod;
}
}
- if( (aWLFindStatus[i] == WLFStatus_Broken) ||
- (aWLFindStatus[i] == WLFStatus_Absent))
- {//Begin and end of WLine must be on boundary point
- //or on seam-edge strictly (if it is possible).
- if(Abs(aU2[i]) <= theTol2D)
- aU2[i] = 0.0;
- else if(Abs(aU2[i] - aPeriod) <= theTol2D)
- aU2[i] = aPeriod;
- else if(Abs(aU2[i] - aUSurf2f) <= theTol2D)
- aU2[i] = aUSurf2f;
- else if(Abs(aU2[i] - aUSurf2l) <= theTol2D)
- aU2[i] = aUSurf2l;
- }
-
- CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs, aV1[i], aV2[i]);
+ ComputationMethods::CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs,
+ aV1[i], aV2[i]);
- if(isFirst)
+ if (isFirst)
{
aV1Prev[i] = aV1[i];
aV2Prev[i] = aV2[i];
//Looking for points into WLine
Standard_Boolean isBroken = Standard_False;
- for(Standard_Integer i = 0; i < aNbWLines; i++)
+ for (Standard_Integer i = 0; i < aNbWLines; i++)
{
- if(!isAddingWLEnabled[i])
+ if (!isAddingWLEnabled[i])
{
Standard_Boolean isBoundIntersect = Standard_False;
- if( (Abs(aV1[i] - aVSurf1f) <= theTol2D) ||
- ((aV1[i]-aVSurf1f)*(aV1Prev[i]-aVSurf1f) < 0.0))
+ if ((Abs(aV1[i] - aVSurf1f) <= aTol2D) ||
+ ((aV1[i] - aVSurf1f)*(aV1Prev[i] - aVSurf1f) < 0.0))
{
isBoundIntersect = Standard_True;
}
- else if( (Abs(aV1[i] - aVSurf1l) <= theTol2D) ||
- ( (aV1[i]-aVSurf1l)*(aV1Prev[i]-aVSurf1l) < 0.0))
+ else if ((Abs(aV1[i] - aVSurf1l) <= aTol2D) ||
+ ((aV1[i] - aVSurf1l)*(aV1Prev[i] - aVSurf1l) < 0.0))
{
isBoundIntersect = Standard_True;
}
- else if( (Abs(aV2[i] - aVSurf2f) <= theTol2D) ||
- ( (aV2[i]-aVSurf2f)*(aV2Prev[i]-aVSurf2f) < 0.0))
+ else if ((Abs(aV2[i] - aVSurf2f) <= aTol2D) ||
+ ((aV2[i] - aVSurf2f)*(aV2Prev[i] - aVSurf2f) < 0.0))
{
isBoundIntersect = Standard_True;
}
- else if( (Abs(aV2[i] - aVSurf2l) <= theTol2D) ||
- ( (aV2[i]-aVSurf2l)*(aV2Prev[i]-aVSurf2l) < 0.0))
+ else if ((Abs(aV2[i] - aVSurf2l) <= aTol2D) ||
+ ((aV2[i] - aVSurf2l)*(aV2Prev[i] - aVSurf2l) < 0.0))
{
isBoundIntersect = Standard_True;
- }
-
- aV1Prev[i] = aV1[i];
- aV2Prev[i] = aV2[i];
+ }
- if(aWLFindStatus[i] == WLFStatus_Broken)
+ if (aWLFindStatus[i] == WLFStatus_Broken)
isBroken = Standard_True;
- if(!isBoundIntersect)
+ if (!isBoundIntersect)
{
continue;
}
}
}
- const Standard_Boolean isInscribe =
- ((aUSurf2f-aU2[i]) <= theTol2D) && ((aU2[i]-aUSurf2l) <= theTol2D) &&
- ((aVSurf1f - aV1[i]) <= theTol2D) && ((aV1[i] - aVSurf1l) <= theTol2D) &&
- ((aVSurf2f - aV2[i]) <= theTol2D) && ((aV2[i] - aVSurf2l) <= theTol2D);
+ // True if the current point already in the domain
+ const Standard_Boolean isInscribe =
+ ((aUSurf2f - aU2[i]) <= aTol2D) && ((aU2[i] - aUSurf2l) <= aTol2D) &&
+ ((aVSurf1f - aV1[i]) <= aTol2D) && ((aV1[i] - aVSurf1l) <= aTol2D) &&
+ ((aVSurf2f - aV2[i]) <= aTol2D) && ((aV2[i] - aVSurf2l) <= aTol2D);
//isVIntersect == TRUE if intersection line intersects two (!)
//V-bounds of cylinder (1st or 2nd - no matter)
const Standard_Boolean isVIntersect =
- ( ((aVSurf1f-aV1[i])*(aVSurf1f-aV1Prev[i]) < RealSmall()) &&
- ((aVSurf1l-aV1[i])*(aVSurf1l-aV1Prev[i]) < RealSmall())) ||
- ( ((aVSurf2f-aV2[i])*(aVSurf2f-aV2Prev[i]) < RealSmall()) &&
- ((aVSurf2l-aV2[i])*(aVSurf2l-aV2Prev[i]) < RealSmall()));
-
+ (((aVSurf1f - aV1[i])*(aVSurf1f - aV1Prev[i]) < RealSmall()) &&
+ ((aVSurf1l - aV1[i])*(aVSurf1l - aV1Prev[i]) < RealSmall())) ||
+ (((aVSurf2f - aV2[i])*(aVSurf2f - aV2Prev[i]) < RealSmall()) &&
+ ((aVSurf2l - aV2[i])*(aVSurf2l - aV2Prev[i]) < RealSmall()));
//isFound1 == TRUE if intersection line intersects V-bounds
// (First or Last - no matter) of the 1st cylynder
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))
+ if (((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1 - aUSurf1l) < aTol2D))
{
isForce = Standard_True;
}
}
- 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);
+ theBW.AddBoundaryPoint(aWLine[i], anU1, aMinCriticalParam, aU2[i],
+ aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], i, isForce,
+ isFound1, isFound2);
const Standard_Boolean isPrevVBound = !isVIntersect &&
- ((Abs(aV1Prev[i] - aVSurf1f) <= theTol2D) ||
- (Abs(aV1Prev[i] - aVSurf1l) <= theTol2D) ||
- (Abs(aV2Prev[i] - aVSurf2f) <= theTol2D) ||
- (Abs(aV2Prev[i] - aVSurf2l) <= theTol2D));
+ ((Abs(aV1Prev[i] - aVSurf1f) <= aTol2D) ||
+ (Abs(aV1Prev[i] - aVSurf1l) <= aTol2D) ||
+ (Abs(aV2Prev[i] - aVSurf2f) <= aTol2D) ||
+ (Abs(aV2Prev[i] - aVSurf2l) <= aTol2D));
+
+ aV1Prev[i] = aV1[i];
+ aV2Prev[i] = aV2[i];
- if((aWLFindStatus[i] == WLFStatus_Exist) && (isFound1 || isFound2) && !isPrevVBound)
+ if ((aWLFindStatus[i] == WLFStatus_Exist) && (isFound1 || isFound2) && !isPrevVBound)
{
aWLFindStatus[i] = WLFStatus_Broken; //start a new line
}
- else if(isInscribe)
+ else if (isInscribe)
{
- if((aWLFindStatus[i] == WLFStatus_Absent) && (isFound1 || isFound2))
+ if ((aWLFindStatus[i] == WLFStatus_Absent) && (isFound1 || isFound2))
{
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(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse, Standard_True,
+ if (aWLine[i]->NbPnts() > 0)
+ {
+ Standard_Real aU2p = 0.0, aV2p = 0.0;
+ if (isReversed)
+ aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
+ else
+ aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
+
+ const Standard_Real aDelta = aU2[i] - aU2p;
+
+ if (2.0 * Abs(aDelta) > aPeriod)
+ {
+ if (aDelta > 0.0)
+ {
+ aU2[i] -= aPeriod;
+ }
+ else
+ {
+ aU2[i] += aPeriod;
+ }
+ }
+ }
+
+ if(AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed, Standard_True,
gp_Pnt2d(anU1, aV1[i]), gp_Pnt2d(aU2[i], aV2[i]),
aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod,
- aWLine[i]->Curve(), i, theTol3D, theTol2D, isForce))
+ aWLine[i]->Curve(), i, aTol3D, aTol2D, isForce))
{
- if(aWLFindStatus[i] == WLFStatus_Absent)
+ if (aWLFindStatus[i] == WLFStatus_Absent)
{
aWLFindStatus[i] = WLFStatus_Exist;
}
}
+ else if (!isFound1 && !isFound2)
+ {//We do not add any point while doing this iteration
+ if (aWLFindStatus[i] == WLFStatus_Exist)
+ {
+ aWLFindStatus[i] = WLFStatus_Broken;
+ }
+ }
+ }
+ }
+ else
+ {//We do not add any point while doing this iteration
+ if (aWLFindStatus[i] == WLFStatus_Exist)
+ {
+ aWLFindStatus[i] = WLFStatus_Broken;
}
}
-
- aV1Prev[i] = aV1[i];
- aV2Prev[i] = aV2[i];
- if(aWLFindStatus[i] == WLFStatus_Broken)
+ if (aWLFindStatus[i] == WLFStatus_Broken)
isBroken = Standard_True;
}//for(Standard_Integer i = 0; i < aNbWLines; i++)
- if(isBroken)
+ if (isBroken)
{//current lines are filled. Go to the next lines
anUf = anU1;
- Standard_Real anULastAdded = RealFirst();
- for(Standard_Integer i = 0; i < aNbWLines; i++)
+ Standard_Boolean isAdded = Standard_True;
+
+ for (Standard_Integer i = 0; i < aNbWLines; i++)
{
- if((aWLFindStatus[i] != WLFStatus_Broken) || aWLine[i]->NbPnts() < 2)
+ if (isAddingWLEnabled[i])
+ {
continue;
+ }
- Standard_Real aU1c = 0.0, aV1c = 0.0;
+ isAdded = Standard_False;
- if(isTheReverse)
- aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS2(aU1c, aV1c);
- else
- aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS1(aU1c, aV1c);
+ Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
+
+ theBW.AddBoundaryPoint(aWLine[i], anU1, aMinCriticalParam, aU2[i],
+ aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], i,
+ Standard_False, isFound1, isFound2);
+
+ if (isFound1 || isFound2)
+ {
+ isAdded = Standard_True;
+ }
+
+ if (aWLine[i]->NbPnts() > 0)
+ {
+ Standard_Real aU2p = 0.0, aV2p = 0.0;
+ if (isReversed)
+ aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
+ else
+ aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
+
+ const Standard_Real aDelta = aU2[i] - aU2p;
+
+ if (2 * Abs(aDelta) > aPeriod)
+ {
+ if (aDelta > 0.0)
+ {
+ aU2[i] -= aPeriod;
+ }
+ else
+ {
+ aU2[i] += aPeriod;
+ }
+ }
+ }
- anULastAdded = Max(anULastAdded, aU1c);
+ if(AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
+ Standard_True, gp_Pnt2d(anU1, aV1[i]),
+ gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
+ aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
+ aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
+ i, aTol3D, aTol2D, Standard_False))
+ {
+ isAdded = Standard_True;
+ }
}
- for(Standard_Integer i = 0; i < aNbWLines; i++)
+ if (!isAdded)
{
- if((aWLFindStatus[i] != WLFStatus_Exist))
- continue;
-
- Standard_Real aU2c = 0.0, aV1c = 0.0, aV2c = 0.0;
- Standard_Boolean isAdded = Standard_False;
+ //Before breaking WL, we must complete it correctly
+ //(e.g. to prolong to the surface boundary).
+ //Therefore, we take the point last added in some WL
+ //(have maximal U1-parameter) and try to add it in
+ //the current WL.
+ Standard_Real anUmaxAdded = RealFirst();
- //Try to add current point
- if(CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2c, aV1c, aV2c))
{
- InscribePoint(aUSurf2f, aUSurf2l, aU2c, theTol2D, aPeriod, Standard_False);
+ Standard_Boolean isChanged = Standard_False;
+ for (Standard_Integer i = 0; i < aNbWLines; i++)
+ {
+ if ((aWLFindStatus[i] == WLFStatus_Absent) || (aWLine[i]->NbPnts() == 0))
+ continue;
- isAdded = AddPointIntoWL( theQuad1, theQuad2, anEquationCoeffs, isTheReverse,
- Standard_True, gp_Pnt2d(anU1, aV1c), gp_Pnt2d(aU2c, aV2c),
- aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
- aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod,
- aWLine[i]->Curve(), i, theTol3D, theTol2D, Standard_False);
- }
+ Standard_Real aU1c = 0.0, aV1c = 0.0;
+ if (isReversed)
+ aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS2(aU1c, aV1c);
+ else
+ aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS1(aU1c, aV1c);
- if(isAdded)
- continue;
+ anUmaxAdded = Max(anUmaxAdded, aU1c);
+ isChanged = Standard_True;
+ }
- if(anULastAdded == anU1)
- {//strictly equal only
- continue;
+ if (!isChanged)
+ { //If anUmaxAdded were not changed in previous cycle then
+ //we would break existing WLines.
+ break;
+ }
}
- //Try to add last added point
-
- if(!CylCylComputeParameters(anULastAdded, i, anEquationCoeffs, aU2c, aV1c, aV2c))
- continue;
+ for (Standard_Integer i = 0; i < aNbWLines; i++)
+ {
+ if (isAddingWLEnabled[i])
+ {
+ continue;
+ }
- InscribePoint(aUSurf2f, aUSurf2l, aU2c, theTol2D, aPeriod, Standard_False);
+ ComputationMethods::CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs,
+ aU2[i], aV1[i], aV2[i]);
- AddPointIntoWL( theQuad1, theQuad2, anEquationCoeffs, isTheReverse, Standard_True,
- gp_Pnt2d(anULastAdded, aV1c), gp_Pnt2d(aU2c, aV2c),
- aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
- aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod,
- aWLine[i]->Curve(), i, theTol3D, theTol2D, Standard_False);
+ AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
+ Standard_True, gp_Pnt2d(anUmaxAdded, aV1[i]),
+ gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
+ aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
+ aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
+ i, aTol3D, aTol2D, Standard_False);
+ }
}
break;
//Step computing
{
- const Standard_Real aDeltaV1 = (aVSurf1l - aVSurf1f)/IntToReal(aNbPoints),
- aDeltaV2 = (aVSurf2l - aVSurf2f)/IntToReal(aNbPoints);
-
+ //Step of aU1-parameter is computed adaptively. The algorithm
+ //aims to provide given aDeltaV1 and aDeltaV2 values (if it is
+ //possible because the intersection line can go along V-isoline)
+ //in every iteration. It allows avoiding "flying" intersection
+ //points too far each from other (see issue #24915).
+
+ const Standard_Real aDeltaV1 = aRangeS1.Delta() / IntToReal(aNbPoints),
+ aDeltaV2 = aRangeS2.Delta() / IntToReal(aNbPoints);
+
math_Matrix aMatr(1, 3, 1, 5);
Standard_Real aMinUexp = RealLast();
-
- for(Standard_Integer i = 0; i < aNbWLines; i++)
+
+ for (Standard_Integer i = 0; i < aNbWLines; i++)
{
- if(theTol2D < (anUexpect[i] - anU1))
+ if (aTol2D < (anUexpect[i] - anU1))
{
continue;
}
- if(aWLFindStatus[i] == WLFStatus_Absent)
+ if (aWLFindStatus[i] == WLFStatus_Absent)
{
anUexpect[i] += aStepMax;
aMinUexp = Min(aMinUexp, anUexpect[i]);
anEquationCoeffs.mVecA2*aCosU2 + anEquationCoeffs.mVecB2*aSinU2 +
anEquationCoeffs.mVecD);
- if(!StepComputing(aMatr, aV1[i], aV2[i], aDeltaV1, aDeltaV2, aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aStepTmp))
+ //The main idea is in solving of linearized system (2)
+ //(see description to ComputationMethods class) in order to find new U1-value
+ //to provide new value V1 or V2, which differs from current one by aDeltaV1 or
+ //aDeltaV2 respectively.
+
+ //While linearizing, following Taylor formulas are used:
+ // cos(x0+dx) = cos(x0) - sin(x0)*dx
+ // sin(x0+dx) = sin(x0) + cos(x0)*dx
+
+ //Consequently cos(U1), cos(U2), sin(U1) and sin(U2) in the system (2)
+ //must be substituted by corresponding values.
+
+ //ATTENTION!!!
+ //The solution is approximate. More over, all requirements to the
+ //linearization must be satisfied in order to obtain quality result.
+
+ if (!StepComputing(aMatr, aV1[i], aV2[i], aDeltaV1, aDeltaV2, aStepTmp))
{
//To avoid cycling-up
anUexpect[i] += aStepMax;
continue;
}
- if(aStepTmp < aStepMin)
+ if (aStepTmp < aStepMin)
aStepTmp = aStepMin;
-
- if(aStepTmp > aStepMax)
+
+ if (aStepTmp > aStepMax)
aStepTmp = aStepMax;
anUexpect[i] = anU1 + aStepTmp;
anU1 = aMinUexp;
}
- if(Precision::PConfusion() >= (anUl - anU1))
+ if (Precision::PConfusion() >= (anUl - anU1))
anU1 = anUl;
anUf = anU1;
- for(Standard_Integer i = 0; i < aNbWLines; i++)
+ for (Standard_Integer i = 0; i < aNbWLines; i++)
{
- if(aWLine[i]->NbPnts() != 1)
+ if (aWLine[i]->NbPnts() != 1)
isAddedIntoWL[i] = Standard_False;
- if(anU1 == anUl)
+ if (anU1 == anUl)
{//strictly equal. Tolerance is considered above.
anUexpect[i] = anUl;
}
}
}
- for(Standard_Integer i = 0; i < aNbWLines; i++)
+ for (Standard_Integer i = 0; i < aNbWLines; i++)
{
- if((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
+ if ((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
{
isTheEmpty = Standard_False;
Standard_Real u1, v1, u2, v2;
IntPatch_Point aP;
aP.SetParameter(u1);
aP.SetParameters(u1, v1, u2, v2);
- aP.SetTolerance(theTol3D);
+ aP.SetTolerance(aTol3D);
aP.SetValue(aWLine[i]->Point(1).Value());
- theSPnt.Append(aP);
+ //Check whether the added point exists.
+ //It is enough to check the last point.
+ if (theSPnt.IsEmpty() ||
+ !theSPnt.Last().PntOn2S().IsSame(aP.PntOn2S(), Precision::Confusion()))
+ {
+ theSPnt.Append(aP);
+ }
}
- else if(aWLine[i]->NbPnts() > 1)
+ else if (aWLine[i]->NbPnts() > 1)
{
Standard_Boolean isGood = Standard_True;
- if(aWLine[i]->NbPnts() == 2)
+ if (aWLine[i]->NbPnts() == 2)
{
const IntSurf_PntOn2S& aPf = aWLine[i]->Point(1);
const IntSurf_PntOn2S& aPl = aWLine[i]->Point(2);
-
- if(aPf.IsSame(aPl, Precision::Confusion()))
+
+ if (aPf.IsSame(aPl, Precision::Confusion()))
isGood = Standard_False;
}
+ else if (aWLine[i]->NbPnts() > 2)
+ {
+ // Sometimes points of the WLine are distributed
+ // linearly and uniformly. However, such position
+ // of the points does not always describe the real intersection
+ // curve. I.e. real tangents at the ends of the intersection
+ // curve can significantly deviate from this "line" direction.
+ // Here we are processing this case by inserting additional points
+ // to the beginning/end of the WLine to make it more precise.
+ // See description to the issue #30082.
+
+ const Standard_Real aSqTol3D = aTol3D*aTol3D;
+ for (Standard_Integer j = 0; j < 2; j++)
+ {
+ // If j == 0 ==> add point at begin of WLine.
+ // If j == 1 ==> add point at end of WLine.
+
+ for (;;)
+ {
+ if (aWLine[i]->NbPnts() >= aNbMaxPoints)
+ {
+ break;
+ }
+
+ // Take 1st and 2nd point to compute the "line" direction.
+ // For our convenience, we make 2nd point be the ends of the WLine
+ // because it will be used for computation of the normals
+ // to the surfaces.
+ const Standard_Integer anIdx1 = j ? aWLine[i]->NbPnts() - 1 : 2;
+ const Standard_Integer anIdx2 = j ? aWLine[i]->NbPnts() : 1;
+
+ const gp_Pnt &aP1 = aWLine[i]->Point(anIdx1).Value();
+ const gp_Pnt &aP2 = aWLine[i]->Point(anIdx2).Value();
+
+ const gp_Vec aDir(aP1, aP2);
+
+ if (aDir.SquareMagnitude() < aSqTol3D)
+ {
+ break;
+ }
+
+ // Compute tangent in first/last point of the WLine.
+ // We do not take into account the flag "isReversed"
+ // because strict direction of the tangent is not
+ // important here (we are interested in the tangent
+ // line itself and nothing to fear if its direction
+ // is reversed).
+ const gp_Vec aN1 = aQuad1.Normale(aP2);
+ const gp_Vec aN2 = aQuad2.Normale(aP2);
+ const gp_Vec aTg(aN1.Crossed(aN2));
+
+ if (aTg.SquareMagnitude() < Precision::SquareConfusion())
+ {
+ // Tangent zone
+ break;
+ }
+
+ // Check of the bending
+ Standard_Real anAngle = aDir.Angle(aTg);
+
+ if (anAngle > M_PI_2)
+ anAngle -= M_PI;
+
+ if (Abs(anAngle) > 0.25) // ~ 14deg.
+ {
+ const Standard_Integer aNbPntsPrev = aWLine[i]->NbPnts();
+ SeekAdditionalPoints(aQuad1, aQuad2, aWLine[i]->Curve(),
+ anEquationCoeffs, i, 3, anIdx1, anIdx2,
+ aTol2D, aPeriod, isReversed);
+
+ if (aWLine[i]->NbPnts() == aNbPntsPrev)
+ {
+ // No points have been added. ==> Exit from a loop.
+ break;
+ }
+ }
+ else
+ {
+ // Good result has been achieved. ==> Exit from a loop.
+ break;
+ }
+ } // for (;;)
+ }
+ }
- if(isGood)
+ if (isGood)
{
isTheEmpty = Standard_False;
isAddedIntoWL[i] = Standard_True;
- SeekAdditionalPoints( theQuad1, theQuad2, aWLine[i]->Curve(),
- anEquationCoeffs, i, aNbPoints, 1,
- aWLine[i]->NbPnts(), aUSurf2f, aUSurf2l,
- theTol2D, aPeriod, isTheReverse);
+ SeekAdditionalPoints(aQuad1, aQuad2, aWLine[i]->Curve(),
+ anEquationCoeffs, i, aNbPoints, 1,
+ aWLine[i]->NbPnts(), aTol2D, aPeriod,
+ isReversed);
- aWLine[i]->ComputeVertexParameters(theTol3D);
+ aWLine[i]->ComputeVertexParameters(aTol3D);
theSlin.Append(aWLine[i]);
}
}
isAddedIntoWL[i] = Standard_False;
}
-#ifdef OCCT_DEBUG
+#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
//aWLine[i]->Dump();
#endif
}
}
-//Delete the points in theSPnt, which
-//lie at least in one of the line in theSlin.
- for(Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
+ //Delete the points in theSPnt, which
+ //lie at least in one of the line in theSlin.
+ for (Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
{
- for(Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); aNbLin++)
+ 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());
const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNbPnt).PntOn2S();
- if( aPntCur.IsSame(aPntFWL1, Precision::Confusion()) ||
- aPntCur.IsSame(aPntLWL1, Precision::Confusion()))
+ if (aPntCur.IsSame(aPntFWL1, aTol3D) ||
+ aPntCur.IsSame(aPntLWL1, aTol3D))
{
theSPnt.Remove(aNbPnt);
aNbPnt--;
}
}
- const Standard_Real aDU = aStepMin + Epsilon(aStepMin);
- //Try to add new points in the neighbourhood of existing point
- for(Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
+ //Try to add new points in the neighborhood of existing point
+ for (Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
{
+ // Standard algorithm (implemented above) could not find any
+ // continuous curve in neighborhood of aPnt2S (e.g. because
+ // this curve is too small; see tests\bugs\modalg_5\bug25292_35 and _36).
+ // Here, we will try to find several new points nearer to aPnt2S.
+
+ // The algorithm below tries to find two points in every
+ // intervals [u1 - aStepMax, u1] and [u1, u1 + aStepMax]
+ // and every new point will be in maximal distance from
+ // u1. If these two points exist they will be joined
+ // by the intersection curve.
+
const IntPatch_Point& aPnt2S = theSPnt.Value(aNbPnt);
Standard_Real u1, v1, u2, v2;
Handle(IntSurf_LineOn2S) aL2S = new IntSurf_LineOn2S();
Handle(IntPatch_WLine) aWLine = new IntPatch_WLine(aL2S, Standard_False);
- aWLine->Curve()->Add(aPnt2S.PntOn2S());
+ aWLine->SetCreatingWayInfo(IntPatch_WLine::IntPatch_WLImpImp);
//Define the index of WLine, which lies the point aPnt2S in.
- Standard_Real anUf = 0.0, anUl = 0.0, aCurU2 = 0.0;
Standard_Integer anIndex = 0;
- if(isTheReverse)
+
+ Standard_Real anUf = 0.0, anUl = 0.0, aCurU2 = 0.0;
+ if (isReversed)
{
anUf = Max(u2 - aStepMax, aUSurf1f);
- anUl = u2;
+ anUl = Min(u2 + aStepMax, aUSurf1l);
aCurU2 = u1;
}
else
{
anUf = Max(u1 - aStepMax, aUSurf1f);
- anUl = u1;
+ anUl = Min(u1 + aStepMax, aUSurf1l);
aCurU2 = u2;
}
- Standard_Real aDelta = RealLast();
- for (Standard_Integer i = 0; i < aNbWLines; i++)
- {
- Standard_Real anU2t = 0.0;
- if(!CylCylComputeParameters(anUl, i, anEquationCoeffs, anU2t))
- continue;
- const Standard_Real aDU = Abs(anU2t - aCurU2);
- if(aDU < aDelta)
+ const Standard_Real anUinf = anUf, anUsup = anUl, anUmid = 0.5*(anUf + anUl);
+
+ {
+ //Find the value of anIndex variable.
+ Standard_Real aDelta = RealLast();
+ for (Standard_Integer i = 0; i < aNbWLines; i++)
{
- aDelta = aDU;
- anIndex = i;
+ Standard_Real anU2t = 0.0;
+ if (!ComputationMethods::CylCylComputeParameters(anUmid, i, anEquationCoeffs, anU2t))
+ continue;
+
+ Standard_Real aDU2 = fmod(Abs(anU2t - aCurU2), aPeriod);
+ aDU2 = Min(aDU2, Abs(aDU2 - aPeriod));
+ if (aDU2 < aDelta)
+ {
+ aDelta = aDU2;
+ anIndex = i;
+ }
}
}
- //Try to fill aWLine by additional points
- while(anUl - anUf > RealSmall())
+ // Bisection method is used in order to find every new point.
+ // I.e. if we need to find intersection point in the interval [anUinf, anUmid]
+ // we check the point anUC = 0.5*(anUinf+anUmid). If it is an intersection point
+ // we try to find another point in the interval [anUinf, anUC] (because we find the point in
+ // maximal distance from anUmid). If it is not then we try to find another point in the
+ // interval [anUC, anUmid]. Next iterations will be made analogically.
+ // When we find intersection point in the interval [anUmid, anUsup] we try to find
+ // another point in the interval [anUC, anUsup] if anUC is intersection point and
+ // in the interval [anUmid, anUC], otherwise.
+
+ Standard_Real anAddedPar[2] = {isReversed ? u2 : u1, isReversed ? u2 : u1};
+
+ for (Standard_Integer aParID = 0; aParID < 2; aParID++)
{
- Standard_Real anU2 = 0.0, anV1 = 0.0, anV2 = 0.0;
- Standard_Boolean isDone =
- CylCylComputeParameters(anUf, anIndex, anEquationCoeffs,
- anU2, anV1, anV2);
- anUf += aDU;
-
- if(!isDone)
+ if (aParID == 0)
{
- continue;
+ anUf = anUinf;
+ anUl = anUmid;
+ }
+ else // if(aParID == 1)
+ {
+ anUf = anUmid;
+ anUl = anUsup;
}
- 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,
- theTol2D, Standard_False, Standard_True))
+ Standard_Real &aPar1 = (aParID == 0) ? anUf : anUl,
+ &aPar2 = (aParID == 0) ? anUl : anUf;
+
+ while (Abs(aPar2 - aPar1) > aStepMin)
{
- break;
+ Standard_Real anUC = 0.5*(anUf + anUl);
+ Standard_Real aU2 = 0.0, aV1 = 0.0, aV2 = 0.0;
+ Standard_Boolean isDone = ComputationMethods::
+ CylCylComputeParameters(anUC, anIndex, anEquationCoeffs, aU2, aV1, aV2);
+
+ if (isDone)
+ {
+ if (Abs(aV1 - aVSurf1f) <= aTol2D)
+ aV1 = aVSurf1f;
+
+ if (Abs(aV1 - aVSurf1l) <= aTol2D)
+ aV1 = aVSurf1l;
+
+ if (Abs(aV2 - aVSurf2f) <= aTol2D)
+ aV2 = aVSurf2f;
+
+ if (Abs(aV2 - aVSurf2l) <= aTol2D)
+ aV2 = aVSurf2l;
+
+ isDone = AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
+ Standard_True, gp_Pnt2d(anUC, aV1), gp_Pnt2d(aU2, aV2),
+ aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
+ aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l,
+ aPeriod, aWLine->Curve(), anIndex, aTol3D,
+ aTol2D, Standard_False, Standard_True);
+ }
+
+ if (isDone)
+ {
+ anAddedPar[0] = Min(anAddedPar[0], anUC);
+ anAddedPar[1] = Max(anAddedPar[1], anUC);
+ aPar2 = anUC;
+ }
+ else
+ {
+ aPar1 = anUC;
+ }
}
}
- if(aWLine->NbPnts() > 1)
+ //Fill aWLine by additional points
+ if (anAddedPar[1] - anAddedPar[0] > aStepMin)
{
- SeekAdditionalPoints( theQuad1, theQuad2, aWLine->Curve(),
+ for (Standard_Integer aParID = 0; aParID < 2; aParID++)
+ {
+ Standard_Real aU2 = 0.0, aV1 = 0.0, aV2 = 0.0;
+ ComputationMethods::CylCylComputeParameters(anAddedPar[aParID], anIndex,
+ anEquationCoeffs, aU2, aV1, aV2);
+
+ AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed, Standard_True,
+ gp_Pnt2d(anAddedPar[aParID], aV1), gp_Pnt2d(aU2, aV2),
+ aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
+ aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod, aWLine->Curve(),
+ anIndex, aTol3D, aTol2D, Standard_False, Standard_False);
+ }
+
+ SeekAdditionalPoints(aQuad1, aQuad2, aWLine->Curve(),
anEquationCoeffs, anIndex, aNbMinPoints,
- 1, aWLine->NbPnts(), aUSurf2f, aUSurf2l,
- theTol2D, aPeriod, isTheReverse);
+ 1, aWLine->NbPnts(), aTol2D, aPeriod,
+ isReversed);
- aWLine->ComputeVertexParameters(theTol3D);
+ aWLine->ComputeVertexParameters(aTol3D);
theSlin.Append(aWLine);
theSPnt.Remove(aNbPnt);
aNbPnt--;
}
}
+
+ return IntPatch_ImpImpIntersection::IntStatus_OK;
+}
+
+//=======================================================================
+//function : IntCyCy
+//purpose :
+//=======================================================================
+IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1,
+ const IntSurf_Quadric& theQuad2,
+ const Standard_Real theTol3D,
+ const Standard_Real theTol2D,
+ const Bnd_Box2d& theUVSurf1,
+ const Bnd_Box2d& theUVSurf2,
+ Standard_Boolean& isTheEmpty,
+ Standard_Boolean& isTheSameSurface,
+ Standard_Boolean& isTheMultiplePoint,
+ IntPatch_SequenceOfLine& theSlin,
+ IntPatch_SequenceOfPoint& theSPnt)
+{
+ isTheEmpty = Standard_True;
+ isTheSameSurface = Standard_False;
+ isTheMultiplePoint = Standard_False;
+ theSlin.Clear();
+ theSPnt.Clear();
+
+ const gp_Cylinder aCyl1 = theQuad1.Cylinder(),
+ aCyl2 = theQuad2.Cylinder();
+
+ IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
+
+ if (!anInter.IsDone())
+ {
+ return IntPatch_ImpImpIntersection::IntStatus_Fail;
+ }
+
+ if(anInter.TypeInter() != IntAna_NoGeometricSolution)
+ {
+ if (CyCyAnalyticalIntersect(theQuad1, theQuad2, anInter,
+ theTol3D, isTheEmpty,
+ isTheSameSurface, isTheMultiplePoint,
+ theSlin, theSPnt))
+ {
+ return IntPatch_ImpImpIntersection::IntStatus_OK;
+ }
+ }
- return Standard_True;
+ //Here, intersection line is not an analytical curve(line, circle, ellipsis etc.)
+
+ Standard_Real aUSBou[2][2], aVSBou[2][2]; //const
+
+ theUVSurf1.Get(aUSBou[0][0], aVSBou[0][0], aUSBou[0][1], aVSBou[0][1]);
+ theUVSurf2.Get(aUSBou[1][0], aVSBou[1][0], aUSBou[1][1], aVSBou[1][1]);
+
+ const Standard_Real aPeriod = 2.0*M_PI;
+ const Standard_Integer aNbWLines = 2;
+
+ const ComputationMethods::stCoeffsValue anEquationCoeffs1(aCyl1, aCyl2);
+ const ComputationMethods::stCoeffsValue anEquationCoeffs2(aCyl2, aCyl1);
+
+ //Boundaries.
+ //Intersection result can include two non-connected regions
+ //(see WorkWithBoundaries::BoundariesComputing(...) method).
+ const Standard_Integer aNbOfBoundaries = 2;
+ Bnd_Range anURange[2][aNbOfBoundaries]; //const
+
+ if (!WorkWithBoundaries::BoundariesComputing(anEquationCoeffs1, aPeriod, anURange[0]))
+ return IntPatch_ImpImpIntersection::IntStatus_OK;
+
+ if (!WorkWithBoundaries::BoundariesComputing(anEquationCoeffs2, aPeriod, anURange[1]))
+ return IntPatch_ImpImpIntersection::IntStatus_OK;
+
+ //anURange[*] can be in different periodic regions in
+ //compare with First-Last surface. E.g. the surface
+ //is full cylinder [0, 2*PI] but anURange is [5, 7].
+ //Trivial common range computation returns [5, 2*PI] and
+ //its summary length is 2*PI-5 == 1.28... only. That is wrong.
+ //This problem can be solved by the following
+ //algorithm:
+ // 1. split anURange[*] by the surface boundary;
+ // 2. shift every new range in order to inscribe it
+ // in [Ufirst, Ulast] of cylinder;
+ // 3. consider only common ranges between [Ufirst, Ulast]
+ // and new ranges.
+ //
+ // In above example, we obtain following:
+ // 1. two ranges: [5, 2*PI] and [2*PI, 7];
+ // 2. after shifting: [5, 2*PI] and [0, 7-2*PI];
+ // 3. Common ranges: ([5, 2*PI] and [0, 2*PI]) == [5, 2*PI],
+ // ([0, 7-2*PI] and [0, 2*PI]) == [0, 7-2*PI];
+ // 4. Their summary length is (2*PI-5)+(7-2*PI-0)==7-5==2 ==> GOOD.
+
+ Standard_Real aSumRange[2] = { 0.0, 0.0 };
+ Handle(NCollection_IncAllocator) anAlloc = new NCollection_IncAllocator;
+ for (Standard_Integer aCID = 0; aCID < 2; aCID++)
+ {
+ anAlloc->Reset();
+ NCollection_List<Bnd_Range> aListOfRng(anAlloc);
+
+ aListOfRng.Append(anURange[aCID][0]);
+ aListOfRng.Append(anURange[aCID][1]);
+
+ const Standard_Real aSplitArr[3] = {aUSBou[aCID][0], aUSBou[aCID][1], 0.0};
+
+ NCollection_List<Bnd_Range>::Iterator anITrRng;
+ for (Standard_Integer aSInd = 0; aSInd < 3; aSInd++)
+ {
+ NCollection_List<Bnd_Range> aLstTemp(aListOfRng);
+ aListOfRng.Clear();
+ for (anITrRng.Init(aLstTemp); anITrRng.More(); anITrRng.Next())
+ {
+ Bnd_Range& aRng = anITrRng.Value();
+ aRng.Split(aSplitArr[aSInd], aListOfRng, aPeriod);
+ }
+ }
+
+ anITrRng.Init(aListOfRng);
+ for (; anITrRng.More(); anITrRng.Next())
+ {
+ Bnd_Range& aCurrRange = anITrRng.Value();
+
+ Bnd_Range aBoundR;
+ aBoundR.Add(aUSBou[aCID][0]);
+ aBoundR.Add(aUSBou[aCID][1]);
+
+ if (!InscribeInterval(aUSBou[aCID][0], aUSBou[aCID][1],
+ aCurrRange, theTol2D, aPeriod))
+ {
+ //If aCurrRange does not have common block with
+ //[Ufirst, Ulast] of cylinder then we will try
+ //to inscribe [Ufirst, Ulast] in the boundaries of aCurrRange.
+ Standard_Real aF = 1.0, aL = 0.0;
+ if (!aCurrRange.GetBounds(aF, aL))
+ continue;
+
+ if ((aL < aUSBou[aCID][0]))
+ {
+ aCurrRange.Shift(aPeriod);
+ }
+ else if (aF > aUSBou[aCID][1])
+ {
+ aCurrRange.Shift(-aPeriod);
+ }
+ }
+
+ aBoundR.Common(aCurrRange);
+
+ const Standard_Real aDelta = aBoundR.Delta();
+
+ if (aDelta > 0.0)
+ {
+ aSumRange[aCID] += aDelta;
+ }
+ }
+ }
+
+ //The bigger range the bigger number 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.
+ const Standard_Boolean isToReverse = (aSumRange[1] > aSumRange[0]);
+
+ if (isToReverse)
+ {
+ const WorkWithBoundaries aBoundWork(theQuad2, theQuad1, anEquationCoeffs2,
+ theUVSurf2, theUVSurf1, aNbWLines,
+ aPeriod, theTol3D, theTol2D, Standard_True);
+
+ return CyCyNoGeometric(aCyl2, aCyl1, aBoundWork, anURange[1], aNbOfBoundaries,
+ isTheEmpty, theSlin, theSPnt);
+ }
+ else
+ {
+ const WorkWithBoundaries aBoundWork(theQuad1, theQuad2, anEquationCoeffs1,
+ theUVSurf1, theUVSurf2, aNbWLines,
+ aPeriod, theTol3D, theTol2D, Standard_False);
+
+ return CyCyNoGeometric(aCyl1, aCyl2, aBoundWork, anURange[0], aNbOfBoundaries,
+ isTheEmpty, theSlin, theSPnt);
+ }
}
//=======================================================================
//curvsol = anaint.Curve(i);
aC=anaint.Curve(i);
aLC.Clear();
- ExploreCurve(Cy, Co, aC, 10.*Tol, aLC);
+ ExploreCurve(Co, aC, 10.*Tol, aLC);
//
aIt.Initialize(aLC);
for (; aIt.More(); aIt.Next()) {
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;
}
//=======================================================================
//function : ExploreCurve
-//purpose :
+//purpose : Splits aC on several curves in the cone apex points.
//=======================================================================
-Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy,
- const gp_Cone& aCo,
- IntAna_Curve& aC,
- const Standard_Real aTol,
- IntAna_ListOfCurve& aLC)
-
+Standard_Boolean ExploreCurve(const gp_Cone& theCo,
+ IntAna_Curve& theCrv,
+ const Standard_Real theTol,
+ IntAna_ListOfCurve& theLC)
{
- Standard_Boolean bFind=Standard_False;
- Standard_Real aTheta, aT1, aT2, aDst;
- gp_Pnt aPapx, aPx;
- //
- //aC.Dump();
- //
- aLC.Clear();
- aLC.Append(aC);
- //
- aPapx=aCo.Apex();
- //
- aC.Domain(aT1, aT2);
- //
- aPx=aC.Value(aT1);
- aDst=aPx.Distance(aPapx);
- if (aDst<aTol) {
- return bFind;
- }
- aPx=aC.Value(aT2);
- aDst=aPx.Distance(aPapx);
- if (aDst<aTol) {
- return bFind;
- }
+ const Standard_Real aSqTol = theTol*theTol;
+ const gp_Pnt aPapx(theCo.Apex());
+
+ Standard_Real aT1, aT2;
+ theCrv.Domain(aT1, aT2);
+
+ theLC.Clear();
//
- bFind=aC.FindParameter(aPapx, aTheta);
- if (!bFind){
- return bFind;
+ TColStd_ListOfReal aLParams;
+ theCrv.FindParameter(aPapx, aLParams);
+ if (aLParams.IsEmpty())
+ {
+ theLC.Append(theCrv);
+ return Standard_False;
}
- //
- aPx=aC.Value(aTheta);
- aDst=aPx.Distance(aPapx);
- if (aDst>aTol) {
- return !bFind;
+
+ for (TColStd_ListIteratorOfListOfReal anItr(aLParams); anItr.More(); anItr.Next())
+ {
+ Standard_Real aPrm = anItr.Value();
+
+ if ((aPrm - aT1) < Precision::PConfusion())
+ continue;
+
+ Standard_Boolean isLast = Standard_False;
+ if ((aT2 - aPrm) < Precision::PConfusion())
+ {
+ aPrm = aT2;
+ isLast = Standard_True;
+ }
+
+ const gp_Pnt aP = theCrv.Value(aPrm);
+ const Standard_Real aSqD = aP.SquareDistance(aPapx);
+ if (aSqD < aSqTol)
+ {
+ IntAna_Curve aC1 = theCrv;
+ aC1.SetDomain(aT1, aPrm);
+ aT1 = aPrm;
+ theLC.Append(aC1);
+ }
+
+ if (isLast)
+ break;
}
- //
- // need to be splitted at aTheta
- IntAna_Curve aC1, aC2;
- //
- aC1=aC;
- aC1.SetDomain(aT1, aTheta);
- aC2=aC;
- aC2.SetDomain(aTheta, aT2);
- //
- aLC.Clear();
- aLC.Append(aC1);
- aLC.Append(aC2);
- //
- 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 (theLC.IsEmpty())
+ {
+ theLC.Append(theCrv);
+ return Standard_False;
}
- if (dR>Tol) {
- bRet=aR1>aR2;
+
+ if ((aT2 - aT1) > Precision::PConfusion())
+ {
+ IntAna_Curve aC1 = theCrv;
+ aC1.SetDomain(aT1, aT2);
+ theLC.Append(aC1);
}
- 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<aSc1) {
- bRet=!bRet;
- }
- }//if (dR<Tol) {
- return bRet;
+
+ return Standard_True;
}
projects / occt.git / blobdiff