src/IntAna/IntAna_Curve.cxx

   1 // Created on: 1992-06-30
   2 // Created by: Laurent BUCHARD
   3 // Copyright (c) 1992-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #ifndef OCCT_DEBUG
  18 #define No_Standard_RangeError
  19 #define No_Standard_OutOfRange
  20 #endif
  21
  22 //----------------------------------------------------------------------
  23 //-- Differents constructeurs sont proposes qui correspondent aux
  24 //-- polynomes en Z :
  25 //--    A(Sin(Theta),Cos(Theta)) Z**2
  26 //--  + B(Sin(Theta),Cos(Theta)) Z
  27 //--  + C(Sin(Theta),Cos(Theta))
  28 //--
  29 //-- Une Courbe est definie sur un domaine
  30 //--
  31 //-- Value retourne le point de parametre U(Theta),V(Theta)
  32 //--       ou V est la solution du polynome A V**2 + B V + C
  33 //--       (Selon les cas, on prend V+ ou V-)
  34 //--
  35 //-- D1u   calcule le vecteur tangent a la courbe
  36 //--       et retourne le booleen Standard_False si ce calcul ne peut
  37 //--       pas etre mene a bien.
  38 //----------------------------------------------------------------------
  39
  40 #include <algorithm>
  41
  42 #include <ElSLib.hxx>
  43 #include <gp_Cone.hxx>
  44 #include <gp_Cylinder.hxx>
  45 #include <gp_Pnt.hxx>
  46 #include <gp_Vec.hxx>
  47 #include <gp_XYZ.hxx>
  48 #include <IntAna_Curve.hxx>
  49 #include <math_DirectPolynomialRoots.hxx>
  50 #include <Precision.hxx>
  51 #include <Standard_DomainError.hxx>
  52
  53 //=======================================================================
  54 //function : IntAna_Curve
  55 //purpose  :
  56 //=======================================================================
  57 IntAna_Curve::IntAna_Curve()
  58 {
  59   typequadric=GeomAbs_OtherSurface;
  60   firstbounded=Standard_False;
  61   lastbounded=Standard_False;
  62 }
  63 //=======================================================================
  64 //function : SetConeQuadValues
  65 //purpose  : Description de l intersection Cone Quadrique
  66 //=======================================================================
  67 void IntAna_Curve::SetConeQuadValues(const gp_Cone& Cone,
  68                                      const Standard_Real Qxx,
  69                                      const Standard_Real Qyy,
  70                                      const Standard_Real Qzz,
  71                                      const Standard_Real Qxy,
  72                                      const Standard_Real Qxz,
  73                                      const Standard_Real Qyz,
  74                                      const Standard_Real Qx,
  75                                      const Standard_Real Qy,
  76                                      const Standard_Real Qz,
  77                                      const Standard_Real Q1,
  78                                      const Standard_Real TOL,
  79                                      const Standard_Real DomInf,
  80                                      const Standard_Real DomSup,
  81                                      const Standard_Boolean twocurves,
  82                                      const Standard_Boolean takezpositive)
  83 {
  84
  85   Ax3        = Cone.Position();
  86   RCyl       = Cone.RefRadius();
  87
  88   Angle      = Cone.SemiAngle();
  89   Standard_Real UnSurTgAngle = 1.0/(Tan(Cone.SemiAngle()));
  90
  91   typequadric= GeomAbs_Cone;
  92
  93   TwoCurves     = twocurves;         //-- deux  Z pour un meme parametre
  94   TakeZPositive = takezpositive;     //-- Prendre sur la courbe le Z Positif
  95                                      //--   ( -B + Sqrt()) et non (-B - Sqrt())
  96
  97
  98   Z0Cte      = Q1;                   //-- Attention On a    Z?Cos Cos(t)
  99   Z0Sin      = 0.0;                  //-- et Non          2 Z?Cos Cos(t) !!!
 100   Z0Cos      = 0.0;                  //-- Ce pour tous les Parametres
 101   Z0CosCos   = 0.0;                  //--  ie pas de Coefficient 2
 102   Z0SinSin   = 0.0;                  //--     devant les termes CS C S
 103   Z0CosSin   = 0.0;
 104
 105   Z1Cte      = 2.0*(UnSurTgAngle)*Qz;
 106   Z1Sin      = Qy+Qy;
 107   Z1Cos      = Qx+Qx;
 108   Z1CosCos   = 0.0;
 109   Z1SinSin   = 0.0;
 110   Z1CosSin   = 0.0;
 111
 112   Z2Cte      = Qzz * UnSurTgAngle*UnSurTgAngle;
 113   Z2Sin      = (UnSurTgAngle+UnSurTgAngle)*Qyz;
 114   Z2Cos      = (UnSurTgAngle+UnSurTgAngle)*Qxz;
 115   Z2CosCos   = Qxx;
 116   Z2SinSin   = Qyy;
 117   Z2CosSin   = Qxy;
 118
 119   Tolerance  = TOL;
 120   DomainInf = DomInf;
 121   DomainSup = DomSup;
 122
 123   RestrictedInf = RestrictedSup = Standard_True;   //-- Le Domaine est Borne
 124   firstbounded = lastbounded = Standard_False;
 125
 126   myFirstParameter = DomainInf;
 127   myLastParameter = (TwoCurves) ? DomainSup + DomainSup - DomainInf :
 128                                   DomainSup;
 129 }
 130
 131 //=======================================================================
 132 //function : SetCylinderQuadValues
 133 //purpose  : Description de l intersection Cylindre Quadrique
 134 //=======================================================================
 135 void IntAna_Curve::SetCylinderQuadValues(const gp_Cylinder& Cyl,
 136                                          const Standard_Real Qxx,
 137                                          const Standard_Real Qyy,
 138                                          const Standard_Real Qzz,
 139                                          const Standard_Real Qxy,
 140                                          const Standard_Real Qxz,
 141                                          const Standard_Real Qyz,
 142                                          const Standard_Real Qx,
 143                                          const Standard_Real Qy,
 144                                          const Standard_Real Qz,
 145                                          const Standard_Real Q1,
 146                                          const Standard_Real TOL,
 147                                          const Standard_Real DomInf,
 148                                          const Standard_Real DomSup,
 149                                          const Standard_Boolean twocurves,
 150                                          const Standard_Boolean takezpositive)
 151 {
 152
 153   Ax3        = Cyl.Position();
 154   RCyl       = Cyl.Radius();
 155   typequadric= GeomAbs_Cylinder;
 156
 157   TwoCurves     = twocurves;         //-- deux  Z pour un meme parametre
 158   TakeZPositive = takezpositive;     //-- Prendre sur la courbe le Z Positif
 159   Standard_Real RCylmul2 = RCyl+RCyl;         //--   ( -B + Sqrt())
 160
 161   Z0Cte      = Q1;
 162   Z0Sin      = RCylmul2*Qy;
 163   Z0Cos      = RCylmul2*Qx;
 164   Z0CosCos   = Qxx*RCyl*RCyl;
 165   Z0SinSin   = Qyy*RCyl*RCyl;
 166   Z0CosSin   = RCyl*RCyl*Qxy;
 167
 168   Z1Cte      = Qz+Qz;
 169   Z1Sin      = RCylmul2*Qyz;
 170   Z1Cos      = RCylmul2*Qxz;
 171   Z1CosCos   = 0.0;
 172   Z1SinSin   = 0.0;
 173   Z1CosSin   = 0.0;
 174
 175   Z2Cte      = Qzz;
 176   Z2Sin      = 0.0;
 177   Z2Cos      = 0.0;
 178   Z2CosCos   = 0.0;
 179   Z2SinSin   = 0.0;
 180   Z2CosSin   = 0.0;
 181
 182   Tolerance  = TOL;
 183   DomainInf = DomInf;
 184   DomainSup = DomSup;
 185
 186   RestrictedInf = RestrictedSup = Standard_True;
 187   firstbounded = lastbounded = Standard_False;
 188
 189   myFirstParameter = DomainInf;
 190   myLastParameter  = (TwoCurves) ? DomainSup + DomainSup - DomainInf :
 191                                    DomainSup;
 192 }
 193
 194 //=======================================================================
 195 //function : IsOpen
 196 //purpose  :
 197 //=======================================================================
 198 Standard_Boolean IntAna_Curve::IsOpen() const
 199 {
 200   return(RestrictedInf && RestrictedSup);
 201 }
 202
 203 //=======================================================================
 204 //function : Domain
 205 //purpose  :
 206 //=======================================================================
 207 void IntAna_Curve::Domain(Standard_Real& theFirst,
 208                           Standard_Real& theLast) const
 209 {
 210   if (RestrictedInf && RestrictedSup)
 211   {
 212     theFirst = myFirstParameter;
 213     theLast = myLastParameter;
 214   }
 215   else
 216   {
 217     throw Standard_DomainError("IntAna_Curve::Domain");
 218   }
 219 }
 220 //=======================================================================
 221 //function : IsConstant
 222 //purpose  :
 223 //=======================================================================
 224 Standard_Boolean IntAna_Curve::IsConstant() const
 225 {
 226   //-- ???  Pas facile de decider a la seule vue des Param.
 227   return(Standard_False);
 228 }
 229
 230 //=======================================================================
 231 //function : IsFirstOpen
 232 //purpose  :
 233 //=======================================================================
 234 Standard_Boolean IntAna_Curve::IsFirstOpen() const
 235 {
 236   return(firstbounded);
 237 }
 238
 239 //=======================================================================
 240 //function : IsLastOpen
 241 //purpose  :
 242 //=======================================================================
 243 Standard_Boolean IntAna_Curve::IsLastOpen() const
 244 {
 245   return(lastbounded);
 246 }
 247 //=======================================================================
 248 //function : SetIsFirstOpen
 249 //purpose  :
 250 //=======================================================================
 251 void IntAna_Curve::SetIsFirstOpen(const Standard_Boolean Flag)
 252 {
 253   firstbounded = Flag;
 254 }
 255
 256 //=======================================================================
 257 //function : SetIsLastOpen
 258 //purpose  :
 259 //=======================================================================
 260 void IntAna_Curve::SetIsLastOpen(const Standard_Boolean Flag)
 261 {
 262   lastbounded = Flag;
 263 }
 264
 265 //=======================================================================
 266 //function : InternalUVValue
 267 //purpose  :
 268 //=======================================================================
 269 void IntAna_Curve::InternalUVValue(const Standard_Real theta,
 270                                    Standard_Real& Param1,
 271                                    Standard_Real& Param2,
 272                                    Standard_Real& A,
 273                                    Standard_Real& B,
 274                                    Standard_Real& C,
 275                                    Standard_Real& cost,
 276                                    Standard_Real& sint,
 277                                    Standard_Real& SigneSqrtDis) const
 278 {
 279   const Standard_Real aRelTolp = 1.0+Epsilon(1.0), aRelTolm = 1.0-Epsilon(1.0);
 280
 281   // Infinitesimal step of increasing curve parameter. See comment below.
 282   const Standard_Real aDT = 100.0*Epsilon(DomainSup + DomainSup - DomainInf);
 283
 284   Standard_Real Theta=theta;
 285   Standard_Boolean SecondSolution=Standard_False;
 286
 287   if ((Theta<DomainInf*aRelTolm) ||
 288       ((Theta>DomainSup*aRelTolp) && (!TwoCurves)) ||
 289       (Theta>(DomainSup + DomainSup - DomainInf)*aRelTolp))
 290   {
 291     SigneSqrtDis = 0.;
 292     throw Standard_DomainError("IntAna_Curve::Domain");
 293   }
 294
 295   if (Abs(Theta - DomainSup) < aDT)
 296   {
 297     // Point of Null-discriminant.
 298     Theta = DomainSup;
 299   }
 300   else if (Theta>DomainSup)
 301   {
 302     Theta = DomainSup + DomainSup - Theta;
 303     SecondSolution=Standard_True;
 304   }
 305
 306   Param1=Theta;
 307
 308   if(!TwoCurves) {
 309     SecondSolution=TakeZPositive;
 310   }
 311   //
 312   cost = Cos(Theta);
 313   sint = Sin(Theta);
 314   const Standard_Real aSin2t = Sin(Theta + Theta);
 315   const Standard_Real aCos2t = Cos(Theta + Theta);
 316
 317   A=Z2Cte+sint*(Z2Sin+sint*Z2SinSin)+cost*(Z2Cos+cost*Z2CosCos)
 318     + Z2CosSin*aSin2t;
 319
 320   const Standard_Real aDA = cost*Z2Sin - sint*Z2Cos +
 321                             aSin2t*(Z2SinSin - Z2CosCos) +
 322                             aCos2t*(Z2CosSin * Z2CosSin);
 323
 324   B=Z1Cte+sint*(Z1Sin+sint*Z1SinSin)+cost*(Z1Cos+cost*Z1CosCos)
 325     + Z1CosSin*aSin2t;
 326
 327   const Standard_Real aDB = Z1Sin*cost - Z1Cos*sint +
 328                             aSin2t*(Z1SinSin - Z1CosCos) +
 329                             aCos2t*(Z1CosSin + Z1CosSin);
 330
 331   C=Z0Cte+sint*(Z0Sin+sint*Z0SinSin)+cost*(Z0Cos+cost*Z0CosCos)
 332     + Z0CosSin*aSin2t;
 333
 334   const Standard_Real aDC = Z0Sin*cost - Z0Cos*sint +
 335                             aSin2t*(Z0SinSin - Z0CosCos) +
 336                             aCos2t*(Z0CosSin + Z0CosSin);
 337
 338   Standard_Real aDiscriminant = B*B-4.0*A*C;
 339
 340   // We consider that infinitesimal dt = aDT.
 341   // Error of discriminant computation is equal to
 342   // (d(Disc)/dt)*dt, where 1st derivative d(Disc)/dt = 2*B*aDB - 4*(A*aDC + C*aDA).
 343
 344   const Standard_Real aTolD = 2.0*aDT*Abs(B*aDB - 2.0*(A*aDC + C*aDA));
 345
 346   if (aDiscriminant < aTolD)
 347     aDiscriminant = 0.0;
 348
 349   if (Abs(A) <= Precision::PConfusion())
 350   {
 351     if (Abs(B) <= Precision::PConfusion())
 352     {
 353       Param2 = 0.0;
 354     }
 355     else
 356     {
 357       Param2 = -C / B;
 358     }
 359   }
 360   else
 361   {
 362     SigneSqrtDis = (SecondSolution) ? Sqrt(aDiscriminant) : -Sqrt(aDiscriminant);
 363     Param2 = (-B + SigneSqrtDis) / (A + A);
 364   }
 365 }
 366
 367 //=======================================================================
 368 //function : Value
 369 //purpose  :
 370 //=======================================================================
 371 gp_Pnt IntAna_Curve::Value(const Standard_Real theta)
 372 {
 373   Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
 374   //
 375   A=0.0;  B=0.0;   C=0.0;
 376   U=0.0;  V=0.0;
 377   sint=0.0;  cost=0.0;
 378   SigneSqrtDis=0.0;
 379   InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
 380   //-- checked the parameter U and Raises Domain Error if Error
 381   return(InternalValue(U,V));
 382 }
 383 //=======================================================================
 384 //function : D1u
 385 //purpose  :
 386 //=======================================================================
 387 Standard_Boolean IntAna_Curve::D1u(const Standard_Real theta,
 388                                    gp_Pnt& Pt,
 389                                    gp_Vec& Vec)
 390 {
 391   //-- Pour detecter le cas ou le calcul est impossible
 392   Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
 393   A=0.0;  B=0.0;   C=0.0;
 394   U=0.0;  V=0.0;
 395   sint=0.0;  cost=0.0;
 396   //
 397   InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
 398   //
 399   Pt = Value(theta);
 400   if (Abs(A)<1.0e-7 || Abs(SigneSqrtDis)<1.0e-10) return(Standard_False);
 401
 402
 403   //-- Approximation de la derivee (mieux que le calcul mathematique!)
 404   Standard_Real dtheta = (DomainSup - DomainInf)*1.0e-6;
 405   Standard_Real theta2 = theta+dtheta;
 406   if ((theta2<DomainInf) || ((theta2>DomainSup) && (!TwoCurves))
 407       || (theta2>(DomainSup + DomainSup - DomainInf + 1.0e-14)))
 408   {
 409     dtheta = -dtheta;
 410     theta2 = theta+dtheta;
 411   }
 412   gp_Pnt P2 = Value(theta2);
 413   dtheta = 1.0/dtheta;
 414   Vec.SetCoord((P2.X()-Pt.X())*dtheta,
 415                (P2.Y()-Pt.Y())*dtheta,
 416                (P2.Z()-Pt.Z())*dtheta);
 417
 418   return(Standard_True);
 419 }
 420 //=======================================================================
 421 //function : FindParameter
 422 //purpose  : Projects P to the ALine. Returns the list of parameters as a results
 423 //            of projection.
 424 //           Sometimes aline can be self-intersected line (see bug #29807 where
 425 //            ALine goes through the cone apex).
 426 //=======================================================================
 427 void IntAna_Curve::FindParameter(const gp_Pnt& theP,
 428                                  TColStd_ListOfReal& theParams) const
 429 {
 430   const Standard_Real aPIpPI = M_PI + M_PI,
 431     anEpsAng = 1.e-8,
 432     InternalPrecision = 1.e-8, //precision of internal algorithm of values computation
 433     aSqTolPrecision = Precision::SquareConfusion(); //for boundary points to check their coincidence with others
 434
 435   Standard_Real aTheta = 0.0;
 436   //
 437   switch (typequadric)
 438   {
 439     case GeomAbs_Cylinder:
 440     {
 441       Standard_Real aZ;
 442       ElSLib::CylinderParameters(Ax3, RCyl, theP, aTheta, aZ);
 443     }
 444     break;
 445
 446     case GeomAbs_Cone:
 447     {
 448       Standard_Real aZ;
 449       ElSLib::ConeParameters(Ax3, RCyl, Angle, theP, aTheta, aZ);
 450     }
 451     break;
 452
 453     default:
 454       return;
 455   }
 456   //
 457   if (!firstbounded && (DomainInf > aTheta) && ((DomainInf - aTheta) <= anEpsAng))
 458   {
 459     aTheta = DomainInf;
 460   }
 461   else if (!lastbounded && (aTheta > DomainSup) && ((aTheta - DomainSup) <= anEpsAng))
 462   {
 463     aTheta = DomainSup;
 464   }
 465   //
 466   if (aTheta < DomainInf)
 467   {
 468     aTheta = aTheta + aPIpPI;
 469   }
 470   else if (aTheta > DomainSup)
 471   {
 472     aTheta = aTheta - aPIpPI;
 473   }
 474
 475   const Standard_Integer aMaxPar = 5;
 476   Standard_Real aParams[aMaxPar] = {DomainInf, DomainSup, aTheta,
 477                                     (TwoCurves)? DomainSup + DomainSup - aTheta : RealLast(),
 478                                     (TwoCurves) ? DomainSup + DomainSup - DomainInf : RealLast()};
 479
 480   std::sort(aParams, aParams + aMaxPar - 1);
 481
 482   for (Standard_Integer i = 0; i < aMaxPar; i++)
 483   {
 484     if (aParams[i] > myLastParameter)
 485       break;
 486
 487     if (aParams[i] < myFirstParameter)
 488       continue;
 489
 490     if (i && (aParams[i] - aParams[i - 1]) < Precision::PConfusion())
 491       continue;
 492
 493     Standard_Real U = 0.0, V= 0.0,
 494                   A = 0.0, B = 0.0, C = 0.0,
 495                   sint = 0.0, cost = 0.0, SigneSqrtDis = 0.0;
 496     InternalUVValue(aParams[i], U, V, A, B, C,
 497                     cost, sint, SigneSqrtDis);
 498     const gp_Pnt aP(InternalValue(U, V));
 499
 500     Standard_Real aSqTol;
 501     if (aParams[i] == aTheta ||
 502         (TwoCurves && aParams[i] == DomainSup + DomainSup - aTheta))
 503       aSqTol = InternalPrecision;
 504     else
 505       aSqTol = aSqTolPrecision;
 506
 507     if (aP.SquareDistance(theP) < aSqTol)
 508     {
 509       theParams.Append(aParams[i]);
 510     }
 511   }
 512 }
 513 //=======================================================================
 514 //function : InternalValue
 515 //purpose  :
 516 //=======================================================================
 517 gp_Pnt IntAna_Curve::InternalValue(const Standard_Real U,
 518                                    const Standard_Real _V) const
 519 {
 520   //-- std::cout<<" ["<<U<<","<<V<<"]";
 521   Standard_Real V = _V;
 522   if(V > 100000.0 )   {   V= 100000.0; }
 523   if(V < -100000.0 )  {   V=-100000.0; }
 524
 525   switch(typequadric) {
 526   case GeomAbs_Cone:
 527     {
 528       //------------------------------------------------
 529       //-- Parametrage : X = V * Cos(U)              ---
 530       //--               Y = V * Sin(U)              ---
 531       //--               Z = (V-RCyl) / Tan(SemiAngle)--
 532       //------------------------------------------------
 533       //-- Angle Vaut Cone.SemiAngle()
 534       return(ElSLib::ConeValue(U,(V-RCyl)/Sin(Angle),Ax3,RCyl,Angle));
 535     }
 536     break;
 537
 538   case GeomAbs_Cylinder:
 539     return(ElSLib::CylinderValue(U,V,Ax3,RCyl));
 540   case GeomAbs_Sphere:
 541     return(ElSLib::SphereValue(U,V,Ax3,RCyl));
 542   default:
 543     return(gp_Pnt(0.0,0.0,0.0));
 544   }
 545 }
 546
 547 //
 548 //=======================================================================
 549 //function : SetDomain
 550 //purpose  :
 551 //=======================================================================
 552 void IntAna_Curve::SetDomain(const Standard_Real theFirst,
 553                              const Standard_Real theLast)
 554 {
 555   if (theLast <= theFirst)
 556   {
 557     throw Standard_DomainError("IntAna_Curve::Domain");
 558   }
 559   //
 560   myFirstParameter = theFirst;
 561   myLastParameter = theLast;
 562 }