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_Surface.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 definitely.
  34 class FuncPreciseSeam: public math_FunctionSetWithDerivatives
  35 {
  36 public:
  37   FuncPreciseSeam(const Handle(Adaptor3d_Surface)& theQSurf, // quadric
  38                   const Handle(Adaptor3d_Surface)& 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 const&)
  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 const&)
 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_Surface)& myQSurf;
 135   const Handle(Adaptor3d_Surface)& 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_Surface)& 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, 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);
 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_Surface)& theQSurf,
 247                                          const Handle(Adaptor3d_Surface)& 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);
 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_Surface)& theQSurf,
 304                                         const Handle(Adaptor3d_Surface)& 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 Analogically 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 really 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     std::cout << "Cannot find UV-coordinate for quadric in the pole."
 464       " See considered comment above. IntPatch_SpecialPoints.cxx,"
 465       " ProcessSphere(...)" << std::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     std::cout << "Cannot find UV-coordinate for quadric in the pole."
 676       " IntPatch_AddSpecialPoints.cxx, ProcessCone(...)" << std::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_Surface)& theQSurf,
 784                                       const Handle(Adaptor3d_Surface)& 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_Surface)& theQSurf,
 955                                             const Handle(Adaptor3d_Surface)& 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