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