projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0033661: Data Exchange, Step Import - Tessellated GDTs are not imported
[occt.git] / src / IntPatch / IntPatch_ImpImpIntersection_4.gxx
CommitLineData
b311480e 1// Created on: 1992-05-07
2// Created by: Jacques GOUSSARD
3// Copyright (c) 1992-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 8// This library is free software; you can redistribute it and/or modify it under
9// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 10// by the Free Software Foundation, with special exception defined in the file
11// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12// distribution for complete text of the license and disclaimer of any warranty.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
7365fad6 17#include <algorithm>
d30895f5 18#include <Bnd_Range.hxx>
d30895f5 19#include <IntAna_ListOfCurve.hxx>
b5ef9d91 20#include <math_Matrix.hxx>
261b7d9e 21#include <NCollection_IncAllocator.hxx>
22#include <Standard_DivideByZero.hxx>
1103eb60 23#include <math_Vector.hxx>
b5ef9d91 24
25//If Abs(a) <= aNulValue then it is considered that a = 0.
26static const Standard_Real aNulValue = 1.0e-11;
27
d30895f5 28static void ShortCosForm( const Standard_Real theCosFactor,
29 const Standard_Real theSinFactor,
30 Standard_Real& theCoeff,
31 Standard_Real& theAngle);
fa0291ff 32//
3306fdd9 33static Standard_Boolean ExploreCurve(const gp_Cone& theCo,
34 IntAna_Curve& aC,
35 const Standard_Real aTol,
36 IntAna_ListOfCurve& aLC);
7fd59977 37
b5ef9d91 38static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
39 const Standard_Real theUlTarget,
40 Standard_Real& theUGiven,
41 const Standard_Real theTol2D,
42 const Standard_Real thePeriod,
43 const Standard_Boolean theFlForce);
44
d30895f5 45
46class ComputationMethods
47{
4dba155d 48 //Every cylinder can be represented by the following equation in parametric form:
49 // S(U,V) = L + R*cos(U)*Xd+R*sin(U)*Yd+V*Zd,
50 //where location L, directions Xd, Yd and Zd have type gp_XYZ.
51
52 //Intersection points between two cylinders can be found from the following system:
53 // S1(U1, V1) = S2(U2, V2)
54 //or
55 // {X01 + R1*cos(U1)*Xx1 + R1*sin(U1)*Yx1 + V1*Zx1 = X02 + R2*cos(U2)*Xx2 + R2*sin(U2)*Yx2 + V2*Zx2
56 // {Y01 + R1*cos(U1)*Xy1 + R1*sin(U1)*Yy1 + V1*Zy1 = Y02 + R2*cos(U2)*Xy2 + R2*sin(U2)*Yy2 + V2*Zy2 (1)
57 // {Z01 + R1*cos(U1)*Xz1 + R1*sin(U1)*Yz1 + V1*Zz1 = Z02 + R2*cos(U2)*Xz2 + R2*sin(U2)*Yz2 + V2*Zz2
58
59 //The formula (1) can be rewritten as follows
60 // {C11*V1+C21*V2=A11*cos(U1)+B11*sin(U1)+A21*cos(U2)+B21*sin(U2)+D1
61 // {C12*V1+C22*V2=A12*cos(U1)+B12*sin(U1)+A22*cos(U2)+B22*sin(U2)+D2 (2)
62 // {C13*V1+C23*V2=A13*cos(U1)+B13*sin(U1)+A23*cos(U2)+B23*sin(U2)+D3
63
64 //Hereafter we consider that in system
65 // {C11*V1+C21*V2=A11*cos(U1)+B11*sin(U1)+A21*cos(U2)+B21*sin(U2)+D1 (3)
66 // {C12*V1+C22*V2=A12*cos(U1)+B12*sin(U1)+A22*cos(U2)+B22*sin(U2)+D2
67 //variables V1 and V2 can be found unambiguously, i.e. determinant
68 // |C11 C21|
69 // | | != 0
70 // |C12 C22|
71 //
72 //In this case, variables V1 and V2 can be found as follows:
73 // {V1 = K11*sin(U1)+K21*sin(U2)+L11*cos(U1)+L21*cos(U2)+M1 = K1*cos(U1-FIV1)+L1*cos(U2-PSIV1)+M1 (4)
74 // {V2 = K12*sin(U1)+K22*sin(U2)+L12*cos(U1)+L22*cos(U2)+M2 = K2*cos(U2-FIV2)+L2*cos(U2-PSIV2)+M2
75
76 //Having substituted result of (4) to the 3rd equation of (2), we will obtain equation
77 // cos(U2-FI2) = B*cos(U1-FI1)+C. (5)
78
79 //I.e. when U1 is taken different given values (from domain), correspond U2 value can be computed
80 //from equation (5). After that, V1 and V2 can be computed from the system (4) (see
81 //CylCylComputeParameters(...) methods).
82
83 //It is important to remark that equation (5) (in general) has two solutions: U2=FI2 +/- f(U1).
84 //Therefore, we are getting here two intersection lines.
85
d30895f5 86public:
87 //Stores equations coefficients
88 struct stCoeffsValue
89 {
90 stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&);
91
92 math_Vector mVecA1;
93 math_Vector mVecA2;
94 math_Vector mVecB1;
95 math_Vector mVecB2;
96 math_Vector mVecC1;
97 math_Vector mVecC2;
98 math_Vector mVecD;
99
100 Standard_Real mK21; //sinU2
101 Standard_Real mK11; //sinU1
102 Standard_Real mL21; //cosU2
103 Standard_Real mL11; //cosU1
104 Standard_Real mM1; //Free member
105
106 Standard_Real mK22; //sinU2
107 Standard_Real mK12; //sinU1
108 Standard_Real mL22; //cosU2
109 Standard_Real mL12; //cosU1
110 Standard_Real mM2; //Free member
111
112 Standard_Real mK1;
113 Standard_Real mL1;
114 Standard_Real mK2;
115 Standard_Real mL2;
116
117 Standard_Real mFIV1;
118 Standard_Real mPSIV1;
119 Standard_Real mFIV2;
120 Standard_Real mPSIV2;
121
122 Standard_Real mB;
123 Standard_Real mC;
124 Standard_Real mFI1;
125 Standard_Real mFI2;
126 };
127
128
129 //! Determines, if U2(U1) function is increasing.
130 static Standard_Boolean CylCylMonotonicity(const Standard_Real theU1par,
131 const Standard_Integer theWLIndex,
132 const stCoeffsValue& theCoeffs,
133 const Standard_Real thePeriod,
134 Standard_Boolean& theIsIncreasing);
135
136 //! Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0,
137 //! esimates the tolerance of U2-computing (estimation result is
138 //! assigned to *theDelta value).
139 static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
140 const Standard_Integer theWLIndex,
141 const stCoeffsValue& theCoeffs,
142 Standard_Real& theU2,
143 Standard_Real* const theDelta = 0);
144
145 static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1,
146 const Standard_Real theU2,
147 const stCoeffsValue& theCoeffs,
148 Standard_Real& theV1,
149 Standard_Real& theV2);
150
151 static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
152 const Standard_Integer theWLIndex,
153 const stCoeffsValue& theCoeffs,
154 Standard_Real& theU2,
155 Standard_Real& theV1,
156 Standard_Real& theV2);
157
158};
159
160ComputationMethods::stCoeffsValue::stCoeffsValue(const gp_Cylinder& theCyl1,
161 const gp_Cylinder& theCyl2):
162 mVecA1(-theCyl1.Radius()*theCyl1.XAxis().Direction().XYZ()),
163 mVecA2(theCyl2.Radius()*theCyl2.XAxis().Direction().XYZ()),
164 mVecB1(-theCyl1.Radius()*theCyl1.YAxis().Direction().XYZ()),
165 mVecB2(theCyl2.Radius()*theCyl2.YAxis().Direction().XYZ()),
166 mVecC1(theCyl1.Axis().Direction().XYZ()),
167 mVecC2(theCyl2.Axis().Direction().XYZ().Reversed()),
168 mVecD(theCyl2.Location().XYZ() - theCyl1.Location().XYZ())
169{
170 enum CoupleOfEquation
171 {
172 COENONE = 0,
173 COE12 = 1,
174 COE23 = 2,
175 COE13 = 3
176 }aFoundCouple = COENONE;
177
178
179 Standard_Real aDetV1V2 = 0.0;
180
181 const Standard_Real aDelta1 = mVecC1(1)*mVecC2(2)-mVecC1(2)*mVecC2(1); //1-2
182 const Standard_Real aDelta2 = mVecC1(2)*mVecC2(3)-mVecC1(3)*mVecC2(2); //2-3
183 const Standard_Real aDelta3 = mVecC1(1)*mVecC2(3)-mVecC1(3)*mVecC2(1); //1-3
184 const Standard_Real anAbsD1 = Abs(aDelta1); //1-2
185 const Standard_Real anAbsD2 = Abs(aDelta2); //2-3
186 const Standard_Real anAbsD3 = Abs(aDelta3); //1-3
187
188 if(anAbsD1 >= anAbsD2)
189 {
190 if(anAbsD3 > anAbsD1)
191 {
192 aFoundCouple = COE13;
193 aDetV1V2 = aDelta3;
194 }
195 else
196 {
197 aFoundCouple = COE12;
198 aDetV1V2 = aDelta1;
199 }
200 }
201 else
202 {
203 if(anAbsD3 > anAbsD2)
204 {
205 aFoundCouple = COE13;
206 aDetV1V2 = aDelta3;
207 }
208 else
209 {
210 aFoundCouple = COE23;
211 aDetV1V2 = aDelta2;
212 }
213 }
214
215 // In point of fact, every determinant (aDelta1, aDelta2 and aDelta3) is
216 // cross-product between directions (i.e. sine of angle).
217 // If sine is too small then sine is (approx.) equal to angle itself.
218 // Therefore, in this case we should compare sine with angular tolerance.
219 // This constant is used for check if axes are parallel (see constructor
220 // AxeOperator::AxeOperator(...) in IntAna_QuadQuadGeo.cxx file).
221 if(Abs(aDetV1V2) < Precision::Angular())
222 {
9775fa61 223 throw Standard_Failure("Error. Exception in divide by zerro (IntCyCyTrim)!!!!");
d30895f5 224 }
225
226 switch(aFoundCouple)
227 {
228 case COE12:
229 break;
230 case COE23:
231 {
232 math_Vector aVTemp(mVecA1);
233 mVecA1(1) = aVTemp(2);
234 mVecA1(2) = aVTemp(3);
235 mVecA1(3) = aVTemp(1);
236
237 aVTemp = mVecA2;
238 mVecA2(1) = aVTemp(2);
239 mVecA2(2) = aVTemp(3);
240 mVecA2(3) = aVTemp(1);
241
242 aVTemp = mVecB1;
243 mVecB1(1) = aVTemp(2);
244 mVecB1(2) = aVTemp(3);
245 mVecB1(3) = aVTemp(1);
246
247 aVTemp = mVecB2;
248 mVecB2(1) = aVTemp(2);
249 mVecB2(2) = aVTemp(3);
250 mVecB2(3) = aVTemp(1);
251
252 aVTemp = mVecC1;
253 mVecC1(1) = aVTemp(2);
254 mVecC1(2) = aVTemp(3);
255 mVecC1(3) = aVTemp(1);
256
257 aVTemp = mVecC2;
258 mVecC2(1) = aVTemp(2);
259 mVecC2(2) = aVTemp(3);
260 mVecC2(3) = aVTemp(1);
261
262 aVTemp = mVecD;
263 mVecD(1) = aVTemp(2);
264 mVecD(2) = aVTemp(3);
265 mVecD(3) = aVTemp(1);
266
267 break;
268 }
269 case COE13:
270 {
271 math_Vector aVTemp = mVecA1;
272 mVecA1(2) = aVTemp(3);
273 mVecA1(3) = aVTemp(2);
274
275 aVTemp = mVecA2;
276 mVecA2(2) = aVTemp(3);
277 mVecA2(3) = aVTemp(2);
278
279 aVTemp = mVecB1;
280 mVecB1(2) = aVTemp(3);
281 mVecB1(3) = aVTemp(2);
282
283 aVTemp = mVecB2;
284 mVecB2(2) = aVTemp(3);
285 mVecB2(3) = aVTemp(2);
286
287 aVTemp = mVecC1;
288 mVecC1(2) = aVTemp(3);
289 mVecC1(3) = aVTemp(2);
290
291 aVTemp = mVecC2;
292 mVecC2(2) = aVTemp(3);
293 mVecC2(3) = aVTemp(2);
294
295 aVTemp = mVecD;
296 mVecD(2) = aVTemp(3);
297 mVecD(3) = aVTemp(2);
298
299 break;
300 }
301 default:
302 break;
303 }
304
305 //------- For V1 (begin)
306 //sinU2
307 mK21 = (mVecC2(2)*mVecB2(1)-mVecC2(1)*mVecB2(2))/aDetV1V2;
308 //sinU1
309 mK11 = (mVecC2(2)*mVecB1(1)-mVecC2(1)*mVecB1(2))/aDetV1V2;
310 //cosU2
311 mL21 = (mVecC2(2)*mVecA2(1)-mVecC2(1)*mVecA2(2))/aDetV1V2;
312 //cosU1
313 mL11 = (mVecC2(2)*mVecA1(1)-mVecC2(1)*mVecA1(2))/aDetV1V2;
314 //Free member
315 mM1 = (mVecC2(2)*mVecD(1)-mVecC2(1)*mVecD(2))/aDetV1V2;
316 //------- For V1 (end)
317
318 //------- For V2 (begin)
319 //sinU2
320 mK22 = (mVecC1(1)*mVecB2(2)-mVecC1(2)*mVecB2(1))/aDetV1V2;
321 //sinU1
322 mK12 = (mVecC1(1)*mVecB1(2)-mVecC1(2)*mVecB1(1))/aDetV1V2;
323 //cosU2
324 mL22 = (mVecC1(1)*mVecA2(2)-mVecC1(2)*mVecA2(1))/aDetV1V2;
325 //cosU1
326 mL12 = (mVecC1(1)*mVecA1(2)-mVecC1(2)*mVecA1(1))/aDetV1V2;
327 //Free member
328 mM2 = (mVecC1(1)*mVecD(2)-mVecC1(2)*mVecD(1))/aDetV1V2;
329 //------- For V1 (end)
330
331 ShortCosForm(mL11, mK11, mK1, mFIV1);
332 ShortCosForm(mL21, mK21, mL1, mPSIV1);
333 ShortCosForm(mL12, mK12, mK2, mFIV2);
334 ShortCosForm(mL22, mK22, mL2, mPSIV2);
335
336 const Standard_Real aA1=mVecC1(3)*mK21+mVecC2(3)*mK22-mVecB2(3), //sinU2
337 aA2=mVecC1(3)*mL21+mVecC2(3)*mL22-mVecA2(3), //cosU2
338 aB1=mVecB1(3)-mVecC1(3)*mK11-mVecC2(3)*mK12, //sinU1
339 aB2=mVecA1(3)-mVecC1(3)*mL11-mVecC2(3)*mL12; //cosU1
340
341 mC =mVecD(3) - mVecC1(3)*mM1 -mVecC2(3)*mM2; //Free
342
343 Standard_Real aA = 0.0;
344
345 ShortCosForm(aB2,aB1,mB,mFI1);
346 ShortCosForm(aA2,aA1,aA,mFI2);
347
348 mB /= aA;
349 mC /= aA;
350}
351
352class WorkWithBoundaries
353{
354public:
355 enum SearchBoundType
356 {
357 SearchNONE = 0,
358 SearchV1 = 1,
359 SearchV2 = 2
360 };
361
362 struct StPInfo
363 {
364 StPInfo()
365 {
366 mySurfID = 0;
367 myU1 = RealLast();
368 myV1 = RealLast();
369 myU2 = RealLast();
370 myV2 = RealLast();
371 }
372
261b7d9e 373 //Equal to 0 for 1st surface non-zero for 2nd one.
d30895f5 374 Standard_Integer mySurfID;
375
376 Standard_Real myU1;
377 Standard_Real myV1;
378 Standard_Real myU2;
379 Standard_Real myV2;
380
381 bool operator>(const StPInfo& theOther) const
382 {
383 return myU1 > theOther.myU1;
384 }
385
386 bool operator<(const StPInfo& theOther) const
387 {
388 return myU1 < theOther.myU1;
389 }
390
391 bool operator==(const StPInfo& theOther) const
392 {
393 return myU1 == theOther.myU1;
394 }
395 };
396
397 WorkWithBoundaries(const IntSurf_Quadric& theQuad1,
398 const IntSurf_Quadric& theQuad2,
399 const ComputationMethods::stCoeffsValue& theCoeffs,
400 const Bnd_Box2d& theUVSurf1,
401 const Bnd_Box2d& theUVSurf2,
402 const Standard_Integer theNbWLines,
403 const Standard_Real thePeriod,
404 const Standard_Real theTol3D,
405 const Standard_Real theTol2D,
406 const Standard_Boolean isTheReverse) :
407 myQuad1(theQuad1), myQuad2(theQuad2), myCoeffs(theCoeffs),
408 myUVSurf1(theUVSurf1), myUVSurf2(theUVSurf2), myNbWLines(theNbWLines),
409 myPeriod(thePeriod), myTol3D(theTol3D), myTol2D(theTol2D),
410 myIsReverse(isTheReverse)
411 {
412 };
413
261b7d9e 414 // Returns parameters of system solved while finding
415 // intersection line
416 const ComputationMethods::stCoeffsValue &SICoeffs() const
417 {
418 return myCoeffs;
419 }
420
421 // Returns quadric correspond to the index theIdx.
422 const IntSurf_Quadric& GetQSurface(const Standard_Integer theIdx) const
423 {
424 if (theIdx <= 1)
425 return myQuad1;
426
427 return myQuad2;
428 }
429
430 // Returns TRUE in case of reverting surfaces
431 Standard_Boolean IsReversed() const
432 {
433 return myIsReverse;
434 }
435
436 // Returns 2D-tolerance
437 Standard_Real Get2dTolerance() const
438 {
439 return myTol2D;
440 }
441
442 // Returns 3D-tolerance
443 Standard_Real Get3dTolerance() const
444 {
445 return myTol3D;
446 }
447
448 // Returns UV-bounds of 1st surface
449 const Bnd_Box2d& UVS1() const
450 {
451 return myUVSurf1;
452 }
453
454 // Returns UV-bounds of 2nd surface
455 const Bnd_Box2d& UVS2() const
456 {
457 return myUVSurf2;
458 }
459
d30895f5 460 void AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL,
461 const Standard_Real theU1,
261b7d9e 462 const Standard_Real theU1Min,
d30895f5 463 const Standard_Real theU2,
464 const Standard_Real theV1,
465 const Standard_Real theV1Prev,
466 const Standard_Real theV2,
467 const Standard_Real theV2Prev,
468 const Standard_Integer theWLIndex,
469 const Standard_Boolean theFlForce,
470 Standard_Boolean& isTheFound1,
471 Standard_Boolean& isTheFound2) const;
472
261b7d9e 473 static Standard_Boolean BoundariesComputing(const ComputationMethods::stCoeffsValue &theCoeffs,
474 const Standard_Real thePeriod,
475 Bnd_Range theURange[]);
d30895f5 476
477 void BoundaryEstimation(const gp_Cylinder& theCy1,
478 const gp_Cylinder& theCy2,
479 Bnd_Range& theOutBoxS1,
480 Bnd_Range& theOutBoxS2) const;
481
482protected:
483
4dba155d 484 //Solves equation (2) (see declaration of ComputationMethods class) in case,
485 //when V1 or V2 (is set by theSBType argument) is known (corresponds to the boundary
486 //and equal to theVzad) but U1 is unknown. Computation is made by numeric methods and
487 //requires initial values (theVInit, theInitU2 and theInitMainVar).
d30895f5 488 Standard_Boolean
489 SearchOnVBounds(const SearchBoundType theSBType,
490 const Standard_Real theVzad,
491 const Standard_Real theVInit,
492 const Standard_Real theInitU2,
493 const Standard_Real theInitMainVar,
494 Standard_Real& theMainVariableValue) const;
495
496 const WorkWithBoundaries& operator=(const WorkWithBoundaries&);
497
498private:
499 friend class ComputationMethods;
500
501 const IntSurf_Quadric& myQuad1;
502 const IntSurf_Quadric& myQuad2;
503 const ComputationMethods::stCoeffsValue& myCoeffs;
504 const Bnd_Box2d& myUVSurf1;
505 const Bnd_Box2d& myUVSurf2;
506 const Standard_Integer myNbWLines;
507 const Standard_Real myPeriod;
508 const Standard_Real myTol3D;
509 const Standard_Real myTol2D;
510 const Standard_Boolean myIsReverse;
511};
512
b5ef9d91 513static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
514 const IntSurf_Quadric& theQuad2,
515 const Handle(IntSurf_LineOn2S)& theLine,
d30895f5 516 const ComputationMethods::stCoeffsValue& theCoeffs,
b5ef9d91 517 const Standard_Integer theWLIndex,
518 const Standard_Integer theMinNbPoints,
519 const Standard_Integer theStartPointOnLine,
520 const Standard_Integer theEndPointOnLine,
b5ef9d91 521 const Standard_Real theTol2D,
522 const Standard_Real thePeriodOfSurf2,
523 const Standard_Boolean isTheReverse);
524
525//=======================================================================
526//function : MinMax
527//purpose : Replaces theParMIN = MIN(theParMIN, theParMAX),
9e20ed57 528// theParMAX = MAX(theParMIN, theParMAX).
b5ef9d91 529//=======================================================================
9e20ed57 530static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX)
531{
532 if(theParMIN > theParMAX)
533 {
534 const Standard_Real aux = theParMAX;
535 theParMAX = theParMIN;
536 theParMIN = aux;
537 }
538}
539
d30895f5 540//=======================================================================
541//function : ExtremaLineLine
542//purpose : Computes extrema between the given lines. Returns parameters
4dba155d 543// on correspond curve (see correspond method for Extrema_ExtElC class).
d30895f5 544//=======================================================================
545static inline void ExtremaLineLine(const gp_Ax1& theC1,
546 const gp_Ax1& theC2,
547 const Standard_Real theCosA,
548 const Standard_Real theSqSinA,
549 Standard_Real& thePar1,
550 Standard_Real& thePar2)
551{
552 const gp_Dir &aD1 = theC1.Direction(),
553 &aD2 = theC2.Direction();
554
555 const gp_XYZ aL1L2 = theC2.Location().XYZ() - theC1.Location().XYZ();
556 const Standard_Real aD1L = aD1.XYZ().Dot(aL1L2),
557 aD2L = aD2.XYZ().Dot(aL1L2);
558
559 thePar1 = (aD1L - theCosA * aD2L) / theSqSinA;
560 thePar2 = (theCosA * aD1L - aD2L) / theSqSinA;
561}
562
b5ef9d91 563//=======================================================================
564//function : VBoundaryPrecise
4dba155d 565//purpose : By default, we shall consider, that V1 and V2 will be increased
566// if U1 is increased. But if it is not, new V1set and/or V2set
b5ef9d91 567// must be computed as [V1current - DeltaV1] (analogically
568// for V2). This function processes this case.
569//=======================================================================
570static void VBoundaryPrecise( const math_Matrix& theMatr,
571 const Standard_Real theV1AfterDecrByDelta,
572 const Standard_Real theV2AfterDecrByDelta,
b5ef9d91 573 Standard_Real& theV1Set,
574 Standard_Real& theV2Set)
575{
576 //Now we are going to define if V1 (and V2) increases
577 //(or decreases) when U1 will increase.
578 const Standard_Integer aNbDim = 3;
579 math_Matrix aSyst(1, aNbDim, 1, aNbDim);
580
581 aSyst.SetCol(1, theMatr.Col(1));
582 aSyst.SetCol(2, theMatr.Col(2));
583 aSyst.SetCol(3, theMatr.Col(4));
584
4dba155d 585 //We have the system (see comment to StepComputing(...) function)
586 // {a11*dV1 + a12*dV2 + a14*dU2 = -a13*dU1
587 // {a21*dV1 + a22*dV2 + a24*dU2 = -a23*dU1
588 // {a31*dV1 + a32*dV2 + a34*dU2 = -a33*dU1
589
b5ef9d91 590 const Standard_Real aDet = aSyst.Determinant();
591
592 aSyst.SetCol(1, theMatr.Col(3));
593 const Standard_Real aDet1 = aSyst.Determinant();
594
595 aSyst.SetCol(1, theMatr.Col(1));
596 aSyst.SetCol(2, theMatr.Col(3));
597
598 const Standard_Real aDet2 = aSyst.Determinant();
599
4dba155d 600 //Now,
601 // dV1 = -dU1*aDet1/aDet
602 // dV2 = -dU1*aDet2/aDet
603
604 //If U1 is increased then dU1 > 0.
605 //If (aDet1/aDet > 0) then dV1 < 0 and
606 //V1 will be decreased after increasing U1.
607
608 //We have analogical situation with V2-parameter.
609
b5ef9d91 610 if(aDet*aDet1 > 0.0)
611 {
7365fad6 612 theV1Set = theV1AfterDecrByDelta;
b5ef9d91 613 }
614
615 if(aDet*aDet2 > 0.0)
616 {
7365fad6 617 theV2Set = theV2AfterDecrByDelta;
b5ef9d91 618 }
619}
620
621//=======================================================================
622//function : DeltaU1Computing
623//purpose : Computes new step for U1 parameter.
624//=======================================================================
625static inline
626 Standard_Boolean DeltaU1Computing(const math_Matrix& theSyst,
627 const math_Vector& theFree,
628 Standard_Real& theDeltaU1Found)
629{
630 Standard_Real aDet = theSyst.Determinant();
631
632 if(Abs(aDet) > aNulValue)
633 {
634 math_Matrix aSyst1(theSyst);
635 aSyst1.SetCol(2, theFree);
636
637 theDeltaU1Found = Abs(aSyst1.Determinant()/aDet);
638 return Standard_True;
639 }
640
641 return Standard_False;
642}
643
644//=======================================================================
645//function : StepComputing
646//purpose :
647//
648//Attention!!!:
649// theMatr must have 3*5-dimension strictly.
650// For system
651// {a11*V1+a12*V2+a13*dU1+a14*dU2=b1;
652// {a21*V1+a22*V2+a23*dU1+a24*dU2=b2;
653// {a31*V1+a32*V2+a33*dU1+a34*dU2=b3;
654// theMatr must be following:
655// (a11 a12 a13 a14 b1)
656// (a21 a22 a23 a24 b2)
657// (a31 a32 a33 a34 b3)
658//=======================================================================
659static Standard_Boolean StepComputing(const math_Matrix& theMatr,
660 const Standard_Real theV1Cur,
661 const Standard_Real theV2Cur,
662 const Standard_Real theDeltaV1,
663 const Standard_Real theDeltaV2,
b5ef9d91 664 Standard_Real& theDeltaU1Found/*,
665 Standard_Real& theDeltaU2Found,
666 Standard_Real& theV1Found,
667 Standard_Real& theV2Found*/)
668{
d30895f5 669#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
7365fad6 670 bool flShow = false;
671
672 if(flShow)
673 {
674 printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
675 theMatr(1,1), theMatr(1,2), theMatr(1,3), theMatr(1,4), theMatr(1,5));
676 printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
677 theMatr(2,1), theMatr(2,2), theMatr(2,3), theMatr(2,4), theMatr(2,5));
678 printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
679 theMatr(3,1), theMatr(3,2), theMatr(3,3), theMatr(3,4), theMatr(3,5));
680 }
681#endif
682
b5ef9d91 683 Standard_Boolean isSuccess = Standard_False;
684 theDeltaU1Found/* = theDeltaU2Found*/ = RealLast();
685 //theV1Found = theV1set;
686 //theV2Found = theV2Set;
687 const Standard_Integer aNbDim = 3;
688
689 math_Matrix aSyst(1, aNbDim, 1, aNbDim);
690 math_Vector aFree(1, aNbDim);
691
692 //By default, increasing V1(U1) and V2(U1) functions is
693 //considered
7365fad6 694 Standard_Real aV1Set = theV1Cur + theDeltaV1,
695 aV2Set = theV2Cur + theDeltaV2;
b5ef9d91 696
7365fad6 697 //However, what is indeed?
698 VBoundaryPrecise( theMatr, theV1Cur - theDeltaV1,
699 theV2Cur - theDeltaV2, aV1Set, aV2Set);
b5ef9d91 700
701 aSyst.SetCol(2, theMatr.Col(3));
702 aSyst.SetCol(3, theMatr.Col(4));
703
704 for(Standard_Integer i = 0; i < 2; i++)
705 {
706 if(i == 0)
707 {//V1 is known
708 aSyst.SetCol(1, theMatr.Col(2));
709 aFree.Set(1, aNbDim, theMatr.Col(5)-aV1Set*theMatr.Col(1));
710 }
711 else
712 {//i==1 => V2 is known
713 aSyst.SetCol(1, theMatr.Col(1));
714 aFree.Set(1, aNbDim, theMatr.Col(5)-aV2Set*theMatr.Col(2));
715 }
716
717 Standard_Real aNewDU = theDeltaU1Found;
718 if(DeltaU1Computing(aSyst, aFree, aNewDU))
719 {
720 isSuccess = Standard_True;
721 if(aNewDU < theDeltaU1Found)
722 {
723 theDeltaU1Found = aNewDU;
724 }
725 }
726 }
727
728 if(!isSuccess)
729 {
730 aFree = theMatr.Col(5) - aV1Set*theMatr.Col(1) - aV2Set*theMatr.Col(2);
731 math_Matrix aSyst1(1, aNbDim, 1, 2);
732 aSyst1.SetCol(1, aSyst.Col(2));
733 aSyst1.SetCol(2, aSyst.Col(3));
734
735 //Now we have overdetermined system.
736
737 const Standard_Real aDet1 = theMatr(1,3)*theMatr(2,4) - theMatr(2,3)*theMatr(1,4);
738 const Standard_Real aDet2 = theMatr(1,3)*theMatr(3,4) - theMatr(3,3)*theMatr(1,4);
739 const Standard_Real aDet3 = theMatr(2,3)*theMatr(3,4) - theMatr(3,3)*theMatr(2,4);
740 const Standard_Real anAbsD1 = Abs(aDet1);
741 const Standard_Real anAbsD2 = Abs(aDet2);
742 const Standard_Real anAbsD3 = Abs(aDet3);
743
744 if(anAbsD1 >= anAbsD2)
745 {
746 if(anAbsD1 >= anAbsD3)
747 {
748 //Det1
749 if(anAbsD1 <= aNulValue)
750 return isSuccess;
751
752 theDeltaU1Found = Abs(aFree(1)*theMatr(2,4) - aFree(2)*theMatr(1,4))/anAbsD1;
753 isSuccess = Standard_True;
754 }
755 else
756 {
757 //Det3
758 if(anAbsD3 <= aNulValue)
759 return isSuccess;
760
761 theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
762 isSuccess = Standard_True;
763 }
764 }
765 else
766 {
767 if(anAbsD2 >= anAbsD3)
768 {
769 //Det2
770 if(anAbsD2 <= aNulValue)
771 return isSuccess;
772
773 theDeltaU1Found = Abs(aFree(1)*theMatr(3,4) - aFree(3)*theMatr(1,4))/anAbsD2;
774 isSuccess = Standard_True;
775 }
776 else
777 {
778 //Det3
779 if(anAbsD3 <= aNulValue)
780 return isSuccess;
781
782 theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
783 isSuccess = Standard_True;
784 }
785 }
786 }
787
788 return isSuccess;
789}
9e20ed57 790
7fd59977 791//=======================================================================
792//function : ProcessBounds
793//purpose :
794//=======================================================================
795void ProcessBounds(const Handle(IntPatch_ALine)& alig, //-- ligne courante
796 const IntPatch_SequenceOfLine& slin,
797 const IntSurf_Quadric& Quad1,
798 const IntSurf_Quadric& Quad2,
799 Standard_Boolean& procf,
800 const gp_Pnt& ptf, //-- Debut Ligne Courante
801 const Standard_Real first, //-- Paramf
802 Standard_Boolean& procl,
803 const gp_Pnt& ptl, //-- Fin Ligne courante
804 const Standard_Real last, //-- Paraml
805 Standard_Boolean& Multpoint,
806 const Standard_Real Tol)
807{
808 Standard_Integer j,k;
809 Standard_Real U1,V1,U2,V2;
810 IntPatch_Point ptsol;
811 Standard_Real d;
812
813 if (procf && procl) {
814 j = slin.Length() + 1;
815 }
816 else {
817 j = 1;
818 }
819
820
821 //-- On prend les lignes deja enregistrees
822
823 while (j <= slin.Length()) {
824 if(slin.Value(j)->ArcType() == IntPatch_Analytic) {
825 const Handle(IntPatch_ALine)& aligold = *((Handle(IntPatch_ALine)*)&slin.Value(j));
826 k = 1;
827
828 //-- On prend les vertex des lignes deja enregistrees
829
830 while (k <= aligold->NbVertex()) {
831 ptsol = aligold->Vertex(k);
832 if (!procf) {
833 d=ptf.Distance(ptsol.Value());
834 if (d <= Tol) {
3306fdd9 835 ptsol.SetTolerance(Tol);
7fd59977 836 if (!ptsol.IsMultiple()) {
837 //-- le point ptsol (de aligold) est declare multiple sur aligold
838 Multpoint = Standard_True;
839 ptsol.SetMultiple(Standard_True);
840 aligold->Replace(k,ptsol);
841 }
842 ptsol.SetParameter(first);
843 alig->AddVertex(ptsol);
844 alig->SetFirstPoint(alig->NbVertex());
845 procf = Standard_True;
846
847 //-- On restore le point avec son parametre sur aligold
848 ptsol = aligold->Vertex(k);
849
850 }
851 }
852 if (!procl) {
853 if (ptl.Distance(ptsol.Value()) <= Tol) {
3306fdd9 854 ptsol.SetTolerance(Tol);
7fd59977 855 if (!ptsol.IsMultiple()) {
856 Multpoint = Standard_True;
857 ptsol.SetMultiple(Standard_True);
858 aligold->Replace(k,ptsol);
859 }
860 ptsol.SetParameter(last);
861 alig->AddVertex(ptsol);
862 alig->SetLastPoint(alig->NbVertex());
863 procl = Standard_True;
864
865 //-- On restore le point avec son parametre sur aligold
866 ptsol = aligold->Vertex(k);
867
868 }
869 }
870 if (procf && procl) {
871 k = aligold->NbVertex()+1;
872 }
873 else {
874 k = k+1;
875 }
876 }
877 if (procf && procl) {
878 j = slin.Length()+1;
879 }
880 else {
881 j = j+1;
882 }
883 }
884 }
3306fdd9 885
886 ptsol.SetTolerance(Tol);
7fd59977 887 if (!procf && !procl) {
888 Quad1.Parameters(ptf,U1,V1);
889 Quad2.Parameters(ptf,U2,V2);
890 ptsol.SetValue(ptf,Tol,Standard_False);
891 ptsol.SetParameters(U1,V1,U2,V2);
892 ptsol.SetParameter(first);
893 if (ptf.Distance(ptl) <= Tol) {
894 ptsol.SetMultiple(Standard_True); // a voir
895 Multpoint = Standard_True; // a voir de meme
896 alig->AddVertex(ptsol);
897 alig->SetFirstPoint(alig->NbVertex());
898
899 ptsol.SetParameter(last);
900 alig->AddVertex(ptsol);
901 alig->SetLastPoint(alig->NbVertex());
902 }
903 else {
904 alig->AddVertex(ptsol);
905 alig->SetFirstPoint(alig->NbVertex());
906 Quad1.Parameters(ptl,U1,V1);
907 Quad2.Parameters(ptl,U2,V2);
908 ptsol.SetValue(ptl,Tol,Standard_False);
909 ptsol.SetParameters(U1,V1,U2,V2);
910 ptsol.SetParameter(last);
911 alig->AddVertex(ptsol);
912 alig->SetLastPoint(alig->NbVertex());
913 }
914 }
915 else if (!procf) {
916 Quad1.Parameters(ptf,U1,V1);
917 Quad2.Parameters(ptf,U2,V2);
918 ptsol.SetValue(ptf,Tol,Standard_False);
919 ptsol.SetParameters(U1,V1,U2,V2);
920 ptsol.SetParameter(first);
921 alig->AddVertex(ptsol);
922 alig->SetFirstPoint(alig->NbVertex());
923 }
924 else if (!procl) {
925 Quad1.Parameters(ptl,U1,V1);
926 Quad2.Parameters(ptl,U2,V2);
927 ptsol.SetValue(ptl,Tol,Standard_False);
928 ptsol.SetParameters(U1,V1,U2,V2);
929 ptsol.SetParameter(last);
930 alig->AddVertex(ptsol);
931 alig->SetLastPoint(alig->NbVertex());
932 }
933}
d30895f5 934
7fd59977 935//=======================================================================
d30895f5 936//function : CyCyAnalyticalIntersect
4dba155d 937//purpose : Checks if intersection curve is analytical (line, circle, ellipse)
938// and returns these curves.
7fd59977 939//=======================================================================
d30895f5 940Standard_Boolean CyCyAnalyticalIntersect( const IntSurf_Quadric& Quad1,
941 const IntSurf_Quadric& Quad2,
942 const IntAna_QuadQuadGeo& theInter,
943 const Standard_Real Tol,
d30895f5 944 Standard_Boolean& Empty,
945 Standard_Boolean& Same,
946 Standard_Boolean& Multpoint,
947 IntPatch_SequenceOfLine& slin,
948 IntPatch_SequenceOfPoint& spnt)
7fd59977 949{
950 IntPatch_Point ptsol;
d30895f5 951
7fd59977 952 Standard_Integer i;
953
954 IntSurf_TypeTrans trans1,trans2;
955 IntAna_ResultType typint;
956
957 gp_Elips elipsol;
958 gp_Lin linsol;
959
960 gp_Cylinder Cy1(Quad1.Cylinder());
961 gp_Cylinder Cy2(Quad2.Cylinder());
962
d30895f5 963 typint = theInter.TypeInter();
964 Standard_Integer NbSol = theInter.NbSolutions();
7fd59977 965 Empty = Standard_False;
966 Same = Standard_False;
967
ecc4f148 968 switch (typint)
d30895f5 969 {
ecc4f148 970 case IntAna_Empty:
7fd59977 971 {
972 Empty = Standard_True;
d30895f5 973 }
7fd59977 974 break;
975
976 case IntAna_Same:
d30895f5 977 {
7fd59977 978 Same = Standard_True;
d30895f5 979 }
7fd59977 980 break;
d30895f5 981
ecc4f148 982 case IntAna_Point:
7fd59977 983 {
d30895f5 984 gp_Pnt psol(theInter.Point(1));
7fd59977 985 ptsol.SetValue(psol,Tol,Standard_True);
d30895f5 986
987 Standard_Real U1,V1,U2,V2;
261b7d9e 988 Quad1.Parameters(psol, U1, V1);
989 Quad2.Parameters(psol, U2, V2);
d30895f5 990
7fd59977 991 ptsol.SetParameters(U1,V1,U2,V2);
992 spnt.Append(ptsol);
993 }
994 break;
995
996 case IntAna_Line:
997 {
998 gp_Pnt ptref;
ecc4f148 999 if (NbSol == 1)
1000 { // Cylinders are tangent to each other by line
d30895f5 1001 linsol = theInter.Line(1);
ecc4f148 1002 ptref = linsol.Location();
d30895f5 1003
1004 //Radius-vectors
ecc4f148 1005 gp_Dir crb1(gp_Vec(ptref,Cy1.Location()));
1006 gp_Dir crb2(gp_Vec(ptref,Cy2.Location()));
d30895f5 1007
1008 //outer normal lines
ecc4f148 1009 gp_Vec norm1(Quad1.Normale(ptref));
1010 gp_Vec norm2(Quad2.Normale(ptref));
1011 IntSurf_Situation situcyl1;
1012 IntSurf_Situation situcyl2;
1013
1014 if (crb1.Dot(crb2) < 0.)
1015 { // centre de courbures "opposes"
d30895f5 1016 //ATTENTION!!!
1017 // Normal and Radius-vector of the 1st(!) cylinder
1018 // is used for judging what the situation of the 2nd(!)
1019 // cylinder is.
1020
ecc4f148 1021 if (norm1.Dot(crb1) > 0.)
1022 {
1023 situcyl2 = IntSurf_Inside;
1024 }
1025 else
1026 {
1027 situcyl2 = IntSurf_Outside;
1028 }
1029
1030 if (norm2.Dot(crb2) > 0.)
1031 {
1032 situcyl1 = IntSurf_Inside;
1033 }
1034 else
1035 {
1036 situcyl1 = IntSurf_Outside;
1037 }
1038 }
1039 else
d30895f5 1040 {
ecc4f148 1041 if (Cy1.Radius() < Cy2.Radius())
1042 {
1043 if (norm1.Dot(crb1) > 0.)
1044 {
1045 situcyl2 = IntSurf_Inside;
1046 }
1047 else
1048 {
1049 situcyl2 = IntSurf_Outside;
1050 }
1051
1052 if (norm2.Dot(crb2) > 0.)
1053 {
1054 situcyl1 = IntSurf_Outside;
1055 }
1056 else
1057 {
1058 situcyl1 = IntSurf_Inside;
1059 }
1060 }
1061 else
1062 {
1063 if (norm1.Dot(crb1) > 0.)
1064 {
1065 situcyl2 = IntSurf_Outside;
1066 }
1067 else
1068 {
1069 situcyl2 = IntSurf_Inside;
1070 }
1071
1072 if (norm2.Dot(crb2) > 0.)
1073 {
1074 situcyl1 = IntSurf_Inside;
1075 }
1076 else
1077 {
1078 situcyl1 = IntSurf_Outside;
1079 }
1080 }
1081 }
1082
1083 Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_True, situcyl1, situcyl2);
1084 slin.Append(glig);
7fd59977 1085 }
ecc4f148 1086 else
1087 {
1088 for (i=1; i <= NbSol; i++)
1089 {
d30895f5 1090 linsol = theInter.Line(i);
ecc4f148 1091 ptref = linsol.Location();
1092 gp_Vec lsd = linsol.Direction();
d30895f5 1093
1094 //Theoretically, qwe = +/- 1.0.
ecc4f148 1095 Standard_Real qwe = lsd.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
1096 if (qwe >0.00000001)
1097 {
1098 trans1 = IntSurf_Out;
1099 trans2 = IntSurf_In;
1100 }
1101 else if (qwe <-0.00000001)
1102 {
1103 trans1 = IntSurf_In;
1104 trans2 = IntSurf_Out;
1105 }
1106 else
1107 {
1108 trans1=trans2=IntSurf_Undecided;
1109 }
1110
1111 Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
1112 slin.Append(glig);
1113 }
7fd59977 1114 }
1115 }
1116 break;
1117
1118 case IntAna_Ellipse:
1119 {
1120 gp_Vec Tgt;
1121 gp_Pnt ptref;
fa0291ff 1122 IntPatch_Point pmult1, pmult2;
d30895f5 1123
1124 elipsol = theInter.Ellipse(1);
1125
fa0291ff 1126 gp_Pnt pttang1(ElCLib::Value(0.5*M_PI, elipsol));
1127 gp_Pnt pttang2(ElCLib::Value(1.5*M_PI, elipsol));
d30895f5 1128
7fd59977 1129 Multpoint = Standard_True;
1130 pmult1.SetValue(pttang1,Tol,Standard_True);
1131 pmult2.SetValue(pttang2,Tol,Standard_True);
1132 pmult1.SetMultiple(Standard_True);
1133 pmult2.SetMultiple(Standard_True);
d30895f5 1134
7fd59977 1135 Standard_Real oU1,oV1,oU2,oV2;
261b7d9e 1136 Quad1.Parameters(pttang1, oU1, oV1);
1137 Quad2.Parameters(pttang1, oU2, oV2);
d30895f5 1138
7fd59977 1139 pmult1.SetParameters(oU1,oV1,oU2,oV2);
261b7d9e 1140 Quad1.Parameters(pttang2,oU1,oV1);
1141 Quad2.Parameters(pttang2,oU2,oV2);
d30895f5 1142
7fd59977 1143 pmult2.SetParameters(oU1,oV1,oU2,oV2);
7fd59977 1144
d30895f5 1145 // on traite la premiere ellipse
1146
7fd59977 1147 //-- Calcul de la Transition de la ligne
1148 ElCLib::D1(0.,elipsol,ptref,Tgt);
d30895f5 1149
1150 //Theoretically, qwe = +/- |Tgt|.
7fd59977 1151 Standard_Real qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
ecc4f148 1152 if (qwe>0.00000001)
1153 {
1154 trans1 = IntSurf_Out;
1155 trans2 = IntSurf_In;
7fd59977 1156 }
ecc4f148 1157 else if (qwe<-0.00000001)
1158 {
1159 trans1 = IntSurf_In;
1160 trans2 = IntSurf_Out;
7fd59977 1161 }
ecc4f148 1162 else
d30895f5 1163 {
ecc4f148 1164 trans1=trans2=IntSurf_Undecided;
7fd59977 1165 }
ecc4f148 1166
7fd59977 1167 //-- Transition calculee au point 0 -> Trans2 , Trans1
1168 //-- car ici, on devarit calculer en PI
1169 Handle(IntPatch_GLine) glig = new IntPatch_GLine(elipsol,Standard_False,trans2,trans1);
fa0291ff 1170 //
1171 {
ecc4f148 1172 Standard_Real aU1, aV1, aU2, aV2;
1173 IntPatch_Point aIP;
1174 gp_Pnt aP (ElCLib::Value(0., elipsol));
1175 //
1176 aIP.SetValue(aP,Tol,Standard_False);
1177 aIP.SetMultiple(Standard_False);
1178 //
261b7d9e 1179 Quad1.Parameters(aP, aU1, aV1);
1180 Quad2.Parameters(aP, aU2, aV2);
d30895f5 1181
ecc4f148 1182 aIP.SetParameters(aU1, aV1, aU2, aV2);
1183 //
1184 aIP.SetParameter(0.);
1185 glig->AddVertex(aIP);
1186 glig->SetFirstPoint(1);
1187 //
1188 aIP.SetParameter(2.*M_PI);
1189 glig->AddVertex(aIP);
1190 glig->SetLastPoint(2);
fa0291ff 1191 }
1192 //
1193 pmult1.SetParameter(0.5*M_PI);
7fd59977 1194 glig->AddVertex(pmult1);
fa0291ff 1195 //
c6541a0c 1196 pmult2.SetParameter(1.5*M_PI);
7fd59977 1197 glig->AddVertex(pmult2);
fa0291ff 1198
1199 //
7fd59977 1200 slin.Append(glig);
1201
1202 //-- Transitions calculee au point 0 OK
fa0291ff 1203 //
7fd59977 1204 // on traite la deuxieme ellipse
d30895f5 1205 elipsol = theInter.Ellipse(2);
7fd59977 1206
1207 Standard_Real param1 = ElCLib::Parameter(elipsol,pttang1);
1208 Standard_Real param2 = ElCLib::Parameter(elipsol,pttang2);
ecc4f148 1209 Standard_Real parampourtransition = 0.0;
1210 if (param1 < param2)
1211 {
1212 pmult1.SetParameter(0.5*M_PI);
1213 pmult2.SetParameter(1.5*M_PI);
1214 parampourtransition = M_PI;
7fd59977 1215 }
1216 else {
ecc4f148 1217 pmult1.SetParameter(1.5*M_PI);
1218 pmult2.SetParameter(0.5*M_PI);
1219 parampourtransition = 0.0;
7fd59977 1220 }
1221
1222 //-- Calcul des transitions de ligne pour la premiere ligne
1223 ElCLib::D1(parampourtransition,elipsol,ptref,Tgt);
d30895f5 1224
1225 //Theoretically, qwe = +/- |Tgt|.
7fd59977 1226 qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
ecc4f148 1227 if(qwe> 0.00000001)
1228 {
1229 trans1 = IntSurf_Out;
1230 trans2 = IntSurf_In;
7fd59977 1231 }
ecc4f148 1232 else if(qwe< -0.00000001)
1233 {
1234 trans1 = IntSurf_In;
1235 trans2 = IntSurf_Out;
7fd59977 1236 }
ecc4f148 1237 else
1238 {
1239 trans1=trans2=IntSurf_Undecided;
7fd59977 1240 }
ecc4f148 1241
7fd59977 1242 //-- La transition a ete calculee sur un point de cette ligne
1243 glig = new IntPatch_GLine(elipsol,Standard_False,trans1,trans2);
fa0291ff 1244 //
1245 {
ecc4f148 1246 Standard_Real aU1, aV1, aU2, aV2;
1247 IntPatch_Point aIP;
1248 gp_Pnt aP (ElCLib::Value(0., elipsol));
1249 //
1250 aIP.SetValue(aP,Tol,Standard_False);
1251 aIP.SetMultiple(Standard_False);
1252 //
7fd59977 1253
261b7d9e 1254 Quad1.Parameters(aP, aU1, aV1);
1255 Quad2.Parameters(aP, aU2, aV2);
ecc4f148 1256
d30895f5 1257 aIP.SetParameters(aU1, aV1, aU2, aV2);
1258 //
1259 aIP.SetParameter(0.);
1260 glig->AddVertex(aIP);
1261 glig->SetFirstPoint(1);
1262 //
1263 aIP.SetParameter(2.*M_PI);
1264 glig->AddVertex(aIP);
1265 glig->SetLastPoint(2);
1266 }
1267 //
1268 glig->AddVertex(pmult1);
1269 glig->AddVertex(pmult2);
1270 //
1271 slin.Append(glig);
ecc4f148 1272 }
ecc4f148 1273 break;
ecc4f148 1274
d30895f5 1275 case IntAna_Parabola:
1276 case IntAna_Hyperbola:
9775fa61 1277 throw Standard_Failure("IntCyCy(): Wrong intersection type!");
ecc4f148 1278
d30895f5 1279 case IntAna_Circle:
1280 // Circle is useful when we will work with trimmed surfaces
1281 // (two cylinders can be tangent by their basises, e.g. circle)
1282 case IntAna_NoGeometricSolution:
1283 default:
1284 return Standard_False;
1285 }
ecc4f148 1286
d30895f5 1287 return Standard_True;
1288}
ecc4f148 1289
d30895f5 1290//=======================================================================
1291//function : ShortCosForm
1292//purpose : Represents theCosFactor*cosA+theSinFactor*sinA as
1293// theCoeff*cos(A-theAngle) if it is possibly (all angles are
1294// in radians).
1295//=======================================================================
1296static void ShortCosForm( const Standard_Real theCosFactor,
1297 const Standard_Real theSinFactor,
1298 Standard_Real& theCoeff,
1299 Standard_Real& theAngle)
1300{
1301 theCoeff = sqrt(theCosFactor*theCosFactor+theSinFactor*theSinFactor);
1302 theAngle = 0.0;
1303 if(IsEqual(theCoeff, 0.0))
1304 {
1305 theAngle = 0.0;
1306 return;
1307 }
ecc4f148 1308
d30895f5 1309 theAngle = acos(Abs(theCosFactor/theCoeff));
ecc4f148 1310
d30895f5 1311 if(theSinFactor > 0.0)
1312 {
1313 if(IsEqual(theCosFactor, 0.0))
1314 {
1315 theAngle = M_PI/2.0;
1316 }
1317 else if(theCosFactor < 0.0)
1318 {
1319 theAngle = M_PI-theAngle;
1320 }
1321 }
1322 else if(IsEqual(theSinFactor, 0.0))
1323 {
1324 if(theCosFactor < 0.0)
1325 {
1326 theAngle = M_PI;
1327 }
1328 }
1329 if(theSinFactor < 0.0)
1330 {
1331 if(theCosFactor > 0.0)
1332 {
1333 theAngle = 2.0*M_PI-theAngle;
1334 }
1335 else if(IsEqual(theCosFactor, 0.0))
1336 {
1337 theAngle = 3.0*M_PI/2.0;
1338 }
1339 else if(theCosFactor < 0.0)
1340 {
1341 theAngle = M_PI+theAngle;
1342 }
1343 }
ecc4f148 1344}
1345
b5ef9d91 1346//=======================================================================
1347//function : CylCylMonotonicity
1348//purpose : Determines, if U2(U1) function is increasing.
1349//=======================================================================
d30895f5 1350Standard_Boolean ComputationMethods::CylCylMonotonicity(const Standard_Real theU1par,
1351 const Standard_Integer theWLIndex,
1352 const stCoeffsValue& theCoeffs,
1353 const Standard_Real thePeriod,
1354 Standard_Boolean& theIsIncreasing)
b5ef9d91 1355{
1356 // U2(U1) = FI2 + (+/-)acos(B*cos(U1 - FI1) + C);
1357 //Accordingly,
1358 //Func. U2(X1) = FI2 + X1;
1359 //Func. X1(X2) = anArccosFactor*X2;
1360 //Func. X2(X3) = acos(X3);
1361 //Func. X3(X4) = B*X4 + C;
1362 //Func. X4(U1) = cos(U1 - FI1).
1363 //
1364 //Consequently,
1365 //U2(X1) is always increasing.
1366 //X1(X2) is increasing if anArccosFactor > 0.0 and
1367 //is decreasing otherwise.
1368 //X2(X3) is always decreasing.
1369 //Therefore, U2(X3) is decreasing if anArccosFactor > 0.0 and
1370 //is increasing otherwise.
1371 //X3(X4) is increasing if B > 0 and is decreasing otherwise.
1372 //X4(U1) is increasing if U1 - FI1 in [PI, 2*PI) and
1373 //is decreasing U1 - FI1 in [0, PI) or if (U1 - FI1 == 2PI).
1374 //After that, we can predict behaviour of U2(U1) function:
1375 //if it is increasing or decreasing.
1376
1377 //For "+/-" sign. If isPlus == TRUE, "+" is chosen, otherwise, we choose "-".
1378 Standard_Boolean isPlus = Standard_False;
1379
1380 switch(theWLIndex)
1381 {
1382 case 0:
1383 isPlus = Standard_True;
1384 break;
1385 case 1:
1386 isPlus = Standard_False;
1387 break;
1388 default:
9775fa61 1389 //throw Standard_Failure("Error. Range Error!!!!");
b5ef9d91 1390 return Standard_False;
1391 }
1392
1393 Standard_Real aU1Temp = theU1par - theCoeffs.mFI1;
1394 InscribePoint(0, thePeriod, aU1Temp, 0.0, thePeriod, Standard_False);
1395
1396 theIsIncreasing = Standard_True;
1397
1398 if(((M_PI - aU1Temp) < RealSmall()) && (aU1Temp < thePeriod))
1399 {
1400 theIsIncreasing = Standard_False;
1401 }
1402
1403 if(theCoeffs.mB < 0.0)
1404 {
1405 theIsIncreasing = !theIsIncreasing;
1406 }
1407
1408 if(!isPlus)
1409 {
1410 theIsIncreasing = !theIsIncreasing;
1411 }
1412
1413 return Standard_True;
1414}
1415
7365fad6 1416//=======================================================================
1417//function : CylCylComputeParameters
1418//purpose : Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0,
261b7d9e 1419// estimates the tolerance of U2-computing (estimation result is
7365fad6 1420// assigned to *theDelta value).
1421//=======================================================================
d30895f5 1422Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
b5ef9d91 1423 const Standard_Integer theWLIndex,
1424 const stCoeffsValue& theCoeffs,
7365fad6 1425 Standard_Real& theU2,
d30895f5 1426 Standard_Real* const theDelta)
b5ef9d91 1427{
7365fad6 1428 //This formula is got from some experience and can be changed.
1429 const Standard_Real aTol0 = Min(10.0*Epsilon(1.0)*theCoeffs.mB, aNulValue);
1430 const Standard_Real aTol = 1.0 - aTol0;
b5ef9d91 1431
7365fad6 1432 if(theWLIndex < 0 || theWLIndex > 1)
b5ef9d91 1433 return Standard_False;
b5ef9d91 1434
7365fad6 1435 const Standard_Real aSign = theWLIndex ? -1.0 : 1.0;
1436
1437 Standard_Real anArg = cos(theU1par - theCoeffs.mFI1);
1438 anArg = theCoeffs.mB*anArg + theCoeffs.mC;
1439
261b7d9e 1440 if(anArg >= aTol)
7365fad6 1441 {
1442 if(theDelta)
1443 *theDelta = 0.0;
b5ef9d91 1444
b5ef9d91 1445 anArg = 1.0;
7365fad6 1446 }
261b7d9e 1447 else if(anArg <= -aTol)
7365fad6 1448 {
1449 if(theDelta)
1450 *theDelta = 0.0;
b5ef9d91 1451
b5ef9d91 1452 anArg = -1.0;
7365fad6 1453 }
1454 else if(theDelta)
1455 {
1456 //There is a case, when
1457 // const double aPar = 0.99999999999721167;
1458 // const double aFI2 = 2.3593296083566181e-006;
1459
1460 //Then
1461 // aPar - cos(aFI2) == -5.10703e-015 ==> cos(aFI2) == aPar.
0177fe26 1462 //Theoretically, in this case
7365fad6 1463 // aFI2 == +/- acos(aPar).
1464 //However,
1465 // acos(aPar) - aFI2 == 2.16362e-009.
1466 //Error is quite big.
1467
1468 //This error should be estimated. Let use following way, which is based
1469 //on expanding to Taylor series.
1470
1471 // acos(p)-acos(p+x) = x/sqrt(1-p*p).
1472
1473 //If p == (1-d) (when p > 0) or p == (-1+d) (when p < 0) then
1474 // acos(p)-acos(p+x) = x/sqrt(d*(2-d)).
1475
1476 //Here always aTol0 <= d <= 1. Max(x) is considered (!) to be equal to aTol0.
1477 //In this case
1478 // 8*aTol0 <= acos(p)-acos(p+x) <= sqrt(2/(2-aTol0)-1),
1479 // because 0 < aTol0 < 1.
1480 //E.g. when aTol0 = 1.0e-11,
1481 // 8.0e-11 <= acos(p)-acos(p+x) < 2.24e-6.
1482
1483 const Standard_Real aDelta = Min(1.0-anArg, 1.0+anArg);
1484 Standard_DivideByZero_Raise_if((aDelta*aDelta < RealSmall()) || (aDelta >= 2.0),
1485 "IntPatch_ImpImpIntersection_4.gxx, CylCylComputeParameters()");
1486 *theDelta = aTol0/sqrt(aDelta*(2.0-aDelta));
1487 }
b5ef9d91 1488
7365fad6 1489 theU2 = acos(anArg);
1490 theU2 = theCoeffs.mFI2 + aSign*theU2;
b5ef9d91 1491
1492 return Standard_True;
1493}
1494
261b7d9e 1495//=======================================================================
1496//function : CylCylComputeParameters
1497//purpose : Computes V1 and V2 (V-parameters of the 1st and 2nd cylinder respectively).
1498//=======================================================================
d30895f5 1499Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1,
b5ef9d91 1500 const Standard_Real theU2,
1501 const stCoeffsValue& theCoeffs,
1502 Standard_Real& theV1,
1503 Standard_Real& theV2)
1504{
1505 theV1 = theCoeffs.mK21 * sin(theU2) +
1506 theCoeffs.mK11 * sin(theU1) +
1507 theCoeffs.mL21 * cos(theU2) +
1508 theCoeffs.mL11 * cos(theU1) + theCoeffs.mM1;
1509
1510 theV2 = theCoeffs.mK22 * sin(theU2) +
1511 theCoeffs.mK12 * sin(theU1) +
1512 theCoeffs.mL22 * cos(theU2) +
1513 theCoeffs.mL12 * cos(theU1) + theCoeffs.mM2;
1514
1515 return Standard_True;
1516}
1517
261b7d9e 1518//=======================================================================
1519//function : CylCylComputeParameters
1520//purpose : Computes U2 (U-parameter of the 2nd cylinder),
1521// V1 and V2 (V-parameters of the 1st and 2nd cylinder respectively).
1522//=======================================================================
d30895f5 1523Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
b5ef9d91 1524 const Standard_Integer theWLIndex,
1525 const stCoeffsValue& theCoeffs,
1526 Standard_Real& theU2,
1527 Standard_Real& theV1,
1528 Standard_Real& theV2)
1529{
1530 if(!CylCylComputeParameters(theU1par, theWLIndex, theCoeffs, theU2))
1531 return Standard_False;
1532
1533 if(!CylCylComputeParameters(theU1par, theU2, theCoeffs, theV1, theV2))
1534 return Standard_False;
1535
1536 return Standard_True;
1537}
1538
ecc4f148 1539//=======================================================================
1540//function : SearchOnVBounds
1541//purpose :
1542//=======================================================================
d30895f5 1543Standard_Boolean WorkWithBoundaries::
1544 SearchOnVBounds(const SearchBoundType theSBType,
ecc4f148 1545 const Standard_Real theVzad,
b5ef9d91 1546 const Standard_Real theVInit,
ecc4f148 1547 const Standard_Real theInitU2,
1548 const Standard_Real theInitMainVar,
d30895f5 1549 Standard_Real& theMainVariableValue) const
ecc4f148 1550{
b5ef9d91 1551 const Standard_Integer aNbDim = 3;
ecc4f148 1552 const Standard_Real aMaxError = 4.0*M_PI; // two periods
1553
b5ef9d91 1554 theMainVariableValue = theInitMainVar;
1555 const Standard_Real aTol2 = 1.0e-18;
1556 Standard_Real aMainVarPrev = theInitMainVar, aU2Prev = theInitU2, anOtherVar = theVInit;
ecc4f148 1557
b5ef9d91 1558 //Structure of aMatr:
1559 // C_{1}*U_{1} & C_{2}*U_{2} & C_{3}*V_{*},
1560 //where C_{1}, C_{2} and C_{3} are math_Vector.
1561 math_Matrix aMatr(1, aNbDim, 1, aNbDim);
1562
ecc4f148 1563 Standard_Real anError = RealLast();
b5ef9d91 1564 Standard_Real anErrorPrev = anError;
ecc4f148 1565 Standard_Integer aNbIter = 0;
1566 do
1567 {
1568 if(++aNbIter > 1000)
1569 return Standard_False;
1570
b5ef9d91 1571 const Standard_Real aSinU1 = sin(aMainVarPrev),
1572 aCosU1 = cos(aMainVarPrev),
1573 aSinU2 = sin(aU2Prev),
1574 aCosU2 = cos(aU2Prev);
ecc4f148 1575
d30895f5 1576 math_Vector aVecFreeMem = (myCoeffs.mVecA2 * aU2Prev +
1577 myCoeffs.mVecB2) * aSinU2 -
1578 (myCoeffs.mVecB2 * aU2Prev -
1579 myCoeffs.mVecA2) * aCosU2 +
1580 (myCoeffs.mVecA1 * aMainVarPrev +
1581 myCoeffs.mVecB1) * aSinU1 -
1582 (myCoeffs.mVecB1 * aMainVarPrev -
1583 myCoeffs.mVecA1) * aCosU1 +
1584 myCoeffs.mVecD;
b5ef9d91 1585
1586 math_Vector aMSum(1, 3);
ecc4f148 1587
1588 switch(theSBType)
1589 {
1590 case SearchV1:
d30895f5 1591 aMatr.SetCol(3, myCoeffs.mVecC2);
1592 aMSum = myCoeffs.mVecC1 * theVzad;
b5ef9d91 1593 aVecFreeMem -= aMSum;
d30895f5 1594 aMSum += myCoeffs.mVecC2*anOtherVar;
ecc4f148 1595 break;
1596
1597 case SearchV2:
d30895f5 1598 aMatr.SetCol(3, myCoeffs.mVecC1);
1599 aMSum = myCoeffs.mVecC2 * theVzad;
b5ef9d91 1600 aVecFreeMem -= aMSum;
d30895f5 1601 aMSum += myCoeffs.mVecC1*anOtherVar;
ecc4f148 1602 break;
1603
1604 default:
1605 return Standard_False;
1606 }
1607
d30895f5 1608 aMatr.SetCol(1, myCoeffs.mVecA1 * aSinU1 - myCoeffs.mVecB1 * aCosU1);
1609 aMatr.SetCol(2, myCoeffs.mVecA2 * aSinU2 - myCoeffs.mVecB2 * aCosU2);
b5ef9d91 1610
1611 Standard_Real aDetMainSyst = aMatr.Determinant();
ecc4f148 1612
b5ef9d91 1613 if(Abs(aDetMainSyst) < aNulValue)
ecc4f148 1614 {
1615 return Standard_False;
1616 }
1617
b5ef9d91 1618 math_Matrix aM1(aMatr), aM2(aMatr), aM3(aMatr);
1619 aM1.SetCol(1, aVecFreeMem);
1620 aM2.SetCol(2, aVecFreeMem);
1621 aM3.SetCol(3, aVecFreeMem);
ecc4f148 1622
b5ef9d91 1623 const Standard_Real aDetMainVar = aM1.Determinant();
1624 const Standard_Real aDetVar1 = aM2.Determinant();
1625 const Standard_Real aDetVar2 = aM3.Determinant();
ecc4f148 1626
1627 Standard_Real aDelta = aDetMainVar/aDetMainSyst-aMainVarPrev;
1628
1629 if(Abs(aDelta) > aMaxError)
1630 return Standard_False;
1631
1632 anError = aDelta*aDelta;
1633 aMainVarPrev += aDelta;
1634
1635 ///
1636 aDelta = aDetVar1/aDetMainSyst-aU2Prev;
1637
1638 if(Abs(aDelta) > aMaxError)
1639 return Standard_False;
1640
1641 anError += aDelta*aDelta;
1642 aU2Prev += aDelta;
1643
1644 ///
1645 aDelta = aDetVar2/aDetMainSyst-anOtherVar;
1646 anError += aDelta*aDelta;
1647 anOtherVar += aDelta;
b5ef9d91 1648
1649 if(anError > anErrorPrev)
1650 {//Method diverges. Keep the best result
51740958 1651 const Standard_Real aSinU1Last = sin(aMainVarPrev),
1652 aCosU1Last = cos(aMainVarPrev),
1653 aSinU2Last = sin(aU2Prev),
1654 aCosU2Last = cos(aU2Prev);
d30895f5 1655 aMSum -= (myCoeffs.mVecA1*aCosU1Last +
1656 myCoeffs.mVecB1*aSinU1Last +
1657 myCoeffs.mVecA2*aCosU2Last +
1658 myCoeffs.mVecB2*aSinU2Last +
1659 myCoeffs.mVecD);
b5ef9d91 1660 const Standard_Real aSQNorm = aMSum.Norm2();
1661 return (aSQNorm < aTol2);
1662 }
1663 else
1664 {
1665 theMainVariableValue = aMainVarPrev;
1666 }
1667
1668 anErrorPrev = anError;
ecc4f148 1669 }
1670 while(anError > aTol2);
1671
1672 theMainVariableValue = aMainVarPrev;
1673
1674 return Standard_True;
1675}
1676
1677//=======================================================================
1678//function : InscribePoint
e002f1ce 1679//purpose : If theFlForce==TRUE theUGiven will be changed forcefully
1680// even if theUGiven is already inscibed in the boundary
1681// (if it is possible; i.e. if new theUGiven is inscribed
1682// in the boundary, too).
ecc4f148 1683//=======================================================================
b5ef9d91 1684Standard_Boolean InscribePoint( const Standard_Real theUfTarget,
1685 const Standard_Real theUlTarget,
1686 Standard_Real& theUGiven,
1687 const Standard_Real theTol2D,
1688 const Standard_Real thePeriod,
1689 const Standard_Boolean theFlForce)
ecc4f148 1690{
7365fad6 1691 if(Precision::IsInfinite(theUGiven))
1692 {
1693 return Standard_False;
1694 }
1695
02effd35 1696 if((theUfTarget - theUGiven <= theTol2D) &&
1697 (theUGiven - theUlTarget <= theTol2D))
ecc4f148 1698 {//it has already inscribed
1699
1700 /*
1701 Utf U Utl
1702 + * +
1703 */
1704
02effd35 1705 if(theFlForce)
1706 {
1707 Standard_Real anUtemp = theUGiven + thePeriod;
1708 if((theUfTarget - anUtemp <= theTol2D) &&
1709 (anUtemp - theUlTarget <= theTol2D))
1710 {
1711 theUGiven = anUtemp;
1712 return Standard_True;
1713 }
1714
1715 anUtemp = theUGiven - thePeriod;
1716 if((theUfTarget - anUtemp <= theTol2D) &&
1717 (anUtemp - theUlTarget <= theTol2D))
1718 {
1719 theUGiven = anUtemp;
1720 }
1721 }
1722
ecc4f148 1723 return Standard_True;
1724 }
1725
d30895f5 1726 const Standard_Real aUf = theUfTarget - theTol2D;
1727 const Standard_Real aUl = aUf + thePeriod;
ecc4f148 1728
d30895f5 1729 theUGiven = ElCLib::InPeriod(theUGiven, aUf, aUl);
ecc4f148 1730
7365fad6 1731 return ((theUfTarget - theUGiven <= theTol2D) &&
1732 (theUGiven - theUlTarget <= theTol2D));
ecc4f148 1733}
1734
1735//=======================================================================
1736//function : InscribeInterval
261b7d9e 1737//purpose : Shifts theRange in order to make at least one of its
1738// boundary in the range [theUfTarget, theUlTarget]
ecc4f148 1739//=======================================================================
1740static Standard_Boolean InscribeInterval(const Standard_Real theUfTarget,
261b7d9e 1741 const Standard_Real theUlTarget,
1742 Bnd_Range &theRange,
1743 const Standard_Real theTol2D,
1744 const Standard_Real thePeriod)
ecc4f148 1745{
261b7d9e 1746 Standard_Real anUpar = 0.0;
1747 if (!theRange.GetMin(anUpar))
1748 {
1749 return Standard_False;
1750 }
1751
1752 const Standard_Real aDelta = theRange.Delta();
1753 if(InscribePoint(theUfTarget, theUlTarget, anUpar,
1754 theTol2D, thePeriod, (Abs(theUlTarget-anUpar) < theTol2D)))
ecc4f148 1755 {
261b7d9e 1756 theRange.SetVoid();
1757 theRange.Add(anUpar);
1758 theRange.Add(anUpar + aDelta);
ecc4f148 1759 }
1760 else
1761 {
8662560e 1762 if (!theRange.GetMax (anUpar))
261b7d9e 1763 {
1764 return Standard_False;
1765 }
1766
02effd35 1767 if(InscribePoint(theUfTarget, theUlTarget, anUpar,
261b7d9e 1768 theTol2D, thePeriod, (Abs(theUfTarget-anUpar) < theTol2D)))
ecc4f148 1769 {
261b7d9e 1770 theRange.SetVoid();
1771 theRange.Add(anUpar);
1772 theRange.Add(anUpar - aDelta);
ecc4f148 1773 }
1774 else
7fd59977 1775 {
1776 return Standard_False;
1777 }
1778 }
ecc4f148 1779
1780 return Standard_True;
1781}
1782
1783//=======================================================================
7365fad6 1784//function : ExcludeNearElements
1785//purpose : Checks if theArr contains two almost equal elements.
1786// If it is true then one of equal elements will be excluded
1787// (made infinite).
1788// Returns TRUE if at least one element of theArr has been changed.
1789//ATTENTION!!!
1790// 1. Every not infinite element of theArr is considered to be
1791// in [0, T] interval (where T is the period);
1792// 2. theArr must be sorted in ascending order.
ecc4f148 1793//=======================================================================
7365fad6 1794static Standard_Boolean ExcludeNearElements(Standard_Real theArr[],
1795 const Standard_Integer theNOfMembers,
79997052 1796 const Standard_Real theUSurf1f,
1797 const Standard_Real theUSurf1l,
7365fad6 1798 const Standard_Real theTol)
ecc4f148 1799{
7365fad6 1800 Standard_Boolean aRetVal = Standard_False;
1801 for(Standard_Integer i = 1; i < theNOfMembers; i++)
ecc4f148 1802 {
7365fad6 1803 Standard_Real &anA = theArr[i],
1804 &anB = theArr[i-1];
ecc4f148 1805
7365fad6 1806 //Here, anA >= anB
ecc4f148 1807
7365fad6 1808 if(Precision::IsInfinite(anA))
1809 break;
ecc4f148 1810
7365fad6 1811 if((anA-anB) < theTol)
ecc4f148 1812 {
79997052 1813 if((anB != 0.0) && (anB != theUSurf1f) && (anB != theUSurf1l))
d30895f5 1814 anA = (anA + anB)/2.0;
79997052 1815 else
1816 anA = anB;
ecc4f148 1817
7365fad6 1818 //Make this element infinite an forget it
1819 //(we will not use it in next iterations).
1820 anB = Precision::Infinite();
1821 aRetVal = Standard_True;
ecc4f148 1822 }
1823 }
ecc4f148 1824
7365fad6 1825 return aRetVal;
1826}
ecc4f148 1827
1828//=======================================================================
1829//function : AddPointIntoWL
1830//purpose : Surf1 is a surface, whose U-par is variable.
261b7d9e 1831// If theFlBefore == TRUE then we enable the U1-parameter
1832// of the added point to be less than U1-parameter of
1833// previously added point (in general U1-parameter is
1834// always increased; therefore, this situation is abnormal).
1835// If theOnlyCheck==TRUE then no point will be added to theLine.
ecc4f148 1836//=======================================================================
1837static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
1838 const IntSurf_Quadric& theQuad2,
d30895f5 1839 const ComputationMethods::stCoeffsValue& theCoeffs,
ecc4f148 1840 const Standard_Boolean isTheReverse,
b5ef9d91 1841 const Standard_Boolean isThePrecise,
ecc4f148 1842 const gp_Pnt2d& thePntOnSurf1,
1843 const gp_Pnt2d& thePntOnSurf2,
1844 const Standard_Real theUfSurf1,
1845 const Standard_Real theUlSurf1,
b5ef9d91 1846 const Standard_Real theUfSurf2,
1847 const Standard_Real theUlSurf2,
1848 const Standard_Real theVfSurf1,
1849 const Standard_Real theVlSurf1,
1850 const Standard_Real theVfSurf2,
1851 const Standard_Real theVlSurf2,
ecc4f148 1852 const Standard_Real thePeriodOfSurf1,
1853 const Handle(IntSurf_LineOn2S)& theLine,
b5ef9d91 1854 const Standard_Integer theWLIndex,
e8feb725 1855 const Standard_Real theTol3D,
02effd35 1856 const Standard_Real theTol2D,
261b7d9e 1857 const Standard_Boolean theFlBefore = Standard_False,
1858 const Standard_Boolean theOnlyCheck = Standard_False)
ecc4f148 1859{
4dba155d 1860 //Check if the point is in the domain or can be inscribed in the domain after adjusting.
1861
ecc4f148 1862 gp_Pnt aPt1(theQuad1.Value(thePntOnSurf1.X(), thePntOnSurf1.Y())),
1863 aPt2(theQuad2.Value(thePntOnSurf2.X(), thePntOnSurf2.Y()));
1864
b5ef9d91 1865 Standard_Real aU1par = thePntOnSurf1.X();
261b7d9e 1866
1867 // aU1par always increases. Therefore, we must reduce its
1868 // value in order to continue creation of WLine.
b5ef9d91 1869 if(!InscribePoint(theUfSurf1, theUlSurf1, aU1par, theTol2D,
261b7d9e 1870 thePeriodOfSurf1, aU1par > 0.5*(theUfSurf1+theUlSurf1)))
ecc4f148 1871 return Standard_False;
1872
261b7d9e 1873 if ((theLine->NbPoints() > 0) &&
dccf8675 1874 ((theUlSurf1 - theUfSurf1) >= (thePeriodOfSurf1 - theTol2D)) &&
261b7d9e 1875 (((aU1par + thePeriodOfSurf1 - theUlSurf1) <= theTol2D) ||
1876 ((aU1par - thePeriodOfSurf1 - theUfSurf1) >= theTol2D)))
1877 {
1878 // aU1par can be adjusted to both theUlSurf1 and theUfSurf1
1879 // with equal possibilities. This code fragment allows choosing
1880 // correct parameter from these two variants.
1881
1882 Standard_Real aU1 = 0.0, aV1 = 0.0;
1883 if (isTheReverse)
1884 {
1885 theLine->Value(theLine->NbPoints()).ParametersOnS2(aU1, aV1);
1886 }
1887 else
1888 {
1889 theLine->Value(theLine->NbPoints()).ParametersOnS1(aU1, aV1);
1890 }
1891
1892 const Standard_Real aDelta = aU1 - aU1par;
1893 if (2.0*Abs(aDelta) > thePeriodOfSurf1)
1894 {
1895 aU1par += Sign(thePeriodOfSurf1, aDelta);
1896 }
1897 }
1898
b5ef9d91 1899 Standard_Real aU2par = thePntOnSurf2.X();
1900 if(!InscribePoint(theUfSurf2, theUlSurf2, aU2par, theTol2D,
1901 thePeriodOfSurf1, Standard_False))
1902 return Standard_False;
1903
1904 Standard_Real aV1par = thePntOnSurf1.Y();
1905 if((aV1par - theVlSurf1 > theTol2D) || (theVfSurf1 - aV1par > theTol2D))
1906 return Standard_False;
1907
1908 Standard_Real aV2par = thePntOnSurf2.Y();
1909 if((aV2par - theVlSurf2 > theTol2D) || (theVfSurf2 - aV2par > theTol2D))
1910 return Standard_False;
1911
4dba155d 1912 //Get intersection point and add it in the WL
ecc4f148 1913 IntSurf_PntOn2S aPnt;
1914
1915 if(isTheReverse)
1916 {
1917 aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
b5ef9d91 1918 aU2par, aV2par,
1919 aU1par, aV1par);
ecc4f148 1920 }
1921 else
1922 {
1923 aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
b5ef9d91 1924 aU1par, aV1par,
1925 aU2par, aV2par);
ecc4f148 1926 }
1927
b5ef9d91 1928 Standard_Integer aNbPnts = theLine->NbPoints();
e8feb725 1929 if(aNbPnts > 0)
1930 {
1931 Standard_Real aUl = 0.0, aVl = 0.0;
1932 const IntSurf_PntOn2S aPlast = theLine->Value(aNbPnts);
1933 if(isTheReverse)
1934 aPlast.ParametersOnS2(aUl, aVl);
1935 else
1936 aPlast.ParametersOnS1(aUl, aVl);
1937
b5ef9d91 1938 if(!theFlBefore && (aU1par <= aUl))
d30895f5 1939 {
1940 //Parameter value must be increased if theFlBefore == FALSE.
261b7d9e 1941
d30895f5 1942 aU1par += thePeriodOfSurf1;
1943
1944 //The condition is as same as in
1945 //InscribePoint(...) function
1946 if((theUfSurf1 - aU1par > theTol2D) ||
1947 (aU1par - theUlSurf1 > theTol2D))
1948 {
1949 //New aU1par is out of target interval.
1950 //Go back to old value.
1951
1952 return Standard_False;
1953 }
e8feb725 1954 }
1955
261b7d9e 1956 if (theOnlyCheck)
1957 return Standard_True;
1958
e8feb725 1959 //theTol2D is minimal step along parameter changed.
1960 //Therefore, if we apply this minimal step two
1961 //neighbour points will be always "same". Consequently,
1962 //we should reduce tolerance for IsSame checking.
1963 const Standard_Real aDTol = 1.0-Epsilon(1.0);
1964 if(aPnt.IsSame(aPlast, theTol3D*aDTol, theTol2D*aDTol))
1965 {
1966 theLine->RemovePoint(aNbPnts);
1967 }
1968 }
1969
261b7d9e 1970 if (theOnlyCheck)
1971 return Standard_True;
1972
ecc4f148 1973 theLine->Add(aPnt);
b5ef9d91 1974
1975 if(!isThePrecise)
1976 return Standard_True;
1977
1978 //Try to precise existing WLine
1979 aNbPnts = theLine->NbPoints();
1980 if(aNbPnts >= 3)
1981 {
1982 Standard_Real aU1 = 0.0, aU2 = 0.0, aU3 = 0.0, aV = 0.0;
1983 if(isTheReverse)
1984 {
1985 theLine->Value(aNbPnts).ParametersOnS2(aU3, aV);
1986 theLine->Value(aNbPnts-1).ParametersOnS2(aU2, aV);
1987 theLine->Value(aNbPnts-2).ParametersOnS2(aU1, aV);
1988 }
1989 else
1990 {
1991 theLine->Value(aNbPnts).ParametersOnS1(aU3, aV);
1992 theLine->Value(aNbPnts-1).ParametersOnS1(aU2, aV);
1993 theLine->Value(aNbPnts-2).ParametersOnS1(aU1, aV);
1994 }
1995
1996 const Standard_Real aStepPrev = aU2-aU1;
1997 const Standard_Real aStep = aU3-aU2;
1998
1999 const Standard_Integer aDeltaStep = RealToInt(aStepPrev/aStep);
2000
2001 if((1 < aDeltaStep) && (aDeltaStep < 2000))
2002 {
4dba155d 2003 //Add new points in case of non-uniform distribution of existing points
b5ef9d91 2004 SeekAdditionalPoints( theQuad1, theQuad2, theLine,
2005 theCoeffs, theWLIndex, aDeltaStep, aNbPnts-2,
d30895f5 2006 aNbPnts-1, theTol2D, thePeriodOfSurf1, isTheReverse);
b5ef9d91 2007 }
2008 }
2009
ecc4f148 2010 return Standard_True;
2011}
2012
2013//=======================================================================
2014//function : AddBoundaryPoint
4dba155d 2015//purpose : Find intersection point on V-boundary.
ecc4f148 2016//=======================================================================
d30895f5 2017void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL,
ecc4f148 2018 const Standard_Real theU1,
261b7d9e 2019 const Standard_Real theU1Min,
ecc4f148 2020 const Standard_Real theU2,
2021 const Standard_Real theV1,
2022 const Standard_Real theV1Prev,
2023 const Standard_Real theV2,
2024 const Standard_Real theV2Prev,
b5ef9d91 2025 const Standard_Integer theWLIndex,
02effd35 2026 const Standard_Boolean theFlForce,
ecc4f148 2027 Standard_Boolean& isTheFound1,
d30895f5 2028 Standard_Boolean& isTheFound2) const
ecc4f148 2029{
2030 Standard_Real aUSurf1f = 0.0, //const
2031 aUSurf1l = 0.0,
2032 aVSurf1f = 0.0,
2033 aVSurf1l = 0.0;
2034 Standard_Real aUSurf2f = 0.0, //const
2035 aUSurf2l = 0.0,
2036 aVSurf2f = 0.0,
2037 aVSurf2l = 0.0;
2038
d30895f5 2039 myUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
2040 myUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
2041
2042 const Standard_Integer aSize = 4;
2043 const Standard_Real anArrVzad[aSize] = {aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l};
ecc4f148 2044
d30895f5 2045 StPInfo aUVPoint[aSize];
b5ef9d91 2046
d30895f5 2047 for(Standard_Integer anIDSurf = 0; anIDSurf < 4; anIDSurf+=2)
ecc4f148 2048 {
d30895f5 2049 const Standard_Real aVf = (anIDSurf == 0) ? theV1 : theV2,
2050 aVl = (anIDSurf == 0) ? theV1Prev : theV2Prev;
ecc4f148 2051
d30895f5 2052 const SearchBoundType aTS = (anIDSurf == 0) ? SearchV1 : SearchV2;
ecc4f148 2053
d30895f5 2054 for(Standard_Integer anIDBound = 0; anIDBound < 2; anIDBound++)
2055 {
d30895f5 2056 const Standard_Integer anIndex = anIDSurf+anIDBound;
b5ef9d91 2057
d30895f5 2058 aUVPoint[anIndex].mySurfID = anIDSurf;
7365fad6 2059
d30895f5 2060 if((Abs(aVf-anArrVzad[anIndex]) > myTol2D) &&
2061 ((aVf-anArrVzad[anIndex])*(aVl-anArrVzad[anIndex]) > 0.0))
ecc4f148 2062 {
d30895f5 2063 continue;
ecc4f148 2064 }
ecc4f148 2065
4dba155d 2066 //Segment [aVf, aVl] intersects at least one V-boundary (first or last)
2067 // (in general, case is possible, when aVf > aVl).
2068
2069 // Precise intersection point
d30895f5 2070 const Standard_Boolean aRes = SearchOnVBounds(aTS, anArrVzad[anIndex],
2071 (anIDSurf == 0) ? theV2 : theV1,
2072 theU2, theU1,
2073 aUVPoint[anIndex].myU1);
2074
261b7d9e 2075 // aUVPoint[anIndex].myU1 is considered to be nearer to theU1 than
2076 // to theU1+/-Period
2077 if (!aRes || (aUVPoint[anIndex].myU1 >= theU1) ||
2078 (aUVPoint[anIndex].myU1 < theU1Min))
ecc4f148 2079 {
4dba155d 2080 //Intersection point is not found or out of the domain
d30895f5 2081 aUVPoint[anIndex].myU1 = RealLast();
2082 continue;
ecc4f148 2083 }
d30895f5 2084 else
ecc4f148 2085 {
4dba155d 2086 //intersection point is found
2087
d30895f5 2088 Standard_Real &aU1 = aUVPoint[anIndex].myU1,
2089 &aU2 = aUVPoint[anIndex].myU2,
2090 &aV1 = aUVPoint[anIndex].myV1,
2091 &aV2 = aUVPoint[anIndex].myV2;
2092 aU2 = theU2;
2093 aV1 = theV1;
2094 aV2 = theV2;
2095
2096 if(!ComputationMethods::
2097 CylCylComputeParameters(aU1, theWLIndex, myCoeffs, aU2, aV1, aV2))
2098 {
4dba155d 2099 // Found point is wrong
d30895f5 2100 aU1 = RealLast();
2101 continue;
2102 }
2103
4dba155d 2104 //Point on true V-boundary.
d30895f5 2105 if(aTS == SearchV1)
2106 aV1 = anArrVzad[anIndex];
2107 else //if(aTS[anIndex] == SearchV2)
2108 aV2 = anArrVzad[anIndex];
ecc4f148 2109 }
2110 }
2111 }
b5ef9d91 2112
4dba155d 2113 // Sort with acceding U1-parameter.
d30895f5 2114 std::sort(aUVPoint, aUVPoint+aSize);
2115
2116 isTheFound1 = isTheFound2 = Standard_False;
b5ef9d91 2117
4dba155d 2118 //Adding found points on boundary in the WLine.
d30895f5 2119 for(Standard_Integer i = 0; i < aSize; i++)
2120 {
2121 if(aUVPoint[i].myU1 == RealLast())
2122 break;
b5ef9d91 2123
d30895f5 2124 if(!AddPointIntoWL(myQuad1, myQuad2, myCoeffs, myIsReverse, Standard_False,
2125 gp_Pnt2d(aUVPoint[i].myU1, aUVPoint[i].myV1),
2126 gp_Pnt2d(aUVPoint[i].myU2, aUVPoint[i].myV2),
2127 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
2128 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, myPeriod,
2129 theWL->Curve(), theWLIndex, myTol3D, myTol2D, theFlForce))
7365fad6 2130 {
d30895f5 2131 continue;
7365fad6 2132 }
ecc4f148 2133
d30895f5 2134 if(aUVPoint[i].mySurfID == 0)
ecc4f148 2135 {
d30895f5 2136 isTheFound1 = Standard_True;
ecc4f148 2137 }
7365fad6 2138 else
2139 {
d30895f5 2140 isTheFound2 = Standard_True;
7365fad6 2141 }
ecc4f148 2142 }
7fd59977 2143}
2144
ecc4f148 2145//=======================================================================
2146//function : SeekAdditionalPoints
4dba155d 2147//purpose : Inserts additional intersection points between neighbor points.
2148// This action is bone with several iterations. During every iteration,
2149// new point is inserted in middle of every interval.
2150// The process will be finished if NbPoints >= theMinNbPoints.
ecc4f148 2151//=======================================================================
2152static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
2153 const IntSurf_Quadric& theQuad2,
b5ef9d91 2154 const Handle(IntSurf_LineOn2S)& theLine,
d30895f5 2155 const ComputationMethods::stCoeffsValue& theCoeffs,
b5ef9d91 2156 const Standard_Integer theWLIndex,
ecc4f148 2157 const Standard_Integer theMinNbPoints,
b5ef9d91 2158 const Standard_Integer theStartPointOnLine,
2159 const Standard_Integer theEndPointOnLine,
ecc4f148 2160 const Standard_Real theTol2D,
2161 const Standard_Real thePeriodOfSurf2,
ecc4f148 2162 const Standard_Boolean isTheReverse)
2163{
b5ef9d91 2164 if(theLine.IsNull())
2165 return;
2166
2167 Standard_Integer aNbPoints = theEndPointOnLine - theStartPointOnLine + 1;
ecc4f148 2168
b5ef9d91 2169 Standard_Real aMinDeltaParam = theTol2D;
2170
2171 {
2172 Standard_Real u1 = 0.0, v1 = 0.0, u2 = 0.0, v2 = 0.0;
2173
2174 if(isTheReverse)
2175 {
2176 theLine->Value(theStartPointOnLine).ParametersOnS2(u1, v1);
2177 theLine->Value(theEndPointOnLine).ParametersOnS2(u2, v2);
2178 }
2179 else
2180 {
2181 theLine->Value(theStartPointOnLine).ParametersOnS1(u1, v1);
2182 theLine->Value(theEndPointOnLine).ParametersOnS1(u2, v2);
2183 }
2184
2185 aMinDeltaParam = Max(Abs(u2 - u1)/IntToReal(theMinNbPoints), aMinDeltaParam);
2186 }
2187
2188 Standard_Integer aLastPointIndex = theEndPointOnLine;
ecc4f148 2189 Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0;
2190
2191 Standard_Integer aNbPointsPrev = 0;
cd0705f6 2192 do
ecc4f148 2193 {
2194 aNbPointsPrev = aNbPoints;
b5ef9d91 2195 for(Standard_Integer fp = theStartPointOnLine, lp = 0; fp < aLastPointIndex; fp = lp + 1)
ecc4f148 2196 {
02effd35 2197 Standard_Real U1f = 0.0, V1f = 0.0; //first point in 1st suraface
2198 Standard_Real U1l = 0.0, V1l = 0.0; //last point in 1st suraface
2199
2200 Standard_Real U2f = 0.0, V2f = 0.0; //first point in 2nd suraface
2201 Standard_Real U2l = 0.0, V2l = 0.0; //last point in 2nd suraface
ecc4f148 2202
2203 lp = fp+1;
2204
2205 if(isTheReverse)
2206 {
b5ef9d91 2207 theLine->Value(fp).ParametersOnS2(U1f, V1f);
2208 theLine->Value(lp).ParametersOnS2(U1l, V1l);
02effd35 2209
b5ef9d91 2210 theLine->Value(fp).ParametersOnS1(U2f, V2f);
2211 theLine->Value(lp).ParametersOnS1(U2l, V2l);
ecc4f148 2212 }
2213 else
2214 {
b5ef9d91 2215 theLine->Value(fp).ParametersOnS1(U1f, V1f);
2216 theLine->Value(lp).ParametersOnS1(U1l, V1l);
02effd35 2217
b5ef9d91 2218 theLine->Value(fp).ParametersOnS2(U2f, V2f);
2219 theLine->Value(lp).ParametersOnS2(U2l, V2l);
02effd35 2220 }
2221
b5ef9d91 2222 if(Abs(U1l - U1f) <= aMinDeltaParam)
02effd35 2223 {
2224 //Step is minimal. It is not necessary to divide it.
2225 continue;
ecc4f148 2226 }
2227
2228 U1prec = 0.5*(U1f+U1l);
2229
d30895f5 2230 if(!ComputationMethods::
2231 CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec))
2232 {
b5ef9d91 2233 continue;
d30895f5 2234 }
ecc4f148 2235
d30895f5 2236 MinMax(U2f, U2l);
2237 if(!InscribePoint(U2f, U2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False))
2238 {
2239 continue;
2240 }
ecc4f148 2241
b5ef9d91 2242 const gp_Pnt aP1(theQuad1.Value(U1prec, V1prec)), aP2(theQuad2.Value(U2prec, V2prec));
2243 const gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
ecc4f148 2244
d30895f5 2245#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
2246 std::cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << std::endl;
7c32c7c4 2247#endif
2248
ecc4f148 2249 IntSurf_PntOn2S anIP;
2250 if(isTheReverse)
2251 {
2252 anIP.SetValue(aPInt, U2prec, V2prec, U1prec, V1prec);
2253 }
2254 else
2255 {
2256 anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec);
2257 }
2258
b5ef9d91 2259 theLine->InsertBefore(lp, anIP);
ecc4f148 2260
b5ef9d91 2261 aNbPoints++;
2262 aLastPointIndex++;
2263 }
2264
2265 if(aNbPoints >= theMinNbPoints)
2266 {
2267 return;
ecc4f148 2268 }
cd0705f6 2269 } while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev));
ecc4f148 2270}
2271
79997052 2272//=======================================================================
2273//function : BoundariesComputing
2274//purpose : Computes true domain of future intersection curve (allows
2275// avoiding excess iterations)
2276//=======================================================================
261b7d9e 2277Standard_Boolean WorkWithBoundaries::
2278 BoundariesComputing(const ComputationMethods::stCoeffsValue &theCoeffs,
2279 const Standard_Real thePeriod,
2280 Bnd_Range theURange[])
79997052 2281{
4dba155d 2282 //All comments to this method is related to the comment
2283 //to ComputationMethods class
2284
2285 //So, we have the equation
2286 // cos(U2-FI2)=B*cos(U1-FI1)+C
2287 //Evidently,
2288 // -1 <= B*cos(U1-FI1)+C <= 1
2289
261b7d9e 2290 if (theCoeffs.mB > 0.0)
79997052 2291 {
4dba155d 2292 // -(1+C)/B <= cos(U1-FI1) <= (1-C)/B
2293
261b7d9e 2294 if (theCoeffs.mB + Abs(theCoeffs.mC) < -1.0)
4dba155d 2295 {
2296 //(1-C)/B < -1 or -(1+C)/B > 1 ==> No solution
2297
79997052 2298 return Standard_False;
2299 }
261b7d9e 2300 else if (theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
4dba155d 2301 {
2302 //(1-C)/B >= 1 and -(1+C)/B <= -1 ==> U=[0;2*PI]+aFI1
261b7d9e 2303 theURange[0].Add(theCoeffs.mFI1);
2304 theURange[0].Add(thePeriod + theCoeffs.mFI1);
79997052 2305 }
261b7d9e 2306 else if ((1 + theCoeffs.mC <= theCoeffs.mB) &&
2307 (theCoeffs.mB <= 1 - theCoeffs.mC))
79997052 2308 {
4dba155d 2309 //(1-C)/B >= 1 and -(1+C)/B >= -1 ==>
2310 //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1),
2311 //where aDAngle = acos(-(myCoeffs.mC + 1) / myCoeffs.mB)
2312
261b7d9e 2313 Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
79997052 2314 if(anArg > 1.0)
2315 anArg = 1.0;
2316 if(anArg < -1.0)
2317 anArg = -1.0;
2318
2319 const Standard_Real aDAngle = acos(anArg);
261b7d9e 2320 theURange[0].Add(theCoeffs.mFI1);
2321 theURange[0].Add(aDAngle + theCoeffs.mFI1);
2322 theURange[1].Add(thePeriod - aDAngle + theCoeffs.mFI1);
2323 theURange[1].Add(thePeriod + theCoeffs.mFI1);
79997052 2324 }
261b7d9e 2325 else if ((1 - theCoeffs.mC <= theCoeffs.mB) &&
2326 (theCoeffs.mB <= 1 + theCoeffs.mC))
79997052 2327 {
4dba155d 2328 //(1-C)/B <= 1 and -(1+C)/B <= -1 ==> U=[aDAngle;2*PI-aDAngle]+aFI1
2329 //where aDAngle = acos((1 - myCoeffs.mC) / myCoeffs.mB)
2330
261b7d9e 2331 Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
79997052 2332 if(anArg > 1.0)
2333 anArg = 1.0;
2334 if(anArg < -1.0)
2335 anArg = -1.0;
2336
2337 const Standard_Real aDAngle = acos(anArg);
261b7d9e 2338 theURange[0].Add(aDAngle + theCoeffs.mFI1);
2339 theURange[0].Add(thePeriod - aDAngle + theCoeffs.mFI1);
79997052 2340 }
261b7d9e 2341 else if (theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
79997052 2342 {
4dba155d 2343 //(1-C)/B <= 1 and -(1+C)/B >= -1 ==>
2344 //(U=[aDAngle1;aDAngle2]+aFI1) ||
2345 //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
2346 //where aDAngle1 = acos((1 - myCoeffs.mC) / myCoeffs.mB),
2347 // aDAngle2 = acos(-(myCoeffs.mC + 1) / myCoeffs.mB).
2348
261b7d9e 2349 Standard_Real anArg1 = (1 - theCoeffs.mC) / theCoeffs.mB,
2350 anArg2 = -(theCoeffs.mC + 1) / theCoeffs.mB;
79997052 2351 if(anArg1 > 1.0)
2352 anArg1 = 1.0;
2353 if(anArg1 < -1.0)
2354 anArg1 = -1.0;
2355
2356 if(anArg2 > 1.0)
2357 anArg2 = 1.0;
2358 if(anArg2 < -1.0)
2359 anArg2 = -1.0;
2360
2361 const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
2362 //(U=[aDAngle1;aDAngle2]+aFI1) ||
2363 //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
261b7d9e 2364 theURange[0].Add(aDAngle1 + theCoeffs.mFI1);
2365 theURange[0].Add(aDAngle2 + theCoeffs.mFI1);
2366 theURange[1].Add(thePeriod - aDAngle2 + theCoeffs.mFI1);
2367 theURange[1].Add(thePeriod - aDAngle1 + theCoeffs.mFI1);
79997052 2368 }
2369 else
2370 {
2371 return Standard_False;
2372 }
2373 }
261b7d9e 2374 else if (theCoeffs.mB < 0.0)
79997052 2375 {
4dba155d 2376 // (1-C)/B <= cos(U1-FI1) <= -(1+C)/B
2377
261b7d9e 2378 if (theCoeffs.mB + Abs(theCoeffs.mC) > 1.0)
4dba155d 2379 {
2380 // -(1+C)/B < -1 or (1-C)/B > 1 ==> No solutions
79997052 2381 return Standard_False;
2382 }
261b7d9e 2383 else if (-theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
4dba155d 2384 {
2385 // -(1+C)/B >= 1 and (1-C)/B <= -1 ==> U=[0;2*PI]+aFI1
261b7d9e 2386 theURange[0].Add(theCoeffs.mFI1);
2387 theURange[0].Add(thePeriod + theCoeffs.mFI1);
79997052 2388 }
261b7d9e 2389 else if ((-theCoeffs.mC - 1 <= theCoeffs.mB) &&
2390 (theCoeffs.mB <= theCoeffs.mC - 1))
79997052 2391 {
4dba155d 2392 // -(1+C)/B >= 1 and (1-C)/B >= -1 ==>
2393 //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1),
2394 //where aDAngle = acos((1 - myCoeffs.mC) / myCoeffs.mB)
2395
261b7d9e 2396 Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
79997052 2397 if(anArg > 1.0)
2398 anArg = 1.0;
2399 if(anArg < -1.0)
2400 anArg = -1.0;
2401
2402 const Standard_Real aDAngle = acos(anArg);
261b7d9e 2403 theURange[0].Add(theCoeffs.mFI1);
2404 theURange[0].Add(aDAngle + theCoeffs.mFI1);
2405 theURange[1].Add(thePeriod - aDAngle + theCoeffs.mFI1);
2406 theURange[1].Add(thePeriod + theCoeffs.mFI1);
79997052 2407 }
261b7d9e 2408 else if ((theCoeffs.mC - 1 <= theCoeffs.mB) &&
2409 (theCoeffs.mB <= -theCoeffs.mB - 1))
79997052 2410 {
4dba155d 2411 // -(1+C)/B <= 1 and (1-C)/B <= -1 ==> U=[aDAngle;2*PI-aDAngle]+aFI1,
2412 //where aDAngle = acos(-(myCoeffs.mC + 1) / myCoeffs.mB).
2413
261b7d9e 2414 Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
79997052 2415 if(anArg > 1.0)
2416 anArg = 1.0;
2417 if(anArg < -1.0)
2418 anArg = -1.0;
2419
2420 const Standard_Real aDAngle = acos(anArg);
261b7d9e 2421 theURange[0].Add(aDAngle + theCoeffs.mFI1);
2422 theURange[0].Add(thePeriod - aDAngle + theCoeffs.mFI1);
79997052 2423 }
261b7d9e 2424 else if (-theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
79997052 2425 {
4dba155d 2426 // -(1+C)/B <= 1 and (1-C)/B >= -1 ==>
2427 //(U=[aDAngle1;aDAngle2]+aFI1) || (U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1),
2428 //where aDAngle1 = acos(-(myCoeffs.mC + 1) / myCoeffs.mB),
2429 // aDAngle2 = acos((1 - myCoeffs.mC) / myCoeffs.mB)
2430
261b7d9e 2431 Standard_Real anArg1 = -(theCoeffs.mC + 1) / theCoeffs.mB,
2432 anArg2 = (1 - theCoeffs.mC) / theCoeffs.mB;
79997052 2433 if(anArg1 > 1.0)
2434 anArg1 = 1.0;
2435 if(anArg1 < -1.0)
2436 anArg1 = -1.0;
2437
2438 if(anArg2 > 1.0)
2439 anArg2 = 1.0;
2440 if(anArg2 < -1.0)
2441 anArg2 = -1.0;
2442
2443 const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
261b7d9e 2444 theURange[0].Add(aDAngle1 + theCoeffs.mFI1);
2445 theURange[0].Add(aDAngle2 + theCoeffs.mFI1);
2446 theURange[1].Add(thePeriod - aDAngle2 + theCoeffs.mFI1);
2447 theURange[1].Add(thePeriod - aDAngle1 + theCoeffs.mFI1);
79997052 2448 }
2449 else
2450 {
2451 return Standard_False;
2452 }
2453 }
2454 else
2455 {
2456 return Standard_False;
2457 }
2458
2459 return Standard_True;
2460}
2461
ecc4f148 2462//=======================================================================
2463//function : CriticalPointsComputing
79997052 2464//purpose : theNbCritPointsMax contains true number of critical points.
2465// It must be initialized correctly before function calling
ecc4f148 2466//=======================================================================
d30895f5 2467static void CriticalPointsComputing(const ComputationMethods::stCoeffsValue& theCoeffs,
ecc4f148 2468 const Standard_Real theUSurf1f,
2469 const Standard_Real theUSurf1l,
2470 const Standard_Real theUSurf2f,
2471 const Standard_Real theUSurf2l,
2472 const Standard_Real thePeriod,
2473 const Standard_Real theTol2D,
7365fad6 2474 Standard_Integer& theNbCritPointsMax,
ecc4f148 2475 Standard_Real theU1crit[])
2476{
d30895f5 2477 //[0...1] - in these points parameter U1 goes through
2478 // the seam-edge of the first cylinder.
2479 //[2...3] - First and last U1 parameter.
2480 //[4...5] - in these points parameter U2 goes through
2481 // the seam-edge of the second cylinder.
2482 //[6...9] - in these points an intersection line goes through
2483 // U-boundaries of the second surface.
79997052 2484 //[10...11] - Boundary of monotonicity interval of U2(U1) function
2485 // (see CylCylMonotonicity() function)
7365fad6 2486
ecc4f148 2487 theU1crit[0] = 0.0;
2488 theU1crit[1] = thePeriod;
2489 theU1crit[2] = theUSurf1f;
2490 theU1crit[3] = theUSurf1l;
2491
2492 const Standard_Real aCOS = cos(theCoeffs.mFI2);
2493 const Standard_Real aBSB = Abs(theCoeffs.mB);
2494 if((theCoeffs.mC - aBSB <= aCOS) && (aCOS <= theCoeffs.mC + aBSB))
2495 {
2496 Standard_Real anArg = (aCOS - theCoeffs.mC) / theCoeffs.mB;
2497 if(anArg > 1.0)
2498 anArg = 1.0;
2499 if(anArg < -1.0)
2500 anArg = -1.0;
2501
2502 theU1crit[4] = -acos(anArg) + theCoeffs.mFI1;
7365fad6 2503 theU1crit[5] = acos(anArg) + theCoeffs.mFI1;
ecc4f148 2504 }
2505
2506 Standard_Real aSf = cos(theUSurf2f - theCoeffs.mFI2);
2507 Standard_Real aSl = cos(theUSurf2l - theCoeffs.mFI2);
2508 MinMax(aSf, aSl);
2509
7365fad6 2510 //In accorance with pure mathematic, theU1crit[6] and [8]
2511 //must be -Precision::Infinite() instead of used +Precision::Infinite()
2512 theU1crit[6] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
2513 -acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
2514 Precision::Infinite();
2515 theU1crit[7] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
2516 -acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
2517 Precision::Infinite();
2518 theU1crit[8] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
2519 acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
2520 Precision::Infinite();
2521 theU1crit[9] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
2522 acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
2523 Precision::Infinite();
2524
79997052 2525 theU1crit[10] = theCoeffs.mFI1;
2526 theU1crit[11] = M_PI+theCoeffs.mFI1;
2527
7365fad6 2528 //preparative treatment of array. This array must have faled to contain negative
2529 //infinity number
2530
2531 for(Standard_Integer i = 0; i < theNbCritPointsMax; i++)
2532 {
2533 if(Precision::IsInfinite(theU1crit[i]))
2534 {
2535 continue;
2536 }
2537
2538 theU1crit[i] = fmod(theU1crit[i], thePeriod);
2539 if(theU1crit[i] < 0.0)
2540 theU1crit[i] += thePeriod;
2541 }
2542
2543 //Here all not infinite elements of theU1crit are in [0, thePeriod) range
ecc4f148 2544
7365fad6 2545 do
ecc4f148 2546 {
7365fad6 2547 std::sort(theU1crit, theU1crit + theNbCritPointsMax);
2548 }
79997052 2549 while(ExcludeNearElements(theU1crit, theNbCritPointsMax,
2550 theUSurf1f, theUSurf1l, theTol2D));
7365fad6 2551
2552 //Here all not infinite elements in theU1crit are different and sorted.
2553
2554 while(theNbCritPointsMax > 0)
2555 {
2556 Standard_Real &anB = theU1crit[theNbCritPointsMax-1];
2557 if(Precision::IsInfinite(anB))
ecc4f148 2558 {
7365fad6 2559 theNbCritPointsMax--;
2560 continue;
ecc4f148 2561 }
7365fad6 2562
2563 //1st not infinte element is found
2564
2565 if(theNbCritPointsMax == 1)
2566 break;
2567
2568 //Here theNbCritPointsMax > 1
2569
2570 Standard_Real &anA = theU1crit[0];
2571
2572 //Compare 1st and last significant elements of theU1crit
2573 //They may still differs by period.
2574
2575 if (Abs(anB - anA - thePeriod) < theTol2D)
2576 {//E.g. anA == 2.0e-17, anB == (thePeriod-1.0e-18)
2577 anA = (anA + anB - thePeriod)/2.0;
2578 anB = Precision::Infinite();
2579 theNbCritPointsMax--;
2580 }
2581
2582 //Out of "while(theNbCritPointsMax > 0)" cycle.
2583 break;
ecc4f148 2584 }
7365fad6 2585
2586 //Attention! Here theU1crit may be unsorted.
ecc4f148 2587}
2588
2589//=======================================================================
d30895f5 2590//function : BoundaryEstimation
2591//purpose : Rough estimation of the parameter range.
2592//=======================================================================
2593void WorkWithBoundaries::BoundaryEstimation(const gp_Cylinder& theCy1,
2594 const gp_Cylinder& theCy2,
2595 Bnd_Range& theOutBoxS1,
2596 Bnd_Range& theOutBoxS2) const
2597{
2598 const gp_Dir &aD1 = theCy1.Axis().Direction(),
2599 &aD2 = theCy2.Axis().Direction();
2600 const Standard_Real aR1 = theCy1.Radius(),
2601 aR2 = theCy2.Radius();
2602
2603 //Let consider a parallelogram. Its edges are parallel to aD1 and aD2.
2604 //Its altitudes are equal to 2*aR1 and 2*aR2 (diameters of the cylinders).
2605 //In fact, this parallelogram is a projection of the cylinders to the plane
2606 //created by the intersected axes aD1 and aD2 (if the axes are skewed then
2607 //one axis can be translated by parallel shifting till intersection).
2608
2609 const Standard_Real aCosA = aD1.Dot(aD2);
95f8c608 2610 const Standard_Real aSqSinA = aD1.XYZ().CrossSquareMagnitude(aD2.XYZ());
d30895f5 2611
2612 //If sine is small then it can be compared with angle.
2613 if (aSqSinA < Precision::Angular()*Precision::Angular())
2614 return;
2615
2616 //Half of delta V. Delta V is a distance between
2617 //projections of two opposite parallelogram vertices
2618 //(joined by the maximal diagonal) to the cylinder axis.
2619 const Standard_Real aSinA = sqrt(aSqSinA);
c9c7286e 2620 const Standard_Real anAbsCosA = Abs(aCosA);
2621 const Standard_Real aHDV1 = (aR1 * anAbsCosA + aR2) / aSinA,
2622 aHDV2 = (aR2 * anAbsCosA + aR1) / aSinA;
d30895f5 2623
2624#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
2625 //The code in this block is created for test only.It is stupidly to create
2626 //OCCT-test for the method, which will be changed possibly never.
2627 std::cout << "Reference values: aHDV1 = " << aHDV1 << "; aHDV2 = " << aHDV2 << std::endl;
2628#endif
2629
2630 //V-parameters of intersection point of the axes (in case of skewed axes,
2631 //see comment above).
2632 Standard_Real aV01 = 0.0, aV02 = 0.0;
2633 ExtremaLineLine(theCy1.Axis(), theCy2.Axis(), aCosA, aSqSinA, aV01, aV02);
2634
2635 theOutBoxS1.Add(aV01 - aHDV1);
2636 theOutBoxS1.Add(aV01 + aHDV1);
2637
2638 theOutBoxS2.Add(aV02 - aHDV2);
2639 theOutBoxS2.Add(aV02 + aHDV2);
2640
2641 theOutBoxS1.Enlarge(Precision::Confusion());
2642 theOutBoxS2.Enlarge(Precision::Confusion());
2643
2644 Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0;
2645 myUVSurf1.Get(aU1, aV1, aU2, aV2);
2646 theOutBoxS1.Common(Bnd_Range(aV1, aV2));
2647
2648 myUVSurf2.Get(aU1, aV1, aU2, aV2);
2649 theOutBoxS2.Common(Bnd_Range(aV1, aV2));
2650}
2651
2652//=======================================================================
261b7d9e 2653//function : CyCyNoGeometric
ecc4f148 2654//purpose :
2655//=======================================================================
261b7d9e 2656static IntPatch_ImpImpIntersection::IntStatus
2657 CyCyNoGeometric(const gp_Cylinder &theCyl1,
2658 const gp_Cylinder &theCyl2,
2659 const WorkWithBoundaries &theBW,
2660 Bnd_Range theRange[],
2661 const Standard_Integer theNbOfRanges /*=2*/,
2662 Standard_Boolean& isTheEmpty,
2663 IntPatch_SequenceOfLine& theSlin,
2664 IntPatch_SequenceOfPoint& theSPnt)
ecc4f148 2665{
261b7d9e 2666 Standard_Real aUSurf1f = 0.0, aUSurf1l = 0.0,
2667 aUSurf2f = 0.0, aUSurf2l = 0.0,
2668 aVSurf1f = 0.0, aVSurf1l = 0.0,
2669 aVSurf2f = 0.0, aVSurf2l = 0.0;
4dba155d 2670
261b7d9e 2671 theBW.UVS1().Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
2672 theBW.UVS2().Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
d30895f5 2673
2674 Bnd_Range aRangeS1, aRangeS2;
261b7d9e 2675 theBW.BoundaryEstimation(theCyl1, theCyl2, aRangeS1, aRangeS2);
d30895f5 2676 if (aRangeS1.IsVoid() || aRangeS2.IsVoid())
e146ef9a 2677 return IntPatch_ImpImpIntersection::IntStatus_OK;
2678
2679 {
261b7d9e 2680 //Quotation of the message from issue #26894 (author MSV):
2681 //"We should return fail status from intersector if the result should be an
2682 //infinite curve of non-analytical type... I propose to define the limit for the
2683 //extent as the radius divided by 1e+2 and multiplied by 1e+7.
2684 //Thus, taking into account the number of valuable digits (15), we provide reliable
2685 //computations with an error not exceeding R/100."
e146ef9a 2686 const Standard_Real aF = 1.0e+5;
261b7d9e 2687 const Standard_Real aMaxV1Range = aF*theCyl1.Radius(), aMaxV2Range = aF*theCyl2.Radius();
e146ef9a 2688 if ((aRangeS1.Delta() > aMaxV1Range) || (aRangeS2.Delta() > aMaxV2Range))
2689 return IntPatch_ImpImpIntersection::IntStatus_InfiniteSectionCurve;
2690 }
64e8b010 2691 //
2692 Standard_Boolean isGoodIntersection = Standard_False;
2693 Standard_Real anOptdu = 0.;
003c363c 2694 for (;;)
64e8b010 2695 {
2696 //Checking parameters of cylinders in order to define "good intersection"
2697 //"Good intersection" means that axes of cylinders are almost perpendicular and
2698 // one radius is much smaller than the other and small cylinder is "inside" big one.
2699 const Standard_Real aToMuchCoeff = 3.;
2700 const Standard_Real aCritAngle = M_PI / 18.; // 10 degree
2701 Standard_Real anR1 = theCyl1.Radius();
2702 Standard_Real anR2 = theCyl2.Radius();
2703 Standard_Real anRmin = 0., anRmax = 0.;
2704 //Radius criterion
2705 if (anR1 > aToMuchCoeff * anR2)
2706 {
2707 anRmax = anR1; anRmin = anR2;
2708 }
2709 else if (anR2 > aToMuchCoeff * anR1)
2710 {
2711 anRmax = anR2; anRmin = anR1;
2712 }
2713 else
2714 {
2715 break;
2716 }
2717 //Angle criterion
2718 const gp_Ax1& anAx1 = theCyl1.Axis();
2719 const gp_Ax1& anAx2 = theCyl2.Axis();
2720 if (!anAx1.IsNormal(anAx2, aCritAngle))
2721 {
2722 break;
2723 }
2724 //Placement criterion
2725 gp_Lin anL1(anAx1), anL2(anAx2);
2726 Standard_Real aDist = anL1.Distance(anL2);
2727 if (aDist > anRmax / 2.)
2728 {
2729 break;
2730 }
ecc4f148 2731
64e8b010 2732 isGoodIntersection = Standard_True;
2733 //Estimation of "optimal" du
2734 //Relative deflection, absolut deflection is Rmin*aDeflection
2735 Standard_Real aDeflection = 0.001;
2736 Standard_Integer aNbP = 3;
2737 if (anRmin * aDeflection > 1.e-3)
2738 {
2739 Standard_Real anAngle = 1.0e0 - aDeflection;
2740 anAngle = 2.0e0 * ACos(anAngle);
2741 aNbP = (Standard_Integer)(2. * M_PI / anAngle) + 1;
2742 }
2743 anOptdu = 2. * M_PI_2 / (Standard_Real)(aNbP - 1);
003c363c 2744 break;
2745 }
64e8b010 2746//
261b7d9e 2747 const ComputationMethods::stCoeffsValue &anEquationCoeffs = theBW.SICoeffs();
2748 const IntSurf_Quadric& aQuad1 = theBW.GetQSurface(1);
2749 const IntSurf_Quadric& aQuad2 = theBW.GetQSurface(2);
2750 const Standard_Boolean isReversed = theBW.IsReversed();
2751 const Standard_Real aTol2D = theBW.Get2dTolerance();
2752 const Standard_Real aTol3D = theBW.Get3dTolerance();
2753 const Standard_Real aPeriod = 2.0*M_PI;
64e8b010 2754 Standard_Integer aNbMaxPoints = 1000;
2755 Standard_Integer aNbMinPoints = 200;
2756 Standard_Real du;
2757 if (isGoodIntersection)
2758 {
2759 du = anOptdu;
2760 aNbMaxPoints = 200;
2761 aNbMinPoints = 50;
2762 }
2763 else
2764 {
2765 du = 2. * M_PI / aNbMaxPoints;
2766 }
2767 Standard_Integer aNbPts = Min(RealToInt((aUSurf1l - aUSurf1f) / du) + 1,
2768 RealToInt(20.0*theCyl1.Radius()));
2769 const Standard_Integer aNbPoints = Min(Max(aNbMinPoints, aNbPts), aNbMaxPoints);
2770 const Standard_Real aStepMin = Max(aTol2D, Precision::PConfusion()),
2771 aStepMax = (aUSurf1l - aUSurf1f > M_PI / 100.0) ?
2772 (aUSurf1l - aUSurf1f) / IntToReal(aNbPoints) : aUSurf1l - aUSurf1f;
2773
2774
4dba155d 2775 //The main idea of the algorithm is to change U1-parameter
261b7d9e 2776 //(U-parameter of theCyl1) from aU1f to aU1l with some step
4dba155d 2777 //(step is adaptive) and to obtain set of intersection points.
2778
261b7d9e 2779 for (Standard_Integer i = 0; i < theNbOfRanges; i++)
ecc4f148 2780 {
261b7d9e 2781 if (theRange[i].IsVoid())
ecc4f148 2782 continue;
2783
261b7d9e 2784 InscribeInterval(aUSurf1f, aUSurf1l, theRange[i], aTol2D, aPeriod);
ecc4f148 2785 }
2786
261b7d9e 2787 if (theRange[0].Union(theRange[1]))
ecc4f148 2788 {
261b7d9e 2789 // Works only if (theNbOfRanges == 2).
2790 theRange[1].SetVoid();
ecc4f148 2791 }
2792
4dba155d 2793 //Critical points are the value of U1-parameter in the points
2794 //where WL must be decomposed.
2795
2796 //When U1 goes through critical points its value is set up to this
2797 //parameter forcefully and the intersection point is added in the line.
2798 //After that, the WL is broken (next U1 value will be correspond to the new WL).
2799
2800 //See CriticalPointsComputing(...) function to get detail information about this array.
79997052 2801 const Standard_Integer aNbCritPointsMax = 12;
261b7d9e 2802 Standard_Real anU1crit[aNbCritPointsMax] = { Precision::Infinite(),
2803 Precision::Infinite(),
2804 Precision::Infinite(),
2805 Precision::Infinite(),
2806 Precision::Infinite(),
2807 Precision::Infinite(),
2808 Precision::Infinite(),
2809 Precision::Infinite(),
2810 Precision::Infinite(),
2811 Precision::Infinite(),
2812 Precision::Infinite(),
2813 Precision::Infinite() };
ecc4f148 2814
4dba155d 2815 //This list of critical points is not full because it does not contain any points
2816 //which intersection line goes through V-bounds of cylinders in.
2817 //They are computed by numerical methods on - line (during algorithm working).
2818 //The moment is caught, when intersection line goes through V-bounds of any cylinder.
2819
7365fad6 2820 Standard_Integer aNbCritPoints = aNbCritPointsMax;
ecc4f148 2821 CriticalPointsComputing(anEquationCoeffs, aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
261b7d9e 2822 aPeriod, aTol2D, aNbCritPoints, anU1crit);
ecc4f148 2823
2824 //Getting Walking-line
2825
b2af2f56 2826 enum WLFStatus
2827 {
4dba155d 2828 // No points have been added in WL
b2af2f56 2829 WLFStatus_Absent = 0,
4dba155d 2830 // WL contains at least one point
261b7d9e 2831 WLFStatus_Exist = 1,
4dba155d 2832 // WL has been finished in some critical point
2833 // We should start new line
b2af2f56 2834 WLFStatus_Broken = 2
2835 };
2836
261b7d9e 2837 const Standard_Integer aNbWLines = 2;
2838 for (Standard_Integer aCurInterval = 0; aCurInterval < theNbOfRanges; aCurInterval++)
ecc4f148 2839 {
4dba155d 2840 //Process every continuous region
7c32c7c4 2841 Standard_Boolean isAddedIntoWL[aNbWLines];
261b7d9e 2842 for (Standard_Integer i = 0; i < aNbWLines; i++)
7c32c7c4 2843 isAddedIntoWL[i] = Standard_False;
ecc4f148 2844
261b7d9e 2845 Standard_Real anUf = 1.0, anUl = 0.0;
2846 if (!theRange[aCurInterval].GetBounds(anUf, anUl))
2847 continue;
2848
2849 const Standard_Boolean isDeltaPeriod = IsEqual(anUl - anUf, aPeriod);
ecc4f148 2850
2851 //Inscribe and sort critical points
261b7d9e 2852 for (Standard_Integer i = 0; i < aNbCritPoints; i++)
7365fad6 2853 {
261b7d9e 2854 InscribePoint(anUf, anUl, anU1crit[i], 0.0, aPeriod, Standard_False);
7365fad6 2855 }
2856
2857 std::sort(anU1crit, anU1crit + aNbCritPoints);
ecc4f148 2858
261b7d9e 2859 while (anUf < anUl)
ecc4f148 2860 {
4dba155d 2861 //Change value of U-parameter on the 1st surface from anUf to anUl
2862 //(anUf will be modified in the cycle body).
2863 //Step is computed adaptively (see comments below).
2864
7c32c7c4 2865 Standard_Real aU2[aNbWLines], aV1[aNbWLines], aV2[aNbWLines];
b2af2f56 2866 WLFStatus aWLFindStatus[aNbWLines];
7c32c7c4 2867 Standard_Real aV1Prev[aNbWLines], aV2Prev[aNbWLines];
b5ef9d91 2868 Standard_Real anUexpect[aNbWLines];
7c32c7c4 2869 Standard_Boolean isAddingWLEnabled[aNbWLines];
2870
2871 Handle(IntSurf_LineOn2S) aL2S[aNbWLines];
2872 Handle(IntPatch_WLine) aWLine[aNbWLines];
261b7d9e 2873 for (Standard_Integer i = 0; i < aNbWLines; i++)
7c32c7c4 2874 {
2875 aL2S[i] = new IntSurf_LineOn2S();
2876 aWLine[i] = new IntPatch_WLine(aL2S[i], Standard_False);
98974dcc 2877 aWLine[i]->SetCreatingWayInfo(IntPatch_WLine::IntPatch_WLImpImp);
b2af2f56 2878 aWLFindStatus[i] = WLFStatus_Absent;
7c32c7c4 2879 isAddingWLEnabled[i] = Standard_True;
2880 aU2[i] = aV1[i] = aV2[i] = 0.0;
2881 aV1Prev[i] = aV2Prev[i] = 0.0;
b5ef9d91 2882 anUexpect[i] = anUf;
7c32c7c4 2883 }
261b7d9e 2884
2885 Standard_Real aCriticalDelta[aNbCritPointsMax] = { 0 };
2886 for (Standard_Integer aCritPID = 0; aCritPID < aNbCritPoints; aCritPID++)
4dba155d 2887 { //We are not interested in elements of aCriticalDelta array
7365fad6 2888 //if their index is greater than or equal to aNbCritPoints
2889
2890 aCriticalDelta[aCritPID] = anUf - anU1crit[aCritPID];
2891 }
ecc4f148 2892
261b7d9e 2893 Standard_Real anU1 = anUf, aMinCriticalParam = anUf;
ecc4f148 2894 Standard_Boolean isFirst = Standard_True;
261b7d9e 2895
2896 while (anU1 <= anUl)
ecc4f148 2897 {
4dba155d 2898 //Change value of U-parameter on the 1st surface from anUf to anUl
2899 //(anUf will be modified in the cycle body). However, this cycle
2900 //can be broken if WL goes though some critical point.
2901 //Step is computed adaptively (see comments below).
2902
261b7d9e 2903 for (Standard_Integer i = 0; i < aNbCritPoints; i++)
ecc4f148 2904 {
261b7d9e 2905 if ((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0)
ecc4f148 2906 {
4dba155d 2907 //WL has gone through i-th critical point
ecc4f148 2908 anU1 = anU1crit[i];
7c32c7c4 2909
261b7d9e 2910 for (Standard_Integer j = 0; j < aNbWLines; j++)
b5ef9d91 2911 {
51740958 2912 aWLFindStatus[j] = WLFStatus_Broken;
2913 anUexpect[j] = anU1;
b5ef9d91 2914 }
7c32c7c4 2915
ecc4f148 2916 break;
2917 }
2918 }
2919
261b7d9e 2920 if (IsEqual(anU1, anUl))
b5ef9d91 2921 {
261b7d9e 2922 for (Standard_Integer i = 0; i < aNbWLines; i++)
b5ef9d91 2923 {
2924 aWLFindStatus[i] = WLFStatus_Broken;
2925 anUexpect[i] = anU1;
2926
261b7d9e 2927 if (isDeltaPeriod)
b5ef9d91 2928 {
7365fad6 2929 //if isAddedIntoWL[i] == TRUE WLine contains only one point
b5ef9d91 2930 //(which was end point of previous WLine). If we will
2931 //add point found on the current step WLine will contain only
2932 //two points. At that both these points will be equal to the
2933 //points found earlier. Therefore, new WLine will repeat
2934 //already existing WLine. Consequently, it is necessary
2935 //to forbid building new line in this case.
2936
2937 isAddingWLEnabled[i] = (!isAddedIntoWL[i]);
2938 }
2939 else
2940 {
261b7d9e 2941 isAddingWLEnabled[i] = ((aTol2D >= (anUexpect[i] - anU1)) ||
b5ef9d91 2942 (aWLFindStatus[i] == WLFStatus_Absent));
2943 }
2944 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2945 }
2946 else
2947 {
261b7d9e 2948 for (Standard_Integer i = 0; i < aNbWLines; i++)
b5ef9d91 2949 {
261b7d9e 2950 isAddingWLEnabled[i] = ((aTol2D >= (anUexpect[i] - anU1)) ||
b5ef9d91 2951 (aWLFindStatus[i] == WLFStatus_Absent));
2952 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2953 }
ecc4f148 2954
261b7d9e 2955 for (Standard_Integer i = 0; i < aNbWLines; i++)
7c32c7c4 2956 {
261b7d9e 2957 const Standard_Integer aNbPntsWL = aWLine[i].IsNull() ? 0 :
2958 aWLine[i]->Curve()->NbPoints();
2959
2960 if ((aWLFindStatus[i] == WLFStatus_Broken) ||
2961 (aWLFindStatus[i] == WLFStatus_Absent))
7365fad6 2962 {//Begin and end of WLine must be on boundary point
2963 //or on seam-edge strictly (if it is possible).
b5ef9d91 2964
261b7d9e 2965 Standard_Real aTol = aTol2D;
d30895f5 2966 ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs,
2967 aU2[i], &aTol);
261b7d9e 2968 InscribePoint(aUSurf2f, aUSurf2l, aU2[i], aTol2D, aPeriod, Standard_False);
02effd35 2969
261b7d9e 2970 aTol = Max(aTol, aTol2D);
02effd35 2971
261b7d9e 2972 if (Abs(aU2[i]) <= aTol)
7365fad6 2973 aU2[i] = 0.0;
261b7d9e 2974 else if (Abs(aU2[i] - aPeriod) <= aTol)
7365fad6 2975 aU2[i] = aPeriod;
261b7d9e 2976 else if (Abs(aU2[i] - aUSurf2f) <= aTol)
7365fad6 2977 aU2[i] = aUSurf2f;
261b7d9e 2978 else if (Abs(aU2[i] - aUSurf2l) <= aTol)
7365fad6 2979 aU2[i] = aUSurf2l;
2980 }
2981 else
2982 {
d30895f5 2983 ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]);
261b7d9e 2984 InscribePoint(aUSurf2f, aUSurf2l, aU2[i], aTol2D, aPeriod, Standard_False);
7365fad6 2985 }
261b7d9e 2986
2987 if (aNbPntsWL == 0)
7c32c7c4 2988 {//the line has not contained any points yet
261b7d9e 2989 if (((aUSurf2f + aPeriod - aUSurf2l) <= 2.0*aTol2D) &&
2990 ((Abs(aU2[i] - aUSurf2f) < aTol2D) ||
2991 (Abs(aU2[i] - aUSurf2l) < aTol2D)))
02effd35 2992 {
e002f1ce 2993 //In this case aU2[i] can have two values: current aU2[i] or
b5ef9d91 2994 //aU2[i]+aPeriod (aU2[i]-aPeriod). It is necessary to choose
2995 //correct value.
e002f1ce 2996
b5ef9d91 2997 Standard_Boolean isIncreasing = Standard_True;
d30895f5 2998 ComputationMethods::CylCylMonotonicity(anU1+aStepMin, i, anEquationCoeffs,
2999 aPeriod, isIncreasing);
7c32c7c4 3000
e002f1ce 3001 //If U2(U1) is increasing and U2 is considered to be equal aUSurf2l
3002 //then after the next step (when U1 will be increased) U2 will be
3003 //increased too. And we will go out of surface boundary.
3004 //Therefore, If U2(U1) is increasing then U2 must be equal aUSurf2f.
3005 //Analogically, if U2(U1) is decreasing.
261b7d9e 3006 if (isIncreasing)
e002f1ce 3007 {
3008 aU2[i] = aUSurf2f;
3009 }
3010 else
3011 {
3012 aU2[i] = aUSurf2l;
7c32c7c4 3013 }
02effd35 3014 }
3015 }
e002f1ce 3016 else
3017 {//aNbPntsWL > 0
261b7d9e 3018 if (((aUSurf2l - aUSurf2f) >= aPeriod) &&
3019 ((Abs(aU2[i] - aUSurf2f) < aTol2D) ||
3020 (Abs(aU2[i] - aUSurf2l) < aTol2D)))
7c32c7c4 3021 {//end of the line
3022 Standard_Real aU2prev = 0.0, aV2prev = 0.0;
261b7d9e 3023 if (isReversed)
7c32c7c4 3024 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS1(aU2prev, aV2prev);
02effd35 3025 else
7c32c7c4 3026 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS2(aU2prev, aV2prev);
3027
261b7d9e 3028 if (2.0*Abs(aU2prev - aU2[i]) > aPeriod)
7c32c7c4 3029 {
261b7d9e 3030 if (aU2prev > aU2[i])
7c32c7c4 3031 aU2[i] += aPeriod;
3032 else
3033 aU2[i] -= aPeriod;
3034 }
02effd35 3035 }
3036 }
02effd35 3037
d30895f5 3038 ComputationMethods::CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs,
3039 aV1[i], aV2[i]);
7c32c7c4 3040
261b7d9e 3041 if (isFirst)
7c32c7c4 3042 {
3043 aV1Prev[i] = aV1[i];
3044 aV2Prev[i] = aV2[i];
02effd35 3045 }
7c32c7c4 3046 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
ecc4f148 3047
7c32c7c4 3048 isFirst = Standard_False;
ecc4f148 3049
7c32c7c4 3050 //Looking for points into WLine
3051 Standard_Boolean isBroken = Standard_False;
261b7d9e 3052 for (Standard_Integer i = 0; i < aNbWLines; i++)
ecc4f148 3053 {
261b7d9e 3054 if (!isAddingWLEnabled[i])
ecc4f148 3055 {
b5ef9d91 3056 Standard_Boolean isBoundIntersect = Standard_False;
261b7d9e 3057 if ((Abs(aV1[i] - aVSurf1f) <= aTol2D) ||
3058 ((aV1[i] - aVSurf1f)*(aV1Prev[i] - aVSurf1f) < 0.0))
b5ef9d91 3059 {
3060 isBoundIntersect = Standard_True;
3061 }
261b7d9e 3062 else if ((Abs(aV1[i] - aVSurf1l) <= aTol2D) ||
3063 ((aV1[i] - aVSurf1l)*(aV1Prev[i] - aVSurf1l) < 0.0))
b5ef9d91 3064 {
3065 isBoundIntersect = Standard_True;
3066 }
261b7d9e 3067 else if ((Abs(aV2[i] - aVSurf2f) <= aTol2D) ||
3068 ((aV2[i] - aVSurf2f)*(aV2Prev[i] - aVSurf2f) < 0.0))
b5ef9d91 3069 {
3070 isBoundIntersect = Standard_True;
3071 }
261b7d9e 3072 else if ((Abs(aV2[i] - aVSurf2l) <= aTol2D) ||
3073 ((aV2[i] - aVSurf2l)*(aV2Prev[i] - aVSurf2l) < 0.0))
b5ef9d91 3074 {
3075 isBoundIntersect = Standard_True;
3076 }
3077
261b7d9e 3078 if (aWLFindStatus[i] == WLFStatus_Broken)
7c32c7c4 3079 isBroken = Standard_True;
02effd35 3080
261b7d9e 3081 if (!isBoundIntersect)
b5ef9d91 3082 {
3083 continue;
3084 }
3085 else
3086 {
3087 anUexpect[i] = anU1;
3088 }
ecc4f148 3089 }
ecc4f148 3090
4dba155d 3091 // True if the current point already in the domain
261b7d9e 3092 const Standard_Boolean isInscribe =
3093 ((aUSurf2f - aU2[i]) <= aTol2D) && ((aU2[i] - aUSurf2l) <= aTol2D) &&
3094 ((aVSurf1f - aV1[i]) <= aTol2D) && ((aV1[i] - aVSurf1l) <= aTol2D) &&
3095 ((aVSurf2f - aV2[i]) <= aTol2D) && ((aV2[i] - aVSurf2l) <= aTol2D);
b5ef9d91 3096
3097 //isVIntersect == TRUE if intersection line intersects two (!)
3098 //V-bounds of cylinder (1st or 2nd - no matter)
3099 const Standard_Boolean isVIntersect =
261b7d9e 3100 (((aVSurf1f - aV1[i])*(aVSurf1f - aV1Prev[i]) < RealSmall()) &&
3101 ((aVSurf1l - aV1[i])*(aVSurf1l - aV1Prev[i]) < RealSmall())) ||
3102 (((aVSurf2f - aV2[i])*(aVSurf2f - aV2Prev[i]) < RealSmall()) &&
3103 ((aVSurf2l - aV2[i])*(aVSurf2l - aV2Prev[i]) < RealSmall()));
b5ef9d91 3104
3105 //isFound1 == TRUE if intersection line intersects V-bounds
3106 // (First or Last - no matter) of the 1st cylynder
3107 //isFound2 == TRUE if intersection line intersects V-bounds
3108 // (First or Last - no matter) of the 2nd cylynder
3109 Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
3110 Standard_Boolean isForce = Standard_False;
3111
5b111128 3112 if (aWLFindStatus[i] == WLFStatus_Absent)
ecc4f148 3113 {
261b7d9e 3114 if (((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1 - aUSurf1l) < aTol2D))
02effd35 3115 {
b5ef9d91 3116 isForce = Standard_True;
3117 }
3118 }
e8feb725 3119
261b7d9e 3120 theBW.AddBoundaryPoint(aWLine[i], anU1, aMinCriticalParam, aU2[i],
3121 aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], i, isForce,
3122 isFound1, isFound2);
02effd35 3123
b5ef9d91 3124 const Standard_Boolean isPrevVBound = !isVIntersect &&
261b7d9e 3125 ((Abs(aV1Prev[i] - aVSurf1f) <= aTol2D) ||
3126 (Abs(aV1Prev[i] - aVSurf1l) <= aTol2D) ||
3127 (Abs(aV2Prev[i] - aVSurf2f) <= aTol2D) ||
3128 (Abs(aV2Prev[i] - aVSurf2l) <= aTol2D));
b5ef9d91 3129
7365fad6 3130 aV1Prev[i] = aV1[i];
3131 aV2Prev[i] = aV2[i];
3132
261b7d9e 3133 if ((aWLFindStatus[i] == WLFStatus_Exist) && (isFound1 || isFound2) && !isPrevVBound)
b5ef9d91 3134 {
3135 aWLFindStatus[i] = WLFStatus_Broken; //start a new line
3136 }
261b7d9e 3137 else if (isInscribe)
b5ef9d91 3138 {
261b7d9e 3139 if ((aWLFindStatus[i] == WLFStatus_Absent) && (isFound1 || isFound2))
b5ef9d91 3140 {
3141 aWLFindStatus[i] = WLFStatus_Exist;
ecc4f148 3142 }
ecc4f148 3143
261b7d9e 3144 if ((aWLFindStatus[i] != WLFStatus_Broken) ||
3145 (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl))
ecc4f148 3146 {
261b7d9e 3147 if (aWLine[i]->NbPnts() > 0)
5b9e1842 3148 {
3149 Standard_Real aU2p = 0.0, aV2p = 0.0;
261b7d9e 3150 if (isReversed)
5b9e1842 3151 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
3152 else
3153 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
3154
3155 const Standard_Real aDelta = aU2[i] - aU2p;
3156
261b7d9e 3157 if (2.0 * Abs(aDelta) > aPeriod)
5b9e1842 3158 {
261b7d9e 3159 if (aDelta > 0.0)
5b9e1842 3160 {
3161 aU2[i] -= aPeriod;
3162 }
3163 else
3164 {
3165 aU2[i] += aPeriod;
3166 }
3167 }
3168 }
3169
261b7d9e 3170 if(AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed, Standard_True,
7c32c7c4 3171 gp_Pnt2d(anU1, aV1[i]), gp_Pnt2d(aU2[i], aV2[i]),
b5ef9d91 3172 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
3173 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod,
261b7d9e 3174 aWLine[i]->Curve(), i, aTol3D, aTol2D, isForce))
ecc4f148 3175 {
261b7d9e 3176 if (aWLFindStatus[i] == WLFStatus_Absent)
e8feb725 3177 {
b2af2f56 3178 aWLFindStatus[i] = WLFStatus_Exist;
e8feb725 3179 }
ecc4f148 3180 }
261b7d9e 3181 else if (!isFound1 && !isFound2)
5b9e1842 3182 {//We do not add any point while doing this iteration
261b7d9e 3183 if (aWLFindStatus[i] == WLFStatus_Exist)
5b9e1842 3184 {
3185 aWLFindStatus[i] = WLFStatus_Broken;
261b7d9e 3186 }
5b9e1842 3187 }
3188 }
3189 }
3190 else
3191 {//We do not add any point while doing this iteration
261b7d9e 3192 if (aWLFindStatus[i] == WLFStatus_Exist)
5b9e1842 3193 {
3194 aWLFindStatus[i] = WLFStatus_Broken;
ecc4f148 3195 }
3196 }
261b7d9e 3197
3198 if (aWLFindStatus[i] == WLFStatus_Broken)
7c32c7c4 3199 isBroken = Standard_True;
3200 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
ecc4f148 3201
261b7d9e 3202 if (isBroken)
ecc4f148 3203 {//current lines are filled. Go to the next lines
3204 anUf = anU1;
b5ef9d91 3205
5b9e1842 3206 Standard_Boolean isAdded = Standard_True;
3207
261b7d9e 3208 for (Standard_Integer i = 0; i < aNbWLines; i++)
b5ef9d91 3209 {
261b7d9e 3210 if (isAddingWLEnabled[i])
5b9e1842 3211 {
b5ef9d91 3212 continue;
5b9e1842 3213 }
b5ef9d91 3214
5b9e1842 3215 isAdded = Standard_False;
b5ef9d91 3216
5b9e1842 3217 Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
b5ef9d91 3218
261b7d9e 3219 theBW.AddBoundaryPoint(aWLine[i], anU1, aMinCriticalParam, aU2[i],
3220 aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], i,
3221 Standard_False, isFound1, isFound2);
b5ef9d91 3222
261b7d9e 3223 if (isFound1 || isFound2)
5b9e1842 3224 {
3225 isAdded = Standard_True;
3226 }
b5ef9d91 3227
261b7d9e 3228 if (aWLine[i]->NbPnts() > 0)
b5ef9d91 3229 {
5b9e1842 3230 Standard_Real aU2p = 0.0, aV2p = 0.0;
261b7d9e 3231 if (isReversed)
5b9e1842 3232 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
3233 else
3234 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
3235
3236 const Standard_Real aDelta = aU2[i] - aU2p;
3237
261b7d9e 3238 if (2 * Abs(aDelta) > aPeriod)
5b9e1842 3239 {
261b7d9e 3240 if (aDelta > 0.0)
5b9e1842 3241 {
3242 aU2[i] -= aPeriod;
3243 }
3244 else
3245 {
3246 aU2[i] += aPeriod;
3247 }
3248 }
3249 }
b5ef9d91 3250
261b7d9e 3251 if(AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
5b9e1842 3252 Standard_True, gp_Pnt2d(anU1, aV1[i]),
3253 gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
3254 aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
3255 aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
261b7d9e 3256 i, aTol3D, aTol2D, Standard_False))
5b9e1842 3257 {
3258 isAdded = Standard_True;
b5ef9d91 3259 }
5b9e1842 3260 }
b5ef9d91 3261
261b7d9e 3262 if (!isAdded)
5b9e1842 3263 {
4dba155d 3264 //Before breaking WL, we must complete it correctly
3265 //(e.g. to prolong to the surface boundary).
3266 //Therefore, we take the point last added in some WL
3267 //(have maximal U1-parameter) and try to add it in
3268 //the current WL.
5b9e1842 3269 Standard_Real anUmaxAdded = RealFirst();
261b7d9e 3270
5b9e1842 3271 {
7365fad6 3272 Standard_Boolean isChanged = Standard_False;
261b7d9e 3273 for (Standard_Integer i = 0; i < aNbWLines; i++)
7365fad6 3274 {
e9a7ec7a 3275 if ((aWLFindStatus[i] == WLFStatus_Absent) || (aWLine[i]->NbPnts() == 0))
7365fad6 3276 continue;
b5ef9d91 3277
7365fad6 3278 Standard_Real aU1c = 0.0, aV1c = 0.0;
261b7d9e 3279 if (isReversed)
7365fad6 3280 aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS2(aU1c, aV1c);
3281 else
3282 aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS1(aU1c, aV1c);
3283
3284 anUmaxAdded = Max(anUmaxAdded, aU1c);
3285 isChanged = Standard_True;
3286 }
3287
261b7d9e 3288 if (!isChanged)
7365fad6 3289 { //If anUmaxAdded were not changed in previous cycle then
3290 //we would break existing WLines.
3291 break;
3292 }
b5ef9d91 3293 }
3294
261b7d9e 3295 for (Standard_Integer i = 0; i < aNbWLines; i++)
5b9e1842 3296 {
261b7d9e 3297 if (isAddingWLEnabled[i])
5b9e1842 3298 {
3299 continue;
3300 }
b5ef9d91 3301
d30895f5 3302 ComputationMethods::CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs,
261b7d9e 3303 aU2[i], aV1[i], aV2[i]);
3304
3305 AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
3306 Standard_True, gp_Pnt2d(anUmaxAdded, aV1[i]),
3307 gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
3308 aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
3309 aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
3310 i, aTol3D, aTol2D, Standard_False);
5b9e1842 3311 }
b5ef9d91 3312 }
3313
ecc4f148 3314 break;
3315 }
3316
7c32c7c4 3317 //Step computing
ecc4f148 3318
7c32c7c4 3319 {
4dba155d 3320 //Step of aU1-parameter is computed adaptively. The algorithm
3321 //aims to provide given aDeltaV1 and aDeltaV2 values (if it is
3322 //possible because the intersection line can go along V-isoline)
3323 //in every iteration. It allows avoiding "flying" intersection
3324 //points too far each from other (see issue #24915).
3325
261b7d9e 3326 const Standard_Real aDeltaV1 = aRangeS1.Delta() / IntToReal(aNbPoints),
3327 aDeltaV2 = aRangeS2.Delta() / IntToReal(aNbPoints);
3328
b5ef9d91 3329 math_Matrix aMatr(1, 3, 1, 5);
ecc4f148 3330
b5ef9d91 3331 Standard_Real aMinUexp = RealLast();
261b7d9e 3332
3333 for (Standard_Integer i = 0; i < aNbWLines; i++)
b5ef9d91 3334 {
261b7d9e 3335 if (aTol2D < (anUexpect[i] - anU1))
b5ef9d91 3336 {
3337 continue;
3338 }
ecc4f148 3339
261b7d9e 3340 if (aWLFindStatus[i] == WLFStatus_Absent)
b5ef9d91 3341 {
3342 anUexpect[i] += aStepMax;
3343 aMinUexp = Min(aMinUexp, anUexpect[i]);
3344 continue;
3345 }
64e8b010 3346 //
3347 if (isGoodIntersection)
3348 {
3349 //Use constant step
3350 anUexpect[i] += aStepMax;
3351 aMinUexp = Min(aMinUexp, anUexpect[i]);
3352
3353 continue;
3354 }
3355 //
7c32c7c4 3356
b5ef9d91 3357 Standard_Real aStepTmp = aStepMax;
ecc4f148 3358
b5ef9d91 3359 const Standard_Real aSinU1 = sin(anU1),
3360 aCosU1 = cos(anU1),
3361 aSinU2 = sin(aU2[i]),
3362 aCosU2 = cos(aU2[i]);
ecc4f148 3363
b5ef9d91 3364 aMatr.SetCol(1, anEquationCoeffs.mVecC1);
3365 aMatr.SetCol(2, anEquationCoeffs.mVecC2);
3366 aMatr.SetCol(3, anEquationCoeffs.mVecA1*aSinU1 - anEquationCoeffs.mVecB1*aCosU1);
3367 aMatr.SetCol(4, anEquationCoeffs.mVecA2*aSinU2 - anEquationCoeffs.mVecB2*aCosU2);
3368 aMatr.SetCol(5, anEquationCoeffs.mVecA1*aCosU1 + anEquationCoeffs.mVecB1*aSinU1 +
3369 anEquationCoeffs.mVecA2*aCosU2 + anEquationCoeffs.mVecB2*aSinU2 +
3370 anEquationCoeffs.mVecD);
ecc4f148 3371
4dba155d 3372 //The main idea is in solving of linearized system (2)
3373 //(see description to ComputationMethods class) in order to find new U1-value
3374 //to provide new value V1 or V2, which differs from current one by aDeltaV1 or
3375 //aDeltaV2 respectively.
3376
3377 //While linearizing, following Taylor formulas are used:
3378 // cos(x0+dx) = cos(x0) - sin(x0)*dx
3379 // sin(x0+dx) = sin(x0) + cos(x0)*dx
3380
3381 //Consequently cos(U1), cos(U2), sin(U1) and sin(U2) in the system (2)
3382 //must be substituted by corresponding values.
3383
3384 //ATTENTION!!!
3385 //The solution is approximate. More over, all requirements to the
3386 //linearization must be satisfied in order to obtain quality result.
3387
261b7d9e 3388 if (!StepComputing(aMatr, aV1[i], aV2[i], aDeltaV1, aDeltaV2, aStepTmp))
b5ef9d91 3389 {
3390 //To avoid cycling-up
3391 anUexpect[i] += aStepMax;
3392 aMinUexp = Min(aMinUexp, anUexpect[i]);
ecc4f148 3393
b5ef9d91 3394 continue;
3395 }
3396
261b7d9e 3397 if (aStepTmp < aStepMin)
b5ef9d91 3398 aStepTmp = aStepMin;
261b7d9e 3399
3400 if (aStepTmp > aStepMax)
b5ef9d91 3401 aStepTmp = aStepMax;
3402
3403 anUexpect[i] = anU1 + aStepTmp;
3404 aMinUexp = Min(aMinUexp, anUexpect[i]);
3405 }
ecc4f148 3406
b5ef9d91 3407 anU1 = aMinUexp;
3408 }
ecc4f148 3409
261b7d9e 3410 if (Precision::PConfusion() >= (anUl - anU1))
ecc4f148 3411 anU1 = anUl;
3412
3413 anUf = anU1;
3414
261b7d9e 3415 for (Standard_Integer i = 0; i < aNbWLines; i++)
baf72cd2 3416 {
261b7d9e 3417 if (aWLine[i]->NbPnts() != 1)
7c32c7c4 3418 isAddedIntoWL[i] = Standard_False;
b5ef9d91 3419
261b7d9e 3420 if (anU1 == anUl)
b5ef9d91 3421 {//strictly equal. Tolerance is considered above.
3422 anUexpect[i] = anUl;
3423 }
baf72cd2 3424 }
ecc4f148 3425 }
ecc4f148 3426
261b7d9e 3427 for (Standard_Integer i = 0; i < aNbWLines; i++)
ecc4f148 3428 {
261b7d9e 3429 if ((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
baf72cd2 3430 {
7c32c7c4 3431 isTheEmpty = Standard_False;
3432 Standard_Real u1, v1, u2, v2;
3433 aWLine[i]->Point(1).Parameters(u1, v1, u2, v2);
3434 IntPatch_Point aP;
3435 aP.SetParameter(u1);
3436 aP.SetParameters(u1, v1, u2, v2);
261b7d9e 3437 aP.SetTolerance(aTol3D);
7c32c7c4 3438 aP.SetValue(aWLine[i]->Point(1).Value());
3439
dcd768a4 3440 //Check whether the added point exists.
3441 //It is enough to check the last point.
3442 if (theSPnt.IsEmpty() ||
3443 !theSPnt.Last().PntOn2S().IsSame(aP.PntOn2S(), Precision::Confusion()))
3444 {
3445 theSPnt.Append(aP);
3446 }
baf72cd2 3447 }
261b7d9e 3448 else if (aWLine[i]->NbPnts() > 1)
baf72cd2 3449 {
7c32c7c4 3450 Standard_Boolean isGood = Standard_True;
ecc4f148 3451
261b7d9e 3452 if (aWLine[i]->NbPnts() == 2)
7c32c7c4 3453 {
3454 const IntSurf_PntOn2S& aPf = aWLine[i]->Point(1);
3455 const IntSurf_PntOn2S& aPl = aWLine[i]->Point(2);
261b7d9e 3456
3457 if (aPf.IsSame(aPl, Precision::Confusion()))
7c32c7c4 3458 isGood = Standard_False;
3459 }
7eb3580b 3460 else if (aWLine[i]->NbPnts() > 2)
3461 {
3462 // Sometimes points of the WLine are distributed
3463 // linearly and uniformly. However, such position
3464 // of the points does not always describe the real intersection
3465 // curve. I.e. real tangents at the ends of the intersection
3466 // curve can significantly deviate from this "line" direction.
3467 // Here we are processing this case by inserting additional points
3468 // to the beginning/end of the WLine to make it more precise.
3469 // See description to the issue #30082.
3470
3471 const Standard_Real aSqTol3D = aTol3D*aTol3D;
3472 for (Standard_Integer j = 0; j < 2; j++)
3473 {
3474 // If j == 0 ==> add point at begin of WLine.
3475 // If j == 1 ==> add point at end of WLine.
3476
3477 for (;;)
3478 {
3479 if (aWLine[i]->NbPnts() >= aNbMaxPoints)
3480 {
3481 break;
3482 }
3483
3484 // Take 1st and 2nd point to compute the "line" direction.
3485 // For our convenience, we make 2nd point be the ends of the WLine
3486 // because it will be used for computation of the normals
3487 // to the surfaces.
3488 const Standard_Integer anIdx1 = j ? aWLine[i]->NbPnts() - 1 : 2;
3489 const Standard_Integer anIdx2 = j ? aWLine[i]->NbPnts() : 1;
3490
3491 const gp_Pnt &aP1 = aWLine[i]->Point(anIdx1).Value();
3492 const gp_Pnt &aP2 = aWLine[i]->Point(anIdx2).Value();
3493
3494 const gp_Vec aDir(aP1, aP2);
3495
3496 if (aDir.SquareMagnitude() < aSqTol3D)
3497 {
3498 break;
3499 }
3500
3501 // Compute tangent in first/last point of the WLine.
3502 // We do not take into account the flag "isReversed"
3503 // because strict direction of the tangent is not
3504 // important here (we are interested in the tangent
3505 // line itself and nothing to fear if its direction
3506 // is reversed).
3507 const gp_Vec aN1 = aQuad1.Normale(aP2);
3508 const gp_Vec aN2 = aQuad2.Normale(aP2);
3509 const gp_Vec aTg(aN1.Crossed(aN2));
3510
3511 if (aTg.SquareMagnitude() < Precision::SquareConfusion())
3512 {
3513 // Tangent zone
3514 break;
3515 }
3516
3517 // Check of the bending
3518 Standard_Real anAngle = aDir.Angle(aTg);
3519
3520 if (anAngle > M_PI_2)
3521 anAngle -= M_PI;
3522
3523 if (Abs(anAngle) > 0.25) // ~ 14deg.
3524 {
3525 const Standard_Integer aNbPntsPrev = aWLine[i]->NbPnts();
3526 SeekAdditionalPoints(aQuad1, aQuad2, aWLine[i]->Curve(),
3527 anEquationCoeffs, i, 3, anIdx1, anIdx2,
3528 aTol2D, aPeriod, isReversed);
3529
3530 if (aWLine[i]->NbPnts() == aNbPntsPrev)
3531 {
3532 // No points have been added. ==> Exit from a loop.
3533 break;
3534 }
3535 }
3536 else
3537 {
3538 // Good result has been achieved. ==> Exit from a loop.
3539 break;
3540 }
3541 } // for (;;)
3542 }
3543 }
ecc4f148 3544
261b7d9e 3545 if (isGood)
7c32c7c4 3546 {
3547 isTheEmpty = Standard_False;
3548 isAddedIntoWL[i] = Standard_True;
261b7d9e 3549 SeekAdditionalPoints(aQuad1, aQuad2, aWLine[i]->Curve(),
3550 anEquationCoeffs, i, aNbPoints, 1,
7eb3580b 3551 aWLine[i]->NbPnts(), aTol2D, aPeriod,
3552 isReversed);
7c32c7c4 3553
261b7d9e 3554 aWLine[i]->ComputeVertexParameters(aTol3D);
7c32c7c4 3555 theSlin.Append(aWLine[i]);
3556 }
3557 }
3558 else
3559 {
3560 isAddedIntoWL[i] = Standard_False;
baf72cd2 3561 }
b5ef9d91 3562
d30895f5 3563#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
cd0705f6 3564 aWLine[i]->Dump(0);
b5ef9d91 3565#endif
3566 }
3567 }
3568 }
3569
3570
7365fad6 3571 //Delete the points in theSPnt, which
3572 //lie at least in one of the line in theSlin.
261b7d9e 3573 for (Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
b5ef9d91 3574 {
261b7d9e 3575 for (Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); aNbLin++)
b5ef9d91 3576 {
261b7d9e 3577 Handle(IntPatch_WLine) aWLine1(Handle(IntPatch_WLine)::
3578 DownCast(theSlin.Value(aNbLin)));
b5ef9d91 3579
3580 const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1);
3581 const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aWLine1->NbPnts());
3582
3583 const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNbPnt).PntOn2S();
590b3f04 3584 if (aPntCur.IsSame(aPntFWL1, aTol3D) ||
3585 aPntCur.IsSame(aPntLWL1, aTol3D))
b5ef9d91 3586 {
3587 theSPnt.Remove(aNbPnt);
3588 aNbPnt--;
3589 break;
ecc4f148 3590 }
3591 }
3592 }
3593
261b7d9e 3594 //Try to add new points in the neighborhood of existing point
3595 for (Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
b5ef9d91 3596 {
261b7d9e 3597 // Standard algorithm (implemented above) could not find any
3598 // continuous curve in neighborhood of aPnt2S (e.g. because
3599 // this curve is too small; see tests\bugs\modalg_5\bug25292_35 and _36).
3600 // Here, we will try to find several new points nearer to aPnt2S.
3601
3602 // The algorithm below tries to find two points in every
3603 // intervals [u1 - aStepMax, u1] and [u1, u1 + aStepMax]
3604 // and every new point will be in maximal distance from
3605 // u1. If these two points exist they will be joined
3606 // by the intersection curve.
3607
b5ef9d91 3608 const IntPatch_Point& aPnt2S = theSPnt.Value(aNbPnt);
3609
3610 Standard_Real u1, v1, u2, v2;
3611 aPnt2S.Parameters(u1, v1, u2, v2);
3612
3613 Handle(IntSurf_LineOn2S) aL2S = new IntSurf_LineOn2S();
3614 Handle(IntPatch_WLine) aWLine = new IntPatch_WLine(aL2S, Standard_False);
98974dcc 3615 aWLine->SetCreatingWayInfo(IntPatch_WLine::IntPatch_WLImpImp);
b5ef9d91 3616
3617 //Define the index of WLine, which lies the point aPnt2S in.
b5ef9d91 3618 Standard_Integer anIndex = 0;
261b7d9e 3619
3620 Standard_Real anUf = 0.0, anUl = 0.0, aCurU2 = 0.0;
3621 if (isReversed)
b5ef9d91 3622 {
3623 anUf = Max(u2 - aStepMax, aUSurf1f);
261b7d9e 3624 anUl = Min(u2 + aStepMax, aUSurf1l);
b5ef9d91 3625 aCurU2 = u1;
3626 }
3627 else
3628 {
3629 anUf = Max(u1 - aStepMax, aUSurf1f);
261b7d9e 3630 anUl = Min(u1 + aStepMax, aUSurf1l);
b5ef9d91 3631 aCurU2 = u2;
3632 }
b5ef9d91 3633
261b7d9e 3634 const Standard_Real anUinf = anUf, anUsup = anUl, anUmid = 0.5*(anUf + anUl);
3635
3636 {
3637 //Find the value of anIndex variable.
3638 Standard_Real aDelta = RealLast();
3639 for (Standard_Integer i = 0; i < aNbWLines; i++)
b5ef9d91 3640 {
261b7d9e 3641 Standard_Real anU2t = 0.0;
3642 if (!ComputationMethods::CylCylComputeParameters(anUmid, i, anEquationCoeffs, anU2t))
3643 continue;
3644
3645 Standard_Real aDU2 = fmod(Abs(anU2t - aCurU2), aPeriod);
3646 aDU2 = Min(aDU2, Abs(aDU2 - aPeriod));
3647 if (aDU2 < aDelta)
3648 {
3649 aDelta = aDU2;
3650 anIndex = i;
3651 }
b5ef9d91 3652 }
3653 }
3654
261b7d9e 3655 // Bisection method is used in order to find every new point.
3656 // I.e. if we need to find intersection point in the interval [anUinf, anUmid]
3657 // we check the point anUC = 0.5*(anUinf+anUmid). If it is an intersection point
3658 // we try to find another point in the interval [anUinf, anUC] (because we find the point in
3659 // maximal distance from anUmid). If it is not then we try to find another point in the
3660 // interval [anUC, anUmid]. Next iterations will be made analogically.
3661 // When we find intersection point in the interval [anUmid, anUsup] we try to find
3662 // another point in the interval [anUC, anUsup] if anUC is intersection point and
3663 // in the interval [anUmid, anUC], otherwise.
3664
dcd768a4 3665 Standard_Real anAddedPar[2] = {isReversed ? u2 : u1, isReversed ? u2 : u1};
5b9e1842 3666
261b7d9e 3667 for (Standard_Integer aParID = 0; aParID < 2; aParID++)
3668 {
3669 if (aParID == 0)
b5ef9d91 3670 {
261b7d9e 3671 anUf = anUinf;
3672 anUl = anUmid;
3673 }
3674 else // if(aParID == 1)
3675 {
3676 anUf = anUmid;
3677 anUl = anUsup;
b5ef9d91 3678 }
3679
261b7d9e 3680 Standard_Real &aPar1 = (aParID == 0) ? anUf : anUl,
3681 &aPar2 = (aParID == 0) ? anUl : anUf;
3682
3683 while (Abs(aPar2 - aPar1) > aStepMin)
b5ef9d91 3684 {
261b7d9e 3685 Standard_Real anUC = 0.5*(anUf + anUl);
3686 Standard_Real aU2 = 0.0, aV1 = 0.0, aV2 = 0.0;
3687 Standard_Boolean isDone = ComputationMethods::
3688 CylCylComputeParameters(anUC, anIndex, anEquationCoeffs, aU2, aV1, aV2);
3689
3690 if (isDone)
3691 {
3692 if (Abs(aV1 - aVSurf1f) <= aTol2D)
3693 aV1 = aVSurf1f;
3694
3695 if (Abs(aV1 - aVSurf1l) <= aTol2D)
3696 aV1 = aVSurf1l;
3697
3698 if (Abs(aV2 - aVSurf2f) <= aTol2D)
3699 aV2 = aVSurf2f;
3700
3701 if (Abs(aV2 - aVSurf2l) <= aTol2D)
3702 aV2 = aVSurf2l;
3703
3704 isDone = AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
3705 Standard_True, gp_Pnt2d(anUC, aV1), gp_Pnt2d(aU2, aV2),
3706 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
3707 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l,
3708 aPeriod, aWLine->Curve(), anIndex, aTol3D,
3709 aTol2D, Standard_False, Standard_True);
3710 }
3711
3712 if (isDone)
3713 {
3714 anAddedPar[0] = Min(anAddedPar[0], anUC);
3715 anAddedPar[1] = Max(anAddedPar[1], anUC);
3716 aPar2 = anUC;
3717 }
3718 else
3719 {
3720 aPar1 = anUC;
3721 }
b5ef9d91 3722 }
3723 }
3724
261b7d9e 3725 //Fill aWLine by additional points
3726 if (anAddedPar[1] - anAddedPar[0] > aStepMin)
b5ef9d91 3727 {
261b7d9e 3728 for (Standard_Integer aParID = 0; aParID < 2; aParID++)
3729 {
3730 Standard_Real aU2 = 0.0, aV1 = 0.0, aV2 = 0.0;
3731 ComputationMethods::CylCylComputeParameters(anAddedPar[aParID], anIndex,
3732 anEquationCoeffs, aU2, aV1, aV2);
3733
3734 AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed, Standard_True,
3735 gp_Pnt2d(anAddedPar[aParID], aV1), gp_Pnt2d(aU2, aV2),
3736 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
3737 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod, aWLine->Curve(),
3738 anIndex, aTol3D, aTol2D, Standard_False, Standard_False);
3739 }
3740
3741 SeekAdditionalPoints(aQuad1, aQuad2, aWLine->Curve(),
b5ef9d91 3742 anEquationCoeffs, anIndex, aNbMinPoints,
261b7d9e 3743 1, aWLine->NbPnts(), aTol2D, aPeriod,
3744 isReversed);
b5ef9d91 3745
261b7d9e 3746 aWLine->ComputeVertexParameters(aTol3D);
b5ef9d91 3747 theSlin.Append(aWLine);
261b7d9e 3748
b5ef9d91 3749 theSPnt.Remove(aNbPnt);
3750 aNbPnt--;
3751 }
3752 }
261b7d9e 3753
e146ef9a 3754 return IntPatch_ImpImpIntersection::IntStatus_OK;
ecc4f148 3755}
7fd59977 3756
3757//=======================================================================
261b7d9e 3758//function : IntCyCy
3759//purpose :
3760//=======================================================================
3761IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1,
3762 const IntSurf_Quadric& theQuad2,
3763 const Standard_Real theTol3D,
3764 const Standard_Real theTol2D,
3765 const Bnd_Box2d& theUVSurf1,
3766 const Bnd_Box2d& theUVSurf2,
3767 Standard_Boolean& isTheEmpty,
3768 Standard_Boolean& isTheSameSurface,
3769 Standard_Boolean& isTheMultiplePoint,
3770 IntPatch_SequenceOfLine& theSlin,
3771 IntPatch_SequenceOfPoint& theSPnt)
3772{
3773 isTheEmpty = Standard_True;
3774 isTheSameSurface = Standard_False;
3775 isTheMultiplePoint = Standard_False;
3776 theSlin.Clear();
3777 theSPnt.Clear();
3778
3779 const gp_Cylinder aCyl1 = theQuad1.Cylinder(),
3780 aCyl2 = theQuad2.Cylinder();
3781
3782 IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
3783
3784 if (!anInter.IsDone())
3785 {
3786 return IntPatch_ImpImpIntersection::IntStatus_Fail;
3787 }
3788
3789 if(anInter.TypeInter() != IntAna_NoGeometricSolution)
3790 {
3791 if (CyCyAnalyticalIntersect(theQuad1, theQuad2, anInter,
3792 theTol3D, isTheEmpty,
3793 isTheSameSurface, isTheMultiplePoint,
3794 theSlin, theSPnt))
3795 {
3796 return IntPatch_ImpImpIntersection::IntStatus_OK;
3797 }
3798 }
3799
3800 //Here, intersection line is not an analytical curve(line, circle, ellipsis etc.)
3801
3802 Standard_Real aUSBou[2][2], aVSBou[2][2]; //const
3803
3804 theUVSurf1.Get(aUSBou[0][0], aVSBou[0][0], aUSBou[0][1], aVSBou[0][1]);
3805 theUVSurf2.Get(aUSBou[1][0], aVSBou[1][0], aUSBou[1][1], aVSBou[1][1]);
3806
3807 const Standard_Real aPeriod = 2.0*M_PI;
3808 const Standard_Integer aNbWLines = 2;
3809
3810 const ComputationMethods::stCoeffsValue anEquationCoeffs1(aCyl1, aCyl2);
3811 const ComputationMethods::stCoeffsValue anEquationCoeffs2(aCyl2, aCyl1);
3812
3813 //Boundaries.
3814 //Intersection result can include two non-connected regions
3815 //(see WorkWithBoundaries::BoundariesComputing(...) method).
3816 const Standard_Integer aNbOfBoundaries = 2;
3817 Bnd_Range anURange[2][aNbOfBoundaries]; //const
3818
3819 if (!WorkWithBoundaries::BoundariesComputing(anEquationCoeffs1, aPeriod, anURange[0]))
3820 return IntPatch_ImpImpIntersection::IntStatus_OK;
3821
3822 if (!WorkWithBoundaries::BoundariesComputing(anEquationCoeffs2, aPeriod, anURange[1]))
3823 return IntPatch_ImpImpIntersection::IntStatus_OK;
3824
3825 //anURange[*] can be in different periodic regions in
3826 //compare with First-Last surface. E.g. the surface
3827 //is full cylinder [0, 2*PI] but anURange is [5, 7].
3828 //Trivial common range computation returns [5, 2*PI] and
3829 //its summary length is 2*PI-5 == 1.28... only. That is wrong.
3830 //This problem can be solved by the following
3831 //algorithm:
3832 // 1. split anURange[*] by the surface boundary;
3833 // 2. shift every new range in order to inscribe it
3834 // in [Ufirst, Ulast] of cylinder;
3835 // 3. consider only common ranges between [Ufirst, Ulast]
3836 // and new ranges.
3837 //
3838 // In above example, we obtain following:
3839 // 1. two ranges: [5, 2*PI] and [2*PI, 7];
3840 // 2. after shifting: [5, 2*PI] and [0, 7-2*PI];
3841 // 3. Common ranges: ([5, 2*PI] and [0, 2*PI]) == [5, 2*PI],
3842 // ([0, 7-2*PI] and [0, 2*PI]) == [0, 7-2*PI];
3843 // 4. Their summary length is (2*PI-5)+(7-2*PI-0)==7-5==2 ==> GOOD.
3844
3845 Standard_Real aSumRange[2] = { 0.0, 0.0 };
3846 Handle(NCollection_IncAllocator) anAlloc = new NCollection_IncAllocator;
3847 for (Standard_Integer aCID = 0; aCID < 2; aCID++)
3848 {
1103eb60 3849 anAlloc->Reset(false);
261b7d9e 3850 NCollection_List<Bnd_Range> aListOfRng(anAlloc);
3851
3852 aListOfRng.Append(anURange[aCID][0]);
3853 aListOfRng.Append(anURange[aCID][1]);
3854
3855 const Standard_Real aSplitArr[3] = {aUSBou[aCID][0], aUSBou[aCID][1], 0.0};
3856
3857 NCollection_List<Bnd_Range>::Iterator anITrRng;
3858 for (Standard_Integer aSInd = 0; aSInd < 3; aSInd++)
3859 {
3860 NCollection_List<Bnd_Range> aLstTemp(aListOfRng);
3861 aListOfRng.Clear();
3862 for (anITrRng.Init(aLstTemp); anITrRng.More(); anITrRng.Next())
3863 {
1103eb60 3864 Bnd_Range& aRng = anITrRng.ChangeValue();
261b7d9e 3865 aRng.Split(aSplitArr[aSInd], aListOfRng, aPeriod);
3866 }
3867 }
3868
3869 anITrRng.Init(aListOfRng);
3870 for (; anITrRng.More(); anITrRng.Next())
3871 {
1103eb60 3872 Bnd_Range& aCurrRange = anITrRng.ChangeValue();
261b7d9e 3873
3874 Bnd_Range aBoundR;
3875 aBoundR.Add(aUSBou[aCID][0]);
3876 aBoundR.Add(aUSBou[aCID][1]);
3877
3878 if (!InscribeInterval(aUSBou[aCID][0], aUSBou[aCID][1],
3879 aCurrRange, theTol2D, aPeriod))
3880 {
3881 //If aCurrRange does not have common block with
3882 //[Ufirst, Ulast] of cylinder then we will try
3883 //to inscribe [Ufirst, Ulast] in the boundaries of aCurrRange.
3884 Standard_Real aF = 1.0, aL = 0.0;
3885 if (!aCurrRange.GetBounds(aF, aL))
3886 continue;
3887
3888 if ((aL < aUSBou[aCID][0]))
3889 {
3890 aCurrRange.Shift(aPeriod);
3891 }
3892 else if (aF > aUSBou[aCID][1])
3893 {
3894 aCurrRange.Shift(-aPeriod);
3895 }
3896 }
3897
3898 aBoundR.Common(aCurrRange);
3899
3900 const Standard_Real aDelta = aBoundR.Delta();
3901
3902 if (aDelta > 0.0)
3903 {
3904 aSumRange[aCID] += aDelta;
3905 }
3906 }
3907 }
3908
3909 //The bigger range the bigger number of points in Walking-line (WLine)
3910 //we will be able to add and consequently we will obtain more
3911 //precise intersection line.
3912 //Every point of WLine is determined as function from U1-parameter,
3913 //where U1 is U-parameter on 1st quadric.
3914 //Therefore, we should use quadric with bigger range as 1st parameter
3915 //in IntCyCy() function.
3916 //On the other hand, there is no point in reversing in case of
3917 //analytical intersection (when result is line, ellipse, point...).
3918 //This result is independent of the arguments order.
3919 const Standard_Boolean isToReverse = (aSumRange[1] > aSumRange[0]);
3920
3921 if (isToReverse)
3922 {
3923 const WorkWithBoundaries aBoundWork(theQuad2, theQuad1, anEquationCoeffs2,
3924 theUVSurf2, theUVSurf1, aNbWLines,
3925 aPeriod, theTol3D, theTol2D, Standard_True);
3926
3927 return CyCyNoGeometric(aCyl2, aCyl1, aBoundWork, anURange[1], aNbOfBoundaries,
3928 isTheEmpty, theSlin, theSPnt);
3929 }
3930 else
3931 {
3932 const WorkWithBoundaries aBoundWork(theQuad1, theQuad2, anEquationCoeffs1,
3933 theUVSurf1, theUVSurf2, aNbWLines,
3934 aPeriod, theTol3D, theTol2D, Standard_False);
3935
3936 return CyCyNoGeometric(aCyl1, aCyl2, aBoundWork, anURange[0], aNbOfBoundaries,
3937 isTheEmpty, theSlin, theSPnt);
3938 }
3939}
3940
3941//=======================================================================
7fd59977 3942//function : IntCySp
3943//purpose :
3944//=======================================================================
3945Standard_Boolean IntCySp(const IntSurf_Quadric& Quad1,
3946 const IntSurf_Quadric& Quad2,
3947 const Standard_Real Tol,
3948 const Standard_Boolean Reversed,
3949 Standard_Boolean& Empty,
3950 Standard_Boolean& Multpoint,
3951 IntPatch_SequenceOfLine& slin,
3952 IntPatch_SequenceOfPoint& spnt)
3953
3954{
3955 Standard_Integer i;
3956
3957 IntSurf_TypeTrans trans1,trans2;
3958 IntAna_ResultType typint;
3959 IntPatch_Point ptsol;
3960 gp_Circ cirsol;
3961
3962 gp_Cylinder Cy;
3963 gp_Sphere Sp;
3964
3965 if (!Reversed) {
3966 Cy = Quad1.Cylinder();
3967 Sp = Quad2.Sphere();
3968 }
3969 else {
3970 Cy = Quad2.Cylinder();
3971 Sp = Quad1.Sphere();
3972 }
3973 IntAna_QuadQuadGeo inter(Cy,Sp,Tol);
3974
3975 if (!inter.IsDone()) {return Standard_False;}
3976
3977 typint = inter.TypeInter();
3978 Standard_Integer NbSol = inter.NbSolutions();
3979 Empty = Standard_False;
3980
3981 switch (typint) {
3982
3983 case IntAna_Empty :
3984 {
3985 Empty = Standard_True;
3986 }
3987 break;
3988
3989 case IntAna_Point :
3990 {
3991 gp_Pnt psol(inter.Point(1));
3992 Standard_Real U1,V1,U2,V2;
3993 Quad1.Parameters(psol,U1,V1);
3994 Quad2.Parameters(psol,U2,V2);
3995 ptsol.SetValue(psol,Tol,Standard_True);
3996 ptsol.SetParameters(U1,V1,U2,V2);
3997 spnt.Append(ptsol);
3998 }
3999 break;
4000
4001 case IntAna_Circle:
4002 {
4003 cirsol = inter.Circle(1);
4004 gp_Vec Tgt;
4005 gp_Pnt ptref;
4006 ElCLib::D1(0.,cirsol,ptref,Tgt);
4007
4008 if (NbSol == 1) {
4009 gp_Vec TestCurvature(ptref,Sp.Location());
4010 gp_Vec Normsp,Normcyl;
4011 if (!Reversed) {
4012 Normcyl = Quad1.Normale(ptref);
4013 Normsp = Quad2.Normale(ptref);
4014 }
4015 else {
4016 Normcyl = Quad2.Normale(ptref);
4017 Normsp = Quad1.Normale(ptref);
4018 }
4019
4020 IntSurf_Situation situcyl;
4021 IntSurf_Situation situsp;
4022
4023 if (Normcyl.Dot(TestCurvature) > 0.) {
4024 situsp = IntSurf_Outside;
4025 if (Normsp.Dot(Normcyl) > 0.) {
4026 situcyl = IntSurf_Inside;
4027 }
4028 else {
4029 situcyl = IntSurf_Outside;
4030 }
4031 }
4032 else {
4033 situsp = IntSurf_Inside;
4034 if (Normsp.Dot(Normcyl) > 0.) {
4035 situcyl = IntSurf_Outside;
4036 }
4037 else {
4038 situcyl = IntSurf_Inside;
4039 }
4040 }
4041 Handle(IntPatch_GLine) glig;
4042 if (!Reversed) {
4043 glig = new IntPatch_GLine(cirsol, Standard_True, situcyl, situsp);
4044 }
4045 else {
4046 glig = new IntPatch_GLine(cirsol, Standard_True, situsp, situcyl);
4047 }
4048 slin.Append(glig);
4049 }
4050 else {
4051 if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
4052 trans1 = IntSurf_Out;
4053 trans2 = IntSurf_In;
4054 }
4055 else {
4056 trans1 = IntSurf_In;
4057 trans2 = IntSurf_Out;
4058 }
4059 Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
4060 slin.Append(glig);
4061
4062 cirsol = inter.Circle(2);
4063 ElCLib::D1(0.,cirsol,ptref,Tgt);
4064 Standard_Real qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
4065 if(qwe> 0.0000001) {
4066 trans1 = IntSurf_Out;
4067 trans2 = IntSurf_In;
4068 }
4069 else if(qwe<-0.0000001) {
4070 trans1 = IntSurf_In;
4071 trans2 = IntSurf_Out;
4072 }
4073 else {
4074 trans1=trans2=IntSurf_Undecided;
4075 }
4076 glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
4077 slin.Append(glig);
4078 }
4079 }
4080 break;
4081
4082 case IntAna_NoGeometricSolution:
4083 {
4084 gp_Pnt psol;
4085 Standard_Real U1,V1,U2,V2;
4086 IntAna_IntQuadQuad anaint(Cy,Sp,Tol);
4087 if (!anaint.IsDone()) {
4088 return Standard_False;
4089 }
4090
4091 if (anaint.NbPnt()==0 && anaint.NbCurve()==0) {
4092 Empty = Standard_True;
4093 }
4094 else {
4095
4096 NbSol = anaint.NbPnt();
4097 for (i = 1; i <= NbSol; i++) {
4098 psol = anaint.Point(i);
4099 Quad1.Parameters(psol,U1,V1);
4100 Quad2.Parameters(psol,U2,V2);
4101 ptsol.SetValue(psol,Tol,Standard_True);
4102 ptsol.SetParameters(U1,V1,U2,V2);
4103 spnt.Append(ptsol);
4104 }
4105
4106 gp_Pnt ptvalid,ptf,ptl;
4107 gp_Vec tgvalid;
4108 Standard_Real first,last,para;
4109 IntAna_Curve curvsol;
4110 Standard_Boolean tgfound;
4111 Standard_Integer kount;
4112
4113 NbSol = anaint.NbCurve();
4114 for (i = 1; i <= NbSol; i++) {
4115 curvsol = anaint.Curve(i);
4116 curvsol.Domain(first,last);
4117 ptf = curvsol.Value(first);
4118 ptl = curvsol.Value(last);
4119
4120 para = last;
4121 kount = 1;
4122 tgfound = Standard_False;
4123
4124 while (!tgfound) {
4125 para = (1.123*first + para)/2.123;
4126 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
4127 if (!tgfound) {
4128 kount ++;
4129 tgfound = kount > 5;
4130 }
4131 }
4132 Handle(IntPatch_ALine) alig;
4133 if (kount <= 5) {
4134 Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
4135 Quad1.Normale(ptvalid));
4136 if(qwe> 0.00000001) {
4137 trans1 = IntSurf_Out;
4138 trans2 = IntSurf_In;
4139 }
4140 else if(qwe<-0.00000001) {
4141 trans1 = IntSurf_In;
4142 trans2 = IntSurf_Out;
4143 }
4144 else {
4145 trans1=trans2=IntSurf_Undecided;
4146 }
4147 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
4148 }
4149 else {
4150 alig = new IntPatch_ALine(curvsol,Standard_False);
4151 }
4152 Standard_Boolean TempFalse1a = Standard_False;
4153 Standard_Boolean TempFalse2a = Standard_False;
4154
4155 //-- ptf et ptl : points debut et fin de alig
4156
4157 ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1a,ptf,first,
4158 TempFalse2a,ptl,last,Multpoint,Tol);
4159 slin.Append(alig);
4160 } //-- boucle sur les lignes
4161 } //-- solution avec au moins une lihne
4162 }
4163 break;
4164
4165 default:
4166 {
4167 return Standard_False;
4168 }
4169 }
4170 return Standard_True;
4171}
4172//=======================================================================
4173//function : IntCyCo
4174//purpose :
4175//=======================================================================
4176Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
4177 const IntSurf_Quadric& Quad2,
4178 const Standard_Real Tol,
4179 const Standard_Boolean Reversed,
4180 Standard_Boolean& Empty,
4181 Standard_Boolean& Multpoint,
4182 IntPatch_SequenceOfLine& slin,
4183 IntPatch_SequenceOfPoint& spnt)
4184
4185{
4186 IntPatch_Point ptsol;
4187
4188 Standard_Integer i;
4189
4190 IntSurf_TypeTrans trans1,trans2;
4191 IntAna_ResultType typint;
4192 gp_Circ cirsol;
4193
4194 gp_Cylinder Cy;
4195 gp_Cone Co;
4196
4197 if (!Reversed) {
4198 Cy = Quad1.Cylinder();
4199 Co = Quad2.Cone();
4200 }
4201 else {
4202 Cy = Quad2.Cylinder();
4203 Co = Quad1.Cone();
4204 }
4205 IntAna_QuadQuadGeo inter(Cy,Co,Tol);
4206
4207 if (!inter.IsDone()) {return Standard_False;}
4208
4209 typint = inter.TypeInter();
4210 Standard_Integer NbSol = inter.NbSolutions();
4211 Empty = Standard_False;
4212
4213 switch (typint) {
4214
4215 case IntAna_Empty : {
4216 Empty = Standard_True;
4217 }
4218 break;
4219
4220 case IntAna_Point :{
4221 gp_Pnt psol(inter.Point(1));
4222 Standard_Real U1,V1,U2,V2;
4223 Quad1.Parameters(psol,U1,V1);
4224 Quad1.Parameters(psol,U2,V2);
4225 ptsol.SetValue(psol,Tol,Standard_True);
4226 ptsol.SetParameters(U1,V1,U2,V2);
4227 spnt.Append(ptsol);
4228 }
4229 break;
4230
4231 case IntAna_Circle: {
4232 gp_Vec Tgt;
4233 gp_Pnt ptref;
4234 Standard_Integer j;
4235 Standard_Real qwe;
4236 //
4237 for(j=1; j<=2; ++j) {
4238 cirsol = inter.Circle(j);
4239 ElCLib::D1(0.,cirsol,ptref,Tgt);
4240 qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
4241 if(qwe> 0.00000001) {
4242 trans1 = IntSurf_Out;
4243 trans2 = IntSurf_In;
4244 }
4245 else if(qwe<-0.00000001) {
4246 trans1 = IntSurf_In;
4247 trans2 = IntSurf_Out;
4248 }
4249 else {
4250 trans1=trans2=IntSurf_Undecided;
4251 }
4252 Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
4253 slin.Append(glig);
4254 }
4255 }
4256 break;
4257
4258 case IntAna_NoGeometricSolution: {
4259 gp_Pnt psol;
4260 Standard_Real U1,V1,U2,V2;
4261 IntAna_IntQuadQuad anaint(Cy,Co,Tol);
4262 if (!anaint.IsDone()) {
4263 return Standard_False;
4264 }
4265
4266 if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) {
4267 Empty = Standard_True;
4268 }
4269 else {
4270 NbSol = anaint.NbPnt();
4271 for (i = 1; i <= NbSol; i++) {
4272 psol = anaint.Point(i);
4273 Quad1.Parameters(psol,U1,V1);
4274 Quad2.Parameters(psol,U2,V2);
4275 ptsol.SetValue(psol,Tol,Standard_True);
4276 ptsol.SetParameters(U1,V1,U2,V2);
4277 spnt.Append(ptsol);
4278 }
4279
4280 gp_Pnt ptvalid, ptf, ptl;
4281 gp_Vec tgvalid;
4282 //
4283 Standard_Real first,last,para;
4284 Standard_Boolean tgfound,firstp,lastp,kept;
4285 Standard_Integer kount;
4286 //
4287 //
4288 //IntAna_Curve curvsol;
4289 IntAna_Curve aC;
4290 IntAna_ListOfCurve aLC;
4291 IntAna_ListIteratorOfListOfCurve aIt;
7fd59977 4292
4293 //
4294 NbSol = anaint.NbCurve();
4295 for (i = 1; i <= NbSol; ++i) {
4296 kept = Standard_False;
4297 //curvsol = anaint.Curve(i);
4298 aC=anaint.Curve(i);
4299 aLC.Clear();
3306fdd9 4300 ExploreCurve(Co, aC, 10.*Tol, aLC);
7fd59977 4301 //
4302 aIt.Initialize(aLC);
4303 for (; aIt.More(); aIt.Next()) {
1103eb60 4304 IntAna_Curve& curvsol=aIt.ChangeValue();
7fd59977 4305 //
4306 curvsol.Domain(first, last);
4307 firstp = !curvsol.IsFirstOpen();
4308 lastp = !curvsol.IsLastOpen();
4309 if (firstp) {
4310 ptf = curvsol.Value(first);
4311 }
4312 if (lastp) {
4313 ptl = curvsol.Value(last);
4314 }
4315 para = last;
4316 kount = 1;
4317 tgfound = Standard_False;
4318
4319 while (!tgfound) {
4320 para = (1.123*first + para)/2.123;
4321 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
4322 if (!tgfound) {
4323 kount ++;
4324 tgfound = kount > 5;
4325 }
4326 }
4327 Handle(IntPatch_ALine) alig;
4328 if (kount <= 5) {
4329 Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
4330 Quad1.Normale(ptvalid));
4331 if(qwe> 0.00000001) {
4332 trans1 = IntSurf_Out;
4333 trans2 = IntSurf_In;
4334 }
4335 else if(qwe<-0.00000001) {
4336 trans1 = IntSurf_In;
4337 trans2 = IntSurf_Out;
4338 }
4339 else {
4340 trans1=trans2=IntSurf_Undecided;
4341 }
4342 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
4343 kept = Standard_True;
4344 }
4345 else {
4346 ptvalid = curvsol.Value(para);
4347 alig = new IntPatch_ALine(curvsol,Standard_False);
4348 kept = Standard_True;
d30895f5 4349 //-- std::cout << "Transition indeterminee" << std::endl;
7fd59977 4350 }
4351 if (kept) {
4352 Standard_Boolean Nfirstp = !firstp;
4353 Standard_Boolean Nlastp = !lastp;
4354 ProcessBounds(alig,slin,Quad1,Quad2,Nfirstp,ptf,first,
4355 Nlastp,ptl,last,Multpoint,Tol);
4356 slin.Append(alig);
4357 }
4358 } // for (; aIt.More(); aIt.Next())
4359 } // for (i = 1; i <= NbSol; ++i)
4360 }
4361 }
4362 break;
4363
4364 default:
4365 return Standard_False;
4366
4367 } // switch (typint)
4368
4369 return Standard_True;
4370}
4371//=======================================================================
4372//function : ExploreCurve
3306fdd9 4373//purpose : Splits aC on several curves in the cone apex points.
7fd59977 4374//=======================================================================
3306fdd9 4375Standard_Boolean ExploreCurve(const gp_Cone& theCo,
4376 IntAna_Curve& theCrv,
4377 const Standard_Real theTol,
4378 IntAna_ListOfCurve& theLC)
7fd59977 4379{
3306fdd9 4380 const Standard_Real aSqTol = theTol*theTol;
4381 const gp_Pnt aPapx(theCo.Apex());
4382
4383 Standard_Real aT1, aT2;
4384 theCrv.Domain(aT1, aT2);
4385
4386 theLC.Clear();
7fd59977 4387 //
3306fdd9 4388 TColStd_ListOfReal aLParams;
4389 theCrv.FindParameter(aPapx, aLParams);
4390 if (aLParams.IsEmpty())
4391 {
4392 theLC.Append(theCrv);
4393 return Standard_False;
7fd59977 4394 }
3306fdd9 4395
4396 for (TColStd_ListIteratorOfListOfReal anItr(aLParams); anItr.More(); anItr.Next())
4397 {
4398 Standard_Real aPrm = anItr.Value();
4399
4400 if ((aPrm - aT1) < Precision::PConfusion())
4401 continue;
4402
4403 Standard_Boolean isLast = Standard_False;
4404 if ((aT2 - aPrm) < Precision::PConfusion())
4405 {
4406 aPrm = aT2;
4407 isLast = Standard_True;
4408 }
4409
4410 const gp_Pnt aP = theCrv.Value(aPrm);
4411 const Standard_Real aSqD = aP.SquareDistance(aPapx);
4412 if (aSqD < aSqTol)
4413 {
4414 IntAna_Curve aC1 = theCrv;
4415 aC1.SetDomain(aT1, aPrm);
4416 aT1 = aPrm;
4417 theLC.Append(aC1);
4418 }
4419
4420 if (isLast)
4421 break;
7fd59977 4422 }
3306fdd9 4423
4424 if (theLC.IsEmpty())
4425 {
4426 theLC.Append(theCrv);
4427 return Standard_False;
7fd59977 4428 }
3306fdd9 4429
4430 if ((aT2 - aT1) > Precision::PConfusion())
4431 {
4432 IntAna_Curve aC1 = theCrv;
4433 aC1.SetDomain(aT1, aT2);
4434 theLC.Append(aC1);
7fd59977 4435 }
3306fdd9 4436
4437 return Standard_True;
7fd59977 4438}
Open CASCADE Technology repository