[ADD] - added stochastic jittering, opacity correction and some other fine tunings

src/modules/directVolumeRendering/WMDirectVolumeRendering.cpp

@@ -100,7 +100,7 @@ void WMDirectVolumeRendering::properties()
     m_propCondition = boost::shared_ptr< WCondition >( new WCondition() );
 
     m_stepCount     = m_properties->addProperty( "Step count",       "The number of steps to walk along the ray during raycasting. A low value "
-                                                                      "may cause artifacts whilst a high value slows down rendering.", 512 );
+                                                                      "may cause artifacts whilst a high value slows down rendering.", 256 );
     m_stepCount->setMin( 1 );
     m_stepCount->setMax( 5000 );
 
@@ -121,6 +121,19 @@ void WMDirectVolumeRendering::properties()
                                                     WPathHelper::getAppPath(), m_propCondition );
     WPropertyHelper::PC_PATHEXISTS::addTo( m_tfLoaderFile );
     m_tfLoaderTrigger = m_tfLoaderGroup->addProperty( "Load", "Triggers loading.", WPVBaseTypes::PV_TRIGGER_READY, m_propCondition );
+
+    // additional artifact removal methods
+    m_improvementGroup = m_properties->addPropertyGroup( "Improvements", "Methods for improving image quality. Most of these methods imply "
+                                                         "additional calculation/texture overhead and therefore slow down rendering." );
+    m_stochasticJitterEnabled = m_improvementGroup->addProperty( "Stochastic Jitter", "With stochastic jitter, wood-grain artifacts can be "
+                                                                                      "removed with the cost of possible noise artifacts.", true,
+                                                                                      m_propCondition );
+
+    m_opacityCorrectionEnabled = m_improvementGroup->addProperty( "Opacity Correction", "If enabled, opacities are assumed to be relative to the "
+                                                                                        "sample count. If disabled, changing the step count "
+                                                                                        "varies brightness of the image.", true,
+                                                                                      m_propCondition );
+
     WModule::properties();
 }
 
@@ -129,6 +142,33 @@ void WMDirectVolumeRendering::requirements()
     m_requirements.push_back( new WGERequirement() );
 }
 
+/**
+ * Generates a white noise texture with given resolution.
+ *
+ * \param resX the resolution
+ *
+ * \return a image with resX*resX resolution.
+ */
+osg::ref_ptr< osg::Image > genWhiteNoise( size_t resX )
+{
+    std::srand( time( 0 ) );
+
+    osg::ref_ptr< osg::Image > randImage = new osg::Image();
+    randImage->allocateImage( resX, resX, 1, GL_LUMINANCE, GL_UNSIGNED_BYTE );
+    unsigned char *randomLuminance = randImage->data();  // should be 4 megs
+    for( unsigned int x = 0; x < resX; x++ )
+    {
+        for( unsigned int y = 0; y < resX; y++ )
+        {
+            // - stylechecker says "use rand_r" but I am not sure about portability.
+            unsigned char r = ( unsigned char )( std::rand() % 255 );  // NOLINT
+            randomLuminance[ ( y * resX ) + x ] = r;
+        }
+    }
+
+    return randImage;
+}
+
 void WMDirectVolumeRendering::moduleMain()
 {
     m_shader = osg::ref_ptr< WShader > ( new WShader( "WMDirectVolumeRendering", m_localPath ) );
@@ -166,7 +206,8 @@ void WMDirectVolumeRendering::moduleMain()
         bool dataUpdated = m_input->updated() || m_gradients->updated();
         boost::shared_ptr< WDataSetScalar > dataSet = m_input->getData();
         bool dataValid   = ( dataSet );
-        bool propUpdated = m_localIlluminationAlgo->changed() || m_tfLoaderEnabled || m_tfLoaderFile->changed() || m_tfLoaderTrigger->changed();
+        bool propUpdated = m_localIlluminationAlgo->changed() || m_tfLoaderEnabled || m_tfLoaderFile->changed() || m_tfLoaderTrigger->changed() ||
+                           m_stochasticJitterEnabled->changed() ||  m_opacityCorrectionEnabled->changed();
 
         // As the data has changed, we need to recreate the texture.
         if ( ( propUpdated || dataUpdated ) && dataValid )
@@ -257,6 +298,37 @@ void WMDirectVolumeRendering::moduleMain()
                 m_tfLoaderTrigger->get( true );
             }
 
+            ////////////////////////////////////////////////////////////////////////////////////////////////////
+            // stochastic jittering texture
+            ////////////////////////////////////////////////////////////////////////////////////////////////////
+
+            // create some random noise
+            if ( m_stochasticJitterEnabled->get( true ) )
+            {
+                const size_t size = 64;
+                osg::ref_ptr< osg::Texture2D > randTexture = new osg::Texture2D();
+                randTexture->setFilter( osg::Texture2D::MIN_FILTER, osg::Texture2D::NEAREST );
+                randTexture->setFilter( osg::Texture2D::MAG_FILTER, osg::Texture2D::NEAREST );
+                randTexture->setWrap( osg::Texture2D::WRAP_S, osg::Texture2D::REPEAT );
+                randTexture->setWrap( osg::Texture2D::WRAP_T, osg::Texture2D::REPEAT );
+                randTexture->setImage( genWhiteNoise( size ) );
+                rootState->setTextureAttributeAndModes( 3, randTexture, osg::StateAttribute::ON );
+                rootState->addUniform( new osg::Uniform( "tex3", 3 ) );
+                m_shader->setDefine( "JITTERTEXTURE_SAMPLER", "tex3" );
+                m_shader->setDefine( "JITTERTEXTURE_SIZEX", size );
+                m_shader->setDefine( "JITTERTEXTURE_ENABLED" );
+            }
+
+            ////////////////////////////////////////////////////////////////////////////////////////////////////
+            // opacity correction
+            ////////////////////////////////////////////////////////////////////////////////////////////////////
+
+            // create some random noise
+            if ( m_opacityCorrectionEnabled->get( true ) )
+            {
+                m_shader->setDefine( "OPACITYCORRECTION_ENABLED" );
+            }
+
             ////////////////////////////////////////////////////////////////////////////////////////////////////
             // setup all those uniforms
             ////////////////////////////////////////////////////////////////////////////////////////////////////

src/modules/directVolumeRendering/WMDirectVolumeRendering.h

@@ -160,6 +160,22 @@ private:
      */
     WPropTrigger m_tfLoaderTrigger;
 
+    /**
+     * All properties for those nice improvement methods.
+     */
+    WPropGroup m_improvementGroup;
+
+    /**
+     * If true, stochastic jittering is used for image quality improvement.
+     */
+    WPropBool m_stochasticJitterEnabled;
+
+    /**
+     * If active, the opacity of the classified fragment gets scaled according to sample count to ensure relative opacities even if sampling
+     * number changes (m_stepCount)
+     */
+    WPropBool m_opacityCorrectionEnabled;
+
     /**
      * A condition used to notify about changes in several properties.
      */

src/modules/directVolumeRendering/shaders/WMDirectVolumeRendering-fragment.glsl

@@ -44,6 +44,7 @@
 uniform sampler3D tex0;
 uniform sampler1D TRANSFERFUNCTION_SAMPLER;
 uniform sampler3D GRADIENTTEXTURE_SAMPLER;
+uniform sampler2D JITTERTEXTURE_SAMPLER;
 
 /////////////////////////////////////////////////////////////////////////////
 // Attributes
@@ -62,14 +63,15 @@ uniform sampler3D GRADIENTTEXTURE_SAMPLER;
  * get a proper maximum distance along the ray.
  *
  * \param d out - this value will contain the maximum distance along the ray untill the end of the cube
+ * \param rayStart in - the start point of the ray in the volume
  *
  * \return the end point
  */
-vec3 findRayEnd( out float d )
+vec3 findRayEnd( in vec3 rayStart, out float d )
 {
     // we need to ensure the vector components are not exactly 0.0 since they are used for division
     vec3 r = v_ray + vec3( 0.000000001 );
-    vec3 p = v_rayStart;
+    vec3 p = rayStart;
 
     // v_ray in cube coordinates is used to check against the unit cube borders
     // when will v_ray reach the front face?
@@ -164,6 +166,20 @@ vec4 localIllumination( in vec3 position, in vec4 color )
 #endif
 }
 
+/**
+ * Calculates a^b. This is a faster, approximate implementation of GLSL's pow().
+ *
+ * \param a base
+ * \param b exponent
+ *
+ * \return a^b.
+ */
+float fastpow( float a, float b )
+{
+    //return exp( log(a)*b );
+    return pow( a, b );
+}
+
 /**
  * Main entry point of the fragment shader.
  */
@@ -173,22 +189,45 @@ void main()
     // when done for each vertex.
     float totalDistance = 0.0;      // the maximal distance along the ray until the BBox ends
     float currentDistance = 0.0;    // accumulated distance along the ray
-    vec3 rayEnd = findRayEnd( totalDistance );
+
+#ifdef JITTERTEXTURE_ENABLED
+    // stochastic jittering can help to void these ugly wood-grain artifacts with larger sampling distances but might
+    // introduce some noise artifacts.
+    float jitter = 0.5 - texture2D( JITTERTEXTURE_SAMPLER, gl_FragCoord.xy / JITTERTEXTURE_SIZEX ).r;
+    vec3 rayStart = v_rayStart + ( v_ray * v_stepDistance * jitter );
+#else
+    vec3 rayStart = v_rayStart;
+#endif
+
+    vec3 rayEnd = findRayEnd( rayStart, totalDistance );
 
     // walk along the ray
     vec4 dst = vec4( 0.0 );
     while ( currentDistance <= totalDistance )
     {
-        // get current value, classify and illuminate
-        vec3 rayPoint = v_rayStart + ( currentDistance * v_ray );
-        vec4 src = localIllumination( rayPoint, transferFunction( texture3D( tex0, rayPoint ).r ) );
+        // do not dynamically branch every cycle for early-ray termination, so do n steps before checking alpha value
+        for ( int i = 0; i < 4; ++i )
+        {
+            // get current value, classify and illuminate
+            vec3 rayPoint = rayStart + ( currentDistance * v_ray );
+            vec4 src = localIllumination( rayPoint, transferFunction( texture3D( tex0, rayPoint ).r ) );
+
+#ifdef OPACITYCORRECTION_ENABLED
+            // opacity correction
+            src.r = 1.0 - fastpow( 1.0 - src.r, v_relativeSteps );
+            src.g = 1.0 - fastpow( 1.0 - src.g, v_relativeSteps );
+            src.b = 1.0 - fastpow( 1.0 - src.b, v_relativeSteps );
+            src.a = 1.0 - fastpow( 1.0 - src.a, v_relativeSteps );
+#endif
 
-        // go to next value
-        currentDistance += v_stepDistance;
+            // apply front-to-back compositing
+            dst = ( 1.0 - dst.a ) * src + dst;
 
-        // apply front-to-back compositing
-        dst = ( 1.0 - dst.a ) * src + dst;
+            // go to next value
+            currentDistance += v_stepDistance;
+        }
 
+        // early ray-termination
         if ( dst.a >= 0.95 )
             break;
     }
@@ -200,9 +239,11 @@ void main()
     //    discard;
     // }
 
-    // set final color
+    // get depth ... currently the frag depth.
     float depth = gl_FragCoord.z;
+
+    // set final color
     gl_FragColor = dst;
-    gl_FragDepth = depth;                    // the depth of the last hit
+    gl_FragDepth = depth;
 }
 

src/modules/directVolumeRendering/shaders/WMDirectVolumeRendering-varyings.glsl

@@ -36,6 +36,9 @@ varying vec3 v_ray;
 // The sampling distance
 varying float v_stepDistance;
 
+// The steps in relation to a default steps of 256.
+varying float v_relativeSteps;
+
 #ifdef LOCALILLUMINATION_PHONG
 // The light source in local coordinates, normalized
 varying vec3 v_lightSource;

src/modules/directVolumeRendering/shaders/WMDirectVolumeRendering-vertex.glsl

@@ -73,6 +73,7 @@ void main()
 
     // to have equidistant sampling for each side of the box, use a fixed step size
     v_stepDistance = 1.0 / float( u_steps );
+    v_relativeSteps = 128.0 /  float( u_steps );
 
 #ifdef LOCALILLUMINATION_PHONG
     // also get the coordinates of the light