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 }