0025442: Visualization, TKOpenGl - prevent inclusion of system header glxext.h
[occt.git] / src / OpenGl / OpenGl_BVHTreeSelector.cxx
1 // Created on: 2013-12-25
2 // Created by: Varvara POSKONINA
3 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 //
5 // This file is part of Open CASCADE Technology software library.
6 //
7 // This library is free software; you can redistribute it and/or modify it under
8 // the terms of the GNU Lesser General Public License version 2.1 as published
9 // by the Free Software Foundation, with special exception defined in the file
10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11 // distribution for complete text of the license and disclaimer of any warranty.
12 //
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
15
16 #include <OpenGl_BVHTreeSelector.hxx>
17 #include <OpenGl_BVHClipPrimitiveSet.hxx>
18
19 #include <vector>
20 #include <algorithm>
21
22 // =======================================================================
23 // function : DotProduct
24 // purpose  : Calculates a dot product of 4-dimensional vectors in homogeneous coordinates
25 // =======================================================================
26 static Standard_ShortReal DotProduct (const OpenGl_Vec4& theA,
27                                       const OpenGl_Vec4& theB)
28 {
29   return theA.x() * theB.x() + theA.y() * theB.y() + theA.z() * theB.z();
30 }
31
32 // =======================================================================
33 // function : BinarySign
34 // purpose  :
35 // =======================================================================
36 static OpenGl_Vec4 BinarySign (const OpenGl_Vec4& theVec)
37 {
38   return OpenGl_Vec4 (theVec.x() > 0.0f ? 1.0f : 0.0f,
39                       theVec.y() > 0.0f ? 1.0f : 0.0f,
40                       theVec.z() > 0.0f ? 1.0f : 0.0f,
41                       theVec.w() > 0.0f ? 1.0f : 0.0f);
42 }
43
44 // =======================================================================
45 // function : InversedBinarySign
46 // purpose  :
47 // =======================================================================
48 static OpenGl_Vec4 InversedBinarySign (const OpenGl_Vec4& theVec)
49 {
50   return OpenGl_Vec4 (theVec.x() > 0.0f ? 0.0f : 1.0f,
51                       theVec.y() > 0.0f ? 0.0f : 1.0f,
52                       theVec.z() > 0.0f ? 0.0f : 1.0f,
53                       theVec.w() > 0.0f ? 0.0f : 1.0f);
54 }
55
56 // =======================================================================
57 // function : OpenGl_BVHTreeSelector
58 // purpose  :
59 // =======================================================================
60 OpenGl_BVHTreeSelector::OpenGl_BVHTreeSelector()
61 : myIsProjectionParallel (Standard_True),
62   myProjectionState      (0),
63   myModelViewState       (0)
64 {
65   //
66 }
67
68 // =======================================================================
69 // function : SetClipVolume
70 // purpose  : Retrieves view volume's planes equations and its vertices from projection and modelview matrices.
71 // =======================================================================
72 void OpenGl_BVHTreeSelector::SetViewVolume (const Handle(Graphic3d_Camera)& theCamera)
73 {
74   myIsProjectionParallel = theCamera->IsOrthographic();
75   const OpenGl_Mat4& aProjMat  = theCamera->ProjectionMatrixF();
76   const OpenGl_Mat4& aModelMat = theCamera->OrientationMatrixF();
77
78   Standard_ShortReal nLeft = 0.0f, nRight = 0.0f, nTop = 0.0f, nBottom = 0.0f;
79   Standard_ShortReal fLeft = 0.0f, fRight = 0.0f, fTop = 0.0f, fBottom = 0.0f;
80   Standard_ShortReal aNear = 0.0f, aFar   = 0.0f;
81   if (!myIsProjectionParallel)
82   {
83     // handle perspective projection
84     aNear   = aProjMat.GetValue (2, 3) / (- 1.0f + aProjMat.GetValue (2, 2));
85     aFar    = aProjMat.GetValue (2, 3) / (  1.0f + aProjMat.GetValue (2, 2));
86     // Near plane
87     nLeft   = aNear * (aProjMat.GetValue (0, 2) - 1.0f) / aProjMat.GetValue (0, 0);
88     nRight  = aNear * (aProjMat.GetValue (0, 2) + 1.0f) / aProjMat.GetValue (0, 0);
89     nTop    = aNear * (aProjMat.GetValue (1, 2) + 1.0f) / aProjMat.GetValue (1, 1);
90     nBottom = aNear * (aProjMat.GetValue (1, 2) - 1.0f) / aProjMat.GetValue (1, 1);
91     // Far plane
92     fLeft   = aFar  * (aProjMat.GetValue (0, 2) - 1.0f) / aProjMat.GetValue (0, 0);
93     fRight  = aFar  * (aProjMat.GetValue (0, 2) + 1.0f) / aProjMat.GetValue (0, 0);
94     fTop    = aFar  * (aProjMat.GetValue (1, 2) + 1.0f) / aProjMat.GetValue (1, 1);
95     fBottom = aFar  * (aProjMat.GetValue (1, 2) - 1.0f) / aProjMat.GetValue (1, 1);
96   }
97   else
98   {
99     // handle orthographic projection
100     aNear   = (1.0f / aProjMat.GetValue (2, 2)) * (aProjMat.GetValue (2, 3) + 1.0f);
101     aFar    = (1.0f / aProjMat.GetValue (2, 2)) * (aProjMat.GetValue (2, 3) - 1.0f);
102     // Near plane
103     nLeft   = ( 1.0f + aProjMat.GetValue (0, 3)) / (-aProjMat.GetValue (0, 0));
104     fLeft   = nLeft;
105     nRight  = ( 1.0f - aProjMat.GetValue (0, 3)) /   aProjMat.GetValue (0, 0);
106     fRight  = nRight;
107     nTop    = ( 1.0f - aProjMat.GetValue (1, 3)) /   aProjMat.GetValue (1, 1);
108     fTop    = nTop;
109     nBottom = (-1.0f - aProjMat.GetValue (1, 3)) /   aProjMat.GetValue (1, 1);
110     fBottom = nBottom;
111   }
112
113   OpenGl_Vec4 aLeftTopNear     (nLeft,  nTop,    -aNear, 1.0f), aRightBottomFar (fRight, fBottom, -aFar, 1.0f);
114   OpenGl_Vec4 aLeftBottomNear  (nLeft,  nBottom, -aNear, 1.0f), aRightTopFar    (fRight, fTop,    -aFar, 1.0f);
115   OpenGl_Vec4 aRightBottomNear (nRight, nBottom, -aNear, 1.0f), aLeftTopFar     (fLeft,  fTop,    -aFar, 1.0f);
116   OpenGl_Vec4 aRightTopNear    (nRight, nTop,    -aNear, 1.0f), aLeftBottomFar  (fLeft,  fBottom, -aFar, 1.0f);
117
118   const OpenGl_Mat4 aViewProj = aModelMat * aProjMat;
119   OpenGl_Mat4 anInvModelView;
120   aModelMat.Inverted(anInvModelView);
121
122   myClipVerts[ClipVert_LeftTopNear]     = anInvModelView * aLeftTopNear;
123   myClipVerts[ClipVert_RightBottomFar]  = anInvModelView * aRightBottomFar;
124   myClipVerts[ClipVert_LeftBottomNear]  = anInvModelView * aLeftBottomNear;
125   myClipVerts[ClipVert_RightTopFar]     = anInvModelView * aRightTopFar;
126   myClipVerts[ClipVert_RightBottomNear] = anInvModelView * aRightBottomNear;
127   myClipVerts[ClipVert_LeftTopFar]      = anInvModelView * aLeftTopFar;
128   myClipVerts[ClipVert_RightTopNear]    = anInvModelView * aRightTopNear;
129   myClipVerts[ClipVert_LeftBottomFar]   = anInvModelView * aLeftBottomFar;
130
131   // UNNORMALIZED!
132   myClipPlanes[Plane_Left]   = aViewProj.GetRow (3) + aViewProj.GetRow (0);
133   myClipPlanes[Plane_Right]  = aViewProj.GetRow (3) - aViewProj.GetRow (0);
134   myClipPlanes[Plane_Top]    = aViewProj.GetRow (3) - aViewProj.GetRow (1);
135   myClipPlanes[Plane_Bottom] = aViewProj.GetRow (3) + aViewProj.GetRow (1);
136   myClipPlanes[Plane_Near]   = aViewProj.GetRow (3) + aViewProj.GetRow (2);
137   myClipPlanes[Plane_Far]    = aViewProj.GetRow (3) - aViewProj.GetRow (2);
138
139   gp_Pnt aPtCenter = theCamera->Center();
140   OpenGl_Vec4 aCenter (static_cast<Standard_ShortReal> (aPtCenter.X()),
141                        static_cast<Standard_ShortReal> (aPtCenter.Y()),
142                        static_cast<Standard_ShortReal> (aPtCenter.Z()),
143                        1.0f);
144
145   for (Standard_Integer aPlaneIter = 0; aPlaneIter < PlanesNB; ++aPlaneIter)
146   {
147     OpenGl_Vec4 anEq = myClipPlanes[aPlaneIter];
148     if (SignedPlanePointDistance (anEq, aCenter) > 0)
149     {
150       anEq *= -1.0f;
151       myClipPlanes[aPlaneIter] = anEq;
152     }
153   }
154 }
155
156 // =======================================================================
157 // function : SignedPlanePointDistance
158 // purpose  :
159 // =======================================================================
160 Standard_ShortReal OpenGl_BVHTreeSelector::SignedPlanePointDistance (const OpenGl_Vec4& theNormal,
161                                                                      const OpenGl_Vec4& thePnt)
162 {
163   const Standard_ShortReal aNormLength = std::sqrt (theNormal.x() * theNormal.x()
164                                                   + theNormal.y() * theNormal.y()
165                                                   + theNormal.z() * theNormal.z());
166
167   if (aNormLength < FLT_EPSILON)
168     return 0.0f;
169
170   const Standard_ShortReal anInvNormLength = 1.0f / aNormLength;
171   const Standard_ShortReal aD  = theNormal.w() * anInvNormLength;
172   const Standard_ShortReal anA = theNormal.x() * anInvNormLength;
173   const Standard_ShortReal aB  = theNormal.y() * anInvNormLength;
174   const Standard_ShortReal aC  = theNormal.z() * anInvNormLength;
175   return aD + (anA * thePnt.x() + aB * thePnt.y() + aC * thePnt.z());
176 }
177
178 // =======================================================================
179 // function : CacheClipPtsProjections
180 // purpose  : Caches view volume's vertices projections along its normals and AABBs dimensions
181 //            Must be called at the beginning of each BVH tree traverse loop
182 // =======================================================================
183 void OpenGl_BVHTreeSelector::CacheClipPtsProjections()
184 {
185   Standard_ShortReal aProjectedVerts[ClipVerticesNB];
186   for (Standard_Integer aPlaneIter = 0; aPlaneIter < PlanesNB; ++aPlaneIter)
187   {
188     const OpenGl_Vec4 aPlane = myClipPlanes[aPlaneIter];
189     for (Standard_Integer aCornerIter = 0; aCornerIter < ClipVerticesNB; ++aCornerIter)
190     {
191       Standard_ShortReal aProjection = DotProduct (aPlane, myClipVerts[aCornerIter]);
192       aProjectedVerts[aCornerIter] = aProjection;
193     }
194     myMaxClipProjectionPts[aPlaneIter] = *std::max_element (aProjectedVerts, aProjectedVerts + ClipVerticesNB);
195     myMinClipProjectionPts[aPlaneIter] = *std::min_element (aProjectedVerts, aProjectedVerts + ClipVerticesNB);
196   }
197
198   OpenGl_Vec4 aDimensions[3] =
199   {
200     OpenGl_Vec4 (1.0f, 0.0f, 0.0f, 1.0f),
201     OpenGl_Vec4 (0.0f, 1.0f, 0.0f, 1.0f),
202     OpenGl_Vec4 (0.0f, 0.0f, 1.0f, 1.0f)
203   };
204
205   for (Standard_Integer aDim = 0; aDim < 3; ++aDim)
206   {
207     for (Standard_Integer aCornerIter = 0; aCornerIter < ClipVerticesNB; ++aCornerIter)
208     {
209       Standard_ShortReal aProjection = DotProduct (aDimensions[aDim], myClipVerts[aCornerIter]);
210       aProjectedVerts[aCornerIter] = aProjection;
211     }
212     myMaxOrthoProjectionPts[aDim] = *std::max_element (aProjectedVerts, aProjectedVerts + ClipVerticesNB);
213     myMinOrthoProjectionPts[aDim] = *std::min_element (aProjectedVerts, aProjectedVerts + ClipVerticesNB);
214   }
215 }
216
217 // =======================================================================
218 // function : Intersect
219 // purpose  : Detects if AABB overlaps view volume using separating axis theorem (SAT)
220 // =======================================================================
221 Standard_Boolean OpenGl_BVHTreeSelector::Intersect (const OpenGl_Vec4& theMinPt,
222                                                     const OpenGl_Vec4& theMaxPt) const
223 {
224   //     E1
225   //    |_ E0
226   //   /
227   //    E2
228   const OpenGl_Vec4  aShiftedBoxMax  = theMaxPt - theMinPt;
229   Standard_ShortReal aBoxProjMax     = 0.0f, aBoxProjMin     = 0.0f;
230   Standard_ShortReal aFrustumProjMax = 0.0f, aFrustumProjMin = 0.0f;
231
232   // E0 test
233   aBoxProjMax = aShiftedBoxMax.x();
234   aFrustumProjMax = myMaxOrthoProjectionPts[0] - DotProduct (OpenGl_Vec4 (1.0f, 0.0f, 0.0f, 1.0f), theMinPt);
235   aFrustumProjMin = myMinOrthoProjectionPts[0] - DotProduct (OpenGl_Vec4 (1.0f, 0.0f, 0.0f, 1.0f), theMinPt);
236   if (aBoxProjMin > aFrustumProjMax
237    || aBoxProjMax < aFrustumProjMin)
238   {
239     return Standard_False;
240   }
241
242   // E1 test
243   aBoxProjMax = aShiftedBoxMax.y();
244   aFrustumProjMax = myMaxOrthoProjectionPts[1] - DotProduct (OpenGl_Vec4 (0.0f, 1.0f, 0.0f, 1.0f), theMinPt);
245   aFrustumProjMin = myMinOrthoProjectionPts[1] - DotProduct (OpenGl_Vec4 (0.0f, 1.0f, 0.0f, 1.0f), theMinPt);
246   if (aBoxProjMin > aFrustumProjMax
247    || aBoxProjMax < aFrustumProjMin)
248   {
249     return Standard_False;
250   }
251
252   // E2 test
253   aBoxProjMax = aShiftedBoxMax.z();
254   aFrustumProjMax = myMaxOrthoProjectionPts[2] - DotProduct (OpenGl_Vec4 (0.0f, 0.0f, 1.0f, 1.0f), theMinPt);
255   aFrustumProjMin = myMinOrthoProjectionPts[2] - DotProduct (OpenGl_Vec4 (0.0f, 0.0f, 1.0f, 1.0f), theMinPt);
256   if (aBoxProjMin > aFrustumProjMax
257    || aBoxProjMax < aFrustumProjMin)
258   {
259     return Standard_False;
260   }
261
262   const Standard_Integer anIncFactor = myIsProjectionParallel ? 2 : 1;
263   for (Standard_Integer aPlaneIter = 0; aPlaneIter < 5; aPlaneIter += anIncFactor)
264   {
265     OpenGl_Vec4 aPlane = myClipPlanes[aPlaneIter];
266     OpenGl_Vec4 aPt1 (0.0f), aPt2 (0.0f);
267     aPt1 = BinarySign (aPlane) * aShiftedBoxMax;
268     aBoxProjMax = DotProduct (aPlane, aPt1);
269     aFrustumProjMax = myMaxClipProjectionPts[aPlaneIter] - DotProduct (aPlane, theMinPt);
270     aFrustumProjMin = myMinClipProjectionPts[aPlaneIter] - DotProduct (aPlane, theMinPt);
271     if (aFrustumProjMin < aBoxProjMax
272      && aBoxProjMax < aFrustumProjMax)
273     {
274       continue;
275     }
276
277     aPt2 = InversedBinarySign (aPlane) * aShiftedBoxMax;
278     aBoxProjMin = DotProduct (aPlane, aPt2);
279     if (aBoxProjMin > aFrustumProjMax
280      || aBoxProjMax < aFrustumProjMin)
281     {
282       return Standard_False;
283     }
284   }
285
286   return Standard_True;
287 }