e2e0498b |
1 | //! Created on: 2016-06-03 |
2 | //! Created by: NIKOLAI BUKHALOV |
3 | //! Copyright (c) 2016 OPEN CASCADE SAS |
4 | // |
5 | // This file is part of Open CASCADE Technology software library. |
6 | // |
7 | // This library is free software; you can redistribute it and/or modify it under |
8 | // the terms of the GNU Lesser General Public License version 2.1 as published |
9 | // by the Free Software Foundation, with special exception defined in the file |
10 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
11 | // distribution for complete text of the license and disclaimer of any warranty. |
12 | // |
13 | // Alternatively, this file may be used under the terms of Open CASCADE |
14 | // commercial license or contractual agreement. |
15 | |
16 | #include <IntPatch_SpecialPoints.hxx> |
17 | |
18 | #include <Adaptor3d_HSurface.hxx> |
19 | #include <Extrema_ExtPS.hxx> |
20 | #include <Extrema_GenLocateExtPS.hxx> |
21 | #include <Geom_ConicalSurface.hxx> |
22 | #include <Geom_SphericalSurface.hxx> |
23 | #include <IntPatch_Point.hxx> |
24 | #include <IntSurf.hxx> |
25 | #include <IntSurf_PntOn2S.hxx> |
26 | #include <Standard_TypeMismatch.hxx> |
27 | #include <math_FunctionSetRoot.hxx> |
28 | #include <math_FunctionSetWithDerivatives.hxx> |
29 | #include <math_Matrix.hxx> |
30 | |
31 | // The function for searching intersection point, which |
32 | // lies in the seam-edge of the quadric definetely. |
33 | class FuncPreciseSeam: public math_FunctionSetWithDerivatives |
34 | { |
35 | public: |
36 | FuncPreciseSeam(const Handle(Adaptor3d_HSurface)& theQSurf, // quadric |
37 | const Handle(Adaptor3d_HSurface)& thePSurf, // another surface |
38 | const Standard_Boolean isTheUSeam): |
39 | myQSurf(theQSurf), |
40 | myPSurf(thePSurf), |
41 | mySeamCoordInd(isTheUSeam? 1 : 0) // Defines, U- or V-seam is used |
42 | { |
43 | }; |
44 | |
45 | virtual Standard_Integer NbVariables() const |
46 | { |
47 | return 3; |
48 | }; |
49 | |
50 | virtual Standard_Integer NbEquations() const |
51 | { |
52 | return 3; |
53 | } |
54 | |
55 | virtual Standard_Boolean Value(const math_Vector& theX, |
56 | math_Vector& theF) |
57 | { |
58 | try |
59 | { |
60 | const Standard_Integer anIndX = theX.Lower(), anIndF = theF.Lower(); |
61 | Standard_Real aUV[] = {0.0, 0.0}; |
62 | aUV[mySeamCoordInd] = theX(anIndX+2); |
63 | const gp_Pnt aP1(myPSurf->Value(theX(anIndX), theX(anIndX+1))); |
64 | const gp_Pnt aP2(myQSurf->Value(aUV[0], aUV[1])); |
65 | |
66 | (aP1.XYZ()-aP2.XYZ()).Coord(theF(anIndF), theF(anIndF+1), theF(anIndF+2)); |
67 | } |
68 | catch(Standard_Failure) |
69 | { |
70 | return Standard_False; |
71 | } |
72 | |
73 | return Standard_True; |
74 | }; |
75 | |
76 | virtual Standard_Boolean Derivatives(const math_Vector& theX, |
77 | math_Matrix& theD) |
78 | { |
79 | try |
80 | { |
81 | const Standard_Integer anIndX = theX.Lower(), |
82 | anIndRD = theD.LowerRow(), |
83 | anIndCD = theD.LowerCol(); |
84 | Standard_Real aUV[] = {0.0, 0.0}; |
85 | aUV[mySeamCoordInd] = theX(anIndX+2); |
86 | |
87 | gp_Pnt aPt; |
88 | |
89 | //0 for U-coordinate, 1 - for V one |
90 | gp_Vec aD1[2], aD2[2]; |
91 | myPSurf->D1(theX(anIndX), theX(anIndX+1), aPt, aD1[0], aD1[1]); |
92 | myQSurf->D1(aUV[0], aUV[1], aPt, aD2[0], aD2[1]); |
93 | |
94 | // d/dX1 |
95 | aD1[0].Coord(theD(anIndRD, anIndCD), |
96 | theD(anIndRD+1, anIndCD), theD(anIndRD+2, anIndCD)); |
97 | |
98 | // d/dX2 |
99 | aD1[1].Coord(theD(anIndRD, anIndCD+1), |
100 | theD(anIndRD+1, anIndCD+1), theD(anIndRD+2, anIndCD+1)); |
101 | |
102 | // d/dX3 |
103 | aD2[mySeamCoordInd].Reversed().Coord(theD(anIndRD, anIndCD+2), |
104 | theD(anIndRD+1, anIndCD+2), theD(anIndRD+2, anIndCD+2)); |
105 | } |
106 | catch(Standard_Failure) |
107 | { |
108 | return Standard_False; |
109 | } |
110 | |
111 | return Standard_True; |
112 | }; |
113 | |
114 | virtual Standard_Boolean Values (const math_Vector& theX, |
115 | math_Vector& theF, |
116 | math_Matrix& theD) |
117 | { |
118 | if(!Value(theX, theF)) |
119 | return Standard_False; |
120 | |
121 | if(!Derivatives(theX, theD)) |
122 | return Standard_False; |
123 | |
124 | return Standard_True; |
125 | } |
126 | |
127 | protected: |
128 | FuncPreciseSeam operator=(FuncPreciseSeam&); |
129 | |
130 | private: |
131 | const Handle(Adaptor3d_HSurface)& myQSurf; |
132 | const Handle(Adaptor3d_HSurface)& myPSurf; |
133 | |
134 | // 1 for U-coordinate, 0 - for V one. |
135 | const Standard_Integer mySeamCoordInd; |
136 | }; |
137 | |
138 | //======================================================================= |
139 | //function : IsPointOnSurface |
140 | //purpose : Checks if thePt is in theSurf (with given tolerance). |
141 | // Returns the foot of projection (theProjPt) and its parameters |
142 | // on theSurf. |
143 | //======================================================================= |
144 | static Standard_Boolean IsPointOnSurface(const Handle(Adaptor3d_HSurface)& theSurf, |
145 | const gp_Pnt& thePt, |
146 | const Standard_Real theTol, |
147 | gp_Pnt& theProjPt, |
148 | Standard_Real& theUpar, |
149 | Standard_Real& theVpar) |
150 | { |
151 | Standard_Boolean aRetVal = Standard_False; |
152 | |
153 | switch(theSurf->GetType()) |
154 | { |
155 | case GeomAbs_Plane: |
156 | case GeomAbs_Cylinder: |
157 | case GeomAbs_Cone: |
158 | case GeomAbs_Sphere: |
159 | case GeomAbs_Torus: |
160 | case GeomAbs_SurfaceOfExtrusion: |
161 | case GeomAbs_SurfaceOfRevolution: |
162 | { |
163 | Extrema_ExtPS anExtr(thePt, theSurf->Surface(), theSurf->UResolution(theTol), |
164 | theSurf->VResolution(theTol), Extrema_ExtFlag_MIN); |
165 | if(!anExtr.IsDone() || (anExtr.NbExt() < 1)) |
166 | { |
167 | aRetVal = Standard_False; |
168 | } |
169 | else |
170 | { |
171 | Standard_Integer anExtrIndex = 1; |
172 | Standard_Real aSqDistMin = anExtr.SquareDistance(anExtrIndex); |
173 | for(Standard_Integer i = anExtrIndex + 1; i <= anExtr.NbExt(); i++) |
174 | { |
175 | const Standard_Real aSqD = anExtr.SquareDistance(i); |
176 | if(aSqD < aSqDistMin) |
177 | { |
178 | aSqDistMin = aSqD; |
179 | anExtrIndex = i; |
180 | } |
181 | } |
182 | |
183 | if(aSqDistMin > theTol*theTol) |
184 | { |
185 | aRetVal = Standard_False; |
186 | } |
187 | else |
188 | { |
189 | theProjPt.SetXYZ(anExtr.Point(anExtrIndex).Value().XYZ()); |
190 | anExtr.Point(anExtrIndex).Parameter(theUpar, theVpar); |
191 | aRetVal = Standard_True; |
192 | } |
193 | } |
194 | } |
195 | break; |
196 | default: |
197 | { |
198 | Extrema_GenLocateExtPS anExtr(theSurf->Surface()); |
199 | anExtr.Perform(thePt, theUpar, theVpar); |
200 | if(!anExtr.IsDone() || (anExtr.SquareDistance() > theTol*theTol)) |
201 | { |
202 | aRetVal = Standard_False; |
203 | } |
204 | else |
205 | { |
206 | anExtr.Point().Parameter(theUpar, theVpar); |
207 | theProjPt.SetXYZ(anExtr.Point().Value().XYZ()); |
208 | aRetVal = Standard_True; |
209 | } |
210 | } |
211 | break; |
212 | } |
213 | |
214 | return aRetVal; |
215 | } |
216 | |
217 | //======================================================================= |
218 | //function : AddCrossUVIsoPoint |
219 | //purpose : theQSurf is the surface possibly containing special point, |
220 | // thePSurf is another surface to intersect. |
221 | //======================================================================= |
222 | Standard_Boolean IntPatch_SpecialPoints:: |
223 | AddCrossUVIsoPoint(const Handle(Adaptor3d_HSurface)& theQSurf, |
224 | const Handle(Adaptor3d_HSurface)& thePSurf, |
225 | const IntSurf_PntOn2S& theRefPt, |
226 | const Standard_Real theTol, |
227 | IntSurf_PntOn2S& theAddedPoint, |
228 | const Standard_Boolean theIsReversed) |
229 | { |
230 | Standard_Real anArrOfPeriod[4] = {0.0, 0.0, 0.0, 0.0}; |
231 | IntSurf::SetPeriod(theIsReversed ? thePSurf : theQSurf, |
232 | theIsReversed ? theQSurf : thePSurf, anArrOfPeriod); |
233 | |
234 | gp_Pnt aPQuad; |
235 | |
236 | //Not quadric point |
237 | Standard_Real aU0 = 0.0, aV0 = 0.0; |
238 | if(theIsReversed) |
239 | theRefPt.ParametersOnS1(aU0, aV0); |
240 | else |
241 | theRefPt.ParametersOnS2(aU0, aV0); |
242 | |
243 | //Quadric point |
244 | Standard_Real aUquad = 0.0, aVquad = 0.0; |
245 | |
246 | theQSurf->D0(aUquad, aVquad, aPQuad); |
247 | |
248 | Extrema_GenLocateExtPS anExtr(thePSurf->Surface()); |
249 | anExtr.Perform(aPQuad, aU0, aV0); |
250 | |
251 | if(!anExtr.IsDone()) |
252 | { |
253 | return Standard_False; |
254 | } |
255 | |
256 | if(anExtr.SquareDistance() > theTol*theTol) |
257 | { |
258 | return Standard_False; |
259 | } |
260 | |
261 | anExtr.Point().Parameter(aU0, aV0); |
262 | gp_Pnt aP0(anExtr.Point().Value()); |
263 | |
264 | if(theIsReversed) |
265 | theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aU0, aV0, aUquad, aVquad); |
266 | else |
267 | theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0); |
268 | |
269 | AdjustPointAndVertex(theRefPt, anArrOfPeriod, theAddedPoint); |
270 | |
271 | return Standard_True; |
272 | } |
273 | |
274 | //======================================================================= |
275 | //function : AddPointOnUorVIso |
276 | //purpose : theQSurf is the surface possibly containing special point, |
277 | // thePSurf is another surface to intersect. |
278 | //======================================================================= |
279 | Standard_Boolean IntPatch_SpecialPoints:: |
280 | AddPointOnUorVIso(const Handle(Adaptor3d_HSurface)& theQSurf, |
281 | const Handle(Adaptor3d_HSurface)& thePSurf, |
282 | const IntSurf_PntOn2S& theRefPt, |
283 | const Standard_Boolean theIsU, |
284 | const math_Vector& theToler, |
285 | const math_Vector& theInitPoint, |
286 | const math_Vector& theInfBound, |
287 | const math_Vector& theSupBound, |
288 | IntSurf_PntOn2S& theAddedPoint, |
289 | const Standard_Boolean theIsReversed) |
290 | { |
291 | Standard_Real anArrOfPeriod[4] = {0.0, 0.0, 0.0, 0.0}; |
292 | IntSurf::SetPeriod(theIsReversed ? thePSurf : theQSurf, |
293 | theIsReversed ? theQSurf : thePSurf, anArrOfPeriod); |
294 | |
295 | FuncPreciseSeam aF(theQSurf, thePSurf, theIsU); |
296 | |
297 | math_FunctionSetRoot aSRF(aF, theToler); |
298 | aSRF.Perform(aF, theInitPoint, theInfBound, theSupBound); |
299 | |
300 | if(!aSRF.IsDone()) |
301 | { |
302 | return Standard_False; |
303 | } |
304 | |
305 | math_Vector aRoots(theInitPoint.Lower(), theInitPoint.Upper()); |
306 | aSRF.Root(aRoots); |
307 | |
308 | //On parametric |
309 | Standard_Real aU0 = aRoots(1), aV0 = aRoots(2); |
310 | |
311 | //On quadric |
312 | Standard_Real aUquad = theIsU ? 0.0 : aRoots(3); |
313 | Standard_Real aVquad = theIsU ? aRoots(3) : 0.0; |
314 | const gp_Pnt aPQuad(theQSurf->Value(aUquad, aVquad)); |
315 | const gp_Pnt aP0(thePSurf->Value(aU0, aV0)); |
316 | |
317 | if(theIsReversed) |
318 | theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aU0, aV0, aUquad, aVquad); |
319 | else |
320 | theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0); |
321 | |
322 | AdjustPointAndVertex(theRefPt, anArrOfPeriod, theAddedPoint); |
323 | return Standard_True; |
324 | } |
325 | |
326 | //======================================================================= |
327 | //function : AddSingularPole |
328 | //purpose : theQSurf is the surface possibly containing special point, |
329 | // thePSurf is another surface to intersect. |
330 | // Returns TRUE, if the pole is an intersection point. |
331 | //======================================================================= |
332 | Standard_Boolean IntPatch_SpecialPoints:: |
333 | AddSingularPole(const Handle(Adaptor3d_HSurface)& theQSurf, |
334 | const Handle(Adaptor3d_HSurface)& thePSurf, |
335 | const IntSurf_PntOn2S& thePtIso, |
336 | const Standard_Real theTol, |
337 | IntPatch_Point& theVertex, |
338 | IntSurf_PntOn2S& theAddedPoint, |
339 | const Standard_Boolean theIsReversed, |
340 | const Standard_Boolean theIsReqRefCheck) |
341 | { |
342 | const Standard_Real aUpPeriod = thePSurf->IsUPeriodic() ? thePSurf->UPeriod() : 0.0; |
343 | const Standard_Real aUqPeriod = theQSurf->IsUPeriodic() ? theQSurf->UPeriod() : 0.0; |
344 | const Standard_Real aVpPeriod = thePSurf->IsVPeriodic() ? thePSurf->VPeriod() : 0.0; |
345 | const Standard_Real aVqPeriod = theQSurf->IsVPeriodic() ? theQSurf->VPeriod() : 0.0; |
346 | |
347 | const Standard_Real anArrOfPeriod[4] = {theIsReversed? aUpPeriod : aUqPeriod, |
348 | theIsReversed? aVpPeriod : aVqPeriod, |
349 | theIsReversed? aUqPeriod : aUpPeriod, |
350 | theIsReversed? aVqPeriod : aVpPeriod}; |
351 | |
352 | //On parametric |
353 | Standard_Real aU0 = 0.0, aV0 = 0.0; |
354 | //aPQuad is Pole |
355 | gp_Pnt aPQuad, aP0; |
356 | Standard_Real aUquad = 0.0, aVquad = 0.0; |
357 | if(theIsReversed) |
358 | theVertex.Parameters(aU0, aV0, aUquad, aVquad); |
359 | else |
360 | theVertex.Parameters(aUquad, aVquad, aU0, aV0); |
361 | |
362 | aUquad = 0.0; |
363 | |
364 | if(theQSurf->GetType() == GeomAbs_Sphere) |
365 | { |
366 | aVquad = Sign(M_PI_2, aVquad); |
367 | } |
368 | else if(theQSurf->GetType() == GeomAbs_Cone) |
369 | { |
370 | const gp_Cone aCo = theQSurf->Cone(); |
371 | const Standard_Real aRadius = aCo.RefRadius(); |
372 | const Standard_Real aSemiAngle = aCo.SemiAngle(); |
373 | aVquad = -aRadius / sin(aSemiAngle); |
374 | } |
375 | else |
376 | { |
9775fa61 |
377 | throw Standard_TypeMismatch( "IntPatch_SpecialPoints::AddSingularPole()," |
e2e0498b |
378 | "Unsupported quadric with Pole"); |
379 | } |
380 | |
381 | theQSurf->D0(aUquad, aVquad, aPQuad); |
382 | |
383 | if (theIsReqRefCheck && (aPQuad.SquareDistance(theVertex.Value()) >= theTol*theTol)) |
384 | { |
385 | return Standard_False; |
386 | } |
387 | |
388 | if(!IsPointOnSurface(thePSurf, aPQuad, theTol, aP0, aU0, aV0)) |
389 | { |
390 | return Standard_False; |
391 | } |
392 | |
393 | //Pole is an intersection point |
394 | //(lies in the quadric and the parametric surface) |
395 | |
396 | if(theIsReversed) |
397 | theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aU0, aV0, aUquad, aVquad); |
398 | else |
399 | theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0); |
400 | |
401 | Standard_Boolean isSame = Standard_False; |
402 | |
403 | if (theAddedPoint.IsSame(theVertex.PntOn2S(), Precision::Confusion())) |
404 | { |
405 | isSame = Standard_True; |
406 | } |
407 | |
408 | //Found pole does not exist in the Walking-line |
409 | //It must be added there (with correct 2D-parameters) |
410 | |
411 | //2D-parameters of theparametric surface have already been found (aU0, aV0). |
412 | //Let find 2D-parameters on the quadric. |
413 | |
414 | //The algorithm depends on the type of the quadric. Here we consider a Sphere only. |
415 | //Analogical result can be made for another types (e.g. cone, but formulas will |
416 | //be different) in case of need. |
417 | |
418 | //First of all, we need in adjusting thePSurf in the coordinate system of the Sphere |
419 | //(in order to make the equation of the sphere maximal simple). However, as it will be |
420 | //shown later, thePSurf is used in algorithm in order to get its derivatives. |
421 | //Therefore, for improving performance, transformation of these vectors is enough |
422 | //(there is no point in transformation of full surface). |
423 | |
424 | gp_Pnt aPtemp; |
425 | gp_Vec aVecDu, aVecDv; |
426 | thePSurf->D1(aU0, aV0, aPtemp, aVecDu, aVecDv); |
427 | |
428 | //Transforms parametric surface in coordinate-system of the quadric |
429 | gp_Trsf aTr; |
430 | aTr.SetTransformation((theQSurf->GetType() == GeomAbs_Sphere) ? |
431 | theQSurf->Sphere().Position() : |
432 | theQSurf->Cone().Position()); |
433 | |
434 | //Derivatives of transformed thePSurf |
435 | aVecDu.Transform(aTr); |
436 | aVecDv.Transform(aTr); |
437 | |
438 | Standard_Boolean isIsoChoosen = Standard_False; |
439 | |
440 | if(theQSurf->GetType() == GeomAbs_Sphere) |
441 | { |
442 | //The intersection point (including the pole) |
443 | //must be satisfied to the following system: |
444 | |
445 | // \left\{\begin{matrix} |
446 | // R*\cos (U_{q})*\cos (V_{q})=S_{x}(U_{s},V_{s}) |
447 | // R*\sin (U_{q})*\cos (V_{q})=S_{y}(U_{s},V_{s}) |
448 | // R*\sin (V_{q})=S_{z}(U_{s},V_{s}) |
449 | // \end{matrix}\right, |
450 | //where |
451 | // R is the radius of the sphere; |
452 | // @S_{x}@, @S_{y}@ and @S_{z}@ are X, Y and Z-coordinates of thePSurf; |
453 | // @U_{s}@ and @V_{s}@ are equal to aU0 and aV0 corespondingly; |
454 | // @U_{q}@ and @V_{q}@ are equal to aUquad and aVquad corespondingly. |
455 | |
456 | //Consequently (from first two equations), |
457 | // \left\{\begin{matrix} |
458 | // \cos (U_{q}) = \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})} |
459 | // \sin (U_{q}) = \frac{S_{y}(U_{s},V_{s})}{R*\cos (V_{q})} |
460 | // \end{matrix}\right. |
461 | |
462 | //For pole, |
463 | // V_{q}=\pm \pi /2 \Rightarrow \cos (V_{q}) = 0 (denominator is equal to 0). |
464 | |
465 | //Therefore, computation U_{q} directly is impossibly. |
466 | // |
467 | //Let @V_{q}@ tends to @\pm \pi /2@. |
468 | //Then (indeterminate form is evaluated in accordance of L'Hospital rule), |
469 | // \cos (U_{q}) = \lim_{V_{q} \to (\pi /2-0)} |
470 | // \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})}= |
471 | // -\lim_{V_{q} \to (\pi /2-0)} |
472 | // \frac{\frac{\partial S_{x}} |
473 | // {\partial U_{s}}*\frac{\mathrm{d} U_{s}} |
474 | // {\mathrm{d} V_{q}}+\frac{\partial S_{x}} |
475 | // {\partial V_{s}}*\frac{\mathrm{d} V_{s}} |
476 | // {\mathrm{d} V_{q}}}{R*\sin (V_{q})} = |
477 | // -\frac{1}{R}*\frac{\mathrm{d} U_{s}} |
478 | // {\mathrm{d} V_{q}}*(\frac{\partial S_{x}} |
479 | // {\partial U_{s}}+\frac{\partial S_{x}} |
480 | // {\partial V_{s}}*\frac{\mathrm{d} V_{s}} |
481 | // {\mathrm{d} U_{s}}) = |
482 | // -\frac{1}{R}*\frac{\mathrm{d} V_{s}} |
483 | // {\mathrm{d} V_{q}}*(\frac{\partial S_{x}} |
484 | // {\partial U_{s}}*\frac{\mathrm{d} U_{s}} |
485 | // {\mathrm{d} V_{s}}+\frac{\partial S_{x}} |
486 | // {\partial V_{s}}). |
487 | |
488 | //Analogicaly for @\sin (U_{q})@ (@S_{x}@ is substituted to @S_{y}@). |
489 | |
490 | //Let mean, that |
491 | // \cos (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \cos (U_{q}) \left | _{V_{q} \to (\pi /2-0)} |
492 | // \sin (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \sin (U_{q}) \left | _{V_{q} \to (\pi /2-0)} |
493 | |
494 | //From the 3rd equation of the system, we obtain |
495 | // \frac{\mathrm{d} (R*\sin (V_{q}))}{\mathrm{d} V_{q}} = |
496 | // \frac{\mathrm{d} S_{z}(U_{s},V_{s})}{\mathrm{d} V_{q}} |
497 | //or |
498 | // R*\cos (V_{q}) = \frac{\partial S_{z}}{\partial U_{s}}* |
499 | // \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}} |
500 | // {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}}. |
501 | |
502 | //If @V_{q}=\pm \pi /2@, then |
503 | // \frac{\partial S_{z}}{\partial U_{s}}* |
504 | // \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}} |
505 | // {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}} = 0. |
506 | |
507 | //Consequently, if @\frac{\partial S_{z}}{\partial U_{s}} \neq 0 @ then |
508 | // \frac{\mathrm{d} U_{s}}{\mathrm{d} V_{s}} = |
509 | // -\frac{\frac{\partial S_{z}}{\partial V_{s}}} |
510 | // {\frac{\partial S_{z}}{\partial U_{s}}}. |
511 | |
512 | //If @ \frac{\partial S_{z}}{\partial V_{s}} \neq 0 @ then |
513 | // \frac{\mathrm{d} V_{s}}{\mathrm{d} U_{s}} = |
514 | // -\frac{\frac{\partial S_{z}}{\partial U_{s}}} |
515 | // {\frac{\partial S_{z}}{\partial V_{s}}} |
516 | |
517 | //Cases, when @ \frac{\partial S_{z}}{\partial U_{s}} = |
518 | //\frac{\partial S_{z}}{\partial V_{s}} = 0 @ are not consider here. |
519 | //The reason is written below. |
520 | |
521 | //Vector with {@ \cos (U_{q}) @, @ \sin (U_{q}) @} coordinates. |
522 | //Ask to pay attention to the fact that this vector is always normalyzed. |
523 | gp_Vec2d aV1; |
524 | |
525 | if( (Abs(aVecDu.Z()) < Precision::PConfusion()) && |
526 | (Abs(aVecDv.Z()) < Precision::PConfusion())) |
527 | { |
528 | //Example of this case is an intersection of a plane with a sphere |
529 | //when the plane tangents the sphere in some pole (i.e. only one |
530 | //intersection point, not line). In this case, U-coordinate of the |
531 | //sphere is undefined (can be realy anything). |
532 | //Another reason is that we have tangent zone around the pole |
533 | //(see bug #26576). |
534 | //Computation of correct value of aUquad is impossible. |
535 | //Therefore, (in oreder to return something) we will consider |
536 | //the intersection line goes along some isoline in neighbourhood |
537 | //of the pole. |
538 | |
539 | #ifdef INTPATCH_ADDSPECIALPOINTS_DEBUG |
540 | cout << "Cannot find UV-coordinate for quadric in the pole." |
541 | " See considered comment above. IntPatch_AddSpecialPoints.cxx," |
542 | " AddSingularPole(...)" << endl; |
543 | #endif |
544 | Standard_Real aUIso = 0.0, aVIso = 0.0; |
545 | if(theIsReversed) |
546 | thePtIso.ParametersOnS2(aUIso, aVIso); |
547 | else |
548 | thePtIso.ParametersOnS1(aUIso, aVIso); |
549 | |
550 | aUquad = aUIso; |
551 | isIsoChoosen = Standard_True; |
552 | } |
553 | else |
554 | { |
555 | if(Abs(aVecDu.Z()) > Abs(aVecDv.Z())) |
556 | { |
557 | const Standard_Real aDusDvs = aVecDv.Z()/aVecDu.Z(); |
558 | aV1.SetCoord( aVecDu.X()*aDusDvs - aVecDv.X(), |
559 | aVecDu.Y()*aDusDvs - aVecDv.Y()); |
560 | } |
561 | else |
562 | { |
563 | const Standard_Real aDvsDus = aVecDu.Z()/aVecDv.Z(); |
564 | aV1.SetCoord( aVecDv.X()*aDvsDus - aVecDu.X(), |
565 | aVecDv.Y()*aDvsDus - aVecDu.Y()); |
566 | } |
567 | |
568 | aV1.Normalize(); |
569 | |
570 | if(Abs(aV1.X()) > Abs(aV1.Y())) |
571 | aUquad = Sign(asin(aV1.Y()), aVquad); |
572 | else |
573 | aUquad = Sign(acos(aV1.X()), aVquad); |
574 | } |
575 | } |
576 | else //if(theQSurf->GetType() == GeomAbs_Cone) |
577 | { |
578 | // This case is not processed. However, |
579 | // it can be done using the same algorithm |
580 | // as for sphere (formulas will be different). |
581 | return Standard_False; |
582 | } |
583 | |
584 | if(theIsReversed) |
585 | theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aU0, aV0, aUquad, aVquad); |
586 | else |
587 | theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0); |
588 | |
589 | if (isSame) |
590 | { |
591 | theVertex.SetValue(theAddedPoint); |
592 | return Standard_True; |
593 | } |
594 | |
595 | if (!isIsoChoosen) |
596 | { |
597 | AdjustPointAndVertex(theVertex.PntOn2S(), anArrOfPeriod, theAddedPoint); |
598 | } |
599 | else |
600 | { |
601 | theVertex.SetValue(theAddedPoint); |
602 | } |
603 | |
604 | return Standard_True; |
605 | } |
606 | |
607 | //======================================================================= |
608 | //function : ContinueAfterSpecialPoint |
609 | //purpose : |
610 | //======================================================================= |
611 | Standard_Boolean IntPatch_SpecialPoints:: |
612 | ContinueAfterSpecialPoint(const Handle(Adaptor3d_HSurface)& theQSurf, |
613 | const Handle(Adaptor3d_HSurface)& thePSurf, |
614 | const IntSurf_PntOn2S& theRefPt, |
615 | const IntPatch_SpecPntType theSPType, |
616 | const Standard_Real theTol2D, |
617 | IntSurf_PntOn2S& theNewPoint, |
618 | const Standard_Boolean theIsReversed) |
619 | { |
620 | if(theSPType == IntPatch_SPntNone) |
621 | return Standard_False; |
622 | |
623 | //If the last point of the line is the pole of the quadric. |
624 | //In this case, Walking-line has been broken in this point. |
625 | //However, new line must start from this point. Here we must |
626 | //find its 2D-coordinates. |
627 | |
628 | //For sphere and cone, some intersection point is satisfied to the system |
629 | // \cos(U_{q}) = S_{x}(U_{s},V_{s})/F(V_{q}) |
630 | // \sin(U_{q}) = S_{y}(U_{s},V_{s})/F(V_{q}) |
631 | |
632 | //where |
633 | // @S_{x}@, @S_{y}@ are X and Y-coordinates of thePSurf; |
634 | // @U_{s}@ and @V_{s}@ are UV-parameters on thePSurf; |
635 | // @U_{q}@ and @V_{q}@ are UV-parameters on theQSurf; |
636 | // @F(V_{q}) @ is some function, which value independs on @U_{q}@ |
637 | // (form of this function depends on the type of the quadric). |
638 | |
639 | //When we go through the pole/apex, the function @F(V_{q}) @ changes sign. |
640 | //Therefore, some cases are possible, when only @\cos(U_{q}) @ or |
641 | //only @ \sin(U_{q}) @ change sign. |
642 | |
643 | //Consequently, when the line goes throug the pole, @U_{q}@ can be |
644 | //changed on @\pi /2 @ (but not less). |
645 | |
646 | if(theNewPoint.IsSame(theRefPt, Precision::Confusion(), theTol2D)) |
647 | { |
648 | return Standard_False; |
649 | } |
650 | |
651 | //Here, in case of pole/apex adding, we forbid "jumping" between two neighbor |
652 | //Walking-point with step greater than pi/4 |
653 | const Standard_Real aPeriod = (theSPType == IntPatch_SPntPole)? M_PI_2 : 2.0*M_PI; |
654 | |
655 | const Standard_Real aUpPeriod = thePSurf->IsUPeriodic() ? thePSurf->UPeriod() : 0.0; |
656 | const Standard_Real aUqPeriod = theQSurf->IsUPeriodic() ? aPeriod : 0.0; |
657 | const Standard_Real aVpPeriod = thePSurf->IsVPeriodic() ? thePSurf->VPeriod() : 0.0; |
658 | const Standard_Real aVqPeriod = theQSurf->IsVPeriodic() ? aPeriod : 0.0; |
659 | |
660 | const Standard_Real anArrOfPeriod[4] = {theIsReversed? aUpPeriod : aUqPeriod, |
661 | theIsReversed? aVpPeriod : aVqPeriod, |
662 | theIsReversed? aUqPeriod : aUpPeriod, |
663 | theIsReversed? aVqPeriod : aVpPeriod}; |
664 | |
665 | AdjustPointAndVertex(theRefPt, anArrOfPeriod, theNewPoint); |
666 | return Standard_True; |
667 | } |
668 | |
669 | //======================================================================= |
670 | //function : AdjustPointAndVertex |
671 | //purpose : |
672 | //======================================================================= |
673 | void IntPatch_SpecialPoints:: |
674 | AdjustPointAndVertex(const IntSurf_PntOn2S &theRefPoint, |
675 | const Standard_Real theArrPeriods[4], |
676 | IntSurf_PntOn2S &theNewPoint, |
677 | IntPatch_Point* const theVertex) |
678 | { |
679 | Standard_Real aRefPar[2] = {0.0, 0.0}; |
680 | Standard_Real aPar[4] = {0.0, 0.0, 0.0, 0.0}; |
681 | theNewPoint.Parameters(aPar[0], aPar[1], aPar[2], aPar[3]); |
682 | |
683 | for(Standard_Integer i = 0; i < 4; i++) |
684 | { |
685 | if(theArrPeriods[i] == 0) |
686 | continue; |
687 | |
688 | const Standard_Real aPeriod = theArrPeriods[i], aHalfPeriod = 0.5*theArrPeriods[i]; |
689 | |
690 | if(i < 2) |
691 | {// 1st surface is used |
692 | theRefPoint.ParametersOnS1(aRefPar[0], aRefPar[1]); |
693 | } |
694 | else |
695 | { |
696 | theRefPoint.ParametersOnS2(aRefPar[0], aRefPar[1]); |
697 | } |
698 | |
699 | const Standard_Integer aRefInd = i%2; |
700 | |
701 | { |
702 | Standard_Real aDeltaPar = aRefPar[aRefInd]-aPar[i]; |
703 | const Standard_Real anIncr = Sign(aPeriod, aDeltaPar); |
704 | while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod)) |
705 | { |
706 | aPar[i] += anIncr; |
707 | aDeltaPar = aRefPar[aRefInd]-aPar[i]; |
708 | } |
709 | } |
710 | } |
711 | |
712 | if(theVertex) |
713 | (*theVertex).SetParameters(aPar[0], aPar[1], aPar[2], aPar[3]); |
714 | |
715 | theNewPoint.SetValue(aPar[0], aPar[1], aPar[2], aPar[3]); |
716 | } |
717 | |