OpenCascade/occt
1
0
Fork 0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-04-04 18:06:22 +03:00
Code Packages Releases Activity
occt/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx
emv 590b3f0416 0031499: Boolean Operations - Custom fuzzy value corrupts the result of CUT 
IntPatch_ImpImpIntersection::CyCyNoGeometric - Use the provided 3D tolerance to compare the points.
2020-04-16 10:01:55 +03:00

4365 lines
142 KiB
Plaintext
Raw Blame History

// Created on: 1992-05-07
// Created by: Jacques GOUSSARD
// Copyright (c) 1992-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.
#include <algorithm>
#include <Bnd_Range.hxx>
#include <IntAna_ListOfCurve.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;
static void ShortCosForm( const Standard_Real theCosFactor,
const Standard_Real theSinFactor,
Standard_Real& theCoeff,
Standard_Real& theAngle);
//
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,
Standard_Real& theUGiven,
const Standard_Real theTol2D,
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 ComputationMethods::stCoeffsValue& theCoeffs,
const Standard_Integer theWLIndex,
const Standard_Integer theMinNbPoints,
const Standard_Integer theStartPointOnLine,
const Standard_Integer theEndPointOnLine,
const Standard_Real theTol2D,
const Standard_Real thePeriodOfSurf2,
const Standard_Boolean isTheReverse);
//=======================================================================
//function : MinMax
//purpose : Replaces theParMIN = MIN(theParMIN, theParMAX),
// theParMAX = MAX(theParMIN, theParMAX).
//=======================================================================
static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX)
{
if(theParMIN > theParMAX)
{
const Standard_Real aux = theParMAX;
theParMAX = theParMIN;
theParMIN = aux;
}
}
//=======================================================================
//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 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,
Standard_Real& theV1Set,
Standard_Real& theV2Set)
{
//Now we are going to define if V1 (and V2) increases
//(or decreases) when U1 will increase.
const Standard_Integer aNbDim = 3;
math_Matrix aSyst(1, aNbDim, 1, aNbDim);
aSyst.SetCol(1, theMatr.Col(1));
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 aDet1 = aSyst.Determinant();
aSyst.SetCol(1, theMatr.Col(1));
aSyst.SetCol(2, theMatr.Col(3));
const Standard_Real aDet2 = aSyst.Determinant();
//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 = theV1AfterDecrByDelta;
}
if(aDet*aDet2 > 0.0)
{
theV2Set = theV2AfterDecrByDelta;
}
}
//=======================================================================
//function : DeltaU1Computing
//purpose : Computes new step for U1 parameter.
//=======================================================================
static inline
Standard_Boolean DeltaU1Computing(const math_Matrix& theSyst,
const math_Vector& theFree,
Standard_Real& theDeltaU1Found)
{
Standard_Real aDet = theSyst.Determinant();
if(Abs(aDet) > aNulValue)
{
math_Matrix aSyst1(theSyst);
aSyst1.SetCol(2, theFree);
theDeltaU1Found = Abs(aSyst1.Determinant()/aDet);
return Standard_True;
}
return Standard_False;
}
//=======================================================================
//function : StepComputing
//purpose :
//
//Attention!!!:
// theMatr must have 3*5-dimension strictly.
// For system
// {a11*V1+a12*V2+a13*dU1+a14*dU2=b1;
// {a21*V1+a22*V2+a23*dU1+a24*dU2=b2;
// {a31*V1+a32*V2+a33*dU1+a34*dU2=b3;
// theMatr must be following:
// (a11 a12 a13 a14 b1)
// (a21 a22 a23 a24 b2)
// (a31 a32 a33 a34 b3)
//=======================================================================
static Standard_Boolean StepComputing(const math_Matrix& theMatr,
const Standard_Real theV1Cur,
const Standard_Real theV2Cur,
const Standard_Real theDeltaV1,
const Standard_Real theDeltaV2,
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;
//theV2Found = theV2Set;
const Standard_Integer aNbDim = 3;
math_Matrix aSyst(1, aNbDim, 1, aNbDim);
math_Vector aFree(1, aNbDim);
//By default, increasing V1(U1) and V2(U1) functions is
//considered
Standard_Real aV1Set = theV1Cur + theDeltaV1,
aV2Set = theV2Cur + theDeltaV2;
//However, what is indeed?
VBoundaryPrecise( theMatr, theV1Cur - theDeltaV1,
theV2Cur - theDeltaV2, aV1Set, aV2Set);
aSyst.SetCol(2, theMatr.Col(3));
aSyst.SetCol(3, theMatr.Col(4));
for(Standard_Integer i = 0; i < 2; i++)
{
if(i == 0)
{//V1 is known
aSyst.SetCol(1, theMatr.Col(2));
aFree.Set(1, aNbDim, theMatr.Col(5)-aV1Set*theMatr.Col(1));
}
else
{//i==1 => V2 is known
aSyst.SetCol(1, theMatr.Col(1));
aFree.Set(1, aNbDim, theMatr.Col(5)-aV2Set*theMatr.Col(2));
}
Standard_Real aNewDU = theDeltaU1Found;
if(DeltaU1Computing(aSyst, aFree, aNewDU))
{
isSuccess = Standard_True;
if(aNewDU < theDeltaU1Found)
{
theDeltaU1Found = aNewDU;
}
}
}
if(!isSuccess)
{
aFree = theMatr.Col(5) - aV1Set*theMatr.Col(1) - aV2Set*theMatr.Col(2);
math_Matrix aSyst1(1, aNbDim, 1, 2);
aSyst1.SetCol(1, aSyst.Col(2));
aSyst1.SetCol(2, aSyst.Col(3));
//Now we have overdetermined system.
const Standard_Real aDet1 = theMatr(1,3)*theMatr(2,4) - theMatr(2,3)*theMatr(1,4);
const Standard_Real aDet2 = theMatr(1,3)*theMatr(3,4) - theMatr(3,3)*theMatr(1,4);
const Standard_Real aDet3 = theMatr(2,3)*theMatr(3,4) - theMatr(3,3)*theMatr(2,4);
const Standard_Real anAbsD1 = Abs(aDet1);
const Standard_Real anAbsD2 = Abs(aDet2);
const Standard_Real anAbsD3 = Abs(aDet3);
if(anAbsD1 >= anAbsD2)
{
if(anAbsD1 >= anAbsD3)
{
//Det1
if(anAbsD1 <= aNulValue)
return isSuccess;
theDeltaU1Found = Abs(aFree(1)*theMatr(2,4) - aFree(2)*theMatr(1,4))/anAbsD1;
isSuccess = Standard_True;
}
else
{
//Det3
if(anAbsD3 <= aNulValue)
return isSuccess;
theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
isSuccess = Standard_True;
}
}
else
{
if(anAbsD2 >= anAbsD3)
{
//Det2
if(anAbsD2 <= aNulValue)
return isSuccess;
theDeltaU1Found = Abs(aFree(1)*theMatr(3,4) - aFree(3)*theMatr(1,4))/anAbsD2;
isSuccess = Standard_True;
}
else
{
//Det3
if(anAbsD3 <= aNulValue)
return isSuccess;
theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
isSuccess = Standard_True;
}
}
}
return isSuccess;
}
//=======================================================================
//function : ProcessBounds
//purpose :
//=======================================================================
void ProcessBounds(const Handle(IntPatch_ALine)& alig, //-- ligne courante
const IntPatch_SequenceOfLine& slin,
const IntSurf_Quadric& Quad1,
const IntSurf_Quadric& Quad2,
Standard_Boolean& procf,
const gp_Pnt& ptf, //-- Debut Ligne Courante
const Standard_Real first, //-- Paramf
Standard_Boolean& procl,
const gp_Pnt& ptl, //-- Fin Ligne courante
const Standard_Real last, //-- Paraml
Standard_Boolean& Multpoint,
const Standard_Real Tol)
{
Standard_Integer j,k;
Standard_Real U1,V1,U2,V2;
IntPatch_Point ptsol;
Standard_Real d;
if (procf && procl) {
j = slin.Length() + 1;
}
else {
j = 1;
}
//-- On prend les lignes deja enregistrees
while (j <= slin.Length()) {
if(slin.Value(j)->ArcType() == IntPatch_Analytic) {
const Handle(IntPatch_ALine)& aligold = *((Handle(IntPatch_ALine)*)&slin.Value(j));
k = 1;
//-- On prend les vertex des lignes deja enregistrees
while (k <= aligold->NbVertex()) {
ptsol = aligold->Vertex(k);
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;
ptsol.SetMultiple(Standard_True);
aligold->Replace(k,ptsol);
}
ptsol.SetParameter(first);
alig->AddVertex(ptsol);
alig->SetFirstPoint(alig->NbVertex());
procf = Standard_True;
//-- On restore le point avec son parametre sur aligold
ptsol = aligold->Vertex(k);
}
}
if (!procl) {
if (ptl.Distance(ptsol.Value()) <= Tol) {
ptsol.SetTolerance(Tol);
if (!ptsol.IsMultiple()) {
Multpoint = Standard_True;
ptsol.SetMultiple(Standard_True);
aligold->Replace(k,ptsol);
}
ptsol.SetParameter(last);
alig->AddVertex(ptsol);
alig->SetLastPoint(alig->NbVertex());
procl = Standard_True;
//-- On restore le point avec son parametre sur aligold
ptsol = aligold->Vertex(k);
}
}
if (procf && procl) {
k = aligold->NbVertex()+1;
}
else {
k = k+1;
}
}
if (procf && procl) {
j = slin.Length()+1;
}
else {
j = j+1;
}
}
}
ptsol.SetTolerance(Tol);
if (!procf && !procl) {
Quad1.Parameters(ptf,U1,V1);
Quad2.Parameters(ptf,U2,V2);
ptsol.SetValue(ptf,Tol,Standard_False);
ptsol.SetParameters(U1,V1,U2,V2);
ptsol.SetParameter(first);
if (ptf.Distance(ptl) <= Tol) {
ptsol.SetMultiple(Standard_True); // a voir
Multpoint = Standard_True; // a voir de meme
alig->AddVertex(ptsol);
alig->SetFirstPoint(alig->NbVertex());
ptsol.SetParameter(last);
alig->AddVertex(ptsol);
alig->SetLastPoint(alig->NbVertex());
}
else {
alig->AddVertex(ptsol);
alig->SetFirstPoint(alig->NbVertex());
Quad1.Parameters(ptl,U1,V1);
Quad2.Parameters(ptl,U2,V2);
ptsol.SetValue(ptl,Tol,Standard_False);
ptsol.SetParameters(U1,V1,U2,V2);
ptsol.SetParameter(last);
alig->AddVertex(ptsol);
alig->SetLastPoint(alig->NbVertex());
}
}
else if (!procf) {
Quad1.Parameters(ptf,U1,V1);
Quad2.Parameters(ptf,U2,V2);
ptsol.SetValue(ptf,Tol,Standard_False);
ptsol.SetParameters(U1,V1,U2,V2);
ptsol.SetParameter(first);
alig->AddVertex(ptsol);
alig->SetFirstPoint(alig->NbVertex());
}
else if (!procl) {
Quad1.Parameters(ptl,U1,V1);
Quad2.Parameters(ptl,U2,V2);
ptsol.SetValue(ptl,Tol,Standard_False);
ptsol.SetParameters(U1,V1,U2,V2);
ptsol.SetParameter(last);
alig->AddVertex(ptsol);
alig->SetLastPoint(alig->NbVertex());
}
}
//=======================================================================
//function : CyCyAnalyticalIntersect
//purpose : Checks if intersection curve is analytical (line, circle, ellipse)
// and returns these curves.
//=======================================================================
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;
IntAna_ResultType typint;
gp_Elips elipsol;
gp_Lin linsol;
gp_Cylinder Cy1(Quad1.Cylinder());
gp_Cylinder Cy2(Quad2.Cylinder());
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(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);
}
break;
case IntAna_Line:
{
gp_Pnt ptref;
if (NbSol == 1)
{ // Cylinders are tangent to each other by line
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;
IntSurf_Situation situcyl2;
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
{
situcyl2 = IntSurf_Outside;
}
if (norm2.Dot(crb2) > 0.)
{
situcyl1 = IntSurf_Inside;
}
else
{
situcyl1 = IntSurf_Outside;
}
}
else
{
if (Cy1.Radius() < Cy2.Radius())
{
if (norm1.Dot(crb1) > 0.)
{
situcyl2 = IntSurf_Inside;
}
else
{
situcyl2 = IntSurf_Outside;
}
if (norm2.Dot(crb2) > 0.)
{
situcyl1 = IntSurf_Outside;
}
else
{
situcyl1 = IntSurf_Inside;
}
}
else
{
if (norm1.Dot(crb1) > 0.)
{
situcyl2 = IntSurf_Outside;
}
else
{
situcyl2 = IntSurf_Inside;
}
if (norm2.Dot(crb2) > 0.)
{
situcyl1 = IntSurf_Inside;
}
else
{
situcyl1 = IntSurf_Outside;
}
}
}
Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_True, situcyl1, situcyl2);
slin.Append(glig);
}
else
{
for (i=1; i <= NbSol; 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)
{
trans1 = IntSurf_Out;
trans2 = IntSurf_In;
}
else if (qwe <-0.00000001)
{
trans1 = IntSurf_In;
trans2 = IntSurf_Out;
}
else
{
trans1=trans2=IntSurf_Undecided;
}
Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
slin.Append(glig);
}
}
}
break;
case IntAna_Ellipse:
{
gp_Vec Tgt;
gp_Pnt ptref;
IntPatch_Point pmult1, pmult2;
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(pttang2,oU1,oV1);
Quad2.Parameters(pttang2,oU2,oV2);
pmult2.SetParameters(oU1,oV1,oU2,oV2);
// on traite la premiere ellipse
//-- Calcul de la Transition de la ligne
ElCLib::D1(0.,elipsol,ptref,Tgt);
//Theoretically, qwe = +/- |Tgt|.
Standard_Real qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
if (qwe>0.00000001)
{
trans1 = IntSurf_Out;
trans2 = IntSurf_In;
}
else if (qwe<-0.00000001)
{
trans1 = IntSurf_In;
trans2 = IntSurf_Out;
}
else
{
trans1=trans2=IntSurf_Undecided;
}
//-- Transition calculee au point 0 -> Trans2 , Trans1
//-- car ici, on devarit calculer en PI
Handle(IntPatch_GLine) glig = new IntPatch_GLine(elipsol,Standard_False,trans2,trans1);
//
{
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);
}
//
pmult1.SetParameter(0.5*M_PI);
glig->AddVertex(pmult1);
//
pmult2.SetParameter(1.5*M_PI);
glig->AddVertex(pmult2);
//
slin.Append(glig);
//-- Transitions calculee au point 0 OK
//
// on traite la deuxieme ellipse
elipsol = theInter.Ellipse(2);
Standard_Real param1 = ElCLib::Parameter(elipsol,pttang1);
Standard_Real param2 = ElCLib::Parameter(elipsol,pttang2);
Standard_Real parampourtransition = 0.0;
if (param1 < param2)
{
pmult1.SetParameter(0.5*M_PI);
pmult2.SetParameter(1.5*M_PI);
parampourtransition = M_PI;
}
else {
pmult1.SetParameter(1.5*M_PI);
pmult2.SetParameter(0.5*M_PI);
parampourtransition = 0.0;
}
//-- 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)
{
trans1 = IntSurf_Out;
trans2 = IntSurf_In;
}
else if(qwe< -0.00000001)
{
trans1 = IntSurf_In;
trans2 = IntSurf_Out;
}
else
{
trans1=trans2=IntSurf_Undecided;
}
//-- La transition a ete calculee sur un point de cette ligne
glig = new IntPatch_GLine(elipsol,Standard_False,trans1,trans2);
//
{
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_Parabola:
case IntAna_Hyperbola:
throw Standard_Failure("IntCyCy(): Wrong intersection type!");
case IntAna_Circle:
// Circle is useful when we will work with trimmed surfaces
// (two cylinders can be tangent by their basises, e.g. circle)
case IntAna_NoGeometricSolution:
default:
return Standard_False;
}
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;
}
}
}
//=======================================================================
//function : CylCylMonotonicity
//purpose : Determines, if U2(U1) function is increasing.
//=======================================================================
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,
//Func. U2(X1) = FI2 + X1;
//Func. X1(X2) = anArccosFactor*X2;
//Func. X2(X3) = acos(X3);
//Func. X3(X4) = B*X4 + C;
//Func. X4(U1) = cos(U1 - FI1).
//
//Consequently,
//U2(X1) is always increasing.
//X1(X2) is increasing if anArccosFactor > 0.0 and
//is decreasing otherwise.
//X2(X3) is always decreasing.
//Therefore, U2(X3) is decreasing if anArccosFactor > 0.0 and
//is increasing otherwise.
//X3(X4) is increasing if B > 0 and is decreasing otherwise.
//X4(U1) is increasing if U1 - FI1 in [PI, 2*PI) and
//is decreasing U1 - FI1 in [0, PI) or if (U1 - FI1 == 2PI).
//After that, we can predict behaviour of U2(U1) function:
//if it is increasing or decreasing.
//For "+/-" sign. If isPlus == TRUE, "+" is chosen, otherwise, we choose "-".
Standard_Boolean isPlus = Standard_False;
switch(theWLIndex)
{
case 0:
isPlus = Standard_True;
break;
case 1:
isPlus = Standard_False;
break;
default:
//throw Standard_Failure("Error. Range Error!!!!");
return Standard_False;
}
Standard_Real aU1Temp = theU1par - theCoeffs.mFI1;
InscribePoint(0, thePeriod, aU1Temp, 0.0, thePeriod, Standard_False);
theIsIncreasing = Standard_True;
if(((M_PI - aU1Temp) < RealSmall()) && (aU1Temp < thePeriod))
{
theIsIncreasing = Standard_False;
}
if(theCoeffs.mB < 0.0)
{
theIsIncreasing = !theIsIncreasing;
}
if(!isPlus)
{
theIsIncreasing = !theIsIncreasing;
}
return Standard_True;
}
//=======================================================================
//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* const theDelta)
{
//This formula is got from some experience and can be changed.
const Standard_Real aTol0 = Min(10.0*Epsilon(1.0)*theCoeffs.mB, aNulValue);
const Standard_Real aTol = 1.0 - aTol0;
if(theWLIndex < 0 || theWLIndex > 1)
return Standard_False;
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;
anArg = 1.0;
}
else if(anArg <= -aTol)
{
if(theDelta)
*theDelta = 0.0;
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 = acos(anArg);
theU2 = theCoeffs.mFI2 + aSign*theU2;
return Standard_True;
}
//=======================================================================
//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,
Standard_Real& theV2)
{
theV1 = theCoeffs.mK21 * sin(theU2) +
theCoeffs.mK11 * sin(theU1) +
theCoeffs.mL21 * cos(theU2) +
theCoeffs.mL11 * cos(theU1) + theCoeffs.mM1;
theV2 = theCoeffs.mK22 * sin(theU2) +
theCoeffs.mK12 * sin(theU1) +
theCoeffs.mL22 * cos(theU2) +
theCoeffs.mL12 * cos(theU1) + theCoeffs.mM2;
return Standard_True;
}
//=======================================================================
//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,
Standard_Real& theV1,
Standard_Real& theV2)
{
if(!CylCylComputeParameters(theU1par, theWLIndex, theCoeffs, theU2))
return Standard_False;
if(!CylCylComputeParameters(theU1par, theU2, theCoeffs, theV1, theV2))
return Standard_False;
return Standard_True;
}
//=======================================================================
//function : SearchOnVBounds
//purpose :
//=======================================================================
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) const
{
const Standard_Integer aNbDim = 3;
const Standard_Real aMaxError = 4.0*M_PI; // two periods
theMainVariableValue = theInitMainVar;
const Standard_Real aTol2 = 1.0e-18;
Standard_Real aMainVarPrev = theInitMainVar, aU2Prev = theInitU2, anOtherVar = theVInit;
//Structure of aMatr:
// C_{1}*U_{1} & C_{2}*U_{2} & C_{3}*V_{*},
//where C_{1}, C_{2} and C_{3} are math_Vector.
math_Matrix aMatr(1, aNbDim, 1, aNbDim);
Standard_Real anError = RealLast();
Standard_Real anErrorPrev = anError;
Standard_Integer aNbIter = 0;
do
{
if(++aNbIter > 1000)
return Standard_False;
const Standard_Real aSinU1 = sin(aMainVarPrev),
aCosU1 = cos(aMainVarPrev),
aSinU2 = sin(aU2Prev),
aCosU2 = cos(aU2Prev);
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, myCoeffs.mVecC2);
aMSum = myCoeffs.mVecC1 * theVzad;
aVecFreeMem -= aMSum;
aMSum += myCoeffs.mVecC2*anOtherVar;
break;
case SearchV2:
aMatr.SetCol(3, myCoeffs.mVecC1);
aMSum = myCoeffs.mVecC2 * theVzad;
aVecFreeMem -= aMSum;
aMSum += myCoeffs.mVecC1*anOtherVar;
break;
default:
return Standard_False;
}
aMatr.SetCol(1, myCoeffs.mVecA1 * aSinU1 - myCoeffs.mVecB1 * aCosU1);
aMatr.SetCol(2, myCoeffs.mVecA2 * aSinU2 - myCoeffs.mVecB2 * aCosU2);
Standard_Real aDetMainSyst = aMatr.Determinant();
if(Abs(aDetMainSyst) < aNulValue)
{
return Standard_False;
}
math_Matrix aM1(aMatr), aM2(aMatr), aM3(aMatr);
aM1.SetCol(1, aVecFreeMem);
aM2.SetCol(2, aVecFreeMem);
aM3.SetCol(3, aVecFreeMem);
const Standard_Real aDetMainVar = aM1.Determinant();
const Standard_Real aDetVar1 = aM2.Determinant();
const Standard_Real aDetVar2 = aM3.Determinant();
Standard_Real aDelta = aDetMainVar/aDetMainSyst-aMainVarPrev;
if(Abs(aDelta) > aMaxError)
return Standard_False;
anError = aDelta*aDelta;
aMainVarPrev += aDelta;
///
aDelta = aDetVar1/aDetMainSyst-aU2Prev;
if(Abs(aDelta) > aMaxError)
return Standard_False;
anError += aDelta*aDelta;
aU2Prev += aDelta;
///
aDelta = aDetVar2/aDetMainSyst-anOtherVar;
anError += aDelta*aDelta;
anOtherVar += aDelta;
if(anError > anErrorPrev)
{//Method diverges. Keep the best result
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);
}
else
{
theMainVariableValue = aMainVarPrev;
}
anErrorPrev = anError;
}
while(anError > aTol2);
theMainVariableValue = aMainVarPrev;
return Standard_True;
}
//=======================================================================
//function : InscribePoint
//purpose : If theFlForce==TRUE theUGiven will be changed forcefully
// even if theUGiven is already inscibed in the boundary
// (if it is possible; i.e. if new theUGiven is inscribed
// in the boundary, too).
//=======================================================================
Standard_Boolean InscribePoint( const Standard_Real theUfTarget,
const Standard_Real theUlTarget,
Standard_Real& theUGiven,
const Standard_Real theTol2D,
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
/*
Utf U Utl
+ * +
*/
if(theFlForce)
{
Standard_Real anUtemp = theUGiven + thePeriod;
if((theUfTarget - anUtemp <= theTol2D) &&
(anUtemp - theUlTarget <= theTol2D))
{
theUGiven = anUtemp;
return Standard_True;
}
anUtemp = theUGiven - thePeriod;
if((theUfTarget - anUtemp <= theTol2D) &&
(anUtemp - theUlTarget <= theTol2D))
{
theUGiven = anUtemp;
}
}
return Standard_True;
}
const Standard_Real aUf = theUfTarget - theTol2D;
const Standard_Real aUl = aUf + thePeriod;
theUGiven = ElCLib::InPeriod(theUGiven, aUf, aUl);
return ((theUfTarget - theUGiven <= theTol2D) &&
(theUGiven - theUlTarget <= theTol2D));
}
//=======================================================================
//function : InscribeInterval
//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,
Bnd_Range &theRange,
const Standard_Real theTol2D,
const Standard_Real thePeriod)
{
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)))
{
theRange.SetVoid();
theRange.Add(anUpar);
theRange.Add(anUpar + aDelta);
}
else
{
if (!theRange.GetMax (anUpar))
{
return Standard_False;
}
if(InscribePoint(theUfTarget, theUlTarget, anUpar,
theTol2D, thePeriod, (Abs(theUfTarget-anUpar) < theTol2D)))
{
theRange.SetVoid();
theRange.Add(anUpar);
theRange.Add(anUpar - aDelta);
}
else
{
return Standard_False;
}
}
return Standard_True;
}
//=======================================================================
//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 Standard_Boolean ExcludeNearElements(Standard_Real theArr[],
const Standard_Integer theNOfMembers,
const Standard_Real theUSurf1f,
const Standard_Real theUSurf1l,
const Standard_Real theTol)
{
Standard_Boolean aRetVal = Standard_False;
for(Standard_Integer i = 1; i < theNOfMembers; i++)
{
Standard_Real &anA = theArr[i],
&anB = theArr[i-1];
//Here, anA >= anB
if(Precision::IsInfinite(anA))
break;
if((anA-anB) < theTol)
{
if((anB != 0.0) && (anB != theUSurf1f) && (anB != theUSurf1l))
anA = (anA + anB)/2.0;
else
anA = anB;
//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 ComputationMethods::stCoeffsValue& theCoeffs,
const Standard_Boolean isTheReverse,
const Standard_Boolean isThePrecise,
const gp_Pnt2d& thePntOnSurf1,
const gp_Pnt2d& thePntOnSurf2,
const Standard_Real theUfSurf1,
const Standard_Real theUlSurf1,
const Standard_Real theUfSurf2,
const Standard_Real theUlSurf2,
const Standard_Real theVfSurf1,
const Standard_Real theVlSurf1,
const Standard_Real theVfSurf2,
const Standard_Real theVlSurf2,
const Standard_Real thePeriodOfSurf1,
const Handle(IntSurf_LineOn2S)& theLine,
const Standard_Integer theWLIndex,
const Standard_Real theTol3D,
const Standard_Real theTol2D,
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, 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))
return Standard_False;
Standard_Real aV1par = thePntOnSurf1.Y();
if((aV1par - theVlSurf1 > theTol2D) || (theVfSurf1 - aV1par > theTol2D))
return Standard_False;
Standard_Real aV2par = thePntOnSurf2.Y();
if((aV2par - theVlSurf2 > theTol2D) || (theVfSurf2 - aV2par > theTol2D))
return Standard_False;
//Get intersection point and add it in the WL
IntSurf_PntOn2S aPnt;
if(isTheReverse)
{
aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
aU2par, aV2par,
aU1par, aV1par);
}
else
{
aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
aU1par, aV1par,
aU2par, aV2par);
}
Standard_Integer aNbPnts = theLine->NbPoints();
if(aNbPnts > 0)
{
Standard_Real aUl = 0.0, aVl = 0.0;
const IntSurf_PntOn2S aPlast = theLine->Value(aNbPnts);
if(isTheReverse)
aPlast.ParametersOnS2(aUl, aVl);
else
aPlast.ParametersOnS1(aUl, aVl);
if(!theFlBefore && (aU1par <= aUl))
{
//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,
//we should reduce tolerance for IsSame checking.
const Standard_Real aDTol = 1.0-Epsilon(1.0);
if(aPnt.IsSame(aPlast, theTol3D*aDTol, theTol2D*aDTol))
{
theLine->RemovePoint(aNbPnts);
}
}
if (theOnlyCheck)
return Standard_True;
theLine->Add(aPnt);
if(!isThePrecise)
return Standard_True;
//Try to precise existing WLine
aNbPnts = theLine->NbPoints();
if(aNbPnts >= 3)
{
Standard_Real aU1 = 0.0, aU2 = 0.0, aU3 = 0.0, aV = 0.0;
if(isTheReverse)
{
theLine->Value(aNbPnts).ParametersOnS2(aU3, aV);
theLine->Value(aNbPnts-1).ParametersOnS2(aU2, aV);
theLine->Value(aNbPnts-2).ParametersOnS2(aU1, aV);
}
else
{
theLine->Value(aNbPnts).ParametersOnS1(aU3, aV);
theLine->Value(aNbPnts-1).ParametersOnS1(aU2, aV);
theLine->Value(aNbPnts-2).ParametersOnS1(aU1, aV);
}
const Standard_Real aStepPrev = aU2-aU1;
const Standard_Real aStep = aU3-aU2;
const Standard_Integer aDeltaStep = RealToInt(aStepPrev/aStep);
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, theTol2D, thePeriodOfSurf1, isTheReverse);
}
}
return Standard_True;
}
//=======================================================================
//function : AddBoundaryPoint
//purpose : Find intersection point on V-boundary.
//=======================================================================
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 theFlForce,
Standard_Boolean& isTheFound1,
Standard_Boolean& isTheFound2) const
{
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;
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};
StPInfo aUVPoint[aSize];
for(Standard_Integer anIDSurf = 0; anIDSurf < 4; anIDSurf+=2)
{
const Standard_Real aVf = (anIDSurf == 0) ? theV1 : theV2,
aVl = (anIDSurf == 0) ? theV1Prev : theV2Prev;
const SearchBoundType aTS = (anIDSurf == 0) ? SearchV1 : SearchV2;
for(Standard_Integer anIDBound = 0; anIDBound < 2; anIDBound++)
{
const Standard_Integer anIndex = anIDSurf+anIDBound;
aUVPoint[anIndex].mySurfID = anIDSurf;
if((Abs(aVf-anArrVzad[anIndex]) > myTol2D) &&
((aVf-anArrVzad[anIndex])*(aVl-anArrVzad[anIndex]) > 0.0))
{
continue;
}
//Segment [aVf, aVl] intersects at least one V-boundary (first or last)
// (in general, case is possible, when aVf > aVl).
// Precise intersection point
const Standard_Boolean aRes = SearchOnVBounds(aTS, anArrVzad[anIndex],
(anIDSurf == 0) ? theV2 : theV1,
theU2, theU1,
aUVPoint[anIndex].myU1);
// aUVPoint[anIndex].myU1 is considered to be nearer to theU1 than
// to theU1+/-Period
if (!aRes || (aUVPoint[anIndex].myU1 >= theU1) ||
(aUVPoint[anIndex].myU1 < theU1Min))
{
//Intersection point is not found or out of the domain
aUVPoint[anIndex].myU1 = RealLast();
continue;
}
else
{
//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;
}
//Point on true V-boundary.
if(aTS == SearchV1)
aV1 = anArrVzad[anIndex];
else //if(aTS[anIndex] == SearchV2)
aV2 = anArrVzad[anIndex];
}
}
}
// Sort with acceding U1-parameter.
std::sort(aUVPoint, aUVPoint+aSize);
isTheFound1 = isTheFound2 = Standard_False;
//Adding found points on boundary in the WLine.
for(Standard_Integer i = 0; i < aSize; i++)
{
if(aUVPoint[i].myU1 == RealLast())
break;
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(aUVPoint[i].mySurfID == 0)
{
isTheFound1 = Standard_True;
}
else
{
isTheFound2 = Standard_True;
}
}
}
//=======================================================================
//function : SeekAdditionalPoints
//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 ComputationMethods::stCoeffsValue& theCoeffs,
const Standard_Integer theWLIndex,
const Standard_Integer theMinNbPoints,
const Standard_Integer theStartPointOnLine,
const Standard_Integer theEndPointOnLine,
const Standard_Real theTol2D,
const Standard_Real thePeriodOfSurf2,
const Standard_Boolean isTheReverse)
{
if(theLine.IsNull())
return;
Standard_Integer aNbPoints = theEndPointOnLine - theStartPointOnLine + 1;
if(aNbPoints >= theMinNbPoints)
{
return;
}
Standard_Real aMinDeltaParam = theTol2D;
{
Standard_Real u1 = 0.0, v1 = 0.0, u2 = 0.0, v2 = 0.0;
if(isTheReverse)
{
theLine->Value(theStartPointOnLine).ParametersOnS2(u1, v1);
theLine->Value(theEndPointOnLine).ParametersOnS2(u2, v2);
}
else
{
theLine->Value(theStartPointOnLine).ParametersOnS1(u1, v1);
theLine->Value(theEndPointOnLine).ParametersOnS1(u2, v2);
}
aMinDeltaParam = Max(Abs(u2 - u1)/IntToReal(theMinNbPoints), aMinDeltaParam);
}
Standard_Integer aLastPointIndex = theEndPointOnLine;
Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0;
Standard_Integer aNbPointsPrev = 0;
while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
{
aNbPointsPrev = aNbPoints;
for(Standard_Integer fp = theStartPointOnLine, lp = 0; fp < aLastPointIndex; fp = lp + 1)
{
Standard_Real U1f = 0.0, V1f = 0.0; //first point in 1st suraface
Standard_Real U1l = 0.0, V1l = 0.0; //last point in 1st suraface
Standard_Real U2f = 0.0, V2f = 0.0; //first point in 2nd suraface
Standard_Real U2l = 0.0, V2l = 0.0; //last point in 2nd suraface
lp = fp+1;
if(isTheReverse)
{
theLine->Value(fp).ParametersOnS2(U1f, V1f);
theLine->Value(lp).ParametersOnS2(U1l, V1l);
theLine->Value(fp).ParametersOnS1(U2f, V2f);
theLine->Value(lp).ParametersOnS1(U2l, V2l);
}
else
{
theLine->Value(fp).ParametersOnS1(U1f, V1f);
theLine->Value(lp).ParametersOnS1(U1l, V1l);
theLine->Value(fp).ParametersOnS2(U2f, V2f);
theLine->Value(lp).ParametersOnS2(U2l, V2l);
}
if(Abs(U1l - U1f) <= aMinDeltaParam)
{
//Step is minimal. It is not necessary to divide it.
continue;
}
U1prec = 0.5*(U1f+U1l);
if(!ComputationMethods::
CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec))
{
continue;
}
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 INTPATCH_IMPIMPINTERSECTION_DEBUG
std::cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << std::endl;
#endif
IntSurf_PntOn2S anIP;
if(isTheReverse)
{
anIP.SetValue(aPInt, U2prec, V2prec, U1prec, V1prec);
}
else
{
anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec);
}
theLine->InsertBefore(lp, anIP);
aNbPoints++;
aLastPointIndex++;
}
if(aNbPoints >= theMinNbPoints)
{
return;
}
}
}
//=======================================================================
//function : BoundariesComputing
//purpose : Computes true domain of future intersection curve (allows
// avoiding excess iterations)
//=======================================================================
Standard_Boolean WorkWithBoundaries::
BoundariesComputing(const ComputationMethods::stCoeffsValue &theCoeffs,
const Standard_Real thePeriod,
Bnd_Range theURange[])
{
//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
if (theCoeffs.mB > 0.0)
{
// -(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 solution
return Standard_False;
}
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 + theCoeffs.mC <= theCoeffs.mB) &&
(theCoeffs.mB <= 1 - theCoeffs.mC))
{
//(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);
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 - theCoeffs.mC <= theCoeffs.mB) &&
(theCoeffs.mB <= 1 + theCoeffs.mC))
{
//(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);
theURange[0].Add(aDAngle + theCoeffs.mFI1);
theURange[0].Add(thePeriod - aDAngle + theCoeffs.mFI1);
}
else if (theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
{
//(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)
anArg1 = -1.0;
if(anArg2 > 1.0)
anArg2 = 1.0;
if(anArg2 < -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)
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
{
return Standard_False;
}
}
else if (theCoeffs.mB < 0.0)
{
// (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 (-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 ((-theCoeffs.mC - 1 <= theCoeffs.mB) &&
(theCoeffs.mB <= theCoeffs.mC - 1))
{
// -(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);
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 ((theCoeffs.mC - 1 <= theCoeffs.mB) &&
(theCoeffs.mB <= -theCoeffs.mB - 1))
{
// -(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);
theURange[0].Add(aDAngle + theCoeffs.mFI1);
theURange[0].Add(thePeriod - aDAngle + theCoeffs.mFI1);
}
else if (-theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
{
// -(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)
anArg1 = -1.0;
if(anArg2 > 1.0)
anArg2 = 1.0;
if(anArg2 < -1.0)
anArg2 = -1.0;
const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
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
{
return Standard_False;
}
}
else
{
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;
}
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(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
do
{
std::sort(theU1crit, theU1crit + theNbCritPointsMax);
}
while(ExcludeNearElements(theU1crit, theNbCritPointsMax,
theUSurf1f, theUSurf1l, theTol2D));
//Here all not infinite elements in theU1crit are different and sorted.
while(theNbCritPointsMax > 0)
{
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
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;
}
//Attention! Here theU1crit may be unsorted.
}
//=======================================================================
//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,
// 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
};
const Standard_Integer aNbWLines = 2;
for (Standard_Integer aCurInterval = 0; aCurInterval < theNbOfRanges; aCurInterval++)
{
//Process every continuous region
Standard_Boolean isAddedIntoWL[aNbWLines];
for (Standard_Integer i = 0; i < aNbWLines; i++)
isAddedIntoWL[i] = Standard_False;
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
for (Standard_Integer i = 0; i < aNbCritPoints; i++)
{
InscribePoint(anUf, anUl, anU1crit[i], 0.0, aPeriod, Standard_False);
}
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];
Standard_Real anUexpect[aNbWLines];
Standard_Boolean isAddingWLEnabled[aNbWLines];
Handle(IntSurf_LineOn2S) aL2S[aNbWLines];
Handle(IntPatch_WLine) aWLine[aNbWLines];
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 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)
{
//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)
{
//WL has gone through i-th critical point
anU1 = anU1crit[i];
for (Standard_Integer j = 0; j < aNbWLines; j++)
{
aWLFindStatus[j] = WLFStatus_Broken;
anUexpect[j] = anU1;
}
break;
}
}
if (IsEqual(anU1, anUl))
{
for (Standard_Integer i = 0; i < aNbWLines; i++)
{
aWLFindStatus[i] = WLFStatus_Broken;
anUexpect[i] = anU1;
if (isDeltaPeriod)
{
//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
//points found earlier. Therefore, new WLine will repeat
//already existing WLine. Consequently, it is necessary
//to forbid building new line in this case.
isAddingWLEnabled[i] = (!isAddedIntoWL[i]);
}
else
{
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++)
{
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++)
{
const Standard_Integer aNbPntsWL = aWLine[i].IsNull() ? 0 :
aWLine[i]->Curve()->NbPoints();
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);
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)
{//the line has not contained any points yet
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;
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)
{
aU2[i] = aUSurf2f;
}
else
{
aU2[i] = aUSurf2l;
}
}
}
else
{//aNbPntsWL > 0
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 (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 (aU2prev > aU2[i])
aU2[i] += aPeriod;
else
aU2[i] -= aPeriod;
}
}
}
ComputationMethods::CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs,
aV1[i], aV2[i]);
if (isFirst)
{
aV1Prev[i] = aV1[i];
aV2Prev[i] = aV2[i];
}
}//for(Standard_Integer i = 0; i < aNbWLines; i++)
isFirst = Standard_False;
//Looking for points into WLine
Standard_Boolean isBroken = Standard_False;
for (Standard_Integer i = 0; i < aNbWLines; i++)
{
if (!isAddingWLEnabled[i])
{
Standard_Boolean isBoundIntersect = Standard_False;
if ((Abs(aV1[i] - aVSurf1f) <= aTol2D) ||
((aV1[i] - aVSurf1f)*(aV1Prev[i] - aVSurf1f) < 0.0))
{
isBoundIntersect = Standard_True;
}
else if ((Abs(aV1[i] - aVSurf1l) <= aTol2D) ||
((aV1[i] - aVSurf1l)*(aV1Prev[i] - aVSurf1l) < 0.0))
{
isBoundIntersect = Standard_True;
}
else if ((Abs(aV2[i] - aVSurf2f) <= aTol2D) ||
((aV2[i] - aVSurf2f)*(aV2Prev[i] - aVSurf2f) < 0.0))
{
isBoundIntersect = Standard_True;
}
else if ((Abs(aV2[i] - aVSurf2l) <= aTol2D) ||
((aV2[i] - aVSurf2l)*(aV2Prev[i] - aVSurf2l) < 0.0))
{
isBoundIntersect = Standard_True;
}
if (aWLFindStatus[i] == WLFStatus_Broken)
isBroken = Standard_True;
if (!isBoundIntersect)
{
continue;
}
else
{
anUexpect[i] = anU1;
}
}
// 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()));
//isFound1 == TRUE if intersection line intersects V-bounds
// (First or Last - no matter) of the 1st cylynder
//isFound2 == TRUE if intersection line intersects V-bounds
// (First or Last - no matter) of the 2nd cylynder
Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
Standard_Boolean isForce = Standard_False;
if (aWLFindStatus[i] == WLFStatus_Absent)
{
if (((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1 - aUSurf1l) < aTol2D))
{
isForce = Standard_True;
}
}
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) <= 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)
{
aWLFindStatus[i] = WLFStatus_Broken; //start a new line
}
else if (isInscribe)
{
if ((aWLFindStatus[i] == WLFStatus_Absent) && (isFound1 || isFound2))
{
aWLFindStatus[i] = WLFStatus_Exist;
}
if ((aWLFindStatus[i] != WLFStatus_Broken) ||
(aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl))
{
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, aTol3D, aTol2D, isForce))
{
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;
}
}
if (aWLFindStatus[i] == WLFStatus_Broken)
isBroken = Standard_True;
}//for(Standard_Integer i = 0; i < aNbWLines; i++)
if (isBroken)
{//current lines are filled. Go to the next lines
anUf = anU1;
Standard_Boolean isAdded = Standard_True;
for (Standard_Integer i = 0; i < aNbWLines; i++)
{
if (isAddingWLEnabled[i])
{
continue;
}
isAdded = Standard_False;
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;
}
}
}
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;
}
}
if (!isAdded)
{
//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();
{
Standard_Boolean isChanged = Standard_False;
for (Standard_Integer i = 0; i < aNbWLines; i++)
{
if ((aWLFindStatus[i] == WLFStatus_Absent) || (aWLine[i]->NbPnts() == 0))
continue;
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);
anUmaxAdded = Max(anUmaxAdded, aU1c);
isChanged = Standard_True;
}
if (!isChanged)
{ //If anUmaxAdded were not changed in previous cycle then
//we would break existing WLines.
break;
}
}
for (Standard_Integer i = 0; i < aNbWLines; i++)
{
if (isAddingWLEnabled[i])
{
continue;
}
ComputationMethods::CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs,
aU2[i], aV1[i], aV2[i]);
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
{
//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++)
{
if (aTol2D < (anUexpect[i] - anU1))
{
continue;
}
if (aWLFindStatus[i] == WLFStatus_Absent)
{
anUexpect[i] += aStepMax;
aMinUexp = Min(aMinUexp, anUexpect[i]);
continue;
}
Standard_Real aStepTmp = aStepMax;
const Standard_Real aSinU1 = sin(anU1),
aCosU1 = cos(anU1),
aSinU2 = sin(aU2[i]),
aCosU2 = cos(aU2[i]);
aMatr.SetCol(1, anEquationCoeffs.mVecC1);
aMatr.SetCol(2, anEquationCoeffs.mVecC2);
aMatr.SetCol(3, anEquationCoeffs.mVecA1*aSinU1 - anEquationCoeffs.mVecB1*aCosU1);
aMatr.SetCol(4, anEquationCoeffs.mVecA2*aSinU2 - anEquationCoeffs.mVecB2*aCosU2);
aMatr.SetCol(5, anEquationCoeffs.mVecA1*aCosU1 + anEquationCoeffs.mVecB1*aSinU1 +
anEquationCoeffs.mVecA2*aCosU2 + anEquationCoeffs.mVecB2*aSinU2 +
anEquationCoeffs.mVecD);
//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;
aMinUexp = Min(aMinUexp, anUexpect[i]);
continue;
}
if (aStepTmp < aStepMin)
aStepTmp = aStepMin;
if (aStepTmp > aStepMax)
aStepTmp = aStepMax;
anUexpect[i] = anU1 + aStepTmp;
aMinUexp = Min(aMinUexp, anUexpect[i]);
}
anU1 = aMinUexp;
}
if (Precision::PConfusion() >= (anUl - anU1))
anU1 = anUl;
anUf = anU1;
for (Standard_Integer i = 0; i < aNbWLines; i++)
{
if (aWLine[i]->NbPnts() != 1)
isAddedIntoWL[i] = Standard_False;
if (anU1 == anUl)
{//strictly equal. Tolerance is considered above.
anUexpect[i] = anUl;
}
}
}
for (Standard_Integer i = 0; i < aNbWLines; i++)
{
if ((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
{
isTheEmpty = Standard_False;
Standard_Real u1, v1, u2, v2;
aWLine[i]->Point(1).Parameters(u1, v1, u2, v2);
IntPatch_Point aP;
aP.SetParameter(u1);
aP.SetParameters(u1, v1, u2, v2);
aP.SetTolerance(aTol3D);
aP.SetValue(aWLine[i]->Point(1).Value());
//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)
{
Standard_Boolean isGood = Standard_True;
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()))
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)
{
isTheEmpty = Standard_False;
isAddedIntoWL[i] = Standard_True;
SeekAdditionalPoints(aQuad1, aQuad2, aWLine[i]->Curve(),
anEquationCoeffs, i, aNbPoints, 1,
aWLine[i]->NbPnts(), aTol2D, aPeriod,
isReversed);
aWLine[i]->ComputeVertexParameters(aTol3D);
theSlin.Append(aWLine[i]);
}
}
else
{
isAddedIntoWL[i] = Standard_False;
}
#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++)
{
for (Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); 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, aTol3D) ||
aPntCur.IsSame(aPntLWL1, aTol3D))
{
theSPnt.Remove(aNbPnt);
aNbPnt--;
break;
}
}
}
//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;
aPnt2S.Parameters(u1, v1, u2, v2);
Handle(IntSurf_LineOn2S) aL2S = new IntSurf_LineOn2S();
Handle(IntPatch_WLine) aWLine = new IntPatch_WLine(aL2S, Standard_False);
aWLine->SetCreatingWayInfo(IntPatch_WLine::IntPatch_WLImpImp);
//Define the index of WLine, which lies the point aPnt2S in.
Standard_Integer anIndex = 0;
Standard_Real anUf = 0.0, anUl = 0.0, aCurU2 = 0.0;
if (isReversed)
{
anUf = Max(u2 - aStepMax, aUSurf1f);
anUl = Min(u2 + aStepMax, aUSurf1l);
aCurU2 = u1;
}
else
{
anUf = Max(u1 - aStepMax, aUSurf1f);
anUl = Min(u1 + aStepMax, aUSurf1l);
aCurU2 = u2;
}
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++)
{
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;
}
}
}
// 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++)
{
if (aParID == 0)
{
anUf = anUinf;
anUl = anUmid;
}
else // if(aParID == 1)
{
anUf = anUmid;
anUl = anUsup;
}
Standard_Real &aPar1 = (aParID == 0) ? anUf : anUl,
&aPar2 = (aParID == 0) ? anUl : anUf;
while (Abs(aPar2 - aPar1) > aStepMin)
{
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;
}
}
}
//Fill aWLine by additional points
if (anAddedPar[1] - anAddedPar[0] > aStepMin)
{
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(), aTol2D, aPeriod,
isReversed);
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;
}
}
//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);
}
}
//=======================================================================
//function : IntCySp
//purpose :
//=======================================================================
Standard_Boolean IntCySp(const IntSurf_Quadric& Quad1,
const IntSurf_Quadric& Quad2,
const Standard_Real Tol,
const Standard_Boolean Reversed,
Standard_Boolean& Empty,
Standard_Boolean& Multpoint,
IntPatch_SequenceOfLine& slin,
IntPatch_SequenceOfPoint& spnt)
{
Standard_Integer i;
IntSurf_TypeTrans trans1,trans2;
IntAna_ResultType typint;
IntPatch_Point ptsol;
gp_Circ cirsol;
gp_Cylinder Cy;
gp_Sphere Sp;
if (!Reversed) {
Cy = Quad1.Cylinder();
Sp = Quad2.Sphere();
}
else {
Cy = Quad2.Cylinder();
Sp = Quad1.Sphere();
}
IntAna_QuadQuadGeo inter(Cy,Sp,Tol);
if (!inter.IsDone()) {return Standard_False;}
typint = inter.TypeInter();
Standard_Integer NbSol = inter.NbSolutions();
Empty = Standard_False;
switch (typint) {
case IntAna_Empty :
{
Empty = 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);
ptsol.SetValue(psol,Tol,Standard_True);
ptsol.SetParameters(U1,V1,U2,V2);
spnt.Append(ptsol);
}
break;
case IntAna_Circle:
{
cirsol = inter.Circle(1);
gp_Vec Tgt;
gp_Pnt ptref;
ElCLib::D1(0.,cirsol,ptref,Tgt);
if (NbSol == 1) {
gp_Vec TestCurvature(ptref,Sp.Location());
gp_Vec Normsp,Normcyl;
if (!Reversed) {
Normcyl = Quad1.Normale(ptref);
Normsp = Quad2.Normale(ptref);
}
else {
Normcyl = Quad2.Normale(ptref);
Normsp = Quad1.Normale(ptref);
}
IntSurf_Situation situcyl;
IntSurf_Situation situsp;
if (Normcyl.Dot(TestCurvature) > 0.) {
situsp = IntSurf_Outside;
if (Normsp.Dot(Normcyl) > 0.) {
situcyl = IntSurf_Inside;
}
else {
situcyl = IntSurf_Outside;
}
}
else {
situsp = IntSurf_Inside;
if (Normsp.Dot(Normcyl) > 0.) {
situcyl = IntSurf_Outside;
}
else {
situcyl = IntSurf_Inside;
}
}
Handle(IntPatch_GLine) glig;
if (!Reversed) {
glig = new IntPatch_GLine(cirsol, Standard_True, situcyl, situsp);
}
else {
glig = new IntPatch_GLine(cirsol, Standard_True, situsp, situcyl);
}
slin.Append(glig);
}
else {
if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
trans1 = IntSurf_Out;
trans2 = IntSurf_In;
}
else {
trans1 = IntSurf_In;
trans2 = IntSurf_Out;
}
Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
slin.Append(glig);
cirsol = inter.Circle(2);
ElCLib::D1(0.,cirsol,ptref,Tgt);
Standard_Real qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
if(qwe> 0.0000001) {
trans1 = IntSurf_Out;
trans2 = IntSurf_In;
}
else if(qwe<-0.0000001) {
trans1 = IntSurf_In;
trans2 = IntSurf_Out;
}
else {
trans1=trans2=IntSurf_Undecided;
}
glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
slin.Append(glig);
}
}
break;
case IntAna_NoGeometricSolution:
{
gp_Pnt psol;
Standard_Real U1,V1,U2,V2;
IntAna_IntQuadQuad anaint(Cy,Sp,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);
}
Standard_Boolean TempFalse1a = Standard_False;
Standard_Boolean TempFalse2a = Standard_False;
//-- ptf et ptl : points debut et fin de alig
ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1a,ptf,first,
TempFalse2a,ptl,last,Multpoint,Tol);
slin.Append(alig);
} //-- boucle sur les lignes
} //-- solution avec au moins une lihne
}
break;
default:
{
return Standard_False;
}
}
return Standard_True;
}
//=======================================================================
//function : IntCyCo
//purpose :
//=======================================================================
Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
const IntSurf_Quadric& Quad2,
const Standard_Real Tol,
const Standard_Boolean Reversed,
Standard_Boolean& Empty,
Standard_Boolean& Multpoint,
IntPatch_SequenceOfLine& slin,
IntPatch_SequenceOfPoint& spnt)
{
IntPatch_Point ptsol;
Standard_Integer i;
IntSurf_TypeTrans trans1,trans2;
IntAna_ResultType typint;
gp_Circ cirsol;
gp_Cylinder Cy;
gp_Cone Co;
if (!Reversed) {
Cy = Quad1.Cylinder();
Co = Quad2.Cone();
}
else {
Cy = Quad2.Cylinder();
Co = Quad1.Cone();
}
IntAna_QuadQuadGeo inter(Cy,Co,Tol);
if (!inter.IsDone()) {return Standard_False;}
typint = inter.TypeInter();
Standard_Integer NbSol = inter.NbSolutions();
Empty = Standard_False;
switch (typint) {
case IntAna_Empty : {
Empty = Standard_True;
}
break;
case IntAna_Point :{
gp_Pnt psol(inter.Point(1));
Standard_Real U1,V1,U2,V2;
Quad1.Parameters(psol,U1,V1);
Quad1.Parameters(psol,U2,V2);
ptsol.SetValue(psol,Tol,Standard_True);
ptsol.SetParameters(U1,V1,U2,V2);
spnt.Append(ptsol);
}
break;
case IntAna_Circle: {
gp_Vec Tgt;
gp_Pnt ptref;
Standard_Integer j;
Standard_Real qwe;
//
for(j=1; j<=2; ++j) {
cirsol = inter.Circle(j);
ElCLib::D1(0.,cirsol,ptref,Tgt);
qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
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;
}
Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
slin.Append(glig);
}
}
break;
case IntAna_NoGeometricSolution: {
gp_Pnt psol;
Standard_Real U1,V1,U2,V2;
IntAna_IntQuadQuad anaint(Cy,Co,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;
Standard_Boolean tgfound,firstp,lastp,kept;
Standard_Integer kount;
//
//
//IntAna_Curve curvsol;
IntAna_Curve aC;
IntAna_ListOfCurve aLC;
IntAna_ListIteratorOfListOfCurve aIt;
//
NbSol = anaint.NbCurve();
for (i = 1; i <= NbSol; ++i) {
kept = Standard_False;
//curvsol = anaint.Curve(i);
aC=anaint.Curve(i);
aLC.Clear();
ExploreCurve(Co, aC, 10.*Tol, aLC);
//
aIt.Initialize(aLC);
for (; aIt.More(); aIt.Next()) {
IntAna_Curve& curvsol=aIt.Value();
//
curvsol.Domain(first, last);
firstp = !curvsol.IsFirstOpen();
lastp = !curvsol.IsLastOpen();
if (firstp) {
ptf = curvsol.Value(first);
}
if (lastp) {
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);
kept = Standard_True;
}
else {
ptvalid = curvsol.Value(para);
alig = new IntPatch_ALine(curvsol,Standard_False);
kept = Standard_True;
//-- std::cout << "Transition indeterminee" << std::endl;
}
if (kept) {
Standard_Boolean Nfirstp = !firstp;
Standard_Boolean Nlastp = !lastp;
ProcessBounds(alig,slin,Quad1,Quad2,Nfirstp,ptf,first,
Nlastp,ptl,last,Multpoint,Tol);
slin.Append(alig);
}
} // for (; aIt.More(); aIt.Next())
} // for (i = 1; i <= NbSol; ++i)
}
}
break;
default:
return Standard_False;
} // switch (typint)
return Standard_True;
}
//=======================================================================
//function : ExploreCurve
//purpose : Splits aC on several curves in the cone apex points.
//=======================================================================
Standard_Boolean ExploreCurve(const gp_Cone& theCo,
IntAna_Curve& theCrv,
const Standard_Real theTol,
IntAna_ListOfCurve& theLC)
{
const Standard_Real aSqTol = theTol*theTol;
const gp_Pnt aPapx(theCo.Apex());
Standard_Real aT1, aT2;
theCrv.Domain(aT1, aT2);
theLC.Clear();
//
TColStd_ListOfReal aLParams;
theCrv.FindParameter(aPapx, aLParams);
if (aLParams.IsEmpty())
{
theLC.Append(theCrv);
return Standard_False;
}
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;
}
if (theLC.IsEmpty())
{
theLC.Append(theCrv);
return Standard_False;
}
if ((aT2 - aT1) > Precision::PConfusion())
{
IntAna_Curve aC1 = theCrv;
aC1.SetDomain(aT1, aT2);
theLC.Append(aC1);
}
return Standard_True;
}
View Git Blame Copy Permalink