b180d9abedd5aacee58977db60a28d560030c379
[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 <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   {
377     throw Standard_TypeMismatch( "IntPatch_SpecialPoints::AddSingularPole(),"
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