* Introducing ContextPreservingRaycaster processor

* Updated ProcessorDecoratorGradient to support controllable LOD for gradient computation
* updated RaycastingProcessor to compute min/max depth values of entry/exit points and set as uniforms if needed by shader

core/glsl/tools/gradient.frag

@@ -27,31 +27,55 @@
 /**
  * Compute the gradient using forward differences on the texture's red channel.
  * \param   tex         3D texture to calculate gradients for
+ * \param   texParams   Texture parameters struct (needs texParams._voxelSize)
  * \param   texCoords   Lookup position in texture coordinates
+ * \param   lod         Explicit LOD in texture to look up
  */
-vec3 computeGradientForwardDifferences(in sampler3D tex, in TextureParameters3D texParams, in vec3 texCoords) {
-    float v = texture(tex, texCoords).r;
-    float dx = textureOffset(tex, texCoords, ivec3(1, 0, 0)).r;
-    float dy = textureOffset(tex, texCoords, ivec3(0, 1, 0)).r;
-    float dz = textureOffset(tex, texCoords, ivec3(0, 0, 1)).r;
+vec3 computeGradientForwardDifferencesLod(in sampler3D tex, in TextureParameters3D texParams, in vec3 texCoords, in float lod) {
+    float v = textureLod(tex, texCoords, lod).r;
+    float dx = textureLodOffset(tex, texCoords, lod, ivec3(1, 0, 0)).r;
+    float dy = textureLodOffset(tex, texCoords, lod, ivec3(0, 1, 0)).r;
+    float dz = textureLodOffset(tex, texCoords, lod, ivec3(0, 0, 1)).r;
     return vec3(v - dx, v - dy, v - dz) * texParams._voxelSize;
 }
 
+/**
+ * Compute the gradient using forward differences on the texture's red channel.
+ * \param   tex         3D texture to calculate gradients for
+ * \param   texParams   Texture parameters struct (needs texParams._voxelSize)
+ * \param   texCoords   Lookup position in texture coordinates
+ */
+vec3 computeGradientForwardDifferences(in sampler3D tex, in TextureParameters3D texParams, in vec3 texCoords) {
+    return computeGradientForwardDifferencesLod(tex, texParams, texCoords, 0.0);
+}
+
 /**
  * Compute the gradient using central differences on the texture's red channel.
  * \param   tex         3D texture to calculate gradients for
+ * \param   texParams   Texture parameters struct (needs texParams._voxelSize)
  * \param   texCoords   Lookup position in texture coordinates
+ * \param   lod         Explicit LOD in texture to look up
  */
-vec3 computeGradientCentralDifferences(in sampler3D tex, in TextureParameters3D texParams, in vec3 texCoords) {
-    float dx = textureOffset(tex, texCoords, ivec3(1, 0, 0)).r;
-    float dy = textureOffset(tex, texCoords, ivec3(0, 1, 0)).r;
-    float dz = textureOffset(tex, texCoords, ivec3(0, 0, 1)).r;
-    float mdx = textureOffset(tex, texCoords, ivec3(-1, 0, 0)).r;
-    float mdy = textureOffset(tex, texCoords, ivec3(0, -1, 0)).r;
-    float mdz = textureOffset(tex, texCoords, ivec3(0, 0, -1)).r;
+vec3 computeGradientCentralDifferencesLod(in sampler3D tex, in TextureParameters3D texParams, in vec3 texCoords, in float lod) {
+    float dx = textureLodOffset(tex, texCoords, lod, ivec3(1, 0, 0)).r;
+    float dy = textureLodOffset(tex, texCoords, lod, ivec3(0, 1, 0)).r;
+    float dz = textureLodOffset(tex, texCoords, lod, ivec3(0, 0, 1)).r;
+    float mdx = textureLodOffset(tex, texCoords, lod, ivec3(-1, 0, 0)).r;
+    float mdy = textureLodOffset(tex, texCoords, lod, ivec3(0, -1, 0)).r;
+    float mdz = textureLodOffset(tex, texCoords, lod, ivec3(0, 0, -1)).r;
     return vec3(mdx - dx, mdy - dy, mdz - dz) * texParams._voxelSize * 0.5;
 }
 
+/**
+ * Compute the gradient using central differences on the texture's red channel.
+ * \param   tex         3D texture to calculate gradients for
+ * \param   texParams   Texture parameters struct (needs texParams._voxelSize)
+ * \param   texCoords   Lookup position in texture coordinates
+ */
+vec3 computeGradientCentralDifferences(in sampler3D tex, in TextureParameters3D texParams, in vec3 texCoords) {
+    return computeGradientCentralDifferencesLod(tex, texParams, texCoords, 0.0);
+}
+
 
 /**
  * Compute the gradient using filtered central differences on the texture's red channel.
@@ -178,3 +202,100 @@ vec3 computeGradientSobel(in sampler3D tex, in vec3 texCoords) {
    
     return -sobelScale * sobel;
 }
+
+
+/**
+ * Compute the Hessian matrix of the input image using central differences twice.
+ * \param   tex         3D gradient texture to calculate Hessian for
+ * \param   texCoords   Lookup position in texture coordinates
+ */
+mat3 computeHessianCentralDifferencesLod(in sampler3D tex, in TextureParameters3D texParams, in vec3 texCoords, in float lod) {
+    vec3 offset = texParams._sizeRCP * pow(2.0, lod);
+
+    // compute first derivatives using central differences
+    vec3 dx  = computeGradientCentralDifferencesLod(tex, texParams, texCoords + vec3(offset.x, 0.0, 0.0), lod);
+    vec3 dy  = computeGradientCentralDifferencesLod(tex, texParams, texCoords + vec3(0.0, offset.y, 0.0), lod);
+    vec3 dz  = computeGradientCentralDifferencesLod(tex, texParams, texCoords + vec3(0.0, 0.0, offset.z), lod);
+    vec3 mdx = computeGradientCentralDifferencesLod(tex, texParams, texCoords - vec3(offset.x, 0.0, 0.0), lod);
+    vec3 mdy = computeGradientCentralDifferencesLod(tex, texParams, texCoords - vec3(0.0, offset.y, 0.0), lod);
+    vec3 mdz = computeGradientCentralDifferencesLod(tex, texParams, texCoords - vec3(0.0, 0.0, offset.z), lod);
+
+    // compute hessian matrix columns (second order derivatives) using central differences of gradients
+    vec3 vdx = (dx - mdx) * texParams._voxelSize * 0.5;
+    vec3 vdy = (dy - mdy) * texParams._voxelSize * 0.5;
+    vec3 vdz = (dz - mdz) * texParams._voxelSize * 0.5;
+
+    return mat3(vdx, vdy, vdz);
+}
+
+/**
+ * Compute the Hessian matrix using central differences on the gradient texture's RGB channel.
+ * \param   tex         3D gradient texture to calculate Hessian for
+ * \param   texCoords   Lookup position in texture coordinates
+ */
+mat3 computeGradientHessianCentralDifferencesLod(in sampler3D tex, in TextureParameters3D texParams, in vec3 texCoords, in float lod) {
+    vec3 offset = texParams._sizeRCP * pow(2.0, lod);
+
+    vec3 dx = textureLod(tex, texCoords + vec3(offset.x, 0.0, 0.0), lod).rgb;
+    vec3 dy = textureLod(tex, texCoords + vec3(0.0, offset.y, 0.0), lod).rgb;
+    vec3 dz = textureLod(tex, texCoords + vec3(0.0, 0.0, offset.z), lod).rgb;
+    vec3 mdx = textureLod(tex, texCoords - vec3(offset.x, 0.0, 0.0), lod).rgb;
+    vec3 mdy = textureLod(tex, texCoords - vec3(0.0, offset.y, 0.0), lod).rgb;
+    vec3 mdz = textureLod(tex, texCoords - vec3(0.0, 0.0, offset.z), lod).rgb;
+
+    float dxx = dx.x - mdx.x;
+    float dxy = dx.y - mdx.y;
+    float dxz = dx.z - mdx.z;
+
+    float dyx = dy.x - mdy.x;
+    float dyy = dy.y - mdy.y;
+    float dyz = dy.z - mdy.z;
+
+    float dzx = dz.x - mdz.x;
+    float dzy = dz.y - mdz.y;
+    float dzz = dz.z - mdz.z;
+
+    mat3 hessian = mat3(dxx, (dyx + dxy) / 2.0, (dzx + dxz) / 2.0,
+                        (dxy + dyx) / 2.0, dyy, (dzy + dyz) / 2.0,
+                        (dxz + dzx) / 2.0, (dyz + dzy) / 2.0, dzz);
+    //hessian *= (texParams._voxelSize * 0.5;
+    return hessian;
+}
+
+
+/**
+ * Compute the Hessian matrix using central differences on the gradient texture's RGB channel.
+ * \param   tex         3D gradient texture to calculate Hessian for
+ * \param   texCoords   Lookup position in texture coordinates
+ */
+mat3 computeHessianCentralDifferences(in sampler3D tex, in TextureParameters3D texParams, in vec3 texCoords) {
+    return computeHessianCentralDifferencesLod(tex, texParams, texCoords, 0.0);
+}
+
+// code inspiured by http://en.wikipedia.org/wiki/Eigenvalue_algorithm#3.C3.973_matrices
+vec3 computeEigenvalues(in mat3 hessian) {
+    const float PI = 3.1415926535897932384626433832795;
+    vec3 toReturn;
+    float p1 = (hessian[1][0] * hessian[1][0]) + (hessian[2][0] * hessian[2][0]) + (hessian[2][1] * hessian[2][1]);
+
+    if (p1 == 0.0) {
+        toReturn.x = hessian[0][0];
+        toReturn.y = hessian[1][1];
+        toReturn.z = hessian[2][2];
+    }
+    else {
+        float q = (hessian[0][0] + hessian[1][1] + hessian[2][2]) / 3.0;
+        float p2 = (hessian[0][0] - q) * (hessian[0][0] - q)   +   (hessian[1][1] - q) * (hessian[1][1] - q)   +   (hessian[2][2] - q) * (hessian[2][2] - q)   +   2.0 * p1;
+        float p = sqrt(p2 / 6.0);
+        mat3 B = (1.0 - p) * (hessian - mat3(q));
+        float r = determinant(B) / 2.0;
+
+        float phi = acos(clamp(r, -1.0, 1.0)) / 3.0;
+
+        toReturn.x = q + (2.0 * p * cos(phi));
+        toReturn.z = q + (2.0 * p * cos(phi + (2.0 * PI / 3.0)));
+        toReturn.y = 3.0 * q - toReturn.x - toReturn.z;
+    }
+
+    return toReturn;
+}
\ No newline at end of file

core/glsl/tools/shading.frag

@@ -48,29 +48,16 @@ float computeAttenuation(in vec3 attenuation, in float d) {
     return min(att, 1.0);
 }
 
-
-/**
- * Computes the ambient term according to the given material color \a ka and ambient intensity \a ia.
- *
- * \param   ka  Material ambient color
- * \param   ia  Ambient light intensity
- */
-vec3 getAmbientTerm(in vec3 ka, in vec3 ia) {
-    return ka * ia;
-}
-
-
 /**
  * Computes the diffuse term according to the given parameters.
  *
- * \param   kd  Material diffuse color
  * \param   id  Diffuse light intensity
  * \param   N   Surface normal
  * \param   L   Normalized light vector
  */
-vec3 getDiffuseTerm(in vec3 kd, in vec3 id, in vec3 N, in vec3 L) {
+vec3 getDiffuseTerm(in vec3 id, in vec3 N, in vec3 L) {
     float NdotL = max(dot(N, L), 0.0);
-    return kd * id * NdotL;
+    return id * NdotL;
 }
 
 
@@ -84,10 +71,10 @@ vec3 getDiffuseTerm(in vec3 kd, in vec3 id, in vec3 N, in vec3 L) {
  * \param   V           View vector
  * \param   shininess   Shininess coefficient
  */
-vec3 getSpecularTerm(in vec3 ks, in vec3 is, in vec3 N, in vec3 L, in vec3 V, in float shininess) {
+vec3 getSpecularTerm(in vec3 is, in vec3 N, in vec3 L, in vec3 V, in float shininess) {
     vec3 H = normalize(V + L);
     float NdotH = pow(max(dot(N, H), 0.0), shininess);
-    return ks * is * NdotH;
+    return is * NdotH;
 }
 
 /**
@@ -110,9 +97,9 @@ vec3 calculatePhongShading(in vec3 position, in LightSource light, in vec3 camer
     float d = length(L);
     L /= d;
 
-    vec3 toReturn = getAmbientTerm(ka, light._ambientColor);
-    toReturn += getDiffuseTerm(kd, light._diffuseColor, N, L);
-    toReturn += getSpecularTerm(ks, light._specularColor, N, L, V, light._shininess);
+    vec3 toReturn = ka * light._ambientColor;   // ambient term
+    toReturn += kd * getDiffuseTerm(light._diffuseColor, N, L);
+    toReturn += ks * getSpecularTerm(light._specularColor, N, L, V, light._shininess);
     #ifdef PHONG_APPLY_ATTENUATION
         toReturn *= computeAttenuation(light._attenuation, d);
     #endif
@@ -137,11 +124,38 @@ vec3 calculatePhongShading(in vec3 position, in LightSource light, in vec3 camer
     float d = length(L);
     L /= d;
 
-    vec3 toReturn = getAmbientTerm(materialColor, light._ambientColor);
-    toReturn += getDiffuseTerm(materialColor, light._diffuseColor, N, L);
-    toReturn += getSpecularTerm(materialColor, light._specularColor, N, L, V, light._shininess);
+    vec3 toReturn = materialColor * light._ambientColor;   // ambient term
+    toReturn += materialColor * getDiffuseTerm(light._diffuseColor, N, L);
+    toReturn += materialColor * getSpecularTerm(light._specularColor, N, L, V, light._shininess);
     #ifdef PHONG_APPLY_ATTENUATION
         toReturn *= computeAttenuation(light._attenuation, d);
     #endif
     return toReturn;
 }
+
+
+/**
+ * Computes the phong shading intensity according to the given parameters.
+ * 
+ * \param   position        sample position
+ * \param   light           LightSource
+ * \param   camera          Camera position
+ * \param   normal          Normal
+ */
+float getPhongShadingIntensity(in vec3 position, in LightSource light, in vec3 camera, in vec3 normal) {
+    vec3 N = normalize(normal);
+    vec3 V = normalize(camera - position);
+
+    // get light source distance for attenuation and normalize light vector
+    vec3 L = light._position - position;
+    float d = length(L);
+    L /= d;
+
+    vec3 toReturn = light._ambientColor;   // ambient term
+    toReturn += getDiffuseTerm(light._diffuseColor, N, L);
+    toReturn += getSpecularTerm(light._specularColor, N, L, V, light._shininess);
+    #ifdef PHONG_APPLY_ATTENUATION
+        toReturn *= computeAttenuation(light._attenuation, d);
+    #endif
+    return (toReturn.x + toReturn.y + toReturn.z) / 3.0;
+}
\ No newline at end of file

core/pipeline/processordecoratorgradient.cpp

@@ -37,26 +37,32 @@ namespace campvis {
 
     ProcessorDecoratorGradient::ProcessorDecoratorGradient()
         : AbstractProcessorDecorator()
-        , _gradientMethod("GradientMethod", "Gradient Computation Method", gradientOptions, 4, AbstractProcessor::INVALID_SHADER | AbstractProcessor::INVALID_RESULT)
+        , p_gradientMethod("GradientMethod", "Gradient Computation Method", gradientOptions, 4, AbstractProcessor::INVALID_SHADER | AbstractProcessor::INVALID_RESULT)
+        , p_lod("GradientLod", "LOD for Gradient Computation", 0.5f, 0.f, 5.f, .1f, 1)
     {
-        _gradientMethod.setValue(1);
+        p_gradientMethod.setValue(1);
+        p_gradientMethod.s_changed.connect(this, &ProcessorDecoratorGradient::onGradientMethodChanged);
     }
 
     ProcessorDecoratorGradient::~ProcessorDecoratorGradient() {
+        p_gradientMethod.s_changed.disconnect(this);
     }
 
     void ProcessorDecoratorGradient::addProperties(HasPropertyCollection* propCollection) {
-        propCollection->addProperty(&_gradientMethod);
+        propCollection->addProperty(&p_gradientMethod);
+        propCollection->addProperty(&p_lod);
     }
 
     std::string ProcessorDecoratorGradient::generateHeader() const {
         std::string toReturn;
-        switch (_gradientMethod.getOptionValue()) {
+        switch (p_gradientMethod.getOptionValue()) {
             case ForwardDifferences:
-                toReturn.append("#define computeGradient(tex, texParams, texCoords) computeGradientForwardDifferences(tex, texParams, texCoords)\n");
+                toReturn.append("#define computeGradient(tex, texParams, texCoords) computeGradientForwardDifferencesLod(tex, texParams, texCoords, _gradientLod)\n");
+                toReturn.append("uniform float _gradientLod = 0.0;\n");
                 break;
             case CentralDifferences:
-                toReturn.append("#define computeGradient(tex, texParams, texCoords) computeGradientCentralDifferences(tex, texParams, texCoords)\n");
+                toReturn.append("#define computeGradient(tex, texParams, texCoords) computeGradientCentralDifferencesLod(tex, texParams, texCoords, _gradientLod)\n");
+                toReturn.append("uniform float _gradientLod = 0.0;\n");
                 break;
             case FilteredCentralDifferences:
                 toReturn.append("#define computeGradient(tex, texParams, texCoords) computeGradientFilteredCentralDifferences(tex, texParams, texCoords)\n");
@@ -72,4 +78,14 @@ namespace campvis {
         return toReturn;
     }
 
+    void ProcessorDecoratorGradient::renderProlog(const DataContainer& dataContainer, tgt::Shader* shader) {
+        if (p_gradientMethod.getOptionValue() == ForwardDifferences || p_gradientMethod.getOptionValue() == CentralDifferences) {
+            shader->setUniform("_gradientLod", p_lod.getValue());
+        }
+    }
+
+    void ProcessorDecoratorGradient::onGradientMethodChanged(const AbstractProperty* prop) {
+        p_lod.setVisible(p_gradientMethod.getOptionValue() == ForwardDifferences || p_gradientMethod.getOptionValue() == CentralDifferences);
+    }
+
 }

core/pipeline/processordecoratorgradient.h

@@ -25,7 +25,9 @@
 #ifndef PROCESSORDECORATORGRADIENT_H__
 #define PROCESSORDECORATORGRADIENT_H__
 
+#include "sigslot/sigslot.h"
 #include "tgt/textureunit.h"
+
 #include "core/pipeline/abstractprocessordecorator.h"
 #include "core/properties/floatingpointproperty.h"
 #include "core/properties/genericproperty.h"
@@ -39,7 +41,7 @@ namespace campvis {
      * generateHeader() to define computeGradient(tex, texParams, texCoords) in GLSL calling
      * the selected function.
      */
-    class CAMPVIS_CORE_API ProcessorDecoratorGradient : public AbstractProcessorDecorator {
+    class CAMPVIS_CORE_API ProcessorDecoratorGradient : public AbstractProcessorDecorator, public sigslot::has_slots<> {
     public:
         /// Method for online-calculating gradients
         enum GradientMethod {
@@ -57,10 +59,17 @@ namespace campvis {
     protected:
         /// \see AbstractProcessorDecorator::addProperties()
         void addProperties(HasPropertyCollection* propCollection);
+        /// \see AbstractProcessorDecorator::renderProlog()
+        virtual void renderProlog(const DataContainer& dataContainer, tgt::Shader* shader);
         /// \see AbstractProcessorDecorator::generateHeader()
         std::string generateHeader() const;
 
-        GenericOptionProperty<GradientMethod> _gradientMethod;  ///< Method for calculating the gradients
+        /// Callback method when p_gradientMethod has changed
+        void onGradientMethodChanged(const AbstractProperty* prop);
+
+        GenericOptionProperty<GradientMethod> p_gradientMethod;  ///< Method for calculating the gradients
+        FloatProperty p_lod;        ///< LOD to use for texture lookup during gradient computation
+
     };
 
 }

core/pipeline/raycastingprocessor.cpp

@@ -30,6 +30,7 @@
 
 #include "core/datastructures/imagedata.h"
 #include "core/datastructures/renderdata.h"
+#include "core/tools/glreduction.h"
 
 #include "core/classification/simpletransferfunction.h"
 
@@ -70,11 +71,18 @@ namespace campvis {
         _shader = ShdrMgr.loadWithCustomGlslVersion("core/glsl/passthrough.vert", "", _fragmentShaderFilename, generateHeader(), _customGlslVersion);
         _shader->setAttributeLocation(0, "in_Position");
         _shader->setAttributeLocation(1, "in_TexCoord");
+
+        _minReduction = new GlReduction(GlReduction::MIN);
+        _maxReduction = new GlReduction(GlReduction::MAX);
     }
 
     void RaycastingProcessor::deinit() {
         ShdrMgr.dispose(_shader);
         _shader = 0;
+
+        delete _minReduction;
+        delete _maxReduction;
+
         VisualizationProcessor::deinit();
     }
 
@@ -85,15 +93,45 @@ namespace campvis {
 
         if (img != 0 && entryPoints != 0 && exitPoints != 0) {
             if (img->getDimensionality() == 3) {
+                // little hack to support LOD texture lookup for the gradients:
+                // if texture does not yet have mipmaps, create them.
+                const tgt::Texture* tex = img->getTexture();
+                if (tex->getFilter() != tgt::Texture::MIPMAP) {
+                    const_cast<tgt::Texture*>(tex)->setFilter(tgt::Texture::MIPMAP);
+                    glGenerateMipmap(GL_TEXTURE_3D);
+                    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
+                    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR);
+                    LGL_ERROR;
+                }
+
                 _shader->activate();
                 _shader->setIgnoreUniformLocationError(true);
+
+                // Compute min/max depth if needed by shader
+                if (_shader->getUniformLocation("_minDepth") != -1) {
+                    _shader->deactivate();
+                    float minDepth = _minReduction->reduce(entryPoints->getDepthTexture()).front();
+                    _shader->activate();
+
+                    _shader->setUniform("_minDepth", minDepth);
+                }
+                if (_shader->getUniformLocation("_maxDepth") != -1) {
+                    _shader->deactivate();
+                    float maxDepth = _maxReduction->reduce(exitPoints->getDepthTexture()).front();
+                    _shader->activate();
+
+                    _shader->setUniform("_maxDepth", maxDepth);
+                }
+
                 decorateRenderProlog(data, _shader);
                 _shader->setUniform("_viewportSizeRCP", 1.f / tgt::vec2(getEffectiveViewportSize()));
                 _shader->setUniform("_jitterStepSizeMultiplier", p_jitterStepSizeMultiplier.getValue());
 
+                // compute sampling step size relative to volume size
                 float samplingStepSize = 1.f / (p_samplingRate.getValue() * tgt::max(img->getSize()));
                 _shader->setUniform("_samplingStepSize", samplingStepSize);
 
+                // compute and set camera parameters
                 const tgt::Camera& cam = p_camera.getValue();
                 float n = cam.getNearDist();
                 float f = cam.getFarDist();
@@ -104,6 +142,7 @@ namespace campvis {
                 _shader->setUniform("const_to_z_w_2", 0.5f*((f+n)/(f-n))+0.5f);
                 _shader->setIgnoreUniformLocationError(false);
 
+                // bind input textures
                 tgt::TextureUnit volumeUnit, entryUnit, exitUnit, tfUnit;
                 img->bind(_shader, volumeUnit, "_volume", "_volumeTextureParams");
                 p_transferFunction.getTF()->bind(_shader, tfUnit);

core/pipeline/raycastingprocessor.h

@@ -42,6 +42,8 @@ namespace tgt {
 }
 
 namespace campvis {
+    class GlReduction;
+
     /**
      * Base class for raycasting processors.
      * Offfers properties various common properties and automatic shader loading/linking.
@@ -138,6 +140,9 @@ namespace campvis {
         tgt::Shader* _shader;                           ///< Shader for raycasting
         bool _bindEntryExitDepthTextures;               ///< Flag whether to also bind the depth textures of the entry-/exit points.
 
+        GlReduction* _minReduction;
+        GlReduction* _maxReduction;
+
         static const std::string loggerCat_;
     };
 

modules/vis/glsl/contextpreservingraycaster.frag

+// ================================================================================================
+// 
+// This file is part of the CAMPVis Software Framework.
+// 
+// If not explicitly stated otherwise: Copyright (C) 2012-2013, all rights reserved,
+//      Christian Schulte zu Berge <christian.szb@in.tum.de>
+//      Chair for Computer Aided Medical Procedures
+//      Technische Universität München
+//      Boltzmannstr. 3, 85748 Garching b. München, Germany
+// 
+// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
+// 
+// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file 
+// except in compliance with the License. You may obtain a copy of the License at
+// 
+// http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software distributed under the 
+// License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, 
+// either express or implied. See the License for the specific language governing permissions 
+// and limitations under the License.
+// 
+// ================================================================================================
+
+layout(location = 0) out vec4 out_Color;     ///< outgoing fragment color
+layout(location = 1) out vec4 out_FHP;       ///< outgoing fragment first hitpoint
+layout(location = 2) out vec4 out_FHN;       ///< outgoing fragment first hit normal
+
+#include "tools/gradient.frag"
+#include "tools/raycasting.frag"
+#include "tools/shading.frag"
+#include "tools/texture2d.frag"
+#include "tools/texture3d.frag"
+#include "tools/transferfunction.frag"
+
+uniform vec2 _viewportSizeRCP;
+uniform float _jitterStepSizeMultiplier;
+
+
+// ray entry points
+uniform sampler2D _entryPoints;
+uniform sampler2D _entryPointsDepth;
+uniform TextureParameters2D _entryParams;
+
+// ray exit points
+uniform sampler2D _exitPoints;
+uniform sampler2D _exitPointsDepth;
+uniform TextureParameters2D _exitParams;
+
+// DRR volume
+uniform sampler3D _volume;
+uniform TextureParameters3D _volumeTextureParams;
+
+// Transfer function
+uniform sampler1D _transferFunction;
+uniform TFParameters1D _transferFunctionParams;
+
+
+uniform LightSource _lightSource;
+uniform vec3 _cameraPosition;
+
+uniform float _samplingStepSize;
+
+uniform float _minDepth;    ///< Minimum depth of entry points
+uniform float _maxDepth;    ///< Maximum depth of exit points
+
+uniform float _kappaS;      ///< k_s parameter from the paper
+uniform float _kappaT;      ///< l_t parameter from the paper
+
+// TODO: copy+paste from Voreen - eliminate or improve.
+const float SAMPLING_BASE_INTERVAL_RCP = 200.0;
+
+/**
+ * Performs the raycasting and returns the final fragment color.
+ */
+vec4 performRaycasting(in vec3 entryPoint, in vec3 exitPoint, in vec2 texCoords) {
+    vec4 result = vec4(0.0);
+    out_FHP = vec4(0.0);
+    float firstHitT = -1.0;
+
+    // calculate ray parameters
+    vec3 direction = exitPoint.rgb - entryPoint.rgb;
+    float t = 0.0;
+    float tend = length(direction);
+    direction = normalize(direction);
+
+    jitterEntryPoint(entryPoint, direction, _samplingStepSize * _jitterStepSizeMultiplier);
+
+    while (t < tend) {
+        // compute sample position
+        vec3 samplePosition = entryPoint.rgb + t * direction;
+
+        // lookup intensity and TF
+        float intensity = texture(_volume, samplePosition).r;
+        vec4 color = lookupTF(_transferFunction, _transferFunctionParams, intensity);
+
+        float depthEntry = texture(_entryPointsDepth, texCoords).z;
+        float depthExit = texture(_exitPointsDepth, texCoords).z;
+        float D = calculateDepthValue(firstHitT/tend, depthEntry, depthExit);
+
+        // perform compositing
+        if (color.a > 0.0) {
+            // compute gradient (needed for shading and normals)
+            vec3 gradient = computeGradient(_volume, _volumeTextureParams, samplePosition) * 10;
+
+            color.rgb = calculatePhongShading(textureToWorld(_volumeTextureParams, samplePosition).xyz, _lightSource, _cameraPosition, gradient, color.rgb, color.rgb, vec3(1.0, 1.0, 1.0));
+            float SI = getPhongShadingIntensity(textureToWorld(_volumeTextureParams, samplePosition).xyz, _lightSource, _cameraPosition, gradient);
+
+            // accomodate for variable sampling rates
+            color.a = 1.0 - pow(1.0 - color.a, _samplingStepSize * SAMPLING_BASE_INTERVAL_RCP);
+
+            // perform context-preserving rendering (the meat is modifying the alpha value of each sample)
+            float pwr = pow(_kappaT * SI * (1 - D) * (1 - result.a), _kappaS);
+            color.a *= pow(length(gradient), pwr);
+
+            result.rgb = mix(color.rgb, result.rgb, result.a);
+            result.a = result.a + (1.0 -result.a) * color.a;
+        }
+
+        // save first hit ray parameter for depth value calculation
+        if (firstHitT < 0.0 && result.a > 0.0) {
+            firstHitT = t;
+            out_FHP = vec4(samplePosition, 1.0);
+            out_FHN = vec4(normalize(computeGradient(_volume, _volumeTextureParams, samplePosition)), 1.0);
+        }