src/Contap/Contap_ContAna.cxx

   1 // Created on: 1993-03-04
   2 // Created by: Jacques GOUSSARD
   3 // Copyright (c) 1993-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
  18 #include <Contap_ContAna.hxx>
  19 #include <gp.hxx>
  20 #include <gp_Cone.hxx>
  21 #include <gp_Cylinder.hxx>
  22 #include <gp_Dir.hxx>
  23 #include <gp_Lin.hxx>
  24 #include <gp_Pnt.hxx>
  25 #include <gp_Sphere.hxx>
  26 #include <gp_XYZ.hxx>
  27 #include <Standard_DomainError.hxx>
  28 #include <Standard_OutOfRange.hxx>
  29 #include <StdFail_NotDone.hxx>
  30
  31 static const Standard_Real Tolpetit = 1.e-8;
  32
  33 Contap_ContAna::Contap_ContAna (): done(Standard_False) {}
  34
  35 void Contap_ContAna::Perform (const gp_Sphere& S,
  36                               const gp_Dir& D)
  37 {
  38   done  = Standard_False;
  39   typL  = GeomAbs_Circle;
  40   pt1   = S.Location();
  41   dir1  = D;
  42   if (Abs(D.Dot(S.XAxis().Direction())) < 0.9999999999999) {
  43     dir2 = D.Crossed(S.XAxis().Direction());
  44   }
  45   else {
  46     dir2 = D.Crossed(S.YAxis().Direction());
  47   }
  48   prm   = S.Radius();
  49   nbSol = 1;
  50   done  = Standard_True;
  51 }
  52
  53 void Contap_ContAna::Perform (const gp_Sphere& S,
  54                               const gp_Dir& D,
  55                               const Standard_Real Angle)
  56 {
  57   done  = Standard_False;
  58   typL  = GeomAbs_Circle;
  59
  60   dir1  = D;
  61   if (Abs(D.Dot(S.XAxis().Direction())) < 0.9999999999999) {
  62     dir2 = D.Crossed(S.XAxis().Direction());
  63   }
  64   else {
  65     dir2 = D.Crossed(S.YAxis().Direction());
  66   }
  67   Standard_Real alpha = (S.Direct() ? Angle : -Angle);
  68   pt1.SetXYZ(S.Location().XYZ() - S.Radius()*sin(alpha)*D.XYZ()) ;
  69   prm   = S.Radius()*cos(alpha);
  70   nbSol = 1;
  71   done  = Standard_True;
  72 }
  73
  74 void Contap_ContAna::Perform (const gp_Sphere& S,
  75                               const gp_Pnt& Eye)
  76 {
  77   done = Standard_False;
  78
  79   Standard_Real radius = S.Radius();
  80   Standard_Real dist = Eye.Distance(S.Location());
  81   if (dist <= radius) {
  82     nbSol = 0;
  83   }
  84   else {
  85     prm = radius*sqrt(1.-radius*radius/(dist*dist));
  86     if (prm < Tolpetit) {
  87       nbSol = 0;
  88     }
  89     else {
  90       gp_XYZ locxyz(S.Location().XYZ());
  91       dir1.SetXYZ(Eye.XYZ()-locxyz);
  92       pt1.SetXYZ(locxyz + (radius*radius/dist)*dir1.XYZ());
  93       if (Abs(dir1.Dot(S.XAxis().Direction())) < 0.9999999999999) {
  94         dir2 = dir1.Crossed(S.XAxis().Direction());
  95       }
  96       else {
  97         dir2 = dir1.Crossed(S.YAxis().Direction());
  98       }
  99       nbSol = 1;
 100       typL = GeomAbs_Circle;
 101     }
 102   }
 103   done = Standard_True;
 104 }
 105
 106 void Contap_ContAna::Perform (const gp_Cylinder& C,
 107                               const gp_Dir& D)
 108 {
 109   done = Standard_False;
 110
 111   gp_XYZ normale(C.Position().Direction().XYZ());
 112   normale.Cross(D.XYZ());
 113   if (normale.Modulus() <= 1e-15) {
 114     nbSol = 0;
 115   }
 116   else {
 117     normale.Normalize();
 118     typL = GeomAbs_Line;
 119     dir1 = C.Position().Direction();
 120     dir2 = dir1;
 121     pt1.SetXYZ(C.Location().XYZ() + C.Radius()*normale);
 122     pt2.SetXYZ(C.Location().XYZ() - C.Radius()*normale);
 123     nbSol = 2;
 124   }
 125
 126   done = Standard_True;
 127 }
 128
 129 void Contap_ContAna::Perform (const gp_Cylinder& C,
 130                               const gp_Dir& D,
 131                               const Standard_Real Angle)
 132 {
 133   done = Standard_False;
 134
 135   Standard_Real Coefcos = D.Dot(C.Position().XDirection());
 136   Standard_Real Coefsin = D.Dot(C.Position().YDirection());
 137   Standard_Real Coefcst = cos(M_PI*0.5 + Angle);
 138
 139   Standard_Real norm1 = Coefcos*Coefcos + Coefsin*Coefsin;
 140   Standard_Real norm2 = sqrt(norm1);
 141
 142   if (Abs(Coefcst) < norm2) {
 143     typL = GeomAbs_Line;
 144     nbSol = 2;
 145     dir1 = dir2 = C.Position().Direction();
 146
 147     if (!C.Direct()) { // The normal is inverted.
 148       Coefcos = -Coefcos;
 149       Coefsin = -Coefsin;
 150     }
 151
 152     // Necessary to solve Coefcos*cos(t) + Coefsin*sin(t) = Coefcst
 153     // and the origins of solution are in the reference of the
 154     // cylinder in (R*cost0, R*sint0,0) and (R*cost1,R*sint1,0)
 155     // By setting cos(phi) = Coefcos/Sqrt(Coefcos**2 + Coefsin**2) and
 156     //           sin(phi) = Coefsin/Sqrt(Coefcos**2 + Coefsin**2)
 157     // and by using trigonometric relations the values of cosinus
 158     // and sinus to the solutions are obtained.
 159
 160     prm = Sqrt(norm1 - Coefcst*Coefcst);
 161     Standard_Real cost0,sint0,cost1,sint1;
 162
 163     cost0 = (Coefcos*Coefcst - Coefsin*prm)/norm1;
 164     cost1 = (Coefcos*Coefcst + Coefsin*prm)/norm1;
 165
 166     sint0 = ( Coefcos*prm + Coefsin*Coefcst)/norm1;
 167     sint1 = (-Coefcos*prm + Coefsin*Coefcst)/norm1;
 168
 169     gp_XYZ Xdir(C.Position().XDirection().XYZ());
 170     gp_XYZ Ydir(C.Position().YDirection().XYZ());
 171     gp_XYZ dirxyz;
 172
 173     dirxyz.SetLinearForm(cost0,Xdir,sint0,Ydir);
 174     dirxyz.Multiply(C.Radius());
 175     pt1.SetXYZ(C.Location().XYZ().Added(dirxyz));
 176
 177     dirxyz.SetLinearForm(cost1,Xdir,sint1,Ydir);
 178     dirxyz.Multiply(C.Radius());
 179     pt2.SetXYZ(C.Location().XYZ().Added(dirxyz));
 180   }
 181   else {
 182     nbSol = 0;
 183   }
 184
 185   done = Standard_True;
 186 }
 187
 188 void Contap_ContAna::Perform (const gp_Cylinder& C,
 189                               const gp_Pnt& Eye)
 190 {
 191   done = Standard_False;
 192
 193   Standard_Real radius = C.Radius();
 194   gp_Lin theaxis(C.Axis());
 195   Standard_Real dist = theaxis.Distance(Eye);
 196   if (dist <= radius) {
 197     nbSol = 0;
 198   }
 199   else {
 200     typL = GeomAbs_Line;
 201     prm = radius*sqrt(1.-radius*radius/(dist*dist));
 202     dir1 = C.Axis().Direction();
 203     dir2 = dir1;
 204     gp_XYZ axeye(theaxis.Normal(Eye).Direction().XYZ()); // orientate the axis to the outside
 205     gp_XYZ normale((theaxis.Direction().Crossed(axeye)).XYZ());
 206 //      normale.Normalize();
 207     pt1.SetXYZ(C.Location().XYZ() + (radius*radius/dist)*axeye);
 208     pt2.SetXYZ(pt1.XYZ() - prm*normale);
 209     pt1.SetXYZ(pt1.XYZ() + prm*normale);
 210     nbSol = 2;
 211   }
 212   done = Standard_True;
 213 }
 214
 215 void Contap_ContAna::Perform (const gp_Cone& C,
 216                               const gp_Dir& D)
 217 {
 218   done = Standard_False;
 219
 220   Standard_Real Tgtalpha = Tan(C.SemiAngle());
 221
 222   Standard_Real Coefcos = D.Dot(C.Position().XDirection());
 223   Standard_Real Coefsin = D.Dot(C.Position().YDirection());
 224   Standard_Real Coefcst = D.Dot(C.Axis().Direction())*Tgtalpha;
 225
 226   Standard_Real norm1 = Coefcos*Coefcos + Coefsin*Coefsin;
 227   Standard_Real norm2 = Sqrt(norm1);
 228 //  if (Abs(Abs(Coefcst)-norm2) <= Tolpetit) { // tol angulaire 1.e-8
 229 //    typL = GeomAbs_Line;
 230 //    nbSol = 1;
 231 //    pt1 = C.Apex();
 232 //    dir1 = D;
 233 //  }
 234 //  else if (Abs(Coefcst) < norm2) {
 235
 236   if (Abs(Coefcst) < norm2) {
 237     typL = GeomAbs_Line;
 238     nbSol = 2;
 239     pt1 = C.Apex();
 240     pt2 = pt1;
 241     // Necessary to solve Coefcos*cos(t) + Coefsin*sin(t) = Coefcst
 242     // and director vectors of solutions are
 243     // cos(t0) * XDirection + sin(t0) * YDirection + ZDirection/Tgtalpha
 244     // By setting cos(phi) = Coefcos/Sqrt(Coefcos**2 + Coefsin**2) and
 245     //           sin(phi) = Coefsin/Sqrt(Coefcos**2 + Coefsin**2)
 246     // and by using trigonometric relations the values of cosinus
 247     // and sinus to the solutions are obtained.
 248
 249     prm = Sqrt(norm1 - Coefcst*Coefcst);
 250     Standard_Real cost0,sint0,cost1,sint1;
 251
 252     cost0 = (Coefcos*Coefcst - Coefsin*prm)/norm1;
 253     cost1 = (Coefcos*Coefcst + Coefsin*prm)/norm1;
 254
 255     sint0 = ( Coefcos*prm + Coefsin*Coefcst)/norm1;
 256     sint1 = (-Coefcos*prm + Coefsin*Coefcst)/norm1;
 257
 258     gp_XYZ Xdir(C.Position().XDirection().XYZ());
 259     gp_XYZ Ydir(C.Position().YDirection().XYZ());
 260     gp_XYZ Zdir(C.Axis().Direction().XYZ());
 261     gp_XYZ dirxyz;
 262     dirxyz.SetLinearForm(cost0,Xdir,sint0,Ydir,1./Tgtalpha,Zdir);
 263     dir1.SetXYZ(dirxyz);
 264     pt1.SetXYZ(pt1.XYZ()+dirxyz);
 265     dirxyz.SetLinearForm(cost1,Xdir,sint1,Ydir,1./Tgtalpha,Zdir);
 266     dir2.SetXYZ(dirxyz);
 267     pt2.SetXYZ(pt2.XYZ()+dirxyz);
 268   }
 269   else {
 270     nbSol = 0;
 271   }
 272   done = Standard_True;
 273 }
 274
 275 void Contap_ContAna::Perform (const gp_Cone& C,
 276                               const gp_Dir& D,
 277                               const Standard_Real Angle)
 278 {
 279   done = Standard_False;
 280   nbSol = 0;
 281
 282   Standard_Real Ang = C.SemiAngle();
 283   Standard_Real Cosa = cos(Ang);
 284   Standard_Real Sina = sin(Ang);
 285
 286   Standard_Real Coefcos = D.Dot(C.Position().XDirection());
 287   Standard_Real Coefsin = D.Dot(C.Position().YDirection());
 288
 289   Standard_Real Coefcst1 = cos(M_PI*0.5 + Angle);
 290
 291   Standard_Real norm1 = Coefcos*Coefcos + Coefsin*Coefsin;
 292   Standard_Real norm2 = Sqrt(norm1);
 293
 294   Standard_Real Coefnz = D.Dot(C.Axis().Direction())*Sina;
 295   Standard_Real Coefcst = (Coefcst1 + Coefnz)/Cosa;
 296
 297   if (Abs(Coefcst) < norm2) {
 298     typL = GeomAbs_Line;
 299     nbSol+= 2;
 300     pt1 = C.Apex();
 301     pt2 = pt1;
 302
 303     // It is requiredto solve Coefcos*cos(t) + Coefsin*sin(t) = Coefcst
 304     // and the director vectors of solutions are
 305     // cos(t0) * XDirection + sin(t0) * YDirection + ZDirection/Tgtalpha
 306     // By setting cos(phi) = Coefcos/Sqrt(Coefcos**2 + Coefsin**2) and
 307     //           sin(phi) = Coefsin/Sqrt(Coefcos**2 + Coefsin**2)
 308     // and by using trigonometric relations the values of cosinus
 309     // and sinus to the solutions are obtained.
 310
 311     prm = Sqrt(norm1 - Coefcst*Coefcst);
 312     Standard_Real cost0,sint0,cost1,sint1;
 313
 314     cost0 = (Coefcos*Coefcst - Coefsin*prm)/norm1;
 315     cost1 = (Coefcos*Coefcst + Coefsin*prm)/norm1;
 316
 317     sint0 = ( Coefcos*prm + Coefsin*Coefcst)/norm1;
 318     sint1 = (-Coefcos*prm + Coefsin*Coefcst)/norm1;
 319
 320     gp_XYZ Xdir(C.Position().XDirection().XYZ());
 321     gp_XYZ Ydir(C.Position().YDirection().XYZ());
 322     gp_XYZ Zdir(C.Axis().Direction().XYZ());
 323     if (!C.Direct()) {
 324       Zdir.Reverse();
 325     }
 326     gp_XYZ dirxyz;
 327     dirxyz.SetLinearForm(cost0,Xdir,sint0,Ydir,Cosa/Sina,Zdir);
 328     dir1.SetXYZ(dirxyz);
 329     pt1.SetXYZ(pt1.XYZ()+dirxyz);
 330     dirxyz.SetLinearForm(cost1,Xdir,sint1,Ydir,Cosa/Sina,Zdir);
 331     dir2.SetXYZ(dirxyz);
 332     pt2.SetXYZ(pt2.XYZ()+dirxyz);
 333   }
 334
 335   Coefcst = (Coefcst1 - Coefnz)/Cosa;
 336
 337   if (Abs(Coefcst) < norm2) {
 338     typL = GeomAbs_Line;
 339     nbSol+= 2;
 340     pt3 = C.Apex();
 341     pt4 = pt3;
 342
 343     prm = Sqrt(norm1 - Coefcst*Coefcst);
 344     Standard_Real cost0,sint0,cost1,sint1;
 345
 346     cost0 = (Coefcos*Coefcst - Coefsin*prm)/norm1;
 347     cost1 = (Coefcos*Coefcst + Coefsin*prm)/norm1;
 348
 349     sint0 = ( Coefcos*prm + Coefsin*Coefcst)/norm1;
 350     sint1 = (-Coefcos*prm + Coefsin*Coefcst)/norm1;
 351
 352     gp_XYZ Xdir(C.Position().XDirection().XYZ());
 353     gp_XYZ Ydir(C.Position().YDirection().XYZ());
 354     gp_XYZ Zdir(C.Axis().Direction().XYZ());
 355     if (!C.Direct()) {
 356       Zdir.Reverse();
 357     }
 358     gp_XYZ dirxyz;
 359     dirxyz.SetLinearForm(cost0,Xdir,sint0,Ydir,-Cosa/Sina,Zdir);
 360     dir3.SetXYZ(dirxyz);
 361     pt3.SetXYZ(pt3.XYZ()+dirxyz);
 362     dirxyz.SetLinearForm(cost1,Xdir,sint1,Ydir,-Cosa/Sina,Zdir);
 363     dir4.SetXYZ(dirxyz);
 364     pt4.SetXYZ(pt4.XYZ()+dirxyz);
 365     if (nbSol == 2) {
 366       pt1 = pt3;
 367       pt2 = pt4;
 368       dir1 = dir3;
 369       dir2 = dir4;
 370     }
 371   }
 372
 373   done = Standard_True;
 374 }
 375
 376 void Contap_ContAna::Perform (const gp_Cone& C,
 377                               const gp_Pnt& Eye)
 378 {
 379   done = Standard_False;
 380
 381   Standard_Real Tgtalpha = Tan(C.SemiAngle());
 382
 383   gp_XYZ apexeye(Eye.XYZ());
 384   apexeye.Subtract(C.Apex().XYZ());
 385
 386   Standard_Real Coefcos = apexeye.Dot(C.Position().XDirection().XYZ());
 387   Standard_Real Coefsin = apexeye.Dot(C.Position().YDirection().XYZ());
 388   Standard_Real Coefcst = apexeye.Dot(C.Axis().Direction().XYZ())*Tgtalpha;
 389
 390   Standard_Real norm1 = Coefcos*Coefcos + Coefsin*Coefsin;
 391   Standard_Real norm2 = Sqrt(Coefcos*Coefcos + Coefsin*Coefsin);
 392 //  if (Abs(Abs(Coefcst)-norm2) <= Tolpetit) { // tol angulaire 1.e-8
 393 //    typL = GeomAbs_Line;
 394 //    nbSol = 1;
 395 //    pt1 = C.Apex();
 396 //    dir1.SetXYZ(apexeye);
 397 //  }
 398 //  else if (Abs(Coefcst) < norm2) {
 399
 400   if (Abs(Coefcst) < norm2) {
 401     typL = GeomAbs_Line;
 402     nbSol = 2;
 403     pt1 = C.Apex();
 404     pt2 = pt1;
 405     // It is required to solve Coefcos*cos(t) + Coefsin*sin(t) = Coefcst
 406     // and the director vectors of solutions are
 407     // cos(t0) * XDirection + sin(t0) * YDirection + ZDirection/Tgtalpha
 408     // By setting cos(phi) = Coefcos/Sqrt(Coefcos**2 + Coefsin**2) and
 409     //           sin(phi) = Coefsin/Sqrt(Coefcos**2 + Coefsin**2)
 410     // and by using trigonometric relations the values of cosinus
 411     // and sinus to the solutions are obtained.
 412
 413     prm = Sqrt(norm1 - Coefcst*Coefcst);
 414     Standard_Real cost0,sint0,cost1,sint1;
 415
 416     cost0 = (Coefcos*Coefcst - Coefsin*prm)/norm1;
 417     cost1 = (Coefcos*Coefcst + Coefsin*prm)/norm1;
 418
 419     sint0 = ( Coefcos*prm + Coefsin*Coefcst)/norm1;
 420     sint1 = (-Coefcos*prm + Coefsin*Coefcst)/norm1;
 421
 422     gp_XYZ Xdir(C.Position().XDirection().XYZ());
 423     gp_XYZ Ydir(C.Position().YDirection().XYZ());
 424     gp_XYZ Zdir(C.Axis().Direction().XYZ());
 425     gp_XYZ dirxyz;
 426     dirxyz.SetLinearForm(cost0,Xdir,sint0,Ydir,1./Tgtalpha,Zdir);
 427     dir1.SetXYZ(dirxyz);
 428     pt1.SetXYZ(pt1.XYZ()+dirxyz);
 429     dirxyz.SetLinearForm(cost1,Xdir,sint1,Ydir,1./Tgtalpha,Zdir);
 430     dir2.SetXYZ(dirxyz);
 431     pt2.SetXYZ(pt2.XYZ()+dirxyz);
 432   }
 433   else {
 434     nbSol = 0;
 435   }
 436   done = Standard_True;
 437 }
 438
 439 gp_Lin Contap_ContAna::Line (const Standard_Integer Index) const
 440 {
 441   if (!done) {StdFail_NotDone::Raise();}
 442   if (typL != GeomAbs_Line || nbSol == 0) {Standard_DomainError::Raise();}
 443   if (Index <=0 || Index > nbSol) {Standard_OutOfRange::Raise();}
 444   switch (Index) {
 445   case 1:
 446     return gp_Lin(pt1,dir1);
 447   case 2:
 448     return gp_Lin(pt2,dir2);
 449   case 3:
 450     return gp_Lin(pt3,dir3);
 451   case 4:
 452     return gp_Lin(pt4,dir4);
 453   }
 454   Standard_OutOfRange::Raise("Program error in Contap_ContAna");
 455   return gp_Lin();
 456 }