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