7cda071df3abc5954722042e067c22e9f6bba7f8
[occt.git] / src / IntPatch / IntPatch_SpecialPoints.cxx
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 <ElCLib.hxx>
20 #include <Extrema_ExtPS.hxx>
21 #include <Extrema_GenLocateExtPS.hxx>
22 #include <Geom_ConicalSurface.hxx>
23 #include <Geom_SphericalSurface.hxx>
24 #include <IntPatch_Point.hxx>
25 #include <IntSurf.hxx>
26 #include <IntSurf_PntOn2S.hxx>
27 #include <Standard_TypeMismatch.hxx>
28 #include <math_FunctionSetRoot.hxx>
29 #include <math_FunctionSetWithDerivatives.hxx>
30 #include <math_Matrix.hxx>
31
32 // The function for searching intersection point, which 
33 // lies in the seam-edge of the quadric definetely.
34 class FuncPreciseSeam: public math_FunctionSetWithDerivatives
35 {
36 public:
37   FuncPreciseSeam(const Handle(Adaptor3d_HSurface)& theQSurf, // quadric
38                   const Handle(Adaptor3d_HSurface)& thePSurf, // another surface
39                   const Standard_Boolean isTheUSeam,
40                   const Standard_Real theIsoParameter): 
41         myQSurf(theQSurf),
42         myPSurf(thePSurf),
43         mySeamCoordInd(isTheUSeam? 1 : 0), // Defines, U- or V-seam is used
44         myIsoParameter(theIsoParameter)
45   {
46   };
47   
48   virtual Standard_Integer NbVariables() const
49   {
50     return 3;
51   };
52
53   virtual Standard_Integer NbEquations() const
54   {
55     return 3;
56   }
57
58   virtual Standard_Boolean Value(const math_Vector& theX,
59                                  math_Vector& theF)
60   {
61     try
62     {
63       const Standard_Integer anIndX = theX.Lower(), anIndF = theF.Lower();
64       Standard_Real aUV[] = {myIsoParameter, myIsoParameter};
65       aUV[mySeamCoordInd] = theX(anIndX+2);
66       const gp_Pnt aP1(myPSurf->Value(theX(anIndX), theX(anIndX+1)));
67       const gp_Pnt aP2(myQSurf->Value(aUV[0], aUV[1]));
68
69       (aP1.XYZ()-aP2.XYZ()).Coord(theF(anIndF), theF(anIndF+1), theF(anIndF+2));
70     }
71     catch(Standard_Failure)
72     {
73       return Standard_False;
74     }
75
76     return Standard_True;
77   };
78
79   virtual Standard_Boolean Derivatives(const math_Vector& theX,
80                                        math_Matrix& theD)
81   {
82     try
83     {
84       const Standard_Integer anIndX = theX.Lower(),
85                              anIndRD = theD.LowerRow(),
86                              anIndCD = theD.LowerCol();
87       Standard_Real aUV[] = {myIsoParameter, myIsoParameter};
88       aUV[mySeamCoordInd] = theX(anIndX+2);
89
90       gp_Pnt aPt;
91
92       //0 for U-coordinate, 1 - for V one
93       gp_Vec aD1[2], aD2[2];
94       myPSurf->D1(theX(anIndX), theX(anIndX+1), aPt, aD1[0], aD1[1]);
95       myQSurf->D1(aUV[0], aUV[1], aPt, aD2[0], aD2[1]);
96
97       // d/dX1
98       aD1[0].Coord(theD(anIndRD, anIndCD),
99                       theD(anIndRD+1, anIndCD), theD(anIndRD+2, anIndCD));
100
101       // d/dX2
102       aD1[1].Coord(theD(anIndRD, anIndCD+1),
103                       theD(anIndRD+1, anIndCD+1), theD(anIndRD+2, anIndCD+1));
104
105       // d/dX3
106       aD2[mySeamCoordInd].Reversed().Coord(theD(anIndRD, anIndCD+2),
107                                 theD(anIndRD+1, anIndCD+2), theD(anIndRD+2, anIndCD+2));
108     }
109     catch(Standard_Failure)
110     {
111       return Standard_False;
112     }
113
114     return Standard_True;
115   };
116
117   virtual Standard_Boolean Values (const math_Vector& theX,
118                                    math_Vector& theF,
119                                    math_Matrix& theD)
120   {
121     if(!Value(theX, theF))
122       return Standard_False;
123
124     if(!Derivatives(theX, theD))
125       return Standard_False;
126
127     return Standard_True;
128   }
129
130 protected:
131   FuncPreciseSeam operator=(FuncPreciseSeam&);
132
133 private:
134   const Handle(Adaptor3d_HSurface)& myQSurf;
135   const Handle(Adaptor3d_HSurface)& myPSurf;
136
137   //! 1 for U-coordinate, 0 - for V one. 
138   const Standard_Integer mySeamCoordInd;
139
140   //! Constant parameter of iso-line
141   const Standard_Real myIsoParameter;
142 };
143
144 //=======================================================================
145 //function : GetTangent
146 //purpose  : Computes tangent having the given parameter.
147 //            See calling method(s) for detailed information
148 //=======================================================================
149 static inline void GetTangent(const Standard_Real theConeSemiAngle,
150                               const Standard_Real theParameter,
151                               gp_XYZ& theResult)
152 {
153   const Standard_Real aW2 = theParameter*theParameter;
154   const Standard_Real aCosUn = (1.0 - aW2) / (1.0 + aW2);
155   const Standard_Real aSinUn = 2.0*theParameter / (1.0 + aW2);
156
157   const Standard_Real aTanA = Tan(theConeSemiAngle);
158   theResult.SetCoord(aTanA*aCosUn, aTanA*aSinUn, 1.0);
159 }
160
161 //=======================================================================
162 //function : IsPointOnSurface
163 //purpose  : Checks if thePt is in theSurf (with given tolerance).
164 //            Returns the foot of projection (theProjPt) and its parameters
165 //           on theSurf.
166 //=======================================================================
167 static Standard_Boolean IsPointOnSurface(const Handle(Adaptor3d_HSurface)& theSurf,
168                                          const gp_Pnt& thePt,
169                                          const Standard_Real theTol,
170                                          gp_Pnt& theProjPt,
171                                          Standard_Real& theUpar,
172                                          Standard_Real& theVpar)
173 {
174   Standard_Boolean aRetVal = Standard_False;
175
176   switch(theSurf->GetType())
177   {
178   case GeomAbs_Plane:
179   case GeomAbs_Cylinder:
180   case GeomAbs_Cone:
181   case GeomAbs_Sphere:
182   case GeomAbs_Torus:
183   case GeomAbs_SurfaceOfExtrusion:
184   case GeomAbs_SurfaceOfRevolution:
185     {
186       Extrema_ExtPS anExtr(thePt, theSurf->Surface(), theSurf->UResolution(theTol),
187                               theSurf->VResolution(theTol), Extrema_ExtFlag_MIN);
188       if(!anExtr.IsDone() || (anExtr.NbExt() < 1))
189       {
190         aRetVal = Standard_False;
191       }
192       else
193       {
194         Standard_Integer anExtrIndex = 1;
195         Standard_Real aSqDistMin = anExtr.SquareDistance(anExtrIndex);
196         for(Standard_Integer i = anExtrIndex + 1; i <= anExtr.NbExt(); i++)
197         {
198           const Standard_Real aSqD = anExtr.SquareDistance(i);
199           if(aSqD < aSqDistMin)
200           {
201             aSqDistMin = aSqD;
202             anExtrIndex = i;
203           }
204         }
205
206         if(aSqDistMin > theTol*theTol)
207         {
208           aRetVal = Standard_False;
209         }
210         else
211         {
212           theProjPt.SetXYZ(anExtr.Point(anExtrIndex).Value().XYZ());
213           anExtr.Point(anExtrIndex).Parameter(theUpar, theVpar);
214           aRetVal = Standard_True;
215         }
216       }
217     }
218     break;
219   default:
220     {
221       Extrema_GenLocateExtPS anExtr(theSurf->Surface());
222       anExtr.Perform(thePt, theUpar, theVpar);
223       if(!anExtr.IsDone() || (anExtr.SquareDistance() > theTol*theTol))
224       {
225         aRetVal = Standard_False;
226       }
227       else
228       {
229         anExtr.Point().Parameter(theUpar, theVpar);
230         theProjPt.SetXYZ(anExtr.Point().Value().XYZ());
231         aRetVal = Standard_True;
232       }
233     }
234     break;
235   }
236
237   return aRetVal;
238 }
239
240 //=======================================================================
241 //function : AddCrossUVIsoPoint
242 //purpose  : theQSurf is the surface possibly containing special point, 
243 //            thePSurf is another surface to intersect.
244 //=======================================================================
245 Standard_Boolean IntPatch_SpecialPoints::
246                       AddCrossUVIsoPoint(const Handle(Adaptor3d_HSurface)& theQSurf,
247                                          const Handle(Adaptor3d_HSurface)& thePSurf,
248                                          const IntSurf_PntOn2S& theRefPt,
249                                          const Standard_Real theTol,
250                                          IntSurf_PntOn2S& theAddedPoint,
251                                          const Standard_Boolean theIsReversed)
252 {
253   Standard_Real anArrOfPeriod[4] = {0.0, 0.0, 0.0, 0.0};
254   IntSurf::SetPeriod(theIsReversed ? thePSurf : theQSurf,
255                      theIsReversed ? theQSurf : thePSurf, anArrOfPeriod);
256
257   gp_Pnt aPQuad;
258
259   //Not quadric point
260   Standard_Real aU0 = 0.0, aV0 = 0.0;
261   if(theIsReversed)
262     theRefPt.ParametersOnS1(aU0, aV0);
263   else
264     theRefPt.ParametersOnS2(aU0, aV0);
265
266   //Quadric point
267   Standard_Real aUquad = 0.0, aVquad = 0.0; 
268
269   theQSurf->D0(aUquad, aVquad, aPQuad);
270
271   Extrema_GenLocateExtPS anExtr(thePSurf->Surface());
272   anExtr.Perform(aPQuad, aU0, aV0);
273
274   if(!anExtr.IsDone())
275   {
276     return Standard_False;
277   }
278
279   if(anExtr.SquareDistance() > theTol*theTol)
280   {
281     return Standard_False;
282   }
283
284   anExtr.Point().Parameter(aU0, aV0);
285   gp_Pnt aP0(anExtr.Point().Value());
286   
287   if(theIsReversed)
288     theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aU0, aV0, aUquad, aVquad);
289   else
290     theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0);
291
292   AdjustPointAndVertex(theRefPt, anArrOfPeriod, theAddedPoint);
293
294   return Standard_True;
295 }
296
297 //=======================================================================
298 //function : AddPointOnUorVIso
299 //purpose  : theQSurf is the surface possibly containing special point, 
300 //            thePSurf is another surface to intersect.
301 //=======================================================================
302 Standard_Boolean IntPatch_SpecialPoints::
303                       AddPointOnUorVIso(const Handle(Adaptor3d_HSurface)& theQSurf,
304                                         const Handle(Adaptor3d_HSurface)& thePSurf,
305                                         const IntSurf_PntOn2S& theRefPt,
306                                         const Standard_Boolean theIsU,
307                                         const Standard_Real theIsoParameter,
308                                         const math_Vector& theToler,
309                                         const math_Vector& theInitPoint,
310                                         const math_Vector& theInfBound,
311                                         const math_Vector& theSupBound,
312                                         IntSurf_PntOn2S& theAddedPoint,
313                                         const Standard_Boolean theIsReversed)
314 {
315   Standard_Real anArrOfPeriod[4] = {0.0, 0.0, 0.0, 0.0};
316   IntSurf::SetPeriod(theIsReversed ? thePSurf : theQSurf,
317                      theIsReversed ? theQSurf : thePSurf, anArrOfPeriod);
318
319   FuncPreciseSeam aF(theQSurf, thePSurf, theIsU, theIsoParameter);
320
321   math_FunctionSetRoot aSRF(aF, theToler);
322   aSRF.Perform(aF, theInitPoint, theInfBound, theSupBound);
323
324   if(!aSRF.IsDone())
325   {
326     return Standard_False;
327   }
328
329   math_Vector aRoots(theInitPoint.Lower(), theInitPoint.Upper());
330   aSRF.Root(aRoots);
331
332   //On parametric
333   Standard_Real aU0 = aRoots(1), aV0 = aRoots(2);
334
335   //On quadric
336   Standard_Real aUquad = theIsU ? 0.0 : aRoots(3);
337   Standard_Real aVquad = theIsU ? aRoots(3) : 0.0; 
338   const gp_Pnt aPQuad(theQSurf->Value(aUquad, aVquad));
339   const gp_Pnt aP0(thePSurf->Value(aU0, aV0));
340
341   if(theIsReversed)
342     theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aU0, aV0, aUquad, aVquad);
343   else
344     theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0);
345
346   AdjustPointAndVertex(theRefPt, anArrOfPeriod, theAddedPoint);
347   return Standard_True;
348 }
349
350 //=======================================================================
351 //function : ProcessSphere
352 //purpose  :
353 /*
354 The intersection point (including the pole)
355 must be satisfied to the following system:
356
357     \left\{\begin{matrix}
358     R*\cos (U_{q})*\cos (V_{q})=S_{x}(U_{s},V_{s})
359     R*\sin (U_{q})*\cos (V_{q})=S_{y}(U_{s},V_{s})
360     R*\sin (V_{q})=S_{z}(U_{s},V_{s})
361     \end{matrix}\right,
362 where
363   R is the radius of the sphere;
364   @S_{x}@, @S_{y}@ and @S_{z}@ are X, Y and Z-coordinates of thePSurf;
365   @U_{s}@ and @V_{s}@ are parameters on the parametric surface;
366   @U_{q}@ and @V_{q}@ are equal to theUquad and theVquad correspondingly.
367
368 Consequently (from first two equations),
369   \left\{\begin{matrix}
370   \cos (U_{q}) = \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})}
371   \sin (U_{q}) = \frac{S_{y}(U_{s},V_{s})}{R*\cos (V_{q})}
372   \end{matrix}\right.
373
374 For pole,
375   V_{q}=\pm \pi /2 \Rightarrow \cos (V_{q}) = 0 (denominator is equal to 0).
376
377 Therefore, computation U_{q} directly is impossibly.
378
379 Let @V_{q}@ tends to @\pm \pi /2@.
380 Then (indeterminate form is evaluated in accordance of L'Hospital rule),
381   \cos (U_{q}) = \lim_{V_{q} \to (\pi /2-0)}
382   \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})}=
383   -\lim_{V_{q} \to (\pi /2-0)}
384   \frac{\frac{\partial S_{x}}
385   {\partial U_{s}}*\frac{\mathrm{d} U_{s}}
386   {\mathrm{d} V_{q}}+\frac{\partial S_{x}}
387   {\partial V_{s}}*\frac{\mathrm{d} V_{s}}
388   {\mathrm{d} V_{q}}}{R*\sin (V_{q})} =
389   -\frac{1}{R}*\frac{\mathrm{d} U_{s}}
390   {\mathrm{d} V_{q}}*(\frac{\partial S_{x}}
391   {\partial U_{s}}+\frac{\partial S_{x}}
392   {\partial V_{s}}*\frac{\mathrm{d} V_{s}}
393   {\mathrm{d} U_{s}}) =
394   -\frac{1}{R}*\frac{\mathrm{d} V_{s}}
395   {\mathrm{d} V_{q}}*(\frac{\partial S_{x}}
396   {\partial U_{s}}*\frac{\mathrm{d} U_{s}}
397   {\mathrm{d} V_{s}}+\frac{\partial S_{x}}
398   {\partial V_{s}}).
399
400 Analogicaly for @\sin (U_{q})@ (@S_{x}@ is substituted to @S_{y}@).
401
402 Let mean, that
403   \cos (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \cos (U_{q}) \left | _{V_{q} \to (\pi /2-0)}
404   \sin (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \sin (U_{q}) \left | _{V_{q} \to (\pi /2-0)}
405
406 From the 3rd equation of the system, we obtain
407   \frac{\mathrm{d} (R*\sin (V_{q}))}{\mathrm{d} V_{q}} =
408   \frac{\mathrm{d} S_{z}(U_{s},V_{s})}{\mathrm{d} V_{q}}
409 or
410   R*\cos (V_{q}) = \frac{\partial S_{z}}{\partial U_{s}}*
411   \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}}
412   {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}}.
413
414 If @V_{q}=\pm \pi /2@, then
415   \frac{\partial S_{z}}{\partial U_{s}}*
416   \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}}
417   {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}} = 0.
418
419 Consequently, if @\frac{\partial S_{z}}{\partial U_{s}} \neq 0 @ then
420   \frac{\mathrm{d} U_{s}}{\mathrm{d} V_{s}} =
421   -\frac{\frac{\partial S_{z}}{\partial V_{s}}}
422   {\frac{\partial S_{z}}{\partial U_{s}}}.
423
424 If @ \frac{\partial S_{z}}{\partial V_{s}} \neq 0 @ then
425   \frac{\mathrm{d} V_{s}}{\mathrm{d} U_{s}} =
426   -\frac{\frac{\partial S_{z}}{\partial U_{s}}}
427   {\frac{\partial S_{z}}{\partial V_{s}}}
428
429 Cases, when @ \frac{\partial S_{z}}{\partial U_{s}} =
430 \frac{\partial S_{z}}{\partial V_{s}} = 0 @ are not consider here.
431 The reason is written below.
432 */
433 //=======================================================================
434 Standard_Boolean IntPatch_SpecialPoints::ProcessSphere(const IntSurf_PntOn2S& thePtIso,
435                                                        const gp_Vec& theDUofPSurf,
436                                                        const gp_Vec& theDVofPSurf,
437                                                        const Standard_Boolean theIsReversed,
438                                                        const Standard_Real theVquad,
439                                                        Standard_Real& theUquad,
440                                                        Standard_Boolean& theIsIsoChoosen)
441 {
442   theIsIsoChoosen = Standard_False;
443
444   //Vector with {@ \cos (U_{q}) @, @ \sin (U_{q}) @} coordinates.
445   //Ask to pay attention to the fact that this vector is always normalized.
446   gp_Vec2d aV1;
447
448   if ((Abs(theDUofPSurf.Z()) < Precision::PConfusion()) &&
449       (Abs(theDVofPSurf.Z()) < Precision::PConfusion()))
450   {
451     //Example of this case is an intersection of a plane with a sphere
452     //when the plane tangents the sphere in some pole (i.e. only one 
453     //intersection point, not line). In this case, U-coordinate of the
454     //sphere is undefined (can be realy anything).
455     //Another reason is that we have tangent zone around the pole
456     //(see bug #26576).
457     //Computation of correct value of theUquad is impossible.
458     //Therefore, (in order to return something) we will consider
459     //the intersection line goes along some isoline in neighborhood
460     //of the pole.
461
462 #ifdef INTPATCH_ADDSPECIALPOINTS_DEBUG
463     cout << "Cannot find UV-coordinate for quadric in the pole."
464       " See considered comment above. IntPatch_SpecialPoints.cxx,"
465       " ProcessSphere(...)" << endl;
466 #endif
467     Standard_Real aUIso = 0.0, aVIso = 0.0;
468     if (theIsReversed)
469       thePtIso.ParametersOnS2(aUIso, aVIso);
470     else
471       thePtIso.ParametersOnS1(aUIso, aVIso);
472
473     theUquad = aUIso;
474     theIsIsoChoosen = Standard_True;
475   }
476   else
477   {
478     if (Abs(theDUofPSurf.Z()) > Abs(theDVofPSurf.Z()))
479     {
480       const Standard_Real aDusDvs = theDVofPSurf.Z() / theDUofPSurf.Z();
481       aV1.SetCoord(theDUofPSurf.X()*aDusDvs - theDVofPSurf.X(),
482                    theDUofPSurf.Y()*aDusDvs - theDVofPSurf.Y());
483     }
484     else
485     {
486       const Standard_Real aDvsDus = theDUofPSurf.Z() / theDVofPSurf.Z();
487       aV1.SetCoord(theDVofPSurf.X()*aDvsDus - theDUofPSurf.X(),
488                    theDVofPSurf.Y()*aDvsDus - theDUofPSurf.Y());
489     }
490
491     aV1.Normalize();
492
493     if (Abs(aV1.X()) > Abs(aV1.Y()))
494       theUquad = Sign(asin(aV1.Y()), theVquad);
495     else
496       theUquad = Sign(acos(aV1.X()), theVquad);
497   }
498
499   return Standard_True;
500 }
501
502 //=======================================================================
503 //function : ProcessCone
504 //purpose  : 
505 /*
506 The intersection point (including the pole)
507 must be satisfied to the following system:
508
509   \left\{\begin {matrix}
510   (V_{q}\sin(a) + R)*\cos(U_{q})) = S_{x}(U_{s}, V_{s})\\
511   (V_{q}\sin(a) + R)*\sin(U_{q})) = S_{y}(U_{s}, V_{s})\\
512   V_{q}\cos(a) = S_{z}(U_{s}, V_{s})
513   \end {matrix}\right,
514 where 
515   R is the radius of the cone;
516   a is its semi-angle;
517   @S_{x}@, @S_{y}@ and @S_{z}@ are X, Y and Z-coordinates of thePSurf;
518   @U_{s}@ and @V_{s}@ are parameters on the parametric surface;
519   @U_{q}@ and @V_{q}@ are equal to theUquad and theVquad correspondingly.
520
521 Consequently (from first two equations), 
522   \left\{\begin{matrix}
523   \cos(U_{q})=\frac{S_{x}(U_{s},V_{s})}{(V_{q}\sin(a)+R)}\\
524   \sin(U_{q})=\frac{S_{y}(U_{s}, V_{s})}{(V_{q}\sin(a)+R)}
525   \end{matrix}\right.
526
527 For pole, the denominator of these two equations is equal to 0.
528 Therefore, computation U_{q} directly is impossibly.
529
530 Let @V_{q}@ tends to @\frac{-R}{\sin(a)})@.
531 Then (indeterminate form is evaluated in accordance of L'Hospital rule),
532
533   \cos (U_{q}) = 
534   \lim_{V_{q} \to (\frac{-R}{\sin(a)})}\frac{S_{x}(U_{s},V_{s})}{(V_{q}\sin(a)+R)}=
535   \frac{1}{\sin(a)}* \lim_{V_{q} \to (\frac{-R}{\sin(a)})}\frac{dU_{s}}{dV_{q}}*
536   (\frac{\partial S_{x}}{\partial U_{s}}+\frac{\partial S_{x}}{\partial V_{s}}*
537   \frac{dV_{s}}{dU_{s}})=
538   \frac{1}{\sin(a)}* \lim_{V_{q} \to (\frac{-R}{\sin(a)})}\frac{dV_{s}}{dV_{q}}*
539   (\frac{\partial S_{x}}{\partial U_{s}}*
540   \frac{dU_{s}}{dV_{s}}+\frac{\partial S_{x}}{\partial V_{s}})
541
542 Analogically for @\sin (U_{q})@ (@S_{x}@ is substituted to @S_{y}@).
543
544 After differentiating 3rd equation of the system, we will obtain
545   \cos(a)=\frac{dS_{z}}{dV_{q}}=\frac{dU_{s}}{dV_{q}}*
546   (\frac{\partial S_{z}}{\partial U_{s}}+\frac{\partial S_{z}}{\partial V_{s}}*
547   \frac{dV_{s}}{dU_{s}})
548 or
549   \frac{dU_{s}}{dV_{q}}=\frac{\cos(a)}{\frac{\partial S_{z}}{\partial U_{s}}+
550   \frac{\partial S_{z}}{\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}
551
552 After substituting we will obtain
553   \cos (U_{q}) = 
554   \cot(a)*\frac{\frac{\partial S_{x}}{\partial U_{s}}+\frac{\partial S_{x}}
555   {\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}{\frac{\partial S_{z}}
556   {\partial U_{s}}+\frac{\partial S_{z}}{\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}
557
558   \sin (U_{q}) =
559   \cot(a)*\frac{\frac{\partial S_{y}}{\partial U_{s}}+\frac{\partial S_{y}}
560   {\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}{\frac{\partial S_{z}}
561   {\partial U_{s}}+\frac{\partial S_{z}}{\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}
562
563 So, we have obtained vector with coordinates {@ \cos (U_{q}) @, @ \sin (U_{q}) @}.
564 Ask to pay attention to the fact that this vector is always normalized.
565 And after normalization this vector will have coordinates
566   {\cos (U_{q}), \sin (U_{q})} = {dS_{x}, dS_{y}}.Normalized().
567
568 It means that we have to compute a tangent to the intersection curve in
569 the cone apex point. After that, just take its X- and Y-coordinates.
570
571 However, we have to compute derivative @\frac{dV_{s}}{dU_{s}}@ in order
572 to compute this vector. In order to find this derivative, we use the
573 information about direction of tangent to the intersection curve.
574 This tangent will be directed along the cone generatrix obtained by intersection
575 of the cone with a plane tangent to 2nd (intersected) surface.
576 */
577 //=======================================================================
578 Standard_Boolean IntPatch_SpecialPoints::ProcessCone(const IntSurf_PntOn2S& thePtIso,
579                                                      const gp_Vec& theDUofPSurf,
580                                                      const gp_Vec& theDVofPSurf,
581                                                      const gp_Cone& theCone,
582                                                      const Standard_Boolean theIsReversed,
583                                                      Standard_Real& theUquad,
584                                                      Standard_Boolean& theIsIsoChoosen)
585 {
586   theIsIsoChoosen = Standard_False;
587
588   // A plane tangent to 2nd (intersected) surface.
589   // Its normal.
590   const gp_XYZ aTgPlaneZ(theDUofPSurf.Crossed(theDVofPSurf).XYZ());
591   const Standard_Real aSqModTg = aTgPlaneZ.SquareModulus();
592   if (aSqModTg < Precision::SquareConfusion())
593   {
594     theIsIsoChoosen = Standard_True;
595   }
596
597   gp_XYZ aTgILine[2];
598   const Standard_Integer aNbTangent = !theIsIsoChoosen? 
599                           GetTangentToIntLineForCone(theCone.SemiAngle(),
600                                                      aTgPlaneZ.Divided(Sqrt(aSqModTg)),
601                                                      aTgILine) : 0;
602
603   if (aNbTangent == 0)
604   {
605     theIsIsoChoosen = Standard_True;
606   }
607   else
608   {
609     const Standard_Real aPeriod = M_PI + M_PI;
610     Standard_Real aUIso = 0.0, aVIso = 0.0;
611     if (theIsReversed)
612       thePtIso.ParametersOnS2(aUIso, aVIso);
613     else
614       thePtIso.ParametersOnS1(aUIso, aVIso);
615
616     aUIso = ElCLib::InPeriod(aUIso, 0.0, aPeriod);
617
618     // Sought U-parameter in the apex point
619
620     // For 2 possible parameters value,
621     // one will be chosen which is nearer
622     // to aUIso. Following variables will help to chose.
623     Standard_Real aMinDelta = RealLast();
624     for (Standard_Integer anIdx = 0; anIdx < aNbTangent; anIdx++)
625     {
626       // Vector {@\cos(a), \sin(a)@}
627       gp_Vec2d aVecCS(aTgILine[anIdx].X(), aTgILine[anIdx].Y());
628       const Standard_Real aSqMod = aVecCS.SquareMagnitude();
629       if (aSqMod < Precision::SquareConfusion())
630       {
631         theIsIsoChoosen = Standard_True;
632         break;
633       }
634
635       // Normalize
636       aVecCS.Divide(Sqrt(aSqMod));
637
638       // Angle in range [0, PI/2]
639       Standard_Real anUq = (Abs(aVecCS.X()) < Abs(aVecCS.Y())) ? ACos(Abs(aVecCS.X())) : ASin(Abs(aVecCS.Y()));
640
641       // Convert angles to the range [0, 2*PI]
642       if (aVecCS.Y() < 0.0)
643       {
644         if (aVecCS.X() > 0.0)
645         {
646           anUq = -anUq;
647         }
648         else
649         {
650           anUq += M_PI;
651         }
652       }
653       else if (aVecCS.X() < 0.0)
654       {
655         anUq = M_PI - anUq;
656       }
657
658       //Select the parameter the nearest to aUIso
659       anUq = ElCLib::InPeriod(anUq, 0.0, aPeriod);
660       Standard_Real aDelta = Abs(anUq - aUIso);
661       if (aDelta > M_PI)
662         aDelta = aPeriod - aDelta;
663
664       if (aDelta < aMinDelta)
665       {
666         aMinDelta = aDelta;
667         theUquad = anUq;
668       }
669     }
670   }
671
672   if (theIsIsoChoosen)
673   {
674 #ifdef INTPATCH_ADDSPECIALPOINTS_DEBUG
675     cout << "Cannot find UV-coordinate for quadric in the pole."
676       " IntPatch_AddSpecialPoints.cxx, ProcessCone(...)" << endl;
677 #endif
678     theIsIsoChoosen = Standard_True;
679
680     Standard_Real aUIso = 0.0, aVIso = 0.0;
681     if (theIsReversed)
682       thePtIso.ParametersOnS2(aUIso, aVIso);
683     else
684       thePtIso.ParametersOnS1(aUIso, aVIso);
685
686     theUquad = aUIso;
687     return Standard_True;
688   }
689   else
690   {
691     return Standard_True;
692   }
693
694   //return Standard_False;
695 }
696
697 //=======================================================================
698 //function : GetTangentToIntLineForCone
699 //purpose  : The following conditions must be satisfied:
700 //1. The cone is represented in its canonical form.
701 //2. The plane goes through the cone apex and has the normal vector thePlnNormal.
702 //3. Vector thePlnNormal has already been normalized
703 /*
704 Let us enter the new coordinate system where the origin will be in the cone apex
705 and axes are the same as in World-Coordinate-System (WCS).
706 There the horizontal plane (which is parallel XY-plane) with the origin
707 (0, 0, 1) will intersect the cone by the circle with center (0, 0, 1),
708 direction {0, 0, 1} and radius tg(a) where a is the cone semi-angle.
709 Its equation will be
710 \left\{\begin{matrix}
711 x(U_{n}) = \tan(a)*\cos(U_{n}) = \tan(a)*\frac{1-\tan^{2}(U_{n}/2)}{1+\tan^{2}(U_{n}/2)}\\
712 y(U_{n}) = \tan(a)*\sin (U_{n}) = \tan(a)*\frac{2*\tan(U_{n}/2)}{1+\tan^{2}(U_{n}/2)}\\
713 z(U_{n}) = 1
714 \end{matrix}\right.
715
716 The given plane has (in this coordinate system) location (0, 0, 0) and
717 the same normal thePlnNormal=={nx,ny,nz}. Its equation is:
718 nx*x+ny*y+nz*z==0
719
720 After substitution circle's equation to the plane's equation
721 we will obtain a quadratic equation
722 aA*w^2 + 2*aB*w + aC = 0.
723 */
724 //=======================================================================
725 Standard_Integer IntPatch_SpecialPoints::GetTangentToIntLineForCone(const Standard_Real theConeSemiAngle,
726                                                                     const gp_XYZ& thePlnNormal,
727                                                                     gp_XYZ theResult[2])
728 {
729   const Standard_Real aNullTol = Epsilon(1.0);
730   const Standard_Real aTanA = Tan(theConeSemiAngle);
731   const Standard_Real aA = thePlnNormal.Z() / aTanA - thePlnNormal.X();
732   const Standard_Real aB = thePlnNormal.Y();
733   const Standard_Real aC = thePlnNormal.Z() / aTanA + thePlnNormal.X();
734
735   if (Abs(aA) < aNullTol)
736   {
737     if (Abs(aB) > aNullTol)
738     {
739       //The plane goes along the cone generatrix.
740       GetTangent(theConeSemiAngle, -aC / (aB + aB), theResult[0]);
741       return 1;
742     }
743
744     //The cone and the plane have only one common point.
745     //It is the cone apex.
746     return 0;
747   }
748
749   //Discriminant of this equation is equal to 
750   Standard_Real aDiscr = thePlnNormal.Z() / Sin(theConeSemiAngle);
751   aDiscr = 1.0 - aDiscr*aDiscr;
752
753   if (Abs(aDiscr) < aNullTol)
754   {
755     //The plane goes along the cone generatrix.
756     // Attention! Mathematically, this cond. is equivalent to 
757     // above processed one (Abs(aA) < aNullTol && (Abs(aB) > aNullTol)).
758     // However, we separate this branch in order to eliminate numerical
759     // instability.
760
761     GetTangent(theConeSemiAngle, -aB / aA, theResult[0]);
762     return 1;
763   }
764   else if (aDiscr > 0.0)
765   {
766     const Standard_Real aRD = Sqrt(aDiscr);
767     GetTangent(theConeSemiAngle, (-aB+aRD)/aA, theResult[0]);
768     GetTangent(theConeSemiAngle, (-aB-aRD)/aA, theResult[1]);
769     return 2;
770   }
771
772   // We will never come here.
773   return 0;
774 }
775
776 //=======================================================================
777 //function : AddSingularPole
778 //purpose  : theQSurf is the surface possibly containing special point, 
779 //            thePSurf is another surface to intersect.
780 //           Returns TRUE, if the pole is an intersection point.
781 //=======================================================================
782 Standard_Boolean IntPatch_SpecialPoints::
783                       AddSingularPole(const Handle(Adaptor3d_HSurface)& theQSurf,
784                                       const Handle(Adaptor3d_HSurface)& thePSurf,
785                                       const IntSurf_PntOn2S& thePtIso,
786                                       IntPatch_Point& theVertex,
787                                       IntSurf_PntOn2S& theAddedPoint,
788                                       const Standard_Boolean theIsReversed,
789                                       const Standard_Boolean theIsReqRefCheck)
790 {
791   //On parametric
792   Standard_Real aU0 = 0.0, aV0 = 0.0;
793   //aPQuad is Pole
794   gp_Pnt aPQuad, aP0;
795   Standard_Real aUquad = 0.0, aVquad = 0.0;
796   if(theIsReversed)
797     theVertex.Parameters(aU0, aV0, aUquad, aVquad);
798   else
799     theVertex.Parameters(aUquad, aVquad, aU0, aV0);
800
801   aUquad = 0.0;
802
803   if(theQSurf->GetType() == GeomAbs_Sphere)
804   {
805     aVquad = Sign(M_PI_2, aVquad);
806   }
807   else if(theQSurf->GetType() == GeomAbs_Cone)
808   {
809     const gp_Cone aCo = theQSurf->Cone();
810     const Standard_Real aRadius = aCo.RefRadius();
811     const Standard_Real aSemiAngle = aCo.SemiAngle();
812     aVquad = -aRadius / sin(aSemiAngle);
813   }
814   else
815   {
816     throw Standard_TypeMismatch( "IntPatch_SpecialPoints::AddSingularPole(),"
817                                   "Unsupported quadric with Pole");
818   }
819
820   theQSurf->D0(aUquad, aVquad, aPQuad);
821   const Standard_Real aTol = theVertex.Tolerance();
822   if (theIsReqRefCheck && (aPQuad.SquareDistance(theVertex.Value()) >= aTol*aTol))
823   {
824     return Standard_False;
825   }
826
827   if (!IsPointOnSurface(thePSurf, aPQuad, aTol, aP0, aU0, aV0))
828   {
829     return Standard_False;
830   }
831
832   //Pole is an intersection point
833   //(lies in the quadric and the parametric surface)
834
835   if(theIsReversed)
836     theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aU0, aV0, aUquad, aVquad);
837   else
838     theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0);
839
840   const Standard_Boolean isSame = theAddedPoint.IsSame(theVertex.PntOn2S(),
841                                                        Precision::Confusion());
842
843   //Found pole does not exist in the Walking-line
844   //It must be added there (with correct 2D-parameters)
845
846   //2D-parameters of thePSurf surface have already been found (aU0, aV0).
847   //Let find 2D-parameters on the quadric.
848
849   //The algorithm depends on the type of the quadric.
850   //Here we consider a Sphere and cone only.
851
852   //First of all, we need in adjusting thePSurf in the coordinate system of the Sphere/Cone
853   //(in order to make its equation maximal simple). However, as it will be
854   //shown later, thePSurf is used in algorithm in order to get its derivatives.
855   //Therefore, for improving performance, transformation of these vectors is enough
856   //(there is no point in transformation of full surface).
857   
858   gp_Pnt aPtemp;
859   gp_Vec aVecDu, aVecDv;
860   thePSurf->D1(aU0, aV0, aPtemp, aVecDu, aVecDv);
861
862   //Transforms parametric surface in coordinate-system of the quadric
863   gp_Trsf aTr;
864   aTr.SetTransformation((theQSurf->GetType() == GeomAbs_Sphere) ?
865                           theQSurf->Sphere().Position() :
866                                       theQSurf->Cone().Position());
867
868   //Derivatives of transformed thePSurf
869   aVecDu.Transform(aTr);
870   aVecDv.Transform(aTr);
871
872   Standard_Boolean isIsoChoosen = Standard_False;
873
874   if(theQSurf->GetType() == GeomAbs_Sphere)
875   {
876     if (!ProcessSphere(thePtIso, aVecDu, aVecDv, theIsReversed, 
877                         aVquad, aUquad, isIsoChoosen))
878     {
879       return Standard_False;
880     }
881   }
882   else //if(theQSurf->GetType() == GeomAbs_Cone)
883   {
884     if (!ProcessCone(thePtIso, aVecDu, aVecDv, theQSurf->Cone(),
885                       theIsReversed, aUquad, isIsoChoosen))
886     {
887       return Standard_False;
888     }
889   }
890
891   if(theIsReversed)
892     theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aU0, aV0, aUquad, aVquad);
893   else
894     theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0);
895
896   if (isSame)
897   {
898     theVertex.SetValue(theAddedPoint);
899     return Standard_True;
900   }
901
902   if (!isIsoChoosen)
903   {
904     Standard_Real anArrOfPeriod[4];
905     if (theIsReversed)
906     {
907       IntSurf::SetPeriod(thePSurf, theQSurf, anArrOfPeriod);
908     }
909     else
910     {
911       IntSurf::SetPeriod(theQSurf, thePSurf, anArrOfPeriod);
912     }
913
914     AdjustPointAndVertex(theVertex.PntOn2S(), anArrOfPeriod, theAddedPoint);
915   }
916   else
917   {
918     theVertex.SetValue(theAddedPoint);
919   }
920
921   return Standard_True;
922 }
923
924 //=======================================================================
925 //function : ContinueAfterSpecialPoint
926 //purpose  : If the last point of the line is the pole of the quadric then
927 //            the Walking-line has been broken in this point.
928 //           However, new line must start from this point. Here we must
929 //            find 2D-coordinates of "this new" point.
930 /*
931 The inters. line in the neighborhood of the Apex/Pole(s) can be
932 approximated by the intersection result of the Cone/Sphere with
933 the plane going through the Apex/Pole and being tangent to the
934 2nd intersected surface. This intersection result is well known.
935
936 In case of sphere, the inters. result is a circle.
937 If we go along this circle and across the Pole then U-parameter of
938 the sphere (@U_{q}@) will change to +/-PI.
939
940 In case of cone, the inters. result is two intersected lines (which
941 can be merged to one in a special case when the plane goes along
942 some generatrix of the cone). The direction of these lines
943 are computed by GetTangentToIntLineForCone(...) method).
944
945 When the real (not lines) inters. curve goes through the cone apex then
946 two variants are possible:
947 a) The tangent line to the inters. curve will be left. In this case
948 U-parameter of the cone (@U_{q}@) will be change to +/-PI.
949 b) Another line (as inters. result of cone + plane) will tangent
950 to the inters. curve. In this case @U_{q}@ must be recomputed.
951 */
952 //=======================================================================
953 Standard_Boolean IntPatch_SpecialPoints::
954                   ContinueAfterSpecialPoint(const Handle(Adaptor3d_HSurface)& theQSurf,
955                                             const Handle(Adaptor3d_HSurface)& thePSurf,
956                                             const IntSurf_PntOn2S& theRefPt,
957                                             const IntPatch_SpecPntType theSPType,
958                                             const Standard_Real theTol2D,
959                                             IntSurf_PntOn2S& theNewPoint,
960                                             const Standard_Boolean theIsReversed)
961 {
962   if(theSPType == IntPatch_SPntNone)
963     return Standard_False;
964
965   if(theNewPoint.IsSame(theRefPt, Precision::Confusion(), theTol2D))
966   {
967     return Standard_False;
968   }
969
970   if ((theSPType == IntPatch_SPntPole) && (theQSurf->GetType() == GeomAbs_Cone))
971   {
972     //Check if the condition b) is satisfied.
973     //Repeat the same steps as in 
974     //IntPatch_SpecialPoints::AddSingularPole(...) method.
975
976     //On parametric
977     Standard_Real aU0 = 0.0, aV0 = 0.0;
978     //On quadric
979     Standard_Real aUquad = 0.0, aVquad = 0.0;
980
981     if (theIsReversed)
982       theNewPoint.Parameters(aU0, aV0, aUquad, aVquad);
983     else
984       theNewPoint.Parameters(aUquad, aVquad, aU0, aV0);
985
986     gp_Pnt aPtemp;
987     gp_Vec aVecDu, aVecDv;
988     thePSurf->D1(aU0, aV0, aPtemp, aVecDu, aVecDv);
989
990     //Transforms parametric surface in coordinate-system of the quadric
991     gp_Trsf aTr;
992     aTr.SetTransformation(theQSurf->Cone().Position());
993
994     //Derivatives of transformed thePSurf
995     aVecDu.Transform(aTr);
996     aVecDv.Transform(aTr);
997
998     Standard_Boolean isIsoChoosen = Standard_False;
999     ProcessCone(theRefPt, aVecDu, aVecDv, theQSurf->Cone(),
1000                 theIsReversed, aUquad, isIsoChoosen);
1001
1002     theNewPoint.SetValue(!theIsReversed, aUquad, aVquad);
1003   }
1004
1005   //As it has already been said, in case of going through the Pole/Apex,
1006   //U-parameter of the quadric surface will change to +/-PI. This rule has some
1007   //exceptions:
1008   //1. When 2nd surface has C0-continuity in the point common with the Apex/Pole.
1009   //    In this case, the tangent line to the intersection curve after the Apex/Pole
1010   //    must be totally recomputed according to the new derivatives of the 2nd surface.
1011   //    Currently, it is not implemented but will be able to be done after the
1012   //    corresponding demand.
1013   //2. The inters. curve has C1 continuity but huge curvature in the point common with 
1014   //    the Apex/Pole. Existing inters. algorithm does not allow putting many points
1015   //    near to the Apex/Pole in order to cover this "sharp" piece of the inters. curve.
1016   //    Therefore, we use adjusting U-parameter of the quadric surface with
1017   //    period PI/2 instead of 2PI. It does not have any mathematical idea
1018   //    but allows creating WLine with more or less uniform distributed points.
1019   //    In other words, we forbid "jumping" between two neighbor Walking-points
1020   //    with step greater than PI/4.
1021
1022   const Standard_Real aPeriod = (theSPType == IntPatch_SPntPole)? M_PI_2 : 2.0*M_PI;
1023
1024   const Standard_Real aUpPeriod = thePSurf->IsUPeriodic() ? thePSurf->UPeriod() : 0.0;
1025   const Standard_Real aUqPeriod = theQSurf->IsUPeriodic() ? aPeriod : 0.0;
1026   const Standard_Real aVpPeriod = thePSurf->IsVPeriodic() ? thePSurf->VPeriod() : 0.0;
1027   const Standard_Real aVqPeriod = theQSurf->IsVPeriodic() ? aPeriod : 0.0;
1028
1029   const Standard_Real anArrOfPeriod[4] = {theIsReversed? aUpPeriod : aUqPeriod,
1030                                           theIsReversed? aVpPeriod : aVqPeriod,
1031                                           theIsReversed? aUqPeriod : aUpPeriod,
1032                                           theIsReversed? aVqPeriod : aVpPeriod};
1033
1034   AdjustPointAndVertex(theRefPt, anArrOfPeriod, theNewPoint);
1035   return Standard_True;
1036 }
1037
1038 //=======================================================================
1039 //function : AdjustPointAndVertex
1040 //purpose  : 
1041 //=======================================================================
1042 void IntPatch_SpecialPoints::
1043                 AdjustPointAndVertex(const IntSurf_PntOn2S &theRefPoint,
1044                                      const Standard_Real theArrPeriods[4],
1045                                      IntSurf_PntOn2S &theNewPoint,
1046                                      IntPatch_Point* const theVertex)
1047 {
1048   Standard_Real aRefPar[2] = {0.0, 0.0};
1049   Standard_Real aPar[4] = {0.0, 0.0, 0.0, 0.0};
1050   theNewPoint.Parameters(aPar[0], aPar[1], aPar[2], aPar[3]);
1051
1052   for(Standard_Integer i = 0; i < 4; i++)
1053   {
1054     if(theArrPeriods[i] == 0)
1055       continue;
1056
1057     const Standard_Real aPeriod = theArrPeriods[i], aHalfPeriod = 0.5*theArrPeriods[i];
1058
1059     if(i < 2)
1060     {// 1st surface is used
1061       theRefPoint.ParametersOnS1(aRefPar[0], aRefPar[1]);
1062     }
1063     else
1064     {
1065       theRefPoint.ParametersOnS2(aRefPar[0], aRefPar[1]);
1066     }
1067
1068     const Standard_Integer aRefInd = i%2;
1069
1070     {
1071       Standard_Real aDeltaPar = aRefPar[aRefInd]-aPar[i];
1072       const Standard_Real anIncr = Sign(aPeriod, aDeltaPar);
1073       while((aDeltaPar > aHalfPeriod) || (aDeltaPar < -aHalfPeriod))
1074       {
1075         aPar[i] += anIncr;
1076         aDeltaPar = aRefPar[aRefInd]-aPar[i];
1077       }
1078     }
1079   }
1080
1081   if(theVertex)
1082     (*theVertex).SetParameters(aPar[0], aPar[1], aPar[2], aPar[3]);
1083
1084   theNewPoint.SetValue(aPar[0], aPar[1], aPar[2], aPar[3]);
1085 }
1086