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