diff --git a/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx b/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx index b2582bd52f..4427763064 100644 --- a/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx +++ b/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx @@ -47,6 +47,44 @@ static Standard_Boolean InscribePoint(const Standard_Real theUfTarget, 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 @@ -397,6 +435,10 @@ public: 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, @@ -459,7 +501,7 @@ static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX) //======================================================================= //function : ExtremaLineLine //purpose : Computes extrema between the given lines. Returns parameters -// on correspond curve. +// on correspond curve (see correspond method for Extrema_ExtElC class). //======================================================================= static inline void ExtremaLineLine(const gp_Ax1& theC1, const gp_Ax1& theC2, @@ -481,8 +523,8 @@ static inline void ExtremaLineLine(const gp_Ax1& theC1, //======================================================================= //function : VBoundaryPrecise -//purpose : By default, we shall consider, that V1 and V2 will increase -// if U1 increases. But if it is not, new V1set and/or V2set +//purpose : By default, we shall consider, that V1 and V2 will be increased +// if U1 is increased. But if it is not, new V1set and/or V2set // must be computed as [V1current - DeltaV1] (analogically // for V2). This function processes this case. //======================================================================= @@ -501,6 +543,11 @@ static void VBoundaryPrecise( const math_Matrix& theMatr, 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)); @@ -511,6 +558,16 @@ static void VBoundaryPrecise( const math_Matrix& theMatr, 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; @@ -834,7 +891,8 @@ void ProcessBounds(const Handle(IntPatch_ALine)& alig, //-- ligne coura //======================================================================= //function : CyCyAnalyticalIntersect -//purpose : +//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, @@ -1670,7 +1728,7 @@ Standard_Boolean InscribePoint( const Standard_Real theUfTarget, //======================================================================= //function : InscribeInterval -//purpose : Adjusts theUfGiven and after that fits theUlGiven to result +//purpose : Adjusts theUfGiven and (after that) adjusts theUlGiven to the result //======================================================================= static Standard_Boolean InscribeInterval(const Standard_Real theUfTarget, const Standard_Real theUlTarget, @@ -1777,6 +1835,8 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1, const Standard_Boolean theFlForce, const Standard_Boolean theFlBefore = 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())); @@ -1798,6 +1858,7 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1, if((aV2par - theVlSurf2 > theTol2D) || (theVfSurf2 - aV2par > theTol2D)) return Standard_False; + //Get intersection point and add it in the WL IntSurf_PntOn2S aPnt; if(isTheReverse) @@ -1882,6 +1943,7 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1, 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); @@ -1893,7 +1955,7 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1, //======================================================================= //function : AddBoundaryPoint -//purpose : +//purpose : Find intersection point on V-boundary. //======================================================================= void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL, const Standard_Real theU1, @@ -1944,6 +2006,10 @@ void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL, 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, @@ -1951,11 +2017,14 @@ void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL, if(!aRes || aUVPoint[anIndex].myU1 >= theU1) { + //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, @@ -1967,10 +2036,12 @@ void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL, 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) @@ -1979,10 +2050,12 @@ void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL, } } + // 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()) @@ -2011,7 +2084,10 @@ void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL, //======================================================================= //function : SeekAdditionalPoints -//purpose : +//purpose : Inserts additional intersection points between neighbor points. +// This action is bone with several iterations. During every iteration, +// new point is inserted in middle of every interval. +// The process will be finished if NbPoints >= theMinNbPoints. //======================================================================= static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, @@ -2145,20 +2221,38 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1, Standard_Boolean WorkWithBoundaries::BoundariesComputing(Standard_Real theU1f[], Standard_Real theU1l[]) const { + //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(myCoeffs.mB > 0.0) { + // -(1+C)/B <= cos(U1-FI1) <= (1-C)/B + if(myCoeffs.mB + Abs(myCoeffs.mC) < -1.0) - {//There is NOT intersection + { + //(1-C)/B < -1 or -(1+C)/B > 1 ==> No solution + return Standard_False; } else if(myCoeffs.mB + Abs(myCoeffs.mC) <= 1.0) - {//U=[0;2*PI]+aFI1 + { + //(1-C)/B >= 1 and -(1+C)/B <= -1 ==> U=[0;2*PI]+aFI1 + theU1f[0] = myCoeffs.mFI1; theU1l[0] = myPeriod + myCoeffs.mFI1; } else if((1 + myCoeffs.mC <= myCoeffs.mB) && (myCoeffs.mB <= 1 - myCoeffs.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 = -(myCoeffs.mC + 1) / myCoeffs.mB; if(anArg > 1.0) anArg = 1.0; @@ -2166,7 +2260,6 @@ Standard_Boolean WorkWithBoundaries::BoundariesComputing(Standard_Real theU1f[], anArg = -1.0; const Standard_Real aDAngle = acos(anArg); - //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1) theU1f[0] = myCoeffs.mFI1; theU1l[0] = aDAngle + myCoeffs.mFI1; theU1f[1] = myPeriod - aDAngle + myCoeffs.mFI1; @@ -2175,6 +2268,9 @@ Standard_Boolean WorkWithBoundaries::BoundariesComputing(Standard_Real theU1f[], else if((1 - myCoeffs.mC <= myCoeffs.mB) && (myCoeffs.mB <= 1 + myCoeffs.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 - myCoeffs.mC) / myCoeffs.mB; if(anArg > 1.0) anArg = 1.0; @@ -2182,13 +2278,17 @@ Standard_Boolean WorkWithBoundaries::BoundariesComputing(Standard_Real theU1f[], anArg = -1.0; const Standard_Real aDAngle = acos(anArg); - //U=[aDAngle;2*PI-aDAngle]+aFI1 - theU1f[0] = aDAngle + myCoeffs.mFI1; theU1l[0] = myPeriod - aDAngle + myCoeffs.mFI1; } else if(myCoeffs.mB - Abs(myCoeffs.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 - myCoeffs.mC) / myCoeffs.mB, anArg2 = -(myCoeffs.mC + 1) / myCoeffs.mB; if(anArg1 > 1.0) @@ -2217,18 +2317,27 @@ Standard_Boolean WorkWithBoundaries::BoundariesComputing(Standard_Real theU1f[], } else if(myCoeffs.mB < 0.0) { + // (1-C)/B <= cos(U1-FI1) <= -(1+C)/B + if(myCoeffs.mB + Abs(myCoeffs.mC) > 1.0) - {//There is NOT intersection + { + // -(1+C)/B < -1 or (1-C)/B > 1 ==> No solutions return Standard_False; } else if(-myCoeffs.mB + Abs(myCoeffs.mC) <= 1.0) - {//U=[0;2*PI]+aFI1 + { + // -(1+C)/B >= 1 and (1-C)/B <= -1 ==> U=[0;2*PI]+aFI1 + theU1f[0] = myCoeffs.mFI1; theU1l[0] = myPeriod + myCoeffs.mFI1; } else if((-myCoeffs.mC - 1 <= myCoeffs.mB) && ( myCoeffs.mB <= myCoeffs.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 - myCoeffs.mC) / myCoeffs.mB; if(anArg > 1.0) anArg = 1.0; @@ -2236,8 +2345,6 @@ Standard_Boolean WorkWithBoundaries::BoundariesComputing(Standard_Real theU1f[], anArg = -1.0; const Standard_Real aDAngle = acos(anArg); - //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1) - theU1f[0] = myCoeffs.mFI1; theU1l[0] = aDAngle + myCoeffs.mFI1; theU1f[1] = myPeriod - aDAngle + myCoeffs.mFI1; @@ -2246,6 +2353,9 @@ Standard_Boolean WorkWithBoundaries::BoundariesComputing(Standard_Real theU1f[], else if((myCoeffs.mC - 1 <= myCoeffs.mB) && (myCoeffs.mB <= -myCoeffs.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 = -(myCoeffs.mC + 1) / myCoeffs.mB; if(anArg > 1.0) anArg = 1.0; @@ -2253,13 +2363,16 @@ Standard_Boolean WorkWithBoundaries::BoundariesComputing(Standard_Real theU1f[], anArg = -1.0; const Standard_Real aDAngle = acos(anArg); - //U=[aDAngle;2*PI-aDAngle]+aFI1 - theU1f[0] = aDAngle + myCoeffs.mFI1; theU1l[0] = myPeriod - aDAngle + myCoeffs.mFI1; } else if(-myCoeffs.mB - Abs(myCoeffs.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 = -(myCoeffs.mC + 1) / myCoeffs.mB, anArg2 = (1 - myCoeffs.mC) / myCoeffs.mB; if(anArg1 > 1.0) @@ -2273,9 +2386,6 @@ Standard_Boolean WorkWithBoundaries::BoundariesComputing(Standard_Real theU1f[], anArg2 = -1.0; const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2); - //(U=[aDAngle1;aDAngle2]+aFI1) || - //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1) - theU1f[0] = aDAngle1 + myCoeffs.mFI1; theU1l[0] = aDAngle2 + myCoeffs.mFI1; theU1f[1] = myPeriod - aDAngle2 + myCoeffs.mFI1; @@ -2528,6 +2638,8 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1, } } + //Here, intersection line is not an analytical curve(line, circle, ellipsis etc.) + Standard_Real aUSurf1f = 0.0, //const aUSurf1l = 0.0, aVSurf1f = 0.0, @@ -2576,7 +2688,9 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1, return IntPatch_ImpImpIntersection::IntStatus_InfiniteSectionCurve; } - //Boundaries + //Boundaries (it is good idea to use Bnd_Range in the future, instead of aU1f and aU1l) + //Intersection result can include two non-connected regions + //(see WorkWithBoundaries::BoundariesComputing(...) method). const Standard_Integer aNbOfBoundaries = 2; Standard_Real aU1f[aNbOfBoundaries] = {-Precision::Infinite(), -Precision::Infinite()}; Standard_Real aU1l[aNbOfBoundaries] = {Precision::Infinite(), Precision::Infinite()}; @@ -2584,6 +2698,10 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1, if(!aBoundWork.BoundariesComputing(aU1f, aU1l)) return IntPatch_ImpImpIntersection::IntStatus_OK; + //The main idea of the algorithm is to change U1-parameter + //(U-parameter of aCyl1) from aU1f to aU1l with some step + //(step is adaptive) and to obtain set of intersection points. + for(Standard_Integer i = 0; i < aNbOfBoundaries; i++) { if(Precision::IsInfinite(aU1f[i]) && Precision::IsInfinite(aU1l[i])) @@ -2606,7 +2724,14 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1, } } - //Critical points + //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(), @@ -2621,6 +2746,11 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1, 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, theTol2D, aNbCritPoints, anU1crit); @@ -2629,13 +2759,19 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1, 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 }; for(Standard_Integer aCurInterval = 0; aCurInterval < aNbOfBoundaries; aCurInterval++) { + //Process every continuous region + if(Precision::IsInfinite(aU1f[aCurInterval]) && Precision::IsInfinite(aU1l[aCurInterval])) continue; @@ -2656,6 +2792,10 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1, 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]; @@ -2675,9 +2815,9 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1, anUexpect[i] = anUf; } - Standard_Real aCriticalDelta[aNbCritPointsMax] = {0}; + Standard_Real aCriticalDelta[aNbCritPointsMax]; for(Standard_Integer aCritPID = 0; aCritPID < aNbCritPoints; aCritPID++) - { //We are not intersted in elements of aCriticalDelta array + { //We are not interested in elements of aCriticalDelta array //if their index is greater than or equal to aNbCritPoints aCriticalDelta[aCritPID] = anUf - anU1crit[aCritPID]; @@ -2685,13 +2825,19 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1, Standard_Real anU1 = 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++) @@ -2875,6 +3021,7 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1, } } + // True if the current point already in the domain const Standard_Boolean isInscribe = ((aUSurf2f-aU2[i]) <= theTol2D) && ((aU2[i]-aUSurf2l) <= theTol2D) && ((aVSurf1f - aV1[i]) <= theTol2D) && ((aV1[i] - aVSurf1l) <= theTol2D) && @@ -3048,6 +3195,11 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1, 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(); { @@ -3099,6 +3251,12 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1, //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); @@ -3135,6 +3293,22 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1, 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 diff --git a/src/gp/gp_Trsf.cxx b/src/gp/gp_Trsf.cxx index 92f773ce5c..56908a19eb 100644 --- a/src/gp/gp_Trsf.cxx +++ b/src/gp/gp_Trsf.cxx @@ -774,9 +774,57 @@ Standard_Boolean gp_Trsf::GetRotation (gp_XYZ& theAxis, //======================================================================= //function : Orthogonalize //purpose : +//ATTENTION!!! +// Orthogonalization is not equivalent transformation. Therefore, +// transformation with source matrix and with orthogonalized matrix can +// lead to different results for one shape. Consequently, source matrix must +// be close to orthogonalized matrix for reducing these differences. //======================================================================= void gp_Trsf::Orthogonalize() { + //Matrix M is called orthogonal if and only if + // M*Transpose(M) == E + //where E is identity matrix. + + //Set of all rows (as of all columns) of matrix M (for gp_Trsf class) is + //orthonormal basis. If this condition is not satisfied then the basis can be + //orthonormalized in accordance with below described algorithm. + + //In 3D-space, we have the linear span of three basis vectors: V1, V2 and V3. + //Correspond orthonormalized basis is formed by vectors Vn1, Vn2 and Vn3. + + //In this case, + // Vn_{i}*Vn_{j} = (i == j)? 1 : 0. + + //The algorithm includes following steps: + + //1. Normalize V1 vector: + // V1n=V1/|V1|; + // + //2. Let + // V2n=V2-m*V1n. + // + //After multiplication two parts of this equation by V1n, + //we will have following equation: + // 0=V2*V1n-m <==> m=V2*V1n. + // + //Consequently, + // V2n=V2-(V2*V1n)*V1n. + + //3. Let + // V3n=V3-m1*V1n-m2*V2n. + // + //After multiplication two parts of this equation by V1n, + //we will have following equation: + // 0=V3*V1n-m1 <==> m1=V3*V1n. + // + //After multiplication two parts of main equation by V2n, + //we will have following equation: + // 0=V3*V2n-m2 <==> m2=V3*V2n. + // + //In conclusion, + // V3n=V3-(V3*V1n)*V1n-(V3*V2n)*V2n. + gp_Mat aTM(matrix); gp_XYZ aV1 = aTM.Column(1); diff --git a/src/gp/gp_Trsf2d.cxx b/src/gp/gp_Trsf2d.cxx index 7ae574abc7..e3b6a4224a 100644 --- a/src/gp/gp_Trsf2d.cxx +++ b/src/gp/gp_Trsf2d.cxx @@ -549,9 +549,17 @@ void gp_Trsf2d::SetValues(const Standard_Real a11, //======================================================================= //function : Orthogonalize //purpose : +//ATTENTION!!! +// Orthogonalization is not equivalent transformation.Therefore, transformation with +// source matrix and with orthogonalized matrix can lead to different results for +// one shape. Consequently, source matrix must be close to orthogonalized +// matrix for reducing these differences. //======================================================================= void gp_Trsf2d::Orthogonalize() { + //See correspond comment in gp_Trsf::Orthogonalize() method in order to make this + //algorithm clear. + gp_Mat2d aTM(matrix); gp_XY aV1 = aTM.Column(1);