src/GeomLib/GeomLib_IsPlanarSurface.cxx

   1 // Created on: 1998-11-23
   2 // Created by: Philippe MANGIN
   3 // Copyright (c) 1998-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 <Geom_BezierCurve.hxx>
  19 #include <Geom_BezierSurface.hxx>
  20 #include <Geom_BSplineCurve.hxx>
  21 #include <Geom_BSplineSurface.hxx>
  22 #include <Geom_Curve.hxx>
  23 #include <Geom_Surface.hxx>
  24 #include <GeomAbs_CurveType.hxx>
  25 #include <GeomAbs_SurfaceType.hxx>
  26 #include <GeomAdaptor_Curve.hxx>
  27 #include <GeomAdaptor_Surface.hxx>
  28 #include <GeomLib.hxx>
  29 #include <GeomLib_IsPlanarSurface.hxx>
  30 #include <gp_Pln.hxx>
  31 #include <StdFail_NotDone.hxx>
  32 #include <TColgp_Array1OfPnt.hxx>
  33 #include <TColgp_HArray1OfPnt.hxx>
  34
  35 static Standard_Boolean Controle(const TColgp_Array1OfPnt& P,
  36                           const gp_Pln& Plan,
  37                           const Standard_Real Tol)
  38 {
  39   Standard_Integer ii;
  40   Standard_Boolean B=Standard_True;
  41
  42   for (ii=1; ii<=P.Length() && B; ii++)
  43       B = (Plan.Distance(P(ii)) < Tol);
  44
  45   return B;
  46 }
  47
  48 static Standard_Boolean Controle(const TColgp_Array1OfPnt& Poles,
  49                                  const Standard_Real Tol,
  50                                  const Handle(Geom_Surface)& S,
  51                                  gp_Pln& Plan)
  52 {
  53   Standard_Boolean IsPlan = Standard_False;
  54   Standard_Boolean Essai = Standard_True;
  55   Standard_Real gx,gy,gz;
  56   Standard_Integer Nb = Poles.Length();
  57     gp_Pnt Bary;
  58   gp_Dir DX, DY;
  59
  60
  61   if (Nb > 10) {
  62     // Test allege (pour une rejection rapide)
  63     TColgp_Array1OfPnt Aux(1,5);
  64     Aux(1) = Poles(1);
  65     Aux(2) = Poles(Nb/3);
  66     Aux(3) = Poles(Nb/2);
  67     Aux(4) = Poles(Nb/2+Nb/3);
  68     Aux(5) =  Poles(Nb);
  69     GeomLib::Inertia(Aux, Bary, DX, DY, gx, gy, gz);
  70     Essai = (gz<Tol);
  71   }
  72
  73   if (Essai) { // Test Grandeur nature...
  74     GeomLib::Inertia(Poles, Bary, DX, DY, gx, gy, gz);
  75     if (gz<Tol && gy>Tol) {
  76       gp_Pnt P;
  77       gp_Vec DU, DV;
  78       Standard_Real umin, umax, vmin, vmax;
  79       S->Bounds(umin, umax, vmin, vmax);
  80       S->D1( (umin+umax)/2, (vmin+vmax)/2, P, DU, DV);
  81       // On prend DX le plus proche possible de DU
  82       gp_Dir du(DU);
  83       Standard_Real Angle1 = du.Angle(DX);
  84       Standard_Real Angle2 = du.Angle(DY);
  85       if (Angle1 > M_PI/2) Angle1 = M_PI-Angle1;
  86       if (Angle2 > M_PI/2) Angle2 = M_PI-Angle2;
  87       if (Angle2 < Angle1) {
  88         du = DY; DY = DX; DX = du;
  89       }
  90       if (DX.Angle(DU) > M_PI/2) DX.Reverse();
  91       if (DY.Angle(DV) > M_PI/2) DY.Reverse();
  92
  93       gp_Ax3 axe(Bary, DX^DY, DX);
  94       Plan.SetPosition(axe);
  95       Plan.SetLocation(Bary);
  96       IsPlan = Standard_True;
  97     }
  98   }
  99   return IsPlan;
 100 }
 101
 102 static Standard_Boolean Controle(const Handle(Geom_Curve)& C,
 103                           const gp_Pln& Plan,
 104                           const Standard_Real Tol)
 105 {
 106   Standard_Boolean B = Standard_True;
 107   Standard_Integer ii, Nb;
 108   GeomAbs_CurveType Type;
 109   GeomAdaptor_Curve AC(C);
 110   Type = AC.GetType();
 111   Handle(TColgp_HArray1OfPnt) TabP;
 112   TabP.Nullify();
 113
 114   switch (Type) {
 115   case GeomAbs_Line :
 116     {
 117       Nb = 2;
 118       break;
 119     }
 120   case GeomAbs_Circle:
 121     {
 122       Nb = 3;
 123       break;
 124     }
 125
 126   case GeomAbs_Ellipse:
 127   case GeomAbs_Hyperbola:
 128   case GeomAbs_Parabola:
 129     {
 130       Nb = 5;
 131       break;
 132     }
 133   case GeomAbs_BezierCurve:
 134     {
 135       Nb =  AC.NbPoles();
 136       Handle (Geom_BezierCurve) BZ = AC.Bezier();
 137       TabP = new (TColgp_HArray1OfPnt) (1, AC.NbPoles());
 138       for (ii=1; ii<=Nb; ii++)
 139         TabP->SetValue(ii, BZ->Pole(ii));
 140       break;
 141     }
 142   case GeomAbs_BSplineCurve:
 143     {
 144       Nb =  AC.NbPoles();
 145       Handle (Geom_BSplineCurve) BZ = AC.BSpline();
 146       TabP = new (TColgp_HArray1OfPnt) (1, AC.NbPoles());
 147       for (ii=1; ii<=Nb; ii++)
 148         TabP->SetValue(ii, BZ->Pole(ii));
 149       break;
 150     }
 151     default :
 152       {
 153         Nb = 8 + 3*AC.NbIntervals(GeomAbs_CN);
 154       }
 155   }
 156
 157   if (TabP.IsNull()) {
 158     Standard_Real u, du, f, l, d;
 159     f = AC.FirstParameter();
 160     l = AC.LastParameter();
 161     du = (l-f)/(Nb-1);
 162     for (ii=1; ii<=Nb && B ; ii++) {
 163       u = (ii-1)*du + f;
 164       d = Plan.Distance(C->Value(u));
 165       B = (d < Tol);
 166     }
 167   }
 168   else {
 169     B = Controle(TabP->Array1(), Plan, Tol);
 170   }
 171
 172   return B;
 173 }
 174
 175
 176
 177 GeomLib_IsPlanarSurface::GeomLib_IsPlanarSurface(const Handle(Geom_Surface)& S,
 178                                                  const Standard_Real Tol)
 179
 180 {
 181   GeomAdaptor_Surface AS(S);
 182   GeomAbs_SurfaceType Type;
 183
 184   Type = AS.GetType();
 185
 186   switch (Type) {
 187   case GeomAbs_Plane :
 188     {
 189       IsPlan = Standard_True;
 190       myPlan = AS.Plane();
 191       break;
 192     }
 193   case GeomAbs_Cylinder :
 194   case GeomAbs_Cone :
 195   case GeomAbs_Sphere :
 196   case GeomAbs_Torus :
 197     {
 198      IsPlan = Standard_False;
 199      break;
 200     }
 201   case GeomAbs_BezierSurface :
 202   case GeomAbs_BSplineSurface :
 203     {
 204       Standard_Integer ii, jj, kk,
 205       NbU = AS.NbUPoles(), NbV = AS.NbVPoles();
 206       TColgp_Array1OfPnt Poles(1, NbU*NbV);
 207       if (Type == GeomAbs_BezierSurface) {
 208         Handle(Geom_BezierSurface) BZ;
 209         BZ = AS.Bezier();
 210         for(ii=1, kk=1; ii<=NbU; ii++)
 211           for(jj=1; jj<=NbV; jj++,kk++)
 212             Poles(kk) = BZ->Pole(ii,jj);
 213       }
 214       else {
 215         Handle(Geom_BSplineSurface) BS;
 216         BS = AS.BSpline();
 217         for(ii=1, kk=1; ii<=NbU; ii++)
 218           for(jj=1; jj<=NbV; jj++,kk++)
 219             Poles(kk) = BS->Pole(ii,jj);
 220       }
 221
 222       IsPlan =  Controle(Poles, Tol, S, myPlan);
 223       break;
 224     }
 225
 226   case GeomAbs_SurfaceOfRevolution :
 227     {
 228       Standard_Boolean Essai = Standard_True;
 229       gp_Pnt P;
 230       gp_Vec DU, DV, Dn;
 231       gp_Dir Dir = AS.AxeOfRevolution().Direction();
 232       Standard_Real Umin, Umax, Vmin, Vmax;
 233       S->Bounds(Umin, Umax, Vmin, Vmax);
 234       S->D1((Umin+Umax)/2, (Vmin+Vmax)/2, P, DU, DV);
 235       if (DU.Magnitude() <= gp::Resolution() ||
 236           DV.Magnitude() <= gp::Resolution())
 237       {
 238         Standard_Real NewU = (Umin+Umax)/2 + (Umax-Umin)*0.1;
 239         Standard_Real NewV = (Vmin+Vmax)/2 + (Vmax-Vmin)*0.1;
 240         S->D1( NewU, NewV, P, DU, DV );
 241       }
 242       Dn = DU^DV;
 243       if (Dn.Magnitude() > 1.e-7) {
 244         Standard_Real angle = Dir.Angle(Dn);
 245         if (angle > M_PI/2) {
 246           angle = M_PI - angle;
 247           Dir.Reverse();
 248         }
 249         Essai = (angle < 0.1);
 250       }
 251
 252       if (Essai) {
 253         gp_Ax3 axe(P, Dir);
 254         axe.SetXDirection(DU);
 255         myPlan.SetPosition(axe);
 256         myPlan.SetLocation(P);
 257         Handle(Geom_Curve) C;
 258         C = S->UIso(Umin);
 259         IsPlan = Controle(C,  myPlan, Tol);
 260       }
 261       else
 262         IsPlan = Standard_False;
 263
 264       break;
 265     }
 266   case GeomAbs_SurfaceOfExtrusion :
 267     {
 268       Standard_Boolean Essai = Standard_False;
 269       Standard_Real Umin, Umax, Vmin, Vmax;
 270       Standard_Real norm;
 271       gp_Vec Du, Dv, Dn;
 272       gp_Pnt P;
 273
 274       S->Bounds(Umin, Umax, Vmin, Vmax);
 275       S->D1((Umin+Umax)/2, (Vmin+Vmax)/2, P, Du, Dv);
 276       if (Du.Magnitude() <= gp::Resolution() ||
 277           Dv.Magnitude() <= gp::Resolution())
 278       {
 279         Standard_Real NewU = (Umin+Umax)/2 + (Umax-Umin)*0.1;
 280         Standard_Real NewV = (Vmin+Vmax)/2 + (Vmax-Vmin)*0.1;
 281         S->D1( NewU, NewV, P, Du, Dv );
 282       }
 283       Dn = Du^Dv;
 284       norm = Dn.Magnitude();
 285       if (norm > 1.e-15) {
 286         Dn /= norm;
 287         Standard_Real angmax = Tol / (Vmax-Vmin);
 288         gp_Dir D(Dn);
 289         Essai = (D.IsNormal(AS.Direction(), angmax));
 290       }
 291       if (Essai) {
 292         gp_Ax3 axe(P, Dn, Du);
 293         myPlan.SetPosition(axe);
 294         myPlan.SetLocation(P);
 295         Handle(Geom_Curve) C;
 296         C = S->VIso((Vmin+Vmax)/2);
 297         IsPlan = Controle(C,  myPlan, Tol);
 298       }
 299       else
 300         IsPlan = Standard_False;
 301       break;
 302     }
 303
 304     default :
 305       {
 306         Standard_Integer NbU,NbV, ii, jj, kk;
 307         NbU = 8 + 3*AS.NbUIntervals(GeomAbs_CN);
 308         NbV = 8 + 3*AS.NbVIntervals(GeomAbs_CN);
 309         Standard_Real Umin, Umax, Vmin, Vmax, du, dv, U, V;
 310         S->Bounds(Umin, Umax, Vmin, Vmax);
 311         du = (Umax-Umin)/(NbU-1);
 312         dv = (Vmax-Vmin)/(NbV-1);
 313         TColgp_Array1OfPnt Pnts(1, NbU*NbV);
 314         for(ii=0, kk=1; ii<NbU; ii++) {
 315           U = Umin + du*ii;
 316           for(jj=0; jj<NbV; jj++,kk++) {
 317             V = Vmin + dv*jj;
 318             S->D0(U,V, Pnts(kk));
 319           }
 320         }
 321
 322         IsPlan =  Controle(Pnts, Tol, S, myPlan);
 323       }
 324     }
 325 }
 326
 327  Standard_Boolean GeomLib_IsPlanarSurface::IsPlanar() const
 328 {
 329   return IsPlan;
 330 }
 331
 332  const gp_Pln& GeomLib_IsPlanarSurface::Plan() const
 333 {
 334   if (!IsPlan) throw StdFail_NotDone(" GeomLib_IsPlanarSurface");
 335   return myPlan;
 336 }