src/IntSurf/IntSurf_Quadric.cxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and/or modify it under
   7 // the terms of the GNU Lesser General Public License version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15
  16 #include <ElCLib.hxx>
  17 #include <ElSLib.hxx>
  18 #include <gp_Cone.hxx>
  19 #include <gp_Cylinder.hxx>
  20 #include <gp_Pln.hxx>
  21 #include <gp_Pnt.hxx>
  22 #include <gp_Sphere.hxx>
  23 #include <gp_Torus.hxx>
  24 #include <gp_Vec.hxx>
  25 #include <IntSurf_Quadric.hxx>
  26 #include <StdFail_NotDone.hxx>
  27
  28 // ============================================================
  29 IntSurf_Quadric::IntSurf_Quadric ():typ(GeomAbs_OtherSurface),
  30    prm1(0.), prm2(0.), prm3(0.), prm4(0.), ax3direc(Standard_False)
  31 {
  32 }
  33 // ============================================================
  34 IntSurf_Quadric::IntSurf_Quadric (const gp_Pln& P):
  35       ax3(P.Position()),typ(GeomAbs_Plane)
  36 {
  37   ax3direc = ax3.Direct();
  38   P.Coefficients(prm1,prm2,prm3,prm4);
  39 }
  40 // ============================================================
  41 IntSurf_Quadric::IntSurf_Quadric (const gp_Cylinder& C):
  42
  43        ax3(C.Position()),lin(ax3.Axis()),typ(GeomAbs_Cylinder)
  44 {
  45   prm2=prm3=prm4=0.0;
  46   ax3direc=ax3.Direct();
  47   prm1=C.Radius();
  48 }
  49 // ============================================================
  50 IntSurf_Quadric::IntSurf_Quadric (const gp_Sphere& S):
  51
  52        ax3(S.Position()),lin(ax3.Axis()),typ(GeomAbs_Sphere)
  53 {
  54   prm2=prm3=prm4=0.0;
  55   ax3direc = ax3.Direct();
  56   prm1=S.Radius();
  57 }
  58 // ============================================================
  59 IntSurf_Quadric::IntSurf_Quadric (const gp_Cone& C):
  60
  61        ax3(C.Position()),typ(GeomAbs_Cone)
  62 {
  63   ax3direc = ax3.Direct();
  64   lin.SetPosition(ax3.Axis());
  65   prm1 = C.RefRadius();
  66   prm2 = C.SemiAngle();
  67   prm3 = Cos(prm2);
  68   prm4 = 0.0;
  69 }
  70 // ============================================================
  71 IntSurf_Quadric::IntSurf_Quadric (const gp_Torus& T):
  72
  73        ax3(T.Position()),typ(GeomAbs_Torus)
  74 {
  75   ax3direc = ax3.Direct();
  76   lin.SetPosition(ax3.Axis());
  77   prm1 = T.MajorRadius();
  78   prm2 = T.MinorRadius();
  79   prm3 = 0.0;
  80   prm4 = 0.0;
  81 }
  82 // ============================================================
  83 void IntSurf_Quadric::SetValue (const gp_Pln& P)
  84 {
  85   typ = GeomAbs_Plane;
  86   ax3 = P.Position();
  87   ax3direc = ax3.Direct();
  88   P.Coefficients(prm1,prm2,prm3,prm4);
  89 }
  90 // ============================================================
  91 void IntSurf_Quadric::SetValue (const gp_Cylinder& C)
  92 {
  93   typ = GeomAbs_Cylinder;
  94   ax3 = C.Position();
  95   ax3direc = ax3.Direct();
  96   lin.SetPosition(ax3.Axis());
  97   prm1 = C.Radius();
  98   prm2=prm3=prm4=0.0;
  99 }
 100 // ============================================================
 101 void IntSurf_Quadric::SetValue (const gp_Sphere& S)
 102 {
 103   typ = GeomAbs_Sphere;
 104   ax3 = S.Position();
 105   ax3direc = ax3.Direct();
 106   lin.SetPosition(ax3.Axis());
 107   prm1 = S.Radius();
 108   prm2=prm3=prm4=0.0;
 109 }
 110 // ============================================================
 111 void IntSurf_Quadric::SetValue (const gp_Cone& C)
 112 {
 113   typ = GeomAbs_Cone;
 114   ax3 = C.Position();
 115   ax3direc = ax3.Direct();
 116   lin.SetPosition(ax3.Axis());
 117   prm1 = C.RefRadius();
 118   prm2 = C.SemiAngle();
 119   prm3 = Cos(prm2);
 120   prm4 = 0.0;
 121 }
 122 // ============================================================
 123 void IntSurf_Quadric::SetValue (const gp_Torus& T)
 124 {
 125   typ = GeomAbs_Torus;
 126   ax3 = T.Position();
 127   ax3direc = ax3.Direct();
 128   lin.SetPosition(ax3.Axis());
 129   prm1 = T.MajorRadius();
 130   prm2 = T.MinorRadius();
 131   prm3 = 0.0;
 132   prm4 = 0.0;
 133 }
 134 // ============================================================
 135 Standard_Real IntSurf_Quadric::Distance (const gp_Pnt& P) const {
 136   switch (typ) {
 137   case GeomAbs_Plane:   // plan
 138     return prm1*P.X() + prm2*P.Y() + prm3*P.Z() + prm4;
 139   case GeomAbs_Cylinder:   // cylindre
 140     return (lin.Distance(P) - prm1);
 141   case GeomAbs_Sphere:   // sphere
 142     return (lin.Location().Distance(P) - prm1);
 143   case GeomAbs_Cone:   // cone
 144     {
 145       Standard_Real dist = lin.Distance(P);
 146       Standard_Real U,V;
 147       ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
 148       gp_Pnt Pp = ElSLib::ConeValue(U,V,ax3,prm1,prm2);
 149       Standard_Real distp = lin.Distance(Pp);
 150       dist = (dist-distp)/prm3;
 151       return(dist);
 152     }
 153   case GeomAbs_Torus:   // torus
 154     {
 155       gp_Pnt O, Pp, PT;
 156       //
 157       O = ax3.Location();
 158       gp_Vec OZ (ax3.Direction());
 159       Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
 160       //
 161       gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
 162         ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
 163       PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
 164       //
 165       Standard_Real dist = P.Distance(PT) - prm2;
 166       return dist;
 167     }
 168   default:
 169     {
 170     }
 171     break;
 172   }
 173   return(0.0);
 174 }
 175 // ============================================================
 176 gp_Vec IntSurf_Quadric::Gradient (const gp_Pnt& P) const {
 177   gp_Vec grad;
 178   switch (typ) {
 179   case GeomAbs_Plane:   // plan
 180     grad.SetCoord(prm1,prm2,prm3);
 181     break;
 182   case GeomAbs_Cylinder:   // cylindre
 183     {
 184       gp_XYZ PP(lin.Location().XYZ());
 185       PP.Add(ElCLib::Parameter(lin,P)*lin.Direction().XYZ());
 186       grad.SetXYZ(P.XYZ()-PP);
 187       Standard_Real N = grad.Magnitude();
 188       if(N>1e-14) { grad.Divide(N); }
 189       else { grad.SetCoord(0.0,0.0,0.0); }
 190     }
 191     break;
 192   case GeomAbs_Sphere:   // sphere
 193     {
 194       gp_XYZ PP(P.XYZ());
 195       grad.SetXYZ((PP-lin.Location().XYZ()));
 196       Standard_Real N = grad.Magnitude();
 197       if(N>1e-14) { grad.Divide(N); }
 198       else { grad.SetCoord(0.0,0.0,0.0); }
 199     }
 200     break;
 201   case GeomAbs_Cone:   // cone
 202     {
 203       Standard_Real U,V;
 204       ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
 205       gp_Pnt Pp = ElSLib::ConeValue(U,V,ax3,prm1,prm2);
 206       gp_Vec D1u,D1v;
 207       ElSLib::ConeD1(U,V,ax3,prm1,prm2,Pp,D1u,D1v);
 208       grad=D1u.Crossed(D1v);
 209       if(ax3direc==Standard_False) {
 210         grad.Reverse();
 211       }
 212       grad.Normalize();
 213     }
 214     break;
 215   case GeomAbs_Torus:   // torus
 216     {
 217       gp_Pnt O, Pp, PT;
 218       //
 219       O = ax3.Location();
 220       gp_Vec OZ (ax3.Direction());
 221       Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
 222       //
 223       gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
 224         ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
 225       PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
 226       //
 227       grad.SetXYZ(P.XYZ() - PT.XYZ());
 228       Standard_Real N = grad.Magnitude();
 229       if(N>1e-14) { grad.Divide(N); }
 230       else { grad.SetCoord(0., 0., 0.); }
 231     }
 232     break;
 233   default:
 234     {}
 235     break;
 236   }
 237   return grad;
 238 }
 239 // ============================================================
 240 void IntSurf_Quadric::ValAndGrad (const gp_Pnt& P,
 241                                   Standard_Real& Dist,
 242                                   gp_Vec& Grad) const
 243 {
 244
 245   switch (typ) {
 246   case GeomAbs_Plane:
 247     {
 248       Dist = prm1*P.X() + prm2*P.Y() + prm3*P.Z() + prm4;
 249       Grad.SetCoord(prm1,prm2,prm3);
 250     }
 251     break;
 252   case GeomAbs_Cylinder:
 253     {
 254       Dist = lin.Distance(P) - prm1;
 255       gp_XYZ PP(lin.Location().XYZ());
 256       PP.Add(ElCLib::Parameter(lin,P)*lin.Direction().XYZ());
 257       Grad.SetXYZ((P.XYZ()-PP));
 258       Standard_Real N = Grad.Magnitude();
 259       if(N>1e-14) { Grad.Divide(N); }
 260       else { Grad.SetCoord(0.0,0.0,0.0); }
 261     }
 262     break;
 263   case GeomAbs_Sphere:
 264     {
 265       Dist = lin.Location().Distance(P) - prm1;
 266       gp_XYZ PP(P.XYZ());
 267       Grad.SetXYZ((PP-lin.Location().XYZ()));
 268       Standard_Real N = Grad.Magnitude();
 269       if(N>1e-14) { Grad.Divide(N); }
 270       else { Grad.SetCoord(0.0,0.0,0.0); }
 271     }
 272     break;
 273   case GeomAbs_Cone:
 274     {
 275       Standard_Real dist = lin.Distance(P);
 276       Standard_Real U,V;
 277       gp_Vec D1u,D1v;
 278       gp_Pnt Pp;
 279       ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
 280       ElSLib::ConeD1(U,V,ax3,prm1,prm2,Pp,D1u,D1v);
 281       Standard_Real distp = lin.Distance(Pp);
 282       dist = (dist-distp)/prm3;
 283       Dist = dist;
 284       Grad=D1u.Crossed(D1v);
 285       if(ax3direc==Standard_False) {
 286         Grad.Reverse();
 287       }
 288       //-- lbr le 7 mars 96
 289       //-- Si le gardient est nul, on est sur l axe
 290       //-- et dans ce cas dist vaut 0
 291       //-- On peut donc renvoyer une valeur quelconque.
 292       if(   Grad.X() > 1e-13
 293          || Grad.Y() > 1e-13
 294          || Grad.Z() > 1e-13) {
 295         Grad.Normalize();
 296       }
 297     }
 298     break;
 299   case GeomAbs_Torus:
 300     {
 301       gp_Pnt O, Pp, PT;
 302       //
 303       O = ax3.Location();
 304       gp_Vec OZ (ax3.Direction());
 305       Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
 306       //
 307       gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
 308         ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
 309       PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
 310       //
 311       Dist = P.Distance(PT) - prm2;
 312       //
 313       Grad.SetXYZ(P.XYZ()-PT.XYZ());
 314       Standard_Real N = Grad.Magnitude();
 315       if(N>1e-14) { Grad.Divide(N); }
 316       else { Grad.SetCoord(0., 0., 0.); }
 317     }
 318     break;
 319   default:
 320     {}
 321     break;
 322   }
 323 }
 324 // ============================================================
 325 gp_Pnt IntSurf_Quadric::Value(const Standard_Real U,
 326                               const Standard_Real V) const
 327 {
 328   switch (typ) {
 329
 330   case GeomAbs_Plane:
 331     return ElSLib::PlaneValue(U,V,ax3);
 332   case GeomAbs_Cylinder:
 333     return ElSLib::CylinderValue(U,V,ax3,prm1);
 334   case GeomAbs_Sphere:
 335     return ElSLib::SphereValue(U,V,ax3,prm1);
 336   case GeomAbs_Cone:
 337     return ElSLib::ConeValue(U,V,ax3,prm1,prm2);
 338   case GeomAbs_Torus:
 339     return ElSLib::TorusValue(U,V,ax3,prm1,prm2);
 340   default:
 341     {
 342       gp_Pnt p(0,0,0);
 343       return(p);
 344     }
 345     //break;
 346   }
 347 // pop : pour NT
 348 //  return gp_Pnt(0,0,0);
 349 }
 350 // ============================================================
 351 void IntSurf_Quadric::D1(const Standard_Real U,
 352                          const Standard_Real V,
 353                          gp_Pnt& P,
 354                          gp_Vec& D1U,
 355                          gp_Vec& D1V) const
 356 {
 357   switch (typ) {
 358   case GeomAbs_Plane:
 359     ElSLib::PlaneD1(U,V,ax3,P,D1U,D1V);
 360     break;
 361   case GeomAbs_Cylinder:
 362     ElSLib::CylinderD1(U,V,ax3,prm1,P,D1U,D1V);
 363     break;
 364   case GeomAbs_Sphere:
 365     ElSLib::SphereD1(U,V,ax3,prm1,P,D1U,D1V);
 366     break;
 367   case GeomAbs_Cone:
 368     ElSLib::ConeD1(U,V,ax3,prm1,prm2,P,D1U,D1V);
 369     break;
 370   case GeomAbs_Torus:
 371     ElSLib::TorusD1(U,V,ax3,prm1,prm2,P,D1U,D1V);
 372     break;
 373   default:
 374     {
 375     }
 376     break;
 377   }
 378 }
 379 // ============================================================
 380 gp_Vec IntSurf_Quadric::DN(const Standard_Real U,
 381                            const Standard_Real V,
 382                            const Standard_Integer Nu,
 383                            const Standard_Integer Nv) const
 384 {
 385   switch (typ) {
 386   case GeomAbs_Plane:
 387     return ElSLib::PlaneDN(U,V,ax3,Nu,Nv);
 388   case GeomAbs_Cylinder:
 389     return ElSLib::CylinderDN(U,V,ax3,prm1,Nu,Nv);
 390   case GeomAbs_Sphere:
 391     return ElSLib::SphereDN(U,V,ax3,prm1,Nu,Nv);
 392   case GeomAbs_Cone:
 393     return ElSLib::ConeDN(U,V,ax3,prm1,prm2,Nu,Nv);
 394   case GeomAbs_Torus:
 395     return ElSLib::TorusDN(U,V,ax3,prm1,prm2,Nu,Nv);
 396   default:
 397     {
 398       gp_Vec v(0,0,0);
 399       return(v);
 400     }
 401     //break;
 402   }
 403 // pop : pour NT
 404 //  return gp_Vec(0,0,0);
 405 }
 406 // ============================================================
 407 gp_Vec IntSurf_Quadric::Normale(const Standard_Real U,
 408                                 const Standard_Real V) const
 409 {
 410   switch (typ) {
 411   case GeomAbs_Plane:
 412     if(ax3direc)
 413       return ax3.Direction();
 414     else
 415       return ax3.Direction().Reversed();
 416   case GeomAbs_Cylinder:
 417     return Normale(Value(U,V));
 418   case GeomAbs_Sphere:
 419     return Normale(Value(U,V));
 420   case GeomAbs_Cone:
 421     {
 422       gp_Pnt P;
 423       gp_Vec D1u,D1v;
 424       ElSLib::ConeD1(U,V,ax3,prm1,prm2,P,D1u,D1v);
 425       if(D1u.Magnitude()<0.0000001) {
 426         gp_Vec Vn(0.0,0.0,0.0);
 427         return(Vn);
 428       }
 429       return(D1u.Crossed(D1v));
 430     }
 431   case GeomAbs_Torus:
 432     return Normale(Value(U,V));
 433   default:
 434     {
 435       gp_Vec v(0,0,0);
 436       return(v);
 437     }
 438   //  break;
 439   }
 440 // pop : pour NT
 441 //  return gp_Vec(0,0,0);
 442 }
 443 // ============================================================
 444 gp_Vec IntSurf_Quadric::Normale (const gp_Pnt& P) const
 445 {
 446   switch (typ) {
 447   case GeomAbs_Plane:
 448     if(ax3direc)
 449       return ax3.Direction();
 450     else
 451       return ax3.Direction().Reversed();
 452   case GeomAbs_Cylinder:
 453     {
 454       if(ax3direc) {
 455         return lin.Normal(P).Direction();
 456       }
 457       else {
 458         gp_Dir D(lin.Normal(P).Direction());
 459         D.Reverse();
 460         return(D);
 461       }
 462     }
 463   case GeomAbs_Sphere:
 464     {
 465       if(ax3direc) {
 466         gp_Vec ax3P(ax3.Location(),P);
 467         return gp_Dir(ax3P);
 468       }
 469       else {
 470         gp_Vec Pax3(P,ax3.Location());
 471         return gp_Dir(Pax3);
 472       }
 473     }
 474   case GeomAbs_Cone:
 475     {
 476       Standard_Real U,V;
 477       ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
 478       return Normale(U,V);
 479     }
 480   case GeomAbs_Torus:
 481     {
 482       gp_Pnt O, Pp, PT;
 483       //
 484       O = ax3.Location();
 485       gp_Vec OZ (ax3.Direction());
 486       Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
 487       //
 488       gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
 489         ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
 490       PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
 491       if (PT.SquareDistance(P) < 1e-14) {
 492         return gp_Dir(OZ);
 493       }
 494       gp_Dir aD(ax3direc ? gp_Vec(PT, P) : gp_Vec(P, PT));
 495       return aD;
 496     }
 497   default:
 498     {
 499       gp_Vec v(0,0,0);
 500       return(v);
 501     } //    break;
 502   }
 503 }
 504 // ============================================================
 505 void IntSurf_Quadric::Parameters (const gp_Pnt& P,
 506                                   Standard_Real& U,
 507                                   Standard_Real& V) const
 508 {
 509   switch (typ) {
 510   case GeomAbs_Plane:
 511     ElSLib::PlaneParameters(ax3,P,U,V);
 512     break;
 513   case GeomAbs_Cylinder:
 514     ElSLib::CylinderParameters(ax3,prm1,P,U,V);
 515     break;
 516   case GeomAbs_Sphere:
 517     ElSLib::SphereParameters(ax3,prm1,P,U,V);
 518     break;
 519   case GeomAbs_Cone:
 520     ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
 521     break;
 522   case GeomAbs_Torus:
 523     ElSLib::TorusParameters(ax3,prm1,prm2,P,U,V);
 524     break;
 525   default:
 526     break;
 527   }
 528 }
 529 // ============================================================
 530
 531
 532
 533
 534