0027607: Visualization - Implement adaptive screen space sampling in path tracing
[occt.git] / src / Shaders / RaytraceBase.fs
index 2356906..6c5f5e4 100644 (file)
@@ -1,3 +1,8 @@
+#ifdef ADAPTIVE_SAMPLING
+  #extension GL_ARB_shader_image_load_store : require
+  #extension GL_NV_shader_atomic_float : require
+#endif
+
 #ifdef USE_TEXTURES
   #extension GL_ARB_bindless_texture : require
 #endif
@@ -66,11 +71,6 @@ uniform samplerBuffer uRaytraceLightSrcTexture;
 //! Environment map texture.
 uniform sampler2D uEnvironmentMapTexture;
 
-//! Input pre-raytracing image rendered by OpenGL.
-uniform sampler2D uOpenGlColorTexture;
-//! Input pre-raytracing depth image rendered by OpenGL.
-uniform sampler2D uOpenGlDepthTexture;
-
 //! Total number of light sources.
 uniform int uLightCount;
 //! Intensity of global ambient light.
@@ -95,6 +95,14 @@ uniform float uSceneEpsilon;
   uniform uvec2 uTextureSamplers[MAX_TEX_NUMBER];
 #endif
 
+#ifdef ADAPTIVE_SAMPLING
+  //! OpenGL image used for accumulating rendering result.
+  volatile restrict layout(size1x32) uniform image2D  uRenderImage;
+
+  //! OpenGL image storing offsets of sampled pixels blocks.
+  coherent restrict layout(size2x32) uniform iimage2D uOffsetImage;
+#endif
+
 //! Top color of gradient background.
 uniform vec4 uBackColorTop = vec4 (0.0);
 //! Bottom color of gradient background.
@@ -240,7 +248,22 @@ vec3 InverseDirection (in vec3 theInput)
 //=======================================================================
 vec4 BackgroundColor()
 {
+#ifdef ADAPTIVE_SAMPLING
+
+  ivec2 aFragCoord = ivec2 (gl_FragCoord.xy);
+
+  ivec2 aTileXY = imageLoad (uOffsetImage, ivec2 (aFragCoord.x / BLOCK_SIZE,
+                                                  aFragCoord.y / BLOCK_SIZE)).xy;
+
+  aTileXY.y += aFragCoord.y % min (uWinSizeY - aTileXY.y, BLOCK_SIZE);
+
+  return mix (uBackColorBot, uBackColorTop, float (aTileXY.y) / uWinSizeY);
+
+#else
+
   return mix (uBackColorBot, uBackColorTop, vPixel.y);
+
+#endif
 }
 
 /////////////////////////////////////////////////////////////////////////////////////////
@@ -263,39 +286,6 @@ SRay GenerateRay (in vec2 thePixel)
   return SRay (mix (aP0, aP1, thePixel.y), aDirection);
 }
 
-// =======================================================================
-// function : ComputeOpenGlDepth
-// purpose  :
-// =======================================================================
-float ComputeOpenGlDepth (in SRay theRay)
-{
-  // a depth in range [0,1]
-  float anOpenGlDepth = texelFetch (uOpenGlDepthTexture, ivec2 (gl_FragCoord.xy), 0).r;
-  // pixel point in NDC-space [-1,1]
-  vec4 aPoint = vec4 (2.0f * vPixel.x - 1.0f,
-                      2.0f * vPixel.y - 1.0f,
-                      2.0f * anOpenGlDepth - 1.0f,
-                      1.0f);
-  vec4 aFinal = uUnviewMat * aPoint;
-  aFinal.xyz *= 1.f / aFinal.w;
-
-  return (anOpenGlDepth < 1.f) ? length (aFinal.xyz - theRay.Origin) : MAXFLOAT;
-}
-
-// =======================================================================
-// function : ComputeOpenGlColor
-// purpose  :
-// =======================================================================
-vec4 ComputeOpenGlColor()
-{
-  vec4 anOpenGlColor = texelFetch (uOpenGlColorTexture, ivec2 (gl_FragCoord.xy), 0);
-  // During blending with factors GL_SRC_ALPHA and GL_ONE_MINUS_SRC_ALPHA (for text and markers)
-  // the alpha channel (written in the color buffer) was squared.
-  anOpenGlColor.a = 1.f - sqrt (anOpenGlColor.a);
-
-  return anOpenGlColor;
-}
-
 // =======================================================================
 // function : IntersectSphere
 // purpose  : Computes ray-sphere intersection
@@ -322,38 +312,44 @@ float IntersectSphere (in SRay theRay, in float theRadius)
 // function : IntersectTriangle
 // purpose  : Computes ray-triangle intersection (branchless version)
 // =======================================================================
-float IntersectTriangle (in SRay theRay,
-                         in vec3 thePnt0,
-                         in vec3 thePnt1,
-                         in vec3 thePnt2,
-                         out vec2 theUV,
-                         out vec3 theNorm)
+void IntersectTriangle (in SRay theRay,
+                        in vec3 thePnt0,
+                        in vec3 thePnt1,
+                        in vec3 thePnt2,
+                        out vec3 theUVT,
+                        out vec3 theNorm)
 {
+  vec3 aToTrg = thePnt0 - theRay.Origin;
+
   vec3 aEdge0 = thePnt1 - thePnt0;
   vec3 aEdge1 = thePnt0 - thePnt2;
 
   theNorm = cross (aEdge1, aEdge0);
 
-  vec3 aEdge2 = (1.0f / dot (theNorm, theRay.Direct)) * (thePnt0 - theRay.Origin);
+  vec3 theVect = cross (theRay.Direct, aToTrg);
 
-  float aTime = dot (theNorm, aEdge2);
+  theUVT = vec3 (dot (theNorm, aToTrg),
+                 dot (theVect, aEdge1),
+                 dot (theVect, aEdge0)) * (1.f / dot (theNorm, theRay.Direct));
 
-  vec3 theVec = cross (theRay.Direct, aEdge2);
+  theUVT.x = any (lessThan (theUVT, ZERO)) || (theUVT.y + theUVT.z) > 1.f ? MAXFLOAT : theUVT.x;
+}
 
-  theUV.x = dot (theVec, aEdge1);
-  theUV.y = dot (theVec, aEdge0);
+#define EMPTY_ROOT ivec4(0)
 
-  return bool (int(aTime >= 0.0f) &
-               int(theUV.x >= 0.0f) &
-               int(theUV.y >= 0.0f) &
-               int(theUV.x + theUV.y <= 1.0f)) ? aTime : MAXFLOAT;
-}
+//! Utility structure containing information about
+//! currently traversing sub-tree of scene's BVH.
+struct SSubTree
+{
+  //! Transformed ray.
+  SRay  TrsfRay;
 
-//! Identifies the absence of intersection.
-#define INALID_HIT ivec4 (-1)
+  //! Inversed ray direction.
+  vec3  Inverse;
 
-//! Global stack shared between traversal functions.
-int Stack[STACK_SIZE];
+  //! Parameters of sub-root node.
+  ivec4 SubData;
+};
 
 #define MATERIAL_AMBN(index) (18 * index + 0)
 #define MATERIAL_DIFF(index) (18 * index + 1)
@@ -366,24 +362,43 @@ int Stack[STACK_SIZE];
 #define MATERIAL_TRS2(index) (18 * index + 8)
 #define MATERIAL_TRS3(index) (18 * index + 9)
 
-struct SSubTree
-{
-  //! Transformed ray.
-  SRay  TrsfRay;
-
-  //! Inversed ray direction.
-  vec3  Inverse;
-
-  //! Parameters of sub-root node.
-  ivec4 SubData;
-};
-
 #define TRS_OFFSET(treelet) treelet.SubData.x
 #define BVH_OFFSET(treelet) treelet.SubData.y
 #define VRT_OFFSET(treelet) treelet.SubData.z
 #define TRG_OFFSET(treelet) treelet.SubData.w
 
-#define EMPTY_ROOT ivec4(0)
+//! Identifies the absence of intersection.
+#define INALID_HIT ivec4 (-1)
+
+//! Global stack shared between traversal functions.
+int Stack[STACK_SIZE];
+
+// =======================================================================
+// function : pop
+// purpose  :
+// =======================================================================
+int pop (inout int theHead)
+{
+  int aData = Stack[theHead];
+
+  int aMask = aData >> 26;
+  int aNode = aMask & 0x3;
+
+  aMask >>= 2;
+
+  if ((aMask & 0x3) == aNode)
+  {
+    --theHead;
+  }
+  else
+  {
+    aMask |= (aMask << 2) & 0x30;
+
+    Stack[theHead] = (aData & 0x03FFFFFF) | (aMask << 26);
+  }
+
+  return (aData & 0x03FFFFFF) + aNode;
+}
 
 // =======================================================================
 // function : SceneNearestHit
@@ -399,55 +414,90 @@ ivec4 SceneNearestHit (in SRay theRay, in vec3 theInverse, inout SIntersect theH
 
   SSubTree aSubTree = SSubTree (theRay, theInverse, EMPTY_ROOT);
 
-  for (bool toContinue = true; toContinue;)
+  for (bool toContinue = true; toContinue; /* none */)
   {
     ivec4 aData = texelFetch (uSceneNodeInfoTexture, aNode);
 
     if (aData.x == 0) // if inner node
     {
-      float aTimeOut;
-      float aTimeLft;
-      float aTimeRgh;
-
       aData.y += BVH_OFFSET (aSubTree);
-      aData.z += BVH_OFFSET (aSubTree);
 
-      vec3 aNodeMinLft = texelFetch (uSceneMinPointTexture, aData.y).xyz;
-      vec3 aNodeMinRgh = texelFetch (uSceneMinPointTexture, aData.z).xyz;
-      vec3 aNodeMaxLft = texelFetch (uSceneMaxPointTexture, aData.y).xyz;
-      vec3 aNodeMaxRgh = texelFetch (uSceneMaxPointTexture, aData.z).xyz;
+      vec4 aHitTimes = vec4 (MAXFLOAT,
+                             MAXFLOAT,
+                             MAXFLOAT,
+                             MAXFLOAT);
+
+      vec3 aRayOriginInverse = -aSubTree.TrsfRay.Origin * aSubTree.Inverse;
+
+      vec3 aNodeMin0 = texelFetch (uSceneMinPointTexture, aData.y +                0).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMin1 = texelFetch (uSceneMinPointTexture, aData.y +                1).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMin2 = texelFetch (uSceneMinPointTexture, aData.y + min (2, aData.z)).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMin3 = texelFetch (uSceneMinPointTexture, aData.y + min (3, aData.z)).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMax0 = texelFetch (uSceneMaxPointTexture, aData.y +                0).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMax1 = texelFetch (uSceneMaxPointTexture, aData.y +                1).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMax2 = texelFetch (uSceneMaxPointTexture, aData.y + min (2, aData.z)).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMax3 = texelFetch (uSceneMaxPointTexture, aData.y + min (3, aData.z)).xyz * aSubTree.Inverse + aRayOriginInverse;
 
-      vec3 aTime0 = (aNodeMinLft - aSubTree.TrsfRay.Origin) * aSubTree.Inverse;
-      vec3 aTime1 = (aNodeMaxLft - aSubTree.TrsfRay.Origin) * aSubTree.Inverse;
+      vec3 aTimeMax = max (aNodeMin0, aNodeMax0);
+      vec3 aTimeMin = min (aNodeMin0, aNodeMax0);
 
-      vec3 aTimeMax = max (aTime0, aTime1);
-      vec3 aTimeMin = min (aTime0, aTime1);
+      float aTimeLeave = min (aTimeMax.x, min (aTimeMax.y, aTimeMax.z));
+      float aTimeEnter = max (aTimeMin.x, max (aTimeMin.y, aTimeMin.z));
 
-      aTime0 = (aNodeMinRgh - aSubTree.TrsfRay.Origin) * aSubTree.Inverse;
-      aTime1 = (aNodeMaxRgh - aSubTree.TrsfRay.Origin) * aSubTree.Inverse;
+      aHitTimes.x = mix (MAXFLOAT, aTimeEnter,
+        aTimeEnter <= aTimeLeave && aTimeEnter <= theHit.Time && aTimeLeave >= 0.f);
 
-      aTimeOut = min (aTimeMax.x, min (aTimeMax.y, aTimeMax.z));
-      aTimeLft = max (aTimeMin.x, max (aTimeMin.y, aTimeMin.z));
+      aTimeMax = max (aNodeMin1, aNodeMax1);
+      aTimeMin = min (aNodeMin1, aNodeMax1);
 
-      int aHitLft = int(aTimeLft <= aTimeOut) & int(aTimeOut >= 0.0f) & int(aTimeLft <= theHit.Time);
+      aTimeLeave = min (aTimeMax.x, min (aTimeMax.y, aTimeMax.z));
+      aTimeEnter = max (aTimeMin.x, max (aTimeMin.y, aTimeMin.z));
 
-      aTimeMax = max (aTime0, aTime1);
-      aTimeMin = min (aTime0, aTime1);
+      aHitTimes.y = mix (MAXFLOAT, aTimeEnter,
+        aTimeEnter <= aTimeLeave && aTimeEnter <= theHit.Time && aTimeLeave >= 0.f);
 
-      aTimeOut = min (aTimeMax.x, min (aTimeMax.y, aTimeMax.z));
-      aTimeRgh = max (aTimeMin.x, max (aTimeMin.y, aTimeMin.z));
+      aTimeMax = max (aNodeMin2, aNodeMax2);
+      aTimeMin = min (aNodeMin2, aNodeMax2);
 
-      int aHitRgh = int(aTimeRgh <= aTimeOut) & int(aTimeOut >= 0.0f) & int(aTimeRgh <= theHit.Time);
+      aTimeLeave = min (aTimeMax.x, min (aTimeMax.y, aTimeMax.z));
+      aTimeEnter = max (aTimeMin.x, max (aTimeMin.y, aTimeMin.z));
 
-      aNode = (aHitLft != 0) ? aData.y : aData.z;
+      aHitTimes.z = mix (MAXFLOAT, aTimeEnter,
+        aTimeEnter <= aTimeLeave && aTimeEnter <= theHit.Time && aTimeLeave >= 0.f && aData.z > 1);
 
-      if (aHitLft + aHitRgh == 2) // hit both children
+      aTimeMax = max (aNodeMin3, aNodeMax3);
+      aTimeMin = min (aNodeMin3, aNodeMax3);
+
+      aTimeLeave = min (aTimeMax.x, min (aTimeMax.y, aTimeMax.z));
+      aTimeEnter = max (aTimeMin.x, max (aTimeMin.y, aTimeMin.z));
+
+      aHitTimes.w = mix (MAXFLOAT, aTimeEnter,
+        aTimeEnter <= aTimeLeave && aTimeEnter <= theHit.Time && aTimeLeave >= 0.f && aData.z > 2);
+
+      ivec4 aChildren = ivec4 (0, 1, 2, 3);
+
+      aChildren.xy = aHitTimes.y < aHitTimes.x ? aChildren.yx : aChildren.xy;
+      aHitTimes.xy = aHitTimes.y < aHitTimes.x ? aHitTimes.yx : aHitTimes.xy;
+      aChildren.zw = aHitTimes.w < aHitTimes.z ? aChildren.wz : aChildren.zw;
+      aHitTimes.zw = aHitTimes.w < aHitTimes.z ? aHitTimes.wz : aHitTimes.zw;
+      aChildren.xz = aHitTimes.z < aHitTimes.x ? aChildren.zx : aChildren.xz;
+      aHitTimes.xz = aHitTimes.z < aHitTimes.x ? aHitTimes.zx : aHitTimes.xz;
+      aChildren.yw = aHitTimes.w < aHitTimes.y ? aChildren.wy : aChildren.yw;
+      aHitTimes.yw = aHitTimes.w < aHitTimes.y ? aHitTimes.wy : aHitTimes.yw;
+      aChildren.yz = aHitTimes.z < aHitTimes.y ? aChildren.zy : aChildren.yz;
+      aHitTimes.yz = aHitTimes.z < aHitTimes.y ? aHitTimes.zy : aHitTimes.yz;
+
+      if (aHitTimes.x != MAXFLOAT)
       {
-        aNode = (aTimeLft < aTimeRgh) ? aData.y : aData.z;
+        int aHitMask = (aHitTimes.w != MAXFLOAT ? aChildren.w : aChildren.z) << 2
+                     | (aHitTimes.z != MAXFLOAT ? aChildren.z : aChildren.y);
+
+        if (aHitTimes.y != MAXFLOAT)
+          Stack[++aHead] = aData.y | (aHitMask << 2 | aChildren.y) << 26;
 
-        Stack[++aHead] = (aTimeLft < aTimeRgh) ? aData.z : aData.y;
+        aNode = aData.y + aChildren.x;
       }
-      else if (aHitLft == aHitRgh) // no hit
+      else
       {
         toContinue = (aHead >= 0);
 
@@ -456,13 +506,14 @@ ivec4 SceneNearestHit (in SRay theRay, in vec3 theInverse, inout SIntersect theH
           aStop = -1; aSubTree = SSubTree (theRay, theInverse, EMPTY_ROOT);
         }
 
-        aNode = Stack[abs (aHead)]; --aHead;
+        if (aHead >= 0)
+          aNode = pop (aHead);
       }
     }
-    else if (aData.x < 0) // leaf node (containg triangles)
+    else if (aData.x < 0) // leaf node (contains triangles)
     {
       vec3 aNormal;
-      vec2 aParams;
+      vec3 aTimeUV;
 
       for (int anIdx = aData.y; anIdx <= aData.z; ++anIdx)
       {
@@ -472,16 +523,15 @@ ivec4 SceneNearestHit (in SRay theRay, in vec3 theInverse, inout SIntersect theH
         vec3 aPoint1 = texelFetch (uGeometryVertexTexture, aTriangle.y += VRT_OFFSET (aSubTree)).xyz;
         vec3 aPoint2 = texelFetch (uGeometryVertexTexture, aTriangle.z += VRT_OFFSET (aSubTree)).xyz;
 
-        float aTime = IntersectTriangle (aSubTree.TrsfRay,
-          aPoint0, aPoint1, aPoint2, aParams, aNormal);
+        IntersectTriangle (aSubTree.TrsfRay, aPoint0, aPoint1, aPoint2, aTimeUV, aNormal);
 
-        if (aTime < theHit.Time)
+        if (aTimeUV.x < theHit.Time)
         {
           aTriIndex = aTriangle;
 
           theTrsfId = TRS_OFFSET (aSubTree);
 
-          theHit = SIntersect (aTime, aParams, aNormal);
+          theHit = SIntersect (aTimeUV.x, aTimeUV.yz, aNormal);
         }
       }
 
@@ -492,7 +542,8 @@ ivec4 SceneNearestHit (in SRay theRay, in vec3 theInverse, inout SIntersect theH
         aStop = -1; aSubTree = SSubTree (theRay, theInverse, EMPTY_ROOT);
       }
 
-      aNode = Stack[abs (aHead)]; --aHead;
+      if (aHead >= 0)
+        aNode = pop (aHead);
     }
     else if (aData.x > 0) // switch node
     {
@@ -540,55 +591,90 @@ float SceneAnyHit (in SRay theRay, in vec3 theInverse, in float theDistance)
 
   SSubTree aSubTree = SSubTree (theRay, theInverse, EMPTY_ROOT);
 
-  for (bool toContinue = true; toContinue;)
+  for (bool toContinue = true; toContinue; /* none */)
   {
     ivec4 aData = texelFetch (uSceneNodeInfoTexture, aNode);
 
     if (aData.x == 0) // if inner node
     {
-      float aTimeOut;
-      float aTimeLft;
-      float aTimeRgh;
-
       aData.y += BVH_OFFSET (aSubTree);
-      aData.z += BVH_OFFSET (aSubTree);
 
-      vec3 aNodeMinLft = texelFetch (uSceneMinPointTexture, aData.y).xyz;
-      vec3 aNodeMinRgh = texelFetch (uSceneMinPointTexture, aData.z).xyz;
-      vec3 aNodeMaxLft = texelFetch (uSceneMaxPointTexture, aData.y).xyz;
-      vec3 aNodeMaxRgh = texelFetch (uSceneMaxPointTexture, aData.z).xyz;
+      vec4 aHitTimes = vec4 (MAXFLOAT,
+                             MAXFLOAT,
+                             MAXFLOAT,
+                             MAXFLOAT);
+
+      vec3 aRayOriginInverse = -aSubTree.TrsfRay.Origin * aSubTree.Inverse;
+
+      vec3 aNodeMin0 = texelFetch (uSceneMinPointTexture, aData.y +                0).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMin1 = texelFetch (uSceneMinPointTexture, aData.y +                1).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMin2 = texelFetch (uSceneMinPointTexture, aData.y + min (2, aData.z)).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMin3 = texelFetch (uSceneMinPointTexture, aData.y + min (3, aData.z)).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMax0 = texelFetch (uSceneMaxPointTexture, aData.y +                0).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMax1 = texelFetch (uSceneMaxPointTexture, aData.y +                1).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMax2 = texelFetch (uSceneMaxPointTexture, aData.y + min (2, aData.z)).xyz * aSubTree.Inverse + aRayOriginInverse;
+      vec3 aNodeMax3 = texelFetch (uSceneMaxPointTexture, aData.y + min (3, aData.z)).xyz * aSubTree.Inverse + aRayOriginInverse;
+
+      vec3 aTimeMax = max (aNodeMin0, aNodeMax0);
+      vec3 aTimeMin = min (aNodeMin0, aNodeMax0);
+
+      float aTimeLeave = min (aTimeMax.x, min (aTimeMax.y, aTimeMax.z));
+      float aTimeEnter = max (aTimeMin.x, max (aTimeMin.y, aTimeMin.z));
+
+      aHitTimes.x = mix (MAXFLOAT, aTimeEnter,
+        aTimeEnter <= aTimeLeave && aTimeEnter <= theDistance && aTimeLeave >= 0.f);
+
+      aTimeMax = max (aNodeMin1, aNodeMax1);
+      aTimeMin = min (aNodeMin1, aNodeMax1);
 
-      vec3 aTime0 = (aNodeMinLft - aSubTree.TrsfRay.Origin) * aSubTree.Inverse;
-      vec3 aTime1 = (aNodeMaxLft - aSubTree.TrsfRay.Origin) * aSubTree.Inverse;
+      aTimeLeave = min (aTimeMax.x, min (aTimeMax.y, aTimeMax.z));
+      aTimeEnter = max (aTimeMin.x, max (aTimeMin.y, aTimeMin.z));
 
-      vec3 aTimeMax = max (aTime0, aTime1);
-      vec3 aTimeMin = min (aTime0, aTime1);
+      aHitTimes.y = mix (MAXFLOAT, aTimeEnter,
+        aTimeEnter <= aTimeLeave && aTimeEnter <= theDistance && aTimeLeave >= 0.f);
 
-      aTime0 = (aNodeMinRgh - aSubTree.TrsfRay.Origin) * aSubTree.Inverse;
-      aTime1 = (aNodeMaxRgh - aSubTree.TrsfRay.Origin) * aSubTree.Inverse;
+      aTimeMax = max (aNodeMin2, aNodeMax2);
+      aTimeMin = min (aNodeMin2, aNodeMax2);
 
-      aTimeOut = min (aTimeMax.x, min (aTimeMax.y, aTimeMax.z));
-      aTimeLft = max (aTimeMin.x, max (aTimeMin.y, aTimeMin.z));
+      aTimeLeave = min (aTimeMax.x, min (aTimeMax.y, aTimeMax.z));
+      aTimeEnter = max (aTimeMin.x, max (aTimeMin.y, aTimeMin.z));
 
-      int aHitLft = int(aTimeLft <= aTimeOut) & int(aTimeOut >= 0.0f) & int(aTimeLft <= theDistance);
+      aHitTimes.z = mix (MAXFLOAT, aTimeEnter,
+        aTimeEnter <= aTimeLeave && aTimeEnter <= theDistance && aTimeLeave >= 0.f && aData.z > 1);
 
-      aTimeMax = max (aTime0, aTime1);
-      aTimeMin = min (aTime0, aTime1);
+      aTimeMax = max (aNodeMin3, aNodeMax3);
+      aTimeMin = min (aNodeMin3, aNodeMax3);
 
-      aTimeOut = min (aTimeMax.x, min (aTimeMax.y, aTimeMax.z));
-      aTimeRgh = max (aTimeMin.x, max (aTimeMin.y, aTimeMin.z));
+      aTimeLeave = min (aTimeMax.x, min (aTimeMax.y, aTimeMax.z));
+      aTimeEnter = max (aTimeMin.x, max (aTimeMin.y, aTimeMin.z));
 
-      int aHitRgh = int(aTimeRgh <= aTimeOut) & int(aTimeOut >= 0.0f) & int(aTimeRgh <= theDistance);
+      aHitTimes.w = mix (MAXFLOAT, aTimeEnter,
+        aTimeEnter <= aTimeLeave && aTimeEnter <= theDistance && aTimeLeave >= 0.f && aData.z > 2);
 
-      aNode = (aHitLft != 0) ? aData.y : aData.z;
+      ivec4 aChildren = ivec4 (0, 1, 2, 3);
 
-      if (aHitLft + aHitRgh == 2) // hit both children
+      aChildren.xy = aHitTimes.y < aHitTimes.x ? aChildren.yx : aChildren.xy;
+      aHitTimes.xy = aHitTimes.y < aHitTimes.x ? aHitTimes.yx : aHitTimes.xy;
+      aChildren.zw = aHitTimes.w < aHitTimes.z ? aChildren.wz : aChildren.zw;
+      aHitTimes.zw = aHitTimes.w < aHitTimes.z ? aHitTimes.wz : aHitTimes.zw;
+      aChildren.xz = aHitTimes.z < aHitTimes.x ? aChildren.zx : aChildren.xz;
+      aHitTimes.xz = aHitTimes.z < aHitTimes.x ? aHitTimes.zx : aHitTimes.xz;
+      aChildren.yw = aHitTimes.w < aHitTimes.y ? aChildren.wy : aChildren.yw;
+      aHitTimes.yw = aHitTimes.w < aHitTimes.y ? aHitTimes.wy : aHitTimes.yw;
+      aChildren.yz = aHitTimes.z < aHitTimes.y ? aChildren.zy : aChildren.yz;
+      aHitTimes.yz = aHitTimes.z < aHitTimes.y ? aHitTimes.zy : aHitTimes.yz;
+
+      if (aHitTimes.x != MAXFLOAT)
       {
-        aNode = (aTimeLft < aTimeRgh) ? aData.y : aData.z;
+        int aHitMask = (aHitTimes.w != MAXFLOAT ? aChildren.w : aChildren.z) << 2
+                     | (aHitTimes.z != MAXFLOAT ? aChildren.z : aChildren.y);
+
+        if (aHitTimes.y != MAXFLOAT)
+          Stack[++aHead] = aData.y | (aHitMask << 2 | aChildren.y) << 26;
 
-        Stack[++aHead] = (aTimeLft < aTimeRgh) ? aData.z : aData.y;
+        aNode = aData.y + aChildren.x;
       }
-      else if (aHitLft == aHitRgh) // no hit
+      else
       {
         toContinue = (aHead >= 0);
 
@@ -597,13 +683,14 @@ float SceneAnyHit (in SRay theRay, in vec3 theInverse, in float theDistance)
           aStop = -1; aSubTree = SSubTree (theRay, theInverse, EMPTY_ROOT);
         }
 
-        aNode = Stack[abs (aHead)]; --aHead;
+        if (aHead >= 0)
+          aNode = pop (aHead);
       }
     }
     else if (aData.x < 0) // leaf node
     {
       vec3 aNormal;
-      vec2 aParams;
+      vec3 aTimeUV;
 
       for (int anIdx = aData.y; anIdx <= aData.z; ++anIdx)
       {
@@ -613,16 +700,15 @@ float SceneAnyHit (in SRay theRay, in vec3 theInverse, in float theDistance)
         vec3 aPoint1 = texelFetch (uGeometryVertexTexture, aTriangle.y += VRT_OFFSET (aSubTree)).xyz;
         vec3 aPoint2 = texelFetch (uGeometryVertexTexture, aTriangle.z += VRT_OFFSET (aSubTree)).xyz;
 
-        float aTime = IntersectTriangle (aSubTree.TrsfRay,
-          aPoint0, aPoint1, aPoint2, aParams, aNormal);
+        IntersectTriangle (aSubTree.TrsfRay, aPoint0, aPoint1, aPoint2, aTimeUV, aNormal);
 
 #ifdef TRANSPARENT_SHADOWS
-        if (aTime < theDistance)
+        if (aTimeUV.x < theDistance)
         {
           aFactor *= 1.f - texelFetch (uRaytraceMaterialTexture, MATERIAL_TRAN (aTriangle.w)).x;
         }
 #else
-        if (aTime < theDistance)
+        if (aTimeUV.x < theDistance)
         {
           aFactor = 0.f;
         }
@@ -636,7 +722,8 @@ float SceneAnyHit (in SRay theRay, in vec3 theInverse, in float theDistance)
         aStop = -1; aSubTree = SSubTree (theRay, theInverse, EMPTY_ROOT);
       }
 
-      aNode = Stack[abs (aHead)]; --aHead;
+      if (aHead >= 0)
+        aNode = pop (aHead);
     }
     else if (aData.x > 0) // switch node
     {
@@ -779,7 +866,6 @@ vec4 Radiance (in SRay theRay, in vec3 theInverse)
 
   int aTrsfId;
 
-  float anOpenGlDepth = ComputeOpenGlDepth (theRay);
   float aRaytraceDepth = MAXFLOAT;
 
   for (int aDepth = 0; aDepth < NB_BOUNCES; ++aDepth)
@@ -801,8 +887,7 @@ vec4 Radiance (in SRay theRay, in vec3 theInverse)
       }
       else
       {
-        vec4 aGlColor = ComputeOpenGlColor();
-        aColor = vec4 (mix (aGlColor.rgb, BackgroundColor().rgb, aGlColor.w), aGlColor.w);
+        aColor = BackgroundColor();
       }
 
       aResult += aWeight.xyz * aColor.xyz; aWeight.w *= aColor.w;
@@ -818,31 +903,22 @@ vec4 Radiance (in SRay theRay, in vec3 theInverse)
                                    dot (aInvTransf1, aHit.Normal),
                                    dot (aInvTransf2, aHit.Normal)));
 
-    // For polygons that are parallel to the screen plane, the depth slope
-    // is equal to 1, resulting in small polygon offset. For polygons that
-    // that are at a large angle to the screen, the depth slope tends to 1,
-    // resulting in a larger polygon offset
-    float aPolygonOffset = uSceneEpsilon * EPS_SCALE /
-      max (abs (dot (theRay.Direct, aHit.Normal)), MIN_SLOPE);
-
-    if (anOpenGlDepth < aHit.Time + aPolygonOffset)
-    {
-      vec4 aGlColor = ComputeOpenGlColor();
-
-      aResult += aWeight.xyz * aGlColor.xyz;
-      aWeight *= aGlColor.w;
-    }
-
     theRay.Origin += theRay.Direct * aHit.Time; // intersection point
 
-    // Evaluate depth
+    // Evaluate depth on first hit
     if (aDepth == 0)
     {
+      // For polygons that are parallel to the screen plane, the depth slope
+      // is equal to 1, resulting in small polygon offset. For polygons that
+      // that are at a large angle to the screen, the depth slope tends to 1,
+      // resulting in a larger polygon offset
+      float aPolygonOffset = uSceneEpsilon * EPS_SCALE /
+        max (abs (dot (theRay.Direct, aHit.Normal)), MIN_SLOPE);
+
       // Hit point in NDC-space [-1,1] (the polygon offset is applied in the world space)
       vec4 aNDCPoint = uViewMat * vec4 (theRay.Origin + theRay.Direct * aPolygonOffset, 1.f);
-      aNDCPoint.xyz *= 1.f / aNDCPoint.w;
 
-      aRaytraceDepth = aNDCPoint.z * 0.5f + 0.5f;
+      aRaytraceDepth = (aNDCPoint.z / aNDCPoint.w) * 0.5f + 0.5f;
     }
 
     vec3 aNormal = SmoothNormal (aHit.UV, aTriIndex);
@@ -945,10 +1021,6 @@ vec4 Radiance (in SRay theRay, in vec3 theInverse)
       {
         theRay.Direct = Refract (theRay.Direct, aNormal, aOpacity.z, aOpacity.w);
       }
-      else
-      {
-        anOpenGlDepth -= aHit.Time + uSceneEpsilon;
-      }
     }
     else
     {
@@ -977,8 +1049,6 @@ vec4 Radiance (in SRay theRay, in vec3 theInverse)
       theInverse = 1.0f / max (abs (theRay.Direct), SMALL);
 
       theInverse = mix (-theInverse, theInverse, step (ZERO, theRay.Direct));
-
-      anOpenGlDepth = MAXFLOAT; // disable combining image with OpenGL output
     }
 
     theRay.Origin += theRay.Direct * uSceneEpsilon;