[译] GLSL 中的视差遮蔽映射（Parallax Occlusion Mapping in GLSL）

关于渲染的中文文章可谓是少之又少，而不少书和中英文技术文章老是大篇幅的进行晦涩难懂的公式推导，这种方式确实表达准确，可苦了数学很差的娃，能找到一篇好的材料进行学习真的是一件很不容易的事情。

我在学习Parallax Mapping的时候无心间找到这篇文章，图文并茂，而且把Bump Mapping系列中的各类技术由简单到复杂逐一介绍给了读者。开心之下就决定把它翻译成中文，以飨后人。

虽然说简单，可是它也不是从零开始的。本文要求读者具有了理解最基本的Normal Mapping的知识，理解切空间。Normal Mapping的文章网上大把抓，切空间的话能够读这篇文章。

下面是正文，原文连接在此。

---

这一课讲如何用GLSL和OpenGL实现各类视差映射技术（一样的技术亦可用DirectX实现）。主要会涵盖以下几种技术：视差映射（Parallax Mapping），带偏移上限的视差映射（Parallax Mapping with Offset Limiting），陡峭视差映射（Steep Parallax Mapping），浮雕视差映射（Relief Parallax Mapping）和视差遮蔽映射(Parallax Occlusion Mapping）。另外本文还会介绍如何实如今视差映射中的自阴影（软阴影）。下面的几个图片展现了几种视差映射技术和简单光照或者法线映射的效果对比。

我就是给原文截了个图，建议进到原文去看图片。

视差映射基础

在计算机图形学中视差映射是法线映射的一个加强版本，它不止改变了光照的做用方式，还在平坦的多边形上建立了3D细节的假象。不会生成任何额外的图元。上面的图片展现了视差映射和发现映射的对比。你可能以为视差映射偏移了原始图元，但其实它只是偏移了用来获取颜色和法线的纹理坐标。

要实现视差映射你须要一张高度贴图。高度图中的每一个像素包含了表面高度的信息。纹理中的高度会被转化成对应的点沉入表面多少的信息。这种状况你得把高度图中读出来的值反过来用。这篇教程中的视差映射会把高度图中的值当深度来用，黑色(0)表明和表面齐平的高度，白色(1)表明最深的凹陷值。

下面的例子会用到3张纹理：高度图，漫反射颜色纹理和法线贴图。一般法线贴图都是从高度图生成出来的。在咱们的栗子中，高度图被当成深度图看待，因此在生成法线贴图以前你得先反转高度图（译者：服了，能不说车轱辘话吗？）。你还能把高度图和法线贴图合并到一张纹理中，把高度存在Alpha通道里，可是为了表述清楚本文仍是把他们分开用了。下面是这三张图：

视差映射技术的主要任务是修改纹理坐标，让平面看起来像是立体的。主要计算都是在Fragment Shader中进行。看看下面的图片。水平线0.0表示彻底没有凹陷的深度，水平线1.0表示凹陷的最大深度。实际的几何体并没改变，其实一直都在0.0水平线上。图中的曲线表明了高度图中存储的高度数据。

设当前点（译者：原文中用的是Fragment，片元。）是图片中用黄色方块高亮出来的那个点，这个点的纹理坐标是T0。向量V是从摄像机到点的方向向量。用坐标T0在高度图上采样，你能获得这个点的高度值H(T0)=0.55。这个值不是0，因此点并非在表面上，而是凹陷下去的。因此你得把向量V继续延长直到与高度图定义出来的表面最近的一个交点。这个交点咱们说它的深度就是H(T1)，它的纹理坐标就是T1。因此咱们就应该用T1的纹理坐标去对颜色和法线贴图进行采样。

因此说，全部视差映射技术的主要目的，就是要精确的计算摄像机的向量V和高度图定义出来的表面的交点。

基本的视差映射的Shader

视差映射的计算是在切空间进行的（跟法线映射同样）。因此指向光源的向量(L)和指向摄像机的向量(V)应该先被变换到切空间。在用视差映射计算出来新的纹理坐标以后，你能够用这个坐标来计算自阴影，能够从漫反射贴图读取颜色以及从发现贴图读取法线。

在这个教程中视差映射的实现是在一个叫parallaxMapping()的函数体中，自阴影是在parallaxSoftShadowMultiplier()中，而后Blinn-Phone光照模型和法线映射的代码是在normalMappingLighting()函数体中。下面的顶点和片元着色器是视差映射和自阴影的基础模板。顶点着色器把光照向量和摄像机向量变换到切空间。片元着色器调用视差映射的相关函数，而后计算自阴影系数，并计算最终光照后的颜色值。

// Basic vertex shader for parallax mapping
#version 330

// attributes
layout(location = 0) in vec3    i_position;    // xyz - position
layout(location = 1) in vec3    i_normal;    // xyz - normal
layout(location = 2) in vec2    i_texcoord0;    // xy - texture coords
layout(location = 3) in vec4 i_tangent; // xyz - tangent, w - handedness

// uniforms
uniform mat4 u_model_mat;
uniform mat4 u_view_mat;
uniform mat4 u_proj_mat;
uniform mat3 u_normal_mat;
uniform vec3 u_light_position;
uniform vec3 u_camera_position;

// data for fragment shader
out vec2    o_texcoords;
out vec3 o_toLightInTangentSpace;
out vec3 o_toCameraInTangentSpace;

///////////////////////////////////////////////////////////////////

void main(void)
{
   // transform to world space
   vec4 worldPosition    = u_model_mat * vec4(i_position, 1);
   vec3 worldNormal    = normalize(u_normal_mat * i_normal);
   vec3 worldTangent    = normalize(u_normal_mat * i_tangent.xyz);

   // calculate vectors to the camera and to the light
   vec3 worldDirectionToLight    = normalize(u_light_position - worldPosition.xyz);
   vec3 worldDirectionToCamera    = normalize(u_camera_position - worldPosition.xyz);

   // calculate bitangent from normal and tangent
   vec3 worldBitangnent    = cross(worldNormal, worldTangent) * i_tangent.w;

   // transform direction to the light to tangent space
   o_toLightInTangentSpace = vec3(
         dot(worldDirectionToLight, worldTangent),
         dot(worldDirectionToLight, worldBitangnent),
         dot(worldDirectionToLight, worldNormal)
      );

   // transform direction to the camera to tangent space
   o_toCameraInTangentSpace= vec3(
         dot(worldDirectionToCamera, worldTangent),
         dot(worldDirectionToCamera, worldBitangnent),
         dot(worldDirectionToCamera, worldNormal)
      );

   // pass texture coordinates to fragment shader
   o_texcoords    = i_texcoord0;

   // calculate screen space position of the vertex
   gl_Position    = u_proj_mat * u_view_mat * worldPosition;
}

// basic fragment shader for Parallax Mapping
#version 330

// data from vertex shader
in vec2    o_texcoords;
in vec3    o_toLightInTangentSpace;
in vec3    o_toCameraInTangentSpace;

// textures
layout(location = 0) uniform sampler2D u_diffuseTexture;
layout(location = 1) uniform sampler2D u_heightTexture;
layout(location = 2) uniform sampler2D u_normalTexture;

// color output to the framebuffer
out vec4    resultingColor;

////////////////////////////////////////

// scale for size of Parallax Mapping effect
uniform float    parallaxScale; // ~0.1

//////////////////////////////////////////////////////
// Implements Parallax Mapping technique
// Returns modified texture coordinates, and last used depth
vec2 parallaxMapping(in vec3 V, in vec2 T, out float parallaxHeight)
{
   // ...
}

//////////////////////////////////////////////////////
// Implements self-shadowing technique - hard or soft shadows
// Returns shadow factor
float parallaxSoftShadowMultiplier(in vec3 L, in vec2 initialTexCoord,
                                       in float initialHeight)
{
   // ...
}

//////////////////////////////////////////////////////
// Calculates lighting by Blinn-Phong model and Normal Mapping
// Returns color of the fragment
vec4 normalMappingLighting(in vec2 T, in vec3 L, in vec3 V, float shadowMultiplier)
{
   // restore normal from normal map
   vec3 N = normalize(texture(u_normalTexture, T).xyz * 2 - 1);
   vec3 D = texture(u_diffuseTexture, T).rgb;

   // ambient lighting
   float iamb = 0.2;
   // diffuse lighting
   float idiff = clamp(dot(N, L), 0, 1);
   // specular lighting
   float ispec = 0;
   if(dot(N, L) > 0.2)
   {
      vec3 R = reflect(-L, N);
      ispec = pow(dot(R, V), 32) / 1.5;
   }

   vec4 resColor;
   resColor.rgb = D * (ambientLighting + (idiff + ispec) * pow(shadowMultiplier, 4));
   resColor.a = 1;

   return resColor;
}

/////////////////////////////////////////////
// Entry point for Parallax Mapping shader
void main(void)
{
   // normalize vectors after vertex shader
   vec3 V = normalize(o_toCameraInTangentSpace);
   vec3 L = normalize(o_toLightInTangentSpace);

   // get new texture coordinates from Parallax Mapping
   float parallaxHeight;
   vec2 T = parallaxMapping(V, o_texcoords, parallaxHeight);

   // get self-shadowing factor for elements of parallax
   float shadowMultiplier = parallaxSoftShadowMultiplier(L, T, parallaxHeight - 0.05);

   // calculate lighting
   resultingColor = normalMappingLighting(T, L, V, shadowMultiplier);
}

视差映射和带偏移上限的视差映射

视差映射中最简单的版本只取一步近似来计算新的纹理坐标，这项技术被简单的称为视差映射。视差映射只有在高度图相对比较平滑，而且不存在复杂的细节时，才能获得相对能够接受的效果。若是摄像机向量和法线向量的夹角过大的话，视差映射的效果会是错误的。视差映射近似计算的核心思想是：

• 从高度图读取纹理坐标T0位置的高度H(T0)

• 根据H(T0)和摄像机向量V，在初始的纹理坐标基础上进行偏移。

偏移纹理坐标的方法以下。由于摄像机向量是在切空间下，而切空间是沿着纹理坐标方向创建的，因此向量V的X和Y份量就能够直接不加换算的用做纹理坐标的偏移量。向量V的Z份量是法向份量，垂直于表面。你能够用Z除X和Y。这就是视差映射技术中对纹理坐标的原始计算。你也能够保留X和Y的值，这样的实现叫带偏移上限的视差映射。带偏移上限的视差映射能够避免在摄像机向量V和法向量N夹角太大时的一些诡异的结果。而后你把V的X和Y份量加到原始纹理坐标上，就获得了沿着V方向的新的纹理坐标。

你得把原纹理位置的深度值H(T0)也算进偏移中，直接把V.xy和H(T0)相乘就行了。

你能够用一个scale变量来控制视差映射效果的幅度。一样，你得把它乘给V.xy。最有意义的scale值在0~0.5之间。更高的值每每会致使视差映射计算出错误的效果（见上图）。你也能够把scale设为负数，这样的话你得把法向量的Z份量反转过来。

下面是偏移后的纹理坐标Tp的最终公式：

下图展现了高度图中的深度值H(T0)是如何影响纹理坐标T0沿着V方向偏移的。此情形下做为结果的Tp是错误的，由于视差映射只是一个近似，而并非找出V和表面的准确交点。

这个方法的主要优势是只须要额外对高度图采样一次，因此性能上杠杠的。下面是shader函数的实现：

vec2 parallaxMapping(in vec3 V, in vec2 T, out float parallaxHeight)
{
   // get depth for this fragment
   float initialHeight = texture(u_heightTexture, o_texcoords).r;

   // calculate amount of offset for Parallax Mapping
   vec2 texCoordOffset = parallaxScale * V.xy / V.z * initialHeight;

   // calculate amount of offset for Parallax Mapping With Offset Limiting
   texCoordOffset = parallaxScale * V.xy * initialHeight;

   // retunr modified texture coordinates
   return o_texcoords - texCoordOffset;
}

陡峭视差映射

陡峭视差映射，不像简单的视差映射近似，并不仅是简单粗暴的对纹理坐标进行偏移而不检查合理性和关联性，会检查结果是否接近于正确值。这种方法的核心思想是把表面的深度切分红等距的若干层。而后从最顶端的一层开始采样高度图，每一次会沿着V的方向偏移纹理坐标。若是点已经低于了表面（当前的层的深度大于采样出的深度），中止检查而且使用最后一次采样的纹理坐标做为结果。

陡峭视差映射的工做方式在下面的图片上举例。深度被分割成8个层，每层的高度值是0.125。每层的纹理坐标偏移是V.xy/V.z * scale/numLayers。从顶层黄色方块的位置开始检查，下面是手动计算步骤：

• 层的深度为0，高度图深度H(T0)大约为0.75。采样到的深度大于层的深度，因此开始下一次迭代。

• 沿着V方向偏移纹理坐标，选定下一层。层深度为0.125，高度图深度H(T1)大约为0.625。采样到的深度大于层的深度，因此开始下一次迭代。

• 沿着V方向偏移纹理坐标，选定下一层。层深度为0.25，高度图深度H(T2)大约为0.4。采样到的深度大于层的深度，因此开始下一次迭代。

• 沿着V方向偏移纹理坐标，选定下一层。层深度为0.375，高度图深度H(T3)大约为0.2。采样到的深度小于层的深度，因此向量V上的当前点在表面之下。咱们找到了纹理坐标Tp=T3是实际交点的近似点。

从上图你能看到，其实纹理坐标T3仍是离交点挺远的。可是这个纹理坐标已经比视差映射要接近正确结果了。若是你想获得更精确的结果，增长层的数量。

陡峭视差映射的主要优点在于它把深度切分红了有限数量的层。若是层数不少，那性能就会低。但若是层数少，就会有明显的锯齿现象产生，就像下面这张图同样。你也能够根据摄像机向量V和多边形法向N之间的夹角来动态的决定层的数量。性能和锯齿的问题在下文的浮雕视差映射和视差遮蔽映射中能够解决。

下面是陡峭视差映射的代码：

vec2 parallaxMapping(in vec3 V, in vec2 T, out float parallaxHeight)
{
   // determine number of layers from angle between V and N
   const float minLayers = 5;
   const float maxLayers = 15;
   float numLayers = mix(maxLayers, minLayers, abs(dot(vec3(0, 0, 1), V)));

   // height of each layer
   float layerHeight = 1.0 / numLayers;
   // depth of current layer
   float currentLayerHeight = 0;
   // shift of texture coordinates for each iteration
   vec2 dtex = parallaxScale * V.xy / V.z / numLayers;

   // current texture coordinates
   vec2 currentTextureCoords = T;

   // get first depth from heightmap
   float heightFromTexture = texture(u_heightTexture, currentTextureCoords).r;

   // while point is above surface
   while(heightFromTexture > currentLayerHeight) 
   {
      // to the next layer
      currentLayerHeight += layerHeight;
      // shift texture coordinates along vector V
      currentTextureCoords -= dtex;
      // get new depth from heightmap
      heightFromTexture = texture(u_heightTexture, currentTextureCoords).r;
   }

   // return results
   parallaxHeight = currentLayerHeight;
   return currentTextureCoords;
}

浮雕视差映射

浮雕视差映射升级了陡峭视差映射，让咱们的shader能找到更精确的纹理坐标。首先你先用浮雕视差映射，而后你能获得交点先后的两个层，和对应的深度值。在下面的图中这两个层分别对应纹理坐标T2和T3。如今你能够用二分法来进一步改进你的结果，每一次搜索迭代可使精确度提高一倍。

下图表达了浮雕视差映射的主要步骤：
（译者：这一段的内容和原文区别较大，由于直接按照原文翻译有不少容易混淆的名词，因此我加入了变量声明。）

• 在陡峭视差映射以后，咱们知道交点确定在T2和T3之间。真实的交点在图上用绿点标出来了。

• 设每次迭代时的纹理坐标变化量ST，它的初始值等于向量V在穿过一个层的深度时的XY份量。

• 设每次迭代时的深度值变化量SH，它的初始值等于一个层的深度。

• 把ST和SH都除以2。

• 把纹理坐标T3沿着反方向偏移ST，把层深度沿反方向偏移SH，获得这次迭代的纹理坐标T4和层深度H(T4)。

• （*）采样高度图，把ST和SH都除以2。

• 若是高度图中的深度值大于当前迭代层的深度H(T4)，则将当前迭代层的深度增长SH，迭代的纹理坐标沿着V的方向增长ST。

• 若是高度图中的深度值小于当前迭代层的深度H(T4)，则将当前迭代层的深度减小SH，迭代的纹理坐标沿着V的相反方向增长ST。

• 从(*)处循环，继续二分搜索，直到规定的次数。

• 最后一步获得的纹理坐标就是浮雕视差映射的结果。

下面是浮雕视差映射的实现：

vec2 parallaxMapping(in vec3 V, in vec2 T, out float parallaxHeight)
{
   // determine required number of layers
   const float minLayers = 10;
   const float maxLayers = 15;
   float numLayers = mix(maxLayers, minLayers, abs(dot(vec3(0, 0, 1), V)));

   // height of each layer
   float layerHeight = 1.0 / numLayers;
   // depth of current layer
   float currentLayerHeight = 0;
   // shift of texture coordinates for each iteration
   vec2 dtex = parallaxScale * V.xy / V.z / numLayers;

   // current texture coordinates
   vec2 currentTextureCoords = T;

   // depth from heightmap
   float heightFromTexture = texture(u_heightTexture, currentTextureCoords).r;

   // while point is above surface
   while(heightFromTexture > currentLayerHeight) 
   {
      // go to the next layer
      currentLayerHeight += layerHeight; 
      // shift texture coordinates along V
      currentTextureCoords -= dtex;
      // new depth from heightmap
      heightFromTexture = texture(u_heightTexture, currentTextureCoords).r;
   }

   ///////////////////////////////////////////////////////////
   // Start of Relief Parallax Mapping

   // decrease shift and height of layer by half
   vec2 deltaTexCoord = dtex / 2;
   float deltaHeight = layerHeight / 2;

   // return to the mid point of previous layer
   currentTextureCoords += deltaTexCoord;
   currentLayerHeight -= deltaHeight;

   // binary search to increase precision of Steep Paralax Mapping
   const int numSearches = 5;
   for(int i=0; i<numSearches; i++)
   {
      // decrease shift and height of layer by half
      deltaTexCoord /= 2;
      deltaHeight /= 2;

      // new depth from heightmap
      heightFromTexture = texture(u_heightTexture, currentTextureCoords).r;

      // shift along or agains vector V
      if(heightFromTexture > currentLayerHeight) // below the surface
      {
         currentTextureCoords -= deltaTexCoord;
         currentLayerHeight += deltaHeight;
      }
      else // above the surface
      {
         currentTextureCoords += deltaTexCoord;
         currentLayerHeight -= deltaHeight;
      }
   }

   // return results
   parallaxHeight = currentLayerHeight;    return currentTextureCoords;
}

视差遮蔽映射

视差遮蔽映射(POM)是陡峭视差映射的另外一个改进版本。浮雕视差映射用了二分搜索法来提高结果精度，可是搜索下降程序性能。视差遮蔽映射旨在比浮雕视差映射更好的性能下获得比陡峭视差映射更好的效果。可是POM的效果要比浮雕视差映射差一些。

视差遮蔽映射简单的对陡峭视差映射的结果进行插值。请看下图，POM使用相交以后的层深度（0.375，陡峭视差映射中止迭代的层），上一个采样深度H(T2)和下一个采样深度H(T3)。从图片中你能看到，视差遮蔽映射的插值结果是在视向量V和H(T2)和H(T3)高度的连线的交点上。这个交点已经足够接近实际交点（标记为绿色的点）了。

图片对应的手动计算步骤：

• nextHeight = H(T3) - currentLayerHeight

• prevHeight = H(T2) - (currentLayerHeight - layerHeight)

• weight = nextHeight / (nextHeight - prevHeight)

• Tp = T(T2) weight + T(T3) (1.0 - weight)

视差遮蔽映射可使用相对较少的采样次数产生很好的结果。但视差遮蔽映射比浮雕视差映射更容易跳太高度图中的小细节，也更容易在高度图数据产生大幅度的变化时获得错误的结果。

这是POM的实现：

vec2 parallaxMapping(in vec3 V, in vec2 T, out float parallaxHeight)
{
   // determine optimal number of layers
   const float minLayers = 10;
   const float maxLayers = 15;
   float numLayers = mix(maxLayers, minLayers, abs(dot(vec3(0, 0, 1), V)));

   // height of each layer
   float layerHeight = 1.0 / numLayers;
   // current depth of the layer
   float curLayerHeight = 0;
   // shift of texture coordinates for each layer
   vec2 dtex = parallaxScale * V.xy / V.z / numLayers;

   // current texture coordinates
   vec2 currentTextureCoords = T;

   // depth from heightmap
   float heightFromTexture = texture(u_heightTexture, currentTextureCoords).r;

   // while point is above the surface
   while(heightFromTexture > curLayerHeight) 
   {
      // to the next layer
      curLayerHeight += layerHeight; 
      // shift of texture coordinates
      currentTextureCoords -= dtex;
      // new depth from heightmap
      heightFromTexture = texture(u_heightTexture, currentTextureCoords).r;
   }

   ///////////////////////////////////////////////////////////

   // previous texture coordinates
   vec2 prevTCoords = currentTextureCoords + texStep;

   // heights for linear interpolation
   float nextH    = heightFromTexture - curLayerHeight;
   float prevH    = texture(u_heightTexture, prevTCoords).r
                           - curLayerHeight + layerHeight;

   // proportions for linear interpolation
   float weight = nextH / (nextH - prevH);

   // interpolation of texture coordinates
   vec2 finalTexCoords = prevTCoords * weight + currentTextureCoords * (1.0-weight);

   // interpolation of depth values
   parallaxHeight = curLayerHeight + prevH * weight + nextH * (1.0 - weight);

   // return result
   return finalTexCoords;
}

视差映射和自阴影

你能够经过和陡峭视差映射很接近的算法来肯定一个点是否处于阴影之中。你要向外搜索，而不是向里。同时纹理坐标的偏移应该从点沿着光的方向，而不是沿着摄像机方向。光源向量L应该在切空间中，跟V同样，它能够直接被用做偏移纹理坐标。自阴影计算的结果是一个阴影系数 - 在[0,1]之间的值。这个数值能够在后面用来调节漫反射和镜面反射的光照强度。

要计算硬边缘的阴影（硬阴影）你要沿着L找到第一个在表面之下的点。若是点在表面之下则阴影系数是0， 不然就是1。好比，在下面的图片上，高度值H(TL1)小于层的高度值Ha，因此这个点在表面如下，阴影系数是0。若是光向量直到水平面0.0也没有找到任何点在表面如下，那咱们的片元就应该是在光照中，阴影系数则为1。阴影的质量极大程度上受到分层数量、scale参数和光向量L和多边形的法向量N之间的角度的影响。若是设置不恰当，阴影会出现锯齿或者更糟。

软阴影会计算沿着光源向量L的多个值，只有在表面如下的点才会包含进来。半阴影的系数根据当前层深度和当前点高度图深度之间的差别来得出。你还得把点到片元的举例计算在内。因此半阴影系数要被乘以(1.0 - stepIndex/numberOfSteps)。要计算最终的阴影系数，你得选出那个最大的半阴影系数。由此获得计算软阴影系数的公式：

这里是软阴影系数的计算步骤（对应于图片）

• 设置shadow factor为0，迭代步数为4。

• 沿着L向前步进到Ha。Ha小于H(TL1)，因此该点在表面之下。计算半阴影系数为Ha-H(TL1)。这是第一次检查，总共的检查次数为4，计算距离影响，将半阴影系数乘以(1.0 - 1.0/4.0)。保存这个半阴影系数。

• 沿着L向前步进到Hb。Hb小于H(TL2)，因此该点在表面之下。计算半阴影系数为Hb-H(TL2)。这事第二次检查，总共的检查次数为4，计算距离影响，将半阴影系数乘以(1.0 - 2.0/4.0)。保存这个半阴影系数。

• 沿着L向前步进，这个点在表面之上。

• 最后一次沿着L向前步进，这个点也在表面之上。

• 迭代的点已经高于了水平线0.0，结束迭代。

• 选取最大的半阴影系数做为最终的阴影系数值。

下面是代码：

float parallaxSoftShadowMultiplier(in vec3 L, in vec2 initialTexCoord,
                                       in float initialHeight)
{
   float shadowMultiplier = 1;

   const float minLayers = 15;
   const float maxLayers = 30;

   // calculate lighting only for surface oriented to the light source
   if(dot(vec3(0, 0, 1), L) > 0)
   {
      // calculate initial parameters
      float numSamplesUnderSurface    = 0;
      shadowMultiplier    = 0;
      float numLayers    = mix(maxLayers, minLayers, abs(dot(vec3(0, 0, 1), L)));
      float layerHeight    = initialHeight / numLayers;
      vec2 texStep    = parallaxScale * L.xy / L.z / numLayers;

      // current parameters
      float currentLayerHeight    = initialHeight - layerHeight;
      vec2 currentTextureCoords    = initialTexCoord + texStep;
      float heightFromTexture    = texture(u_heightTexture, currentTextureCoords).r;
      int stepIndex    = 1;

      // while point is below depth 0.0 )
      while(currentLayerHeight > 0)
      {
         // if point is under the surface
         if(heightFromTexture < currentLayerHeight)
         {
            // calculate partial shadowing factor
            numSamplesUnderSurface    += 1;
            float newShadowMultiplier    = (currentLayerHeight - heightFromTexture) *
                                             (1.0 - stepIndex / numLayers);
            shadowMultiplier    = max(shadowMultiplier, newShadowMultiplier);
         }

         // offset to the next layer
         stepIndex    += 1;
         currentLayerHeight    -= layerHeight;
         currentTextureCoords    += texStep;
         heightFromTexture    = texture(u_heightTexture, currentTextureCoords).r;
      }

      // Shadowing factor should be 1 if there were no points under the surface
      if(numSamplesUnderSurface < 1)
      {
         shadowMultiplier = 1;
      }
      else
      {
         shadowMultiplier = 1.0 - shadowMultiplier;
      }
   }
   return shadowMultiplier;
}