0027719: HLRBrep_Algo incorrect output
[occt.git] / src / HLRBRep / HLRBRep_Surface.cxx
1 // Created on: 1992-03-13
2 // Created by: Christophe MARION
3 // Copyright (c) 1992-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 <BRepAdaptor_Surface.hxx>
19 #include <BRepClass_FaceClassifier.hxx>
20 #include <Geom_BezierSurface.hxx>
21 #include <Geom_BSplineSurface.hxx>
22 #include <gp_Pln.hxx>
23 #include <gp_Pnt.hxx>
24 #include <gp_Vec.hxx>
25 #include <GProp_PEquation.hxx>
26 #include <HLRAlgo_Projector.hxx>
27 #include <HLRBRep_BSurfaceTool.hxx>
28 #include <HLRBRep_Curve.hxx>
29 #include <HLRBRep_Surface.hxx>
30 #include <Standard_DomainError.hxx>
31 #include <Standard_NoSuchObject.hxx>
32 #include <Standard_OutOfRange.hxx>
33 #include <TColgp_Array1OfPnt.hxx>
34 #include <TColStd_Array2OfReal.hxx>
35 #include <TopoDS_Face.hxx>
36
37 //=======================================================================
38 //function : HLRBRep_Surface
39 //purpose  : 
40 //=======================================================================
41 HLRBRep_Surface::HLRBRep_Surface ()
42 {
43 }
44
45 //=======================================================================
46 //function : Surface
47 //purpose  : 
48 //=======================================================================
49
50 void HLRBRep_Surface::Surface (const TopoDS_Face& F)
51 {
52   //mySurf.Initialize(F,Standard_False);
53   mySurf.Initialize(F,Standard_True);
54   GeomAbs_SurfaceType typ = HLRBRep_BSurfaceTool::GetType(mySurf);
55   switch (typ) {
56
57   case GeomAbs_Plane :
58   case GeomAbs_Cylinder :
59   case GeomAbs_Cone :
60   case GeomAbs_Sphere :
61   case GeomAbs_Torus :
62     // unchanged type
63     myType = typ;
64     break;
65     
66   case GeomAbs_BezierSurface :
67     if (HLRBRep_BSurfaceTool::UDegree(mySurf) == 1 &&
68         HLRBRep_BSurfaceTool::VDegree(mySurf) == 1) {
69       myType = GeomAbs_Plane;
70     }
71     else
72       myType = typ;
73     break;
74     
75     default :
76     myType = GeomAbs_OtherSurface;
77     break;
78   }
79 }
80
81 //=======================================================================
82 //function : SideRowsOfPoles
83 //purpose  : 
84 //=======================================================================
85
86 Standard_Boolean 
87 HLRBRep_Surface::SideRowsOfPoles (const Standard_Real tol,
88                                   const Standard_Integer nbuPoles,
89                                   const Standard_Integer nbvPoles,
90                                   TColgp_Array2OfPnt& Pnt) const
91 {
92   Standard_Integer iu,iv;
93   Standard_Real x0,y0,x,y,z;
94   Standard_Boolean result;
95   Standard_Real tole = (Standard_Real)tol;
96   const gp_Trsf& T = ((HLRAlgo_Projector*) myProj)->Transformation();
97
98   for (iu = 1; iu <= nbuPoles; iu++) {
99     
100     for (iv = 1; iv <= nbvPoles; iv++)
101       Pnt(iu,iv).Transform(T);
102   }
103   result = Standard_True;
104
105   for (iu = 1; iu <= nbuPoles && result; iu++) {         // Side iso u ?
106     Pnt(iu,1).Coord(x0,y0,z);
107     
108     for (iv = 2; iv <= nbvPoles && result; iv++) {
109       Pnt(iu,iv).Coord(x,y,z);
110       result = Abs(x-x0) < tole && Abs(y-y0) < tole;
111     }
112   }
113   if (result) return result;
114   result = Standard_True;
115   
116   for (iv = 1; iv <= nbvPoles && result; iv++) {         // Side iso v ?
117     Pnt(1,iv).Coord(x0,y0,z);
118     
119     for (iu = 2; iu <= nbuPoles && result; iu++) {
120       Pnt(iu,iv).Coord(x,y,z);
121       result = Abs(x-x0) < tole && Abs(y-y0) < tole;
122     }
123   }
124   if (result) return result;
125
126   // Are the Poles in a Side Plane ?
127   TColgp_Array1OfPnt p(1,nbuPoles*nbvPoles);
128   Standard_Integer i = 0;
129
130   for (iu = 1; iu <= nbuPoles; iu++) {
131     
132     for (iv = 1; iv <= nbvPoles; iv++) { 
133       i++;
134       p(i) = Pnt(iu,iv);
135     }
136   }
137
138   GProp_PEquation Pl(p,(Standard_Real)tol);
139   if (Pl.IsPlanar())
140     result = Abs(Pl.Plane().Axis().Direction().Z()) < 0.0001;
141
142   return result;
143 }
144
145 //=======================================================================
146 //function : IsSide
147 //purpose  : 
148 //=======================================================================
149
150 Standard_Boolean 
151 HLRBRep_Surface::IsSide (const Standard_Real tolF,
152                          const Standard_Real toler) const
153 {
154   gp_Pnt Pt;
155   gp_Vec D;
156   Standard_Real r;
157
158   if (myType == GeomAbs_Plane) {
159     gp_Pln Pl = Plane();
160     gp_Ax1 A  = Pl.Axis();
161     Pt = A.Location();
162     D  = A.Direction();
163     Pt.Transform(((HLRAlgo_Projector*) myProj)->Transformation());
164     D .Transform(((HLRAlgo_Projector*) myProj)->Transformation());
165     if (((HLRAlgo_Projector*) myProj)->Perspective()) {
166       r = D.Z() * ((HLRAlgo_Projector*) myProj)->Focus() - 
167         ( D.X() * Pt.X() + D.Y() * Pt.Y() + D.Z() * Pt.Z() );
168     }
169     else r= D.Z();
170     return Abs(r) < toler;
171   }
172   else if (myType == GeomAbs_Cylinder) {
173     if (((HLRAlgo_Projector*) myProj)->Perspective()) return Standard_False;
174     gp_Cylinder Cyl = HLRBRep_BSurfaceTool::Cylinder(mySurf);
175     gp_Ax1 A = Cyl.Axis();
176     D  = A.Direction();
177     D .Transform(((HLRAlgo_Projector*) myProj)->Transformation());
178     r = Sqrt(D.X() * D.X() + D.Y() * D.Y());
179     return r < toler;
180   }
181   else if (myType == GeomAbs_Cone) {
182     if (!((HLRAlgo_Projector*) myProj)->Perspective()) return Standard_False;
183     gp_Cone Con = HLRBRep_BSurfaceTool::Cone(mySurf);
184     Pt = Con.Apex();
185     Pt.Transform(((HLRAlgo_Projector*) myProj)->Transformation());
186     Standard_Real tol = 0.001;
187     return Pt.IsEqual(gp_Pnt(0,0,((HLRAlgo_Projector*) myProj)->Focus()),tol);
188   }
189   else if (myType == GeomAbs_BezierSurface) {
190     if (((HLRAlgo_Projector*) myProj)->Perspective()) return Standard_False;
191     Standard_Integer nu = HLRBRep_BSurfaceTool::NbUPoles(mySurf);
192     Standard_Integer nv = HLRBRep_BSurfaceTool::NbVPoles(mySurf);
193     TColgp_Array2OfPnt Pnt(1,nu,1,nv);
194     HLRBRep_BSurfaceTool::Bezier(mySurf)->Poles(Pnt);
195     return SideRowsOfPoles (tolF,nu,nv,Pnt);
196   }
197   else if (myType == GeomAbs_BSplineSurface) {
198     if (((HLRAlgo_Projector*) myProj)->Perspective()) return Standard_False;
199     Standard_Integer nu = HLRBRep_BSurfaceTool::NbUPoles(mySurf);
200     Standard_Integer nv = HLRBRep_BSurfaceTool::NbVPoles(mySurf);
201     TColgp_Array2OfPnt Pnt(1,nu,1,nv);
202     TColStd_Array2OfReal W(1,nu,1,nv);
203     HLRBRep_BSurfaceTool::BSpline(mySurf)->Poles(Pnt);
204     HLRBRep_BSurfaceTool::BSpline(mySurf)->Weights(W);
205     return SideRowsOfPoles (tolF,nu,nv,Pnt);
206   }
207   else return Standard_False;
208 }
209
210 //=======================================================================
211 //function : IsAbove
212 //purpose  : 
213 //=======================================================================
214
215 Standard_Boolean  
216 HLRBRep_Surface::IsAbove (const Standard_Boolean back,
217                           const Standard_Address A,
218                           const Standard_Real tol) const
219
220   Standard_Boolean planar = (myType == GeomAbs_Plane);
221   if (planar) {
222     gp_Pln Pl = Plane();
223     Standard_Real a,b,c,d;
224     Pl.Coefficients(a,b,c,d);
225     Standard_Real u,u1,u2,dd,x,y,z;
226     gp_Pnt P;
227     u1 = ((HLRBRep_Curve*)A)->Parameter3d
228       (((HLRBRep_Curve*)A)->FirstParameter());
229     u2 = ((HLRBRep_Curve*)A)->Parameter3d
230       (((HLRBRep_Curve*)A)->LastParameter());
231     u=u1;
232     ((HLRBRep_Curve*)A)->D0(u,P);
233     P.Coord(x,y,z);
234     dd = a*x + b*y + c*z + d;
235     if (back) dd = -dd;
236     if (dd < -tol) return Standard_False;
237     if (((HLRBRep_Curve*)A)->GetType() != GeomAbs_Line) {
238       Standard_Integer nbPnt = 30;
239       Standard_Real step = (u2-u1)/(nbPnt+1);
240       for (Standard_Integer i = 1; i <= nbPnt; i++) {
241         u += step;
242         ((HLRBRep_Curve*)A)->D0(u,P);
243         P.Coord(x,y,z);
244         dd = a*x + b*y + c*z + d;
245         if (back) dd = -dd;
246         if (dd < -tol) return Standard_False;
247       }
248     }
249     u = u2;
250     ((HLRBRep_Curve*)A)->D0(u,P);
251     P.Coord(x,y,z);
252     dd = a*x + b*y + c*z + d;
253     if (back) dd = -dd;
254     if (dd < -tol) return Standard_False;
255     return Standard_True;
256   }
257   else return Standard_False; 
258 }
259
260 //=======================================================================
261 //function : Value
262 //purpose  : 
263 //=======================================================================
264
265 gp_Pnt HLRBRep_Surface::Value (const Standard_Real U,
266                                const Standard_Real V) const
267
268   gp_Pnt P;
269   D0(U,V,P);
270   return P;
271 }
272
273 //=======================================================================
274 //function : Plane
275 //purpose  : 
276 //=======================================================================
277
278 gp_Pln  HLRBRep_Surface::Plane () const 
279 {
280   GeomAbs_SurfaceType typ = HLRBRep_BSurfaceTool::GetType(mySurf);
281   switch (typ) {
282   case GeomAbs_BezierSurface :
283     {
284       gp_Pnt P;
285       gp_Vec D1U;
286       gp_Vec D1V;
287       D1(0.5,0.5,P,D1U,D1V);
288       return gp_Pln(P,gp_Dir(D1U.Crossed(D1V)));
289     }
290     
291     default :
292     return HLRBRep_BSurfaceTool::Plane(mySurf);
293   }
294 }