Pixel perfect volume shader

Now the volume should render with pixel perfect accuracy.
The drawback is that it is significantly slower than the old
approximate solution.

Change-Id: I4df7e9a28a27b56150c71a85a4c1fb69a215dabc
Reviewed-by: Tomi Korpipää <tomi.korpipaa@digia.com>

3 files changed, 128 insertions, 75 deletions

diff --git a/src/datavisualization/engine/abstract3drenderer.cpp b/src/datavisualization/engine/abstract3drenderer.cpp
index 090a833a..78f4b600 100644
--- a/src/datavisualization/engine/abstract3drenderer.cpp
+++ b/src/datavisualization/engine/abstract3drenderer.cpp
@@ -1320,15 +1320,15 @@ void Abstract3DRenderer::drawCustomItems(RenderingState state,
                     } else {
                         // Precalculate texture dimensions so we can optimize
                         // ray stepping to hit every texture layer.
-                        QVector3D textureDimensions(float(item->textureWidth()),
-                                                    float(item->textureHeight()),
-                                                    float(item->textureDepth()));
+                        QVector3D textureDimensions(1.0f / float(item->textureWidth()),
+                                                    1.0f / float(item->textureHeight()),
+                                                    1.0f / float(item->textureDepth()));
                         shader->setUniformValue(shader->textureDimensions(), textureDimensions);
 
-                        int sampleCount = qMax(item->textureWidth(), item->textureHeight());
-                        sampleCount = qMax(sampleCount, item->textureDepth());
+                        // Worst case scenario sample count
+                        int sampleCount = item->textureWidth() + item->textureHeight()
+                                + item->textureDepth();
                         shader->setUniformValue(shader->sampleCount(), sampleCount);
-
                     }
                     m_drawer->drawObject(shader, item->mesh(), 0, 0, item->texture());
                 } else
diff --git a/src/datavisualization/engine/shaders/texture3d.frag b/src/datavisualization/engine/shaders/texture3d.frag
index 1192ae85..6dd7f78c 100644
--- a/src/datavisualization/engine/shaders/texture3d.frag
+++ b/src/datavisualization/engine/shaders/texture3d.frag
@@ -11,7 +11,10 @@ uniform highp int sampleCount; // This is the maximum sample count
 uniform highp float alphaMultiplier;
 uniform highp int preserveOpacity;
 
-const highp float alphaThreshold = 0.0001;
+// Ray traveling straight through a single 'alpha thickness' applies 100% of the encountered alpha.
+// Rays traveling shorter distances apply a fraction. This is used to normalize the alpha over
+// entire volume, regardless of texture dimensions
+const highp float alphaThicknesses = 32.0;
 
 void main() {
     highp vec3 rayDir = -(cameraPositionRelativeToModel - pos);
@@ -41,69 +44,122 @@ void main() {
         if (rayDir.z < 0)
             boxBounds.z = -1.0;
         highp vec3 t = (boxBounds - rayStart) * invRayDir;
-        highp float minT = min(t.x, t.y);
-        minT = min(minT, t.z);
+        highp float minT = min(t.x, min(t.y, t.z));
         rayStop = rayStart + minT * rayDir;
     }
 
+    // Convert intersections to texture coords
+    rayStart = 0.5 * (rayStart + 1.0);
+    rayStop = 0.5 * (rayStop + 1.0);
+
+    highp vec3 ray = rayStop - rayStart;
+
+    // Avoid artifacts from divisions by zero
+    if (ray.x == 0)
+        ray.x = 0.000000001;
+    if (ray.y == 0)
+        ray.y = 0.000000001;
+    if (ray.z == 0)
+        ray.z = 0.000000001;
+
+    highp vec3 absRay = abs(ray);
+    highp vec3 invAbsRay = 1.0 / absRay;
+    highp float fullDist = length(ray);
+    highp vec3 curPos = rayStart;
+
+    highp vec4 curColor = vec4(0, 0, 0, 0);
+    highp float curAlpha = 0.0;
+    highp float curLen = 0.0;
+    highp vec3 curRgb = vec3(0, 0, 0);
+
     highp vec4 destColor = vec4(0, 0, 0, 0);
-    highp float totalAlpha = 0.0;
-
-    if (rayStart != rayStop) {
-        // Convert intersections to texture coords
-        rayStart = 0.5 * (rayStart + 1.0);
-        rayStop = 0.5 * (rayStop + 1.0);
-
-        highp vec3 ray = rayStop - rayStart;
-        highp float fullDist = abs(length(ray));
-        highp float rayX = abs(ray.x) * textureDimensions.x;
-        highp float rayY = abs(ray.y) * textureDimensions.y;
-        highp float rayZ = abs(ray.z) * textureDimensions.z;
-        highp float maxRayDim = max(rayX, rayY);
-        maxRayDim = max(maxRayDim, rayZ);
-        int maxCount = int(floor(maxRayDim));
-
-        highp vec3 step = ray / maxRayDim;
-        highp float stepSize = fullDist / maxRayDim;
-
-        rayStart += (step * 0.001);
-        highp vec3 curPos = rayStart;
-
-        // Adjust alpha multiplier according to the step size to get uniform alpha effect
-        // regardless of the ray angle.
-        highp float totalAlphaMultiplier = (stepSize / (1.0 / sampleCount)) * alphaMultiplier;
-
-        highp vec4 curColor = vec4(0, 0, 0, 0);
-        highp vec3 curRgb = vec3(0, 0, 0);
-        highp float curAlpha = 0.0;
-
-        // Raytrace into volume, need to sample pixels along the eye ray until we hit opacity 1
-        for (int i = 0; i < sampleCount; i++) {
-            curColor = texture3D(textureSampler, curPos);
-            if (color8Bit != 0)
-                curColor = colorIndex[int(curColor.r * 255.0)];
-
-            // Unless we have explicit alpha multiplier, we want to preserve opacity anyway
-            if (curColor.a == 1.0 && (preserveOpacity != 0 || alphaMultiplier == 1.0))
+    highp float totalOpacity = 1.0;
+    highp float nextOpacity = 1.0;
+
+    highp float extraAlphaMultiplier = fullDist * alphaThicknesses;
+
+    // nextEdges vector indicates the next edges of the texel boundaries along each axis that
+    // the ray is about to cross. The first edges are offset by a fraction of a texel to
+    // avoid artifacts from rounding errors later.
+    highp vec3 nextEdges = vec3(floor(curPos.x / textureDimensions.x) * textureDimensions.x,
+                                floor(curPos.y / textureDimensions.y) * textureDimensions.y,
+                                floor(curPos.z / textureDimensions.z) * textureDimensions.z);
+
+    highp vec3 textureSteps = textureDimensions;
+    highp vec3 textureOffset = textureDimensions * 0.001;
+    if (ray.x > 0) {
+        nextEdges.x += textureDimensions.x + textureOffset.x;
+    } else {
+        nextEdges.x -= textureOffset.x;
+        textureSteps.x = -textureDimensions.x;
+    }
+    if (ray.y > 0) {
+        nextEdges.y += textureDimensions.y + textureOffset.y;
+    } else {
+        nextEdges.y -= textureOffset.y;
+        textureSteps.y = -textureDimensions.y;
+    }
+    if (ray.z > 0) {
+        nextEdges.z += textureDimensions.z + textureOffset.z;
+    } else {
+        nextEdges.z -= textureOffset.z;
+        textureSteps.z = -textureDimensions.z;
+    }
+
+    // Raytrace into volume, need to sample pixels along the eye ray until we hit opacity 1
+    for (int i = 0; i < sampleCount; i++) {
+        curColor = texture3D(textureSampler, curPos);
+        if (color8Bit != 0)
+            curColor = colorIndex[int(curColor.r * 255.0)];
+
+        // Find which dimension has least to go to figure out the next step distance
+        highp vec3 delta = abs(nextEdges - curPos);
+        highp vec3 modDelta = delta * invAbsRay;
+        highp float minDelta = min(modDelta.x, min(modDelta.y, modDelta.z));
+        highp float stepSize;
+        if (minDelta == modDelta.x) {
+            nextEdges.x += textureSteps.x;
+            stepSize = delta.x * invAbsRay.x;
+        }
+        if (minDelta == modDelta.y) {
+            nextEdges.y += textureSteps.y;
+            stepSize = delta.y * invAbsRay.y;
+        }
+        if (minDelta == modDelta.z) {
+            nextEdges.z += textureSteps.z;
+            stepSize = delta.z * invAbsRay.z;
+        }
+
+        curPos += stepSize * ray;
+        curLen += stepSize;
+
+        if (curColor.a >= 0.0) {
+            curAlpha = alphaMultiplier * curColor.a;
+            if (curAlpha >= 1.0 || (curColor.a == 1.0 && preserveOpacity == 1))
                 curAlpha = 1.0;
             else
-                curAlpha = clamp(curColor.a * totalAlphaMultiplier, 0.0, 1.0);
-
-            if (curAlpha > alphaThreshold) {
-                curRgb = curColor.rgb * curAlpha * (1.0 - totalAlpha);
-                destColor.rgb += curRgb;
-                totalAlpha += curAlpha;
+                curAlpha = curAlpha * extraAlphaMultiplier * stepSize;
+            highp float nextOpacity = totalOpacity - curAlpha;
+            // If opacity goes beyond full opacity, we need to adjust current alpha according
+            // to the fraction of the distance the material is visible, so that we don't get
+            // box artifacts around texels.
+            if (nextOpacity < 0.0) {
+                curAlpha *= totalOpacity / curAlpha;
+                nextOpacity = 0.0;
             }
-            if (i == maxCount || totalAlpha >= 1.0)
-                break;
-            curPos += step;
+            curRgb = curColor.rgb * curAlpha * (totalOpacity + nextOpacity);
+            totalOpacity = nextOpacity;
+            destColor.rgb += curRgb;
         }
+
+        if (curLen >= 1.0 || totalOpacity <= 0.0)
+            break;
     }
 
     // Brighten up the final color if there is some transparency left
-    if (totalAlpha > alphaThreshold && totalAlpha < 1.0)
-        destColor *= 1.0 / totalAlpha;
+    if (totalOpacity >= 0.0 && totalOpacity < 1.0)
+        destColor *= (1.0 - (totalOpacity * 0.5)) / (1.0 - totalOpacity);
 
-    destColor.a = totalAlpha;
+    destColor.a = (1.0 - totalOpacity);
     gl_FragColor = clamp(destColor, 0.0, 1.0);
 }
diff --git a/src/datavisualization/engine/shaders/texture3dslice.frag b/src/datavisualization/engine/shaders/texture3dslice.frag
index 3d4c9030..00584744 100644
--- a/src/datavisualization/engine/shaders/texture3dslice.frag
+++ b/src/datavisualization/engine/shaders/texture3dslice.frag
@@ -14,8 +14,6 @@ const highp vec3 xPlaneNormal = vec3(1.0, 0, 0);
 const highp vec3 yPlaneNormal = vec3(0, 1.0, 0);
 const highp vec3 zPlaneNormal = vec3(0, 0, 1.0);
 
-const highp float alphaThreshold = 0.0001;
-
 void main() {
     // Find out where ray intersects the slice planes
     highp vec3 rayDir = -(cameraPositionRelativeToModel - pos);
@@ -24,7 +22,7 @@ void main() {
     // Flip Y and Z so QImage bits work directly for texture and first image is in the front
     rayStart.yz = -rayStart.yz;
     rayDir.yz = -rayDir.yz;
-    highp float tFar = 2.0f;
+    highp float minT = 2.0f;
     if (rayDir.x != 0.0 && rayDir.y != 0.0 && rayDir.z != 0.0) {
         highp vec3 boxBounds = vec3(1.0, 1.0, 1.0);
         highp vec3 invRayDir = 1.0 / rayDir;
@@ -35,24 +33,23 @@ void main() {
         if (rayDir.z < 0)
             boxBounds.z = -1.0;
         highp vec3 t = (boxBounds - rayStart) * invRayDir;
-        tFar = min(t.x, t.y);
-        tFar = min(tFar, t.z);
+        minT = min(t.x, min(t.y, t.z));
     }
 
     highp vec3 xPoint = vec3(volumeSliceIndices.x, 0, 0);
     highp vec3 yPoint = vec3(0, volumeSliceIndices.y, 0);
     highp vec3 zPoint = vec3(0, 0, volumeSliceIndices.z);
-    highp float firstD = tFar + 1.0;
+    highp float firstD = minT + 1.0;
     highp float secondD = firstD;
     highp float thirdD = firstD;
     if (volumeSliceIndices.x >= -1.0) {
         highp float dx = dot(xPoint - rayStart, xPlaneNormal) / dot(rayDir, xPlaneNormal);
-        if (dx >= 0.0 && dx <= tFar)
+        if (dx >= 0.0 && dx <= minT)
             firstD = min(dx, firstD);
     }
     if (volumeSliceIndices.y >= -1.0) {
         highp float dy = dot(yPoint - rayStart, yPlaneNormal) / dot(rayDir, yPlaneNormal);
-        if (dy >= 0.0 && dy <= tFar) {
+        if (dy >= 0.0 && dy <= minT) {
             if (dy < firstD) {
                 secondD = firstD;
                 firstD = dy;
@@ -64,7 +61,7 @@ void main() {
     if (volumeSliceIndices.z >= -1.0) {
         highp float dz = dot(zPoint - rayStart, zPlaneNormal) / dot(rayDir, zPlaneNormal);
         if (dz >= 0.0) {
-            if (dz < firstD && dz <= tFar) {
+            if (dz < firstD && dz <= minT) {
                 thirdD = secondD;
                 secondD = firstD;
                 firstD = dz;
@@ -85,14 +82,14 @@ void main() {
 
     // Convert intersection to texture coords
 
-    if (firstD <= tFar) {
+    if (firstD <= minT) {
         highp vec3 firstTex = rayStart + rayDir * firstD;
         firstTex = 0.5 * (firstTex + 1.0);
         curColor = texture3D(textureSampler, firstTex);
         if (color8Bit != 0)
             curColor = colorIndex[int(curColor.r * 255.0)];
 
-        if (curColor.a > alphaThreshold) {
+        if (curColor.a > 0.0) {
             curAlpha = curColor.a;
             if (curColor.a == 1.0 && preserveOpacity != 0)
                 curAlpha = 1.0;
@@ -101,13 +98,13 @@ void main() {
             destColor.rgb = curColor.rgb * curAlpha;
             totalAlpha = curAlpha;
         }
-        if (secondD <= tFar && totalAlpha < 1.0) {
+        if (secondD <= minT && totalAlpha < 1.0) {
             highp vec3 secondTex = rayStart + rayDir * secondD;
             secondTex = 0.5 * (secondTex + 1.0);
             curColor = texture3D(textureSampler, secondTex);
             if (color8Bit != 0)
                 curColor = colorIndex[int(curColor.r * 255.0)];
-            if (curColor.a > alphaThreshold) {
+            if (curColor.a > 0.0) {
                 if (curColor.a == 1.0 && preserveOpacity != 0)
                     curAlpha = 1.0;
                 else
@@ -116,11 +113,11 @@ void main() {
                 destColor.rgb += curRgb;
                 totalAlpha += curAlpha;
             }
-            if (thirdD <= tFar && totalAlpha < 1.0) {
+            if (thirdD <= minT && totalAlpha < 1.0) {
                 highp vec3 thirdTex = rayStart + rayDir * thirdD;
                 thirdTex = 0.5 * (thirdTex + 1.0);
                 curColor = texture3D(textureSampler, thirdTex);
-                if (curColor.a > alphaThreshold) {
+                if (curColor.a > 0.0) {
                     if (color8Bit != 0)
                         curColor = colorIndex[int(curColor.r * 255.0)];
                     if (curColor.a == 1.0 && preserveOpacity != 0)
@@ -136,7 +133,7 @@ void main() {
     }
 
     // Brighten up the final color if there is some transparency left
-    if (totalAlpha > alphaThreshold && totalAlpha < 1.0)
+    if (totalAlpha > 0.0 && totalAlpha < 1.0)
         destColor *= 1.0 / totalAlpha;
 
     destColor.a = totalAlpha;