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