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