0031313: Foundation Classes - Dump improvement for classes
[occt.git] / 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 }