[occt.git] / src / OpenGl / OpenGl_BVHTreeSelector.cxx

// Created on: 2013-12-25
// Created by: Varvara POSKONINA
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#include <OpenGl_BVHTreeSelector.hxx>
#include <OpenGl_BVHClipPrimitiveSet.hxx>

#include <vector>
#include <algorithm>

// =======================================================================
// function : DotProduct
// purpose  : Calculates a dot product of 4-dimensional vectors in homogeneous coordinates
// =======================================================================
static Standard_ShortReal DotProduct (const OpenGl_Vec4& theA,
                                      const OpenGl_Vec4& theB)
{
  return theA.x() * theB.x() + theA.y() * theB.y() + theA.z() * theB.z();
}

// =======================================================================
// function : BinarySign
// purpose  :
// =======================================================================
static OpenGl_Vec4 BinarySign (const OpenGl_Vec4& theVec)
{
  return OpenGl_Vec4 (theVec.x() > 0.0f ? 1.0f : 0.0f,
                      theVec.y() > 0.0f ? 1.0f : 0.0f,
                      theVec.z() > 0.0f ? 1.0f : 0.0f,
                      theVec.w() > 0.0f ? 1.0f : 0.0f);
}

// =======================================================================
// function : InversedBinarySign
// purpose  :
// =======================================================================
static OpenGl_Vec4 InversedBinarySign (const OpenGl_Vec4& theVec)
{
  return OpenGl_Vec4 (theVec.x() > 0.0f ? 0.0f : 1.0f,
                      theVec.y() > 0.0f ? 0.0f : 1.0f,
                      theVec.z() > 0.0f ? 0.0f : 1.0f,
                      theVec.w() > 0.0f ? 0.0f : 1.0f);
}

// =======================================================================
// function : OpenGl_BVHTreeSelector
// purpose  :
// =======================================================================
OpenGl_BVHTreeSelector::OpenGl_BVHTreeSelector()
: myIsProjectionParallel (Standard_True),
  myProjectionState      (0),
  myModelViewState       (0)
{
  //
}

// =======================================================================
// function : SetClipVolume
// purpose  : Retrieves view volume's planes equations and its vertices from projection and modelview matrices.
// =======================================================================
void OpenGl_BVHTreeSelector::SetViewVolume (const Handle(Graphic3d_Camera)& theCamera)
{
  myIsProjectionParallel = theCamera->IsOrthographic();
  const OpenGl_Mat4& aProjMat  = theCamera->ProjectionMatrixF();
  const OpenGl_Mat4& aModelMat = theCamera->OrientationMatrixF();

  Standard_ShortReal nLeft = 0.0f, nRight = 0.0f, nTop = 0.0f, nBottom = 0.0f;
  Standard_ShortReal fLeft = 0.0f, fRight = 0.0f, fTop = 0.0f, fBottom = 0.0f;
  Standard_ShortReal aNear = 0.0f, aFar   = 0.0f;
  if (!myIsProjectionParallel)
  {
    // handle perspective projection
    aNear   = aProjMat.GetValue (2, 3) / (- 1.0f + aProjMat.GetValue (2, 2));
    aFar    = aProjMat.GetValue (2, 3) / (  1.0f + aProjMat.GetValue (2, 2));
    // Near plane
    nLeft   = aNear * (aProjMat.GetValue (0, 2) - 1.0f) / aProjMat.GetValue (0, 0);
    nRight  = aNear * (aProjMat.GetValue (0, 2) + 1.0f) / aProjMat.GetValue (0, 0);
    nTop    = aNear * (aProjMat.GetValue (1, 2) + 1.0f) / aProjMat.GetValue (1, 1);
    nBottom = aNear * (aProjMat.GetValue (1, 2) - 1.0f) / aProjMat.GetValue (1, 1);
    // Far plane
    fLeft   = aFar  * (aProjMat.GetValue (0, 2) - 1.0f) / aProjMat.GetValue (0, 0);
    fRight  = aFar  * (aProjMat.GetValue (0, 2) + 1.0f) / aProjMat.GetValue (0, 0);
    fTop    = aFar  * (aProjMat.GetValue (1, 2) + 1.0f) / aProjMat.GetValue (1, 1);
    fBottom = aFar  * (aProjMat.GetValue (1, 2) - 1.0f) / aProjMat.GetValue (1, 1);
  }
  else
  {
    // handle orthographic projection
    aNear   = (1.0f / aProjMat.GetValue (2, 2)) * (aProjMat.GetValue (2, 3) + 1.0f);
    aFar    = (1.0f / aProjMat.GetValue (2, 2)) * (aProjMat.GetValue (2, 3) - 1.0f);
    // Near plane
    nLeft   = ( 1.0f + aProjMat.GetValue (0, 3)) / (-aProjMat.GetValue (0, 0));
    fLeft   = nLeft;
    nRight  = ( 1.0f - aProjMat.GetValue (0, 3)) /   aProjMat.GetValue (0, 0);
    fRight  = nRight;
    nTop    = ( 1.0f - aProjMat.GetValue (1, 3)) /   aProjMat.GetValue (1, 1);
    fTop    = nTop;
    nBottom = (-1.0f - aProjMat.GetValue (1, 3)) /   aProjMat.GetValue (1, 1);
    fBottom = nBottom;
  }

  OpenGl_Vec4 aLeftTopNear     (nLeft,  nTop,    -aNear, 1.0f), aRightBottomFar (fRight, fBottom, -aFar, 1.0f);
  OpenGl_Vec4 aLeftBottomNear  (nLeft,  nBottom, -aNear, 1.0f), aRightTopFar    (fRight, fTop,    -aFar, 1.0f);
  OpenGl_Vec4 aRightBottomNear (nRight, nBottom, -aNear, 1.0f), aLeftTopFar     (fLeft,  fTop,    -aFar, 1.0f);
  OpenGl_Vec4 aRightTopNear    (nRight, nTop,    -aNear, 1.0f), aLeftBottomFar  (fLeft,  fBottom, -aFar, 1.0f);

  const OpenGl_Mat4 aViewProj = aModelMat * aProjMat;
  OpenGl_Mat4 anInvModelView;
  aModelMat.Inverted(anInvModelView);

  myClipVerts[ClipVert_LeftTopNear]     = anInvModelView * aLeftTopNear;
  myClipVerts[ClipVert_RightBottomFar]  = anInvModelView * aRightBottomFar;
  myClipVerts[ClipVert_LeftBottomNear]  = anInvModelView * aLeftBottomNear;
  myClipVerts[ClipVert_RightTopFar]     = anInvModelView * aRightTopFar;
  myClipVerts[ClipVert_RightBottomNear] = anInvModelView * aRightBottomNear;
  myClipVerts[ClipVert_LeftTopFar]      = anInvModelView * aLeftTopFar;
  myClipVerts[ClipVert_RightTopNear]    = anInvModelView * aRightTopNear;
  myClipVerts[ClipVert_LeftBottomFar]   = anInvModelView * aLeftBottomFar;

  // UNNORMALIZED!
  myClipPlanes[Plane_Left]   = aViewProj.GetRow (3) + aViewProj.GetRow (0);
  myClipPlanes[Plane_Right]  = aViewProj.GetRow (3) - aViewProj.GetRow (0);
  myClipPlanes[Plane_Top]    = aViewProj.GetRow (3) - aViewProj.GetRow (1);
  myClipPlanes[Plane_Bottom] = aViewProj.GetRow (3) + aViewProj.GetRow (1);
  myClipPlanes[Plane_Near]   = aViewProj.GetRow (3) + aViewProj.GetRow (2);
  myClipPlanes[Plane_Far]    = aViewProj.GetRow (3) - aViewProj.GetRow (2);

  gp_Pnt aPtCenter = theCamera->Center();
  OpenGl_Vec4 aCenter (static_cast<Standard_ShortReal> (aPtCenter.X()),
                       static_cast<Standard_ShortReal> (aPtCenter.Y()),
                       static_cast<Standard_ShortReal> (aPtCenter.Z()),
                       1.0f);

  for (Standard_Integer aPlaneIter = 0; aPlaneIter < PlanesNB; ++aPlaneIter)
  {
    OpenGl_Vec4 anEq = myClipPlanes[aPlaneIter];
    if (SignedPlanePointDistance (anEq, aCenter) > 0)
    {
      anEq *= -1.0f;
      myClipPlanes[aPlaneIter] = anEq;
    }
  }
}

// =======================================================================
// function : SignedPlanePointDistance
// purpose  :
// =======================================================================
Standard_ShortReal OpenGl_BVHTreeSelector::SignedPlanePointDistance (const OpenGl_Vec4& theNormal,
                                                                     const OpenGl_Vec4& thePnt)
{
  const Standard_ShortReal aNormLength = std::sqrt (theNormal.x() * theNormal.x()
                                                  + theNormal.y() * theNormal.y()
                                                  + theNormal.z() * theNormal.z());
  const Standard_ShortReal anInvNormLength = 1.0f / aNormLength;
  const Standard_ShortReal aD  = theNormal.w() * anInvNormLength;
  const Standard_ShortReal anA = theNormal.x() * anInvNormLength;
  const Standard_ShortReal aB  = theNormal.y() * anInvNormLength;
  const Standard_ShortReal aC  = theNormal.z() * anInvNormLength;
  return aD + (anA * thePnt.x() + aB * thePnt.y() + aC * thePnt.z());
}

// =======================================================================
// function : CacheClipPtsProjections
// purpose  : Caches view volume's vertices projections along its normals and AABBs dimensions
//            Must be called at the beginning of each BVH tree traverse loop
// =======================================================================
void OpenGl_BVHTreeSelector::CacheClipPtsProjections()
{
  Standard_ShortReal aProjectedVerts[ClipVerticesNB];
  for (Standard_Integer aPlaneIter = 0; aPlaneIter < PlanesNB; ++aPlaneIter)
  {
    const OpenGl_Vec4 aPlane = myClipPlanes[aPlaneIter];
    for (Standard_Integer aCornerIter = 0; aCornerIter < ClipVerticesNB; ++aCornerIter)
    {
      Standard_ShortReal aProjection = DotProduct (aPlane, myClipVerts[aCornerIter]);
      aProjectedVerts[aCornerIter] = aProjection;
    }
    myMaxClipProjectionPts[aPlaneIter] = *std::max_element (aProjectedVerts, aProjectedVerts + ClipVerticesNB);
    myMinClipProjectionPts[aPlaneIter] = *std::min_element (aProjectedVerts, aProjectedVerts + ClipVerticesNB);
  }

  OpenGl_Vec4 aDimensions[3] =
  {
    OpenGl_Vec4 (1.0f, 0.0f, 0.0f, 1.0f),
    OpenGl_Vec4 (0.0f, 1.0f, 0.0f, 1.0f),
    OpenGl_Vec4 (0.0f, 0.0f, 1.0f, 1.0f)
  };

  for (Standard_Integer aDim = 0; aDim < 3; ++aDim)
  {
    for (Standard_Integer aCornerIter = 0; aCornerIter < ClipVerticesNB; ++aCornerIter)
    {
      Standard_ShortReal aProjection = DotProduct (aDimensions[aDim], myClipVerts[aCornerIter]);
      aProjectedVerts[aCornerIter] = aProjection;
    }
    myMaxOrthoProjectionPts[aDim] = *std::max_element (aProjectedVerts, aProjectedVerts + ClipVerticesNB);
    myMinOrthoProjectionPts[aDim] = *std::min_element (aProjectedVerts, aProjectedVerts + ClipVerticesNB);
  }
}

// =======================================================================
// function : Intersect
// purpose  : Detects if AABB overlaps view volume using separating axis theorem (SAT)
// =======================================================================
Standard_Boolean OpenGl_BVHTreeSelector::Intersect (const OpenGl_Vec4& theMinPt,
                                                    const OpenGl_Vec4& theMaxPt) const
{
  //     E1
  //    |_ E0
  //   /
  //    E2
  const OpenGl_Vec4  aShiftedBoxMax  = theMaxPt - theMinPt;
  Standard_ShortReal aBoxProjMax     = 0.0f, aBoxProjMin     = 0.0f;
  Standard_ShortReal aFrustumProjMax = 0.0f, aFrustumProjMin = 0.0f;

  // E0 test
  aBoxProjMax = aShiftedBoxMax.x();
  aFrustumProjMax = myMaxOrthoProjectionPts[0] - DotProduct (OpenGl_Vec4 (1.0f, 0.0f, 0.0f, 1.0f), theMinPt);
  aFrustumProjMin = myMinOrthoProjectionPts[0] - DotProduct (OpenGl_Vec4 (1.0f, 0.0f, 0.0f, 1.0f), theMinPt);
  if (aBoxProjMin > aFrustumProjMax
   || aBoxProjMax < aFrustumProjMin)
  {
    return Standard_False;
  }

  // E1 test
  aBoxProjMax = aShiftedBoxMax.y();
  aFrustumProjMax = myMaxOrthoProjectionPts[1] - DotProduct (OpenGl_Vec4 (0.0f, 1.0f, 0.0f, 1.0f), theMinPt);
  aFrustumProjMin = myMinOrthoProjectionPts[1] - DotProduct (OpenGl_Vec4 (0.0f, 1.0f, 0.0f, 1.0f), theMinPt);
  if (aBoxProjMin > aFrustumProjMax
   || aBoxProjMax < aFrustumProjMin)
  {
    return Standard_False;
  }

  // E2 test
  aBoxProjMax = aShiftedBoxMax.z();
  aFrustumProjMax = myMaxOrthoProjectionPts[2] - DotProduct (OpenGl_Vec4 (0.0f, 0.0f, 1.0f, 1.0f), theMinPt);
  aFrustumProjMin = myMinOrthoProjectionPts[2] - DotProduct (OpenGl_Vec4 (0.0f, 0.0f, 1.0f, 1.0f), theMinPt);
  if (aBoxProjMin > aFrustumProjMax
   || aBoxProjMax < aFrustumProjMin)
  {
    return Standard_False;
  }

  const Standard_Integer anIncFactor = myIsProjectionParallel ? 2 : 1;
  for (Standard_Integer aPlaneIter = 0; aPlaneIter < 5; aPlaneIter += anIncFactor)
  {
    OpenGl_Vec4 aPlane = myClipPlanes[aPlaneIter];
    OpenGl_Vec4 aPt1 (0.0f), aPt2 (0.0f);
    aPt1 = BinarySign (aPlane) * aShiftedBoxMax;
    aBoxProjMax = DotProduct (aPlane, aPt1);
    aFrustumProjMax = myMaxClipProjectionPts[aPlaneIter] - DotProduct (aPlane, theMinPt);
    aFrustumProjMin = myMinClipProjectionPts[aPlaneIter] - DotProduct (aPlane, theMinPt);
    if (aFrustumProjMin < aBoxProjMax
     && aBoxProjMax < aFrustumProjMax)
    {
      continue;
    }

    aPt2 = InversedBinarySign (aPlane) * aShiftedBoxMax;
    aBoxProjMin = DotProduct (aPlane, aPt2);
    if (aBoxProjMin > aFrustumProjMax
     || aBoxProjMax < aFrustumProjMin)
    {
      return Standard_False;
    }
  }

  return Standard_True;
}
Commit	Line	Data
b7cd4ba7	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	const Standard_ShortReal anInvNormLength = 1.0f / aNormLength;
167	const Standard_ShortReal aD = theNormal.w() * anInvNormLength;
168	const Standard_ShortReal anA = theNormal.x() * anInvNormLength;
169	const Standard_ShortReal aB = theNormal.y() * anInvNormLength;
170	const Standard_ShortReal aC = theNormal.z() * anInvNormLength;
171	return aD + (anA * thePnt.x() + aB * thePnt.y() + aC * thePnt.z());
172	}
173
174	// =======================================================================
175	// function : CacheClipPtsProjections
176	// purpose : Caches view volume's vertices projections along its normals and AABBs dimensions
177	// Must be called at the beginning of each BVH tree traverse loop
178	// =======================================================================
179	void OpenGl_BVHTreeSelector::CacheClipPtsProjections()
180	{
181	Standard_ShortReal aProjectedVerts[ClipVerticesNB];
182	for (Standard_Integer aPlaneIter = 0; aPlaneIter < PlanesNB; ++aPlaneIter)
183	{
184	const OpenGl_Vec4 aPlane = myClipPlanes[aPlaneIter];
185	for (Standard_Integer aCornerIter = 0; aCornerIter < ClipVerticesNB; ++aCornerIter)
186	{
187	Standard_ShortReal aProjection = DotProduct (aPlane, myClipVerts[aCornerIter]);
188	aProjectedVerts[aCornerIter] = aProjection;
189	}
190	myMaxClipProjectionPts[aPlaneIter] = *std::max_element (aProjectedVerts, aProjectedVerts + ClipVerticesNB);
191	myMinClipProjectionPts[aPlaneIter] = *std::min_element (aProjectedVerts, aProjectedVerts + ClipVerticesNB);
192	}
193
194	OpenGl_Vec4 aDimensions[3] =
195	{
196	OpenGl_Vec4 (1.0f, 0.0f, 0.0f, 1.0f),
197	OpenGl_Vec4 (0.0f, 1.0f, 0.0f, 1.0f),
198	OpenGl_Vec4 (0.0f, 0.0f, 1.0f, 1.0f)
199	};
200
201	for (Standard_Integer aDim = 0; aDim < 3; ++aDim)
202	{
203	for (Standard_Integer aCornerIter = 0; aCornerIter < ClipVerticesNB; ++aCornerIter)
204	{
205	Standard_ShortReal aProjection = DotProduct (aDimensions[aDim], myClipVerts[aCornerIter]);
206	aProjectedVerts[aCornerIter] = aProjection;
207	}
208	myMaxOrthoProjectionPts[aDim] = *std::max_element (aProjectedVerts, aProjectedVerts + ClipVerticesNB);
209	myMinOrthoProjectionPts[aDim] = *std::min_element (aProjectedVerts, aProjectedVerts + ClipVerticesNB);
210	}
211	}
212
213	// =======================================================================
214	// function : Intersect
215	// purpose : Detects if AABB overlaps view volume using separating axis theorem (SAT)
216	// =======================================================================
217	Standard_Boolean OpenGl_BVHTreeSelector::Intersect (const OpenGl_Vec4& theMinPt,
218	const OpenGl_Vec4& theMaxPt) const
219	{
220	// E1
221	// \|_ E0
222	// /
223	// E2
224	const OpenGl_Vec4 aShiftedBoxMax = theMaxPt - theMinPt;
225	Standard_ShortReal aBoxProjMax = 0.0f, aBoxProjMin = 0.0f;
226	Standard_ShortReal aFrustumProjMax = 0.0f, aFrustumProjMin = 0.0f;
227
228	// E0 test
229	aBoxProjMax = aShiftedBoxMax.x();
230	aFrustumProjMax = myMaxOrthoProjectionPts[0] - DotProduct (OpenGl_Vec4 (1.0f, 0.0f, 0.0f, 1.0f), theMinPt);
231	aFrustumProjMin = myMinOrthoProjectionPts[0] - DotProduct (OpenGl_Vec4 (1.0f, 0.0f, 0.0f, 1.0f), theMinPt);
232	if (aBoxProjMin > aFrustumProjMax
233	\|\| aBoxProjMax < aFrustumProjMin)
234	{
235	return Standard_False;
236	}
237
238	// E1 test
239	aBoxProjMax = aShiftedBoxMax.y();
240	aFrustumProjMax = myMaxOrthoProjectionPts[1] - DotProduct (OpenGl_Vec4 (0.0f, 1.0f, 0.0f, 1.0f), theMinPt);
241	aFrustumProjMin = myMinOrthoProjectionPts[1] - DotProduct (OpenGl_Vec4 (0.0f, 1.0f, 0.0f, 1.0f), theMinPt);
242	if (aBoxProjMin > aFrustumProjMax
243	\|\| aBoxProjMax < aFrustumProjMin)
244	{
245	return Standard_False;
246	}
247
248	// E2 test
249	aBoxProjMax = aShiftedBoxMax.z();
250	aFrustumProjMax = myMaxOrthoProjectionPts[2] - DotProduct (OpenGl_Vec4 (0.0f, 0.0f, 1.0f, 1.0f), theMinPt);
251	aFrustumProjMin = myMinOrthoProjectionPts[2] - DotProduct (OpenGl_Vec4 (0.0f, 0.0f, 1.0f, 1.0f), theMinPt);
252	if (aBoxProjMin > aFrustumProjMax
253	\|\| aBoxProjMax < aFrustumProjMin)
254	{
255	return Standard_False;
256	}
257
258	const Standard_Integer anIncFactor = myIsProjectionParallel ? 2 : 1;
259	for (Standard_Integer aPlaneIter = 0; aPlaneIter < 5; aPlaneIter += anIncFactor)
260	{
261	OpenGl_Vec4 aPlane = myClipPlanes[aPlaneIter];
262	OpenGl_Vec4 aPt1 (0.0f), aPt2 (0.0f);
263	aPt1 = BinarySign (aPlane) * aShiftedBoxMax;
264	aBoxProjMax = DotProduct (aPlane, aPt1);
265	aFrustumProjMax = myMaxClipProjectionPts[aPlaneIter] - DotProduct (aPlane, theMinPt);
266	aFrustumProjMin = myMinClipProjectionPts[aPlaneIter] - DotProduct (aPlane, theMinPt);
267	if (aFrustumProjMin < aBoxProjMax
268	&& aBoxProjMax < aFrustumProjMax)
269	{
270	continue;
271	}
272
273	aPt2 = InversedBinarySign (aPlane) * aShiftedBoxMax;
274	aBoxProjMin = DotProduct (aPlane, aPt2);
275	if (aBoxProjMin > aFrustumProjMax
276	\|\| aBoxProjMax < aFrustumProjMin)
277	{
278	return Standard_False;
279	}
280	}
281
282	return Standard_True;
283	}