/*========================================================================= Program: Visualization Toolkit Module: vtkGPUVolumeRayCastMapper_HeaderFS.glsl Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen All rights reserved. See Copyright.txt or http://www.kitware.com/Copyright.htm for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notice for more information. =========================================================================*/ // High quality volume renderer for axis-aligned hexahedra // Implementation by Stephane Marchesin (stephane.marchesin@gmail.com) // CEA/DIF - Commissariat a l'Energie Atomique, Centre DAM Ile-De-France // BP12, F-91297 Arpajon, France. // // This file implements the paper // "High-Quality, Semi-Analytical Volume Rendering for AMR Data", // Stephane Marchesin and Guillaume Colin de Verdiere, IEEE Vis 2009. // inputs of FS // texcoord0.xyz : position in the brick space [0,1]^3 // texcoord0.w : cell size // texcoord1.xyz : position in object space // texcoord2.xyzw: node data 0,1,2,3 // texcoord3.xyzw: node data 4,5,6,7 uniform sampler3D preintegration_table; uniform vec3 observer; uniform float length_max; vec4 sample(float sample0, float sample1, float length) { float corrected_length = length * gl_TexCoord[0].w / length_max ; return texture3D(preintegration_table, vec3(sample0,sample1,corrected_length)); } void main() { vec3 pos = gl_TexCoord[0].xyz; vec3 progression; vec3 dist1,dist2,dist; vec3 l=vec3(1.0,0.0,0.0); float length; float cell_length = gl_TexCoord[0].w; progression.xyz = gl_TexCoord[1].xyz - observer.xyz; progression = normalize(progression); dist1.xyz = abs((1.0-pos.xyz)/progression.xyz); dist2.xyz = abs((pos.xyz)/progression.xyz); if (progression.x>0.0) dist.x=dist1.x; else dist.x=dist2.x; if (progression.y>0.0) dist.y=dist1.y; else dist.y=dist2.y; if (progression.z>0.0) dist.z=dist1.z; else dist.z=dist2.z; length = min(dist.x,min(dist.y,dist.z)); vec3 p1 = pos, p2 = pos + vec3(length) * progression; float s0 = gl_TexCoord[2].x; float s1 = gl_TexCoord[2].y; float s2 = gl_TexCoord[2].z; float s3 = gl_TexCoord[2].w; float s4 = gl_TexCoord[3].x; float s5 = gl_TexCoord[3].y; float s6 = gl_TexCoord[3].z; float s7 = gl_TexCoord[3].w; float x0 = p1.x, x1 = p2.x - p1.x, y0 = p1.y, y1 = p2.y - p1.y, z0 = p1.z, z1 = p2.z - p1.z; float a = (s3 - s0 + s1 + s4 + s6 - s2 - s5 - s7) *x1*y1*z1; float b = (-x0*y1*z1 - x1*y0*z1 - x1*y1*z0 + x1*z1)*s7 + (x0*y1*z1 + x1*y0*z1 + x1*y1*z0)*s6 + (y1*z1 - x0*y1*z1 - x1*y0*z1 - x1*y1*z0)*s5 + (-x1*z1 + x1*y1*z0 - y1*z1 + x0*y1*z1 + x1*y0*z1)*s4 + (-x1*z1 + x1*y0*z1 + x1*y1*z0 - x1*y1 + x0*y1*z1)*s3 + (-x1*y0*z1 - x0*y1*z1 + x1*y1 - x1*y1*z0)*s2 + (x1*y1 + y1*z1 - x1*y1*z0 + x1*z1 - x0*y1*z1 - x1*y0*z1)*s0 + (x1*y1*z0 - y1*z1 - x1*y1 + x0*y1*z1 + x1*y0*z1)*s1; float c = (-x0*y0*z1 + x0*z1 + x1*z0 - x1*y0*z0 - x0*y1*z0)*s7 + (x1*y0*z0 + x0*y1*z0 + x0*y0*z1)*s6 + (y0*z1 - x0*y1*z0 - x1*y0*z0 + y1*z0 - x0*y0*z1)*s5 + (x0*y0*z1 + z1 - y0*z1 - y1*z0 - x0*z1 - x1*z0 + x0*y1*z0 + x1*y0*z0)*s4 + (x1*y0*z0 + x0*y0*z1 + x0*y1*z0 - x1*z0 - x1*y0 - x0*z1 - x0*y1 + x1)*s3 + (x0*y1 + x1*y0 - x0*y0*z1 - x0*y1*z0 - x1*y0*z0)*s2 + (-x1*y0 + x0*y1*z0 - y0*z1 - x0*y1 + x0*y0*z1 + y1 + x1*y0*z0 - y1*z0)*s1 + (-x0*y1*z0 - z1 + x1*y0 - x0*y0*z1 - x1*y0*z0 - y1 + y0*z1 + x1*z0 + y1*z0 + x0*y1 - x1 + x0*z1)*s0; float d = (x0*z0 - x0*y0*z0)*s7 + (y0*z0 - x0*y0*z0)*s5 + (-x0*z0 - y0*z0 + x0*y0*z0 + z0)*s4 + (-x0*z0 + x0 + x0*y0*z0 - x0*y0)*s3 + (x0*y0 - x0*y0*z0)*s2 + (-y0*z0 - x0*y0 + y0 + x0*y0*z0)*s1 + (-y0 - z0 - x0*y0*z0 + x0*z0 + y0*z0 - x0 + x0*y0 + 1.0)*s0 + s6*x0*y0*z0; float r[4]; r[0] = 0.0; r[1] = 0.0; r[2] = 0.0; r[3] = 0.0; int numsteps = 0; // at this point P(t) = a.t^3 + b.t^2 + c.t + d if ( (abs(a)<=0.00001) && (abs(b)<=0.00001) ) { // P(t) is linear numsteps = 0; } else if (abs(a)<=0.00001) { // P(t) is quadratic r[0] = -c/(2.0*b); if ((r[0] <= 0.0) || (r[0] >= 1.0)) { numsteps = 0; } else { numsteps = 1; } } else { // P(t) is cubic // small optimization here : we divide delta by 4, // and simplify r[0]/r[1] by 2 float delta = b*b - 3.0*a*c; if (delta < 0.0) { numsteps = 0; } else { numsteps = 2; r[0] = (-b - sqrt(delta))/(3.0*a); r[1] = (-b + sqrt(delta))/(3.0*a); if ((r[1] <= 0.0) || (r[1] >= 1.0)) { numsteps--; } if ((r[0] <= 0.0) || (r[0] >= 1.0)) { numsteps--; r[0] = r[1]; } } } #if 0 // handle light extrema as well if (abs(e)>0.00001) { // Q(t) is quadratic if ((-f/(2.0*e) > 0.0) && (-f/(2.0*e) < 1.0)) { r[numsteps] = -f/(2.0*e); numsteps++; } } #endif vec4 result, val0, val1, val2, val3; float sample0,sample1,sample2,sample3,sample4; if (numsteps==0) { // single preintegration over [0,1] // evaluate the scalar value at the 2 points : // sample0 at t = 0.0; // sample1 at t = 1.0; sample0 = d; sample1 = d + c + b + a; // preintegrate over [0,1.0] -> [sample0,sample1] val0 = sample(sample0,sample1,length); // blend values result.rgba = val0.rgba; } else if (numsteps==1) { // double preintegration over [0,r[0]] and [r[0],1.0] // evaluate the scalar value at the 3 points : // sample0 at t = 0.0; // sample1 at t = r[0]; // sample2 at t = 1.0; sample0 = d; sample1 = d + r[0]* (c + r[0]* (b + r[0]*a)); sample2 = d + c + b + a; // preintegrate over [0,r[0]] -> [sample0,sample1] val0 = sample(sample0,sample1,r[0]*length); // preintegrate over [r[0],1] -> [sample1,sample2] val1 = sample(sample1,sample2,(1.0 - r[0])*length); // blend values result.rgba = val0.rgba + vec4(1.0 - val0.a) * val1.rgba; } else if (numsteps==2) { // numsteps==2 // triple preintegration over [0,r[0]], [r[0],r[1]] and [r[1],1.0] if (r[1] [sample0,sample1] val0 = sample(sample0,sample1,r[0]*length); // preintegrate over [r[0],r[1]] -> [sample1,sample2] val1 = sample(sample1,sample2,(r[1] - r[0])*length); // preintegrate over [r[1],1] -> [sample2,sample3] val2 = sample(sample2,sample3,(1.0 - r[1])*length); // blend values result.rgba = val0.rgba + vec4(1.0 - val0.a) * (val1.rgba + vec4(1.0 - val1.a) * val2.rgba); } else { // numsteps==3 // triple preintegration over [0,r[0]], [r[0],r[1]], [r[1],r[2]] // and [r[2],1.0] if (r[0]>r[1]) { float tmp = r[0]; r[0] = r[1]; r[1] = tmp; } if (r[1]>r[2]) { float tmp = r[2]; r[2] = r[1]; r[1] = tmp; } if (r[0]>r[1]) { float tmp = r[0]; r[0] = r[1]; r[1] = tmp; } // evaluate the scalar value at the 4 points : // sample0 at t = 0.0; // sample1 at t = r[0]; // sample2 at t = r[1]; // sample3 at t = 1.0; sample0 = d; sample1 = d + r[0]* (c + r[0]* (b + r[0]*a)); sample2 = d + r[1]* (c + r[1]* (b + r[1]*a)); sample3 = d + r[2]* (c + r[2]* (b + r[2]*a)); sample4 = d + c + b + a; // preintegrate over [0,r[0]] -> [sample0,sample1] val0 = sample(sample0,sample1,r[0]*length); // preintegrate over [r[0],r[1]] -> [sample1,sample2] val1 = sample(sample1,sample2,(r[1] - r[0])*length); // preintegrate over [r[1],r[2]] -> [sample2,sample3] val2 = sample(sample2,sample3,(r[2] - r[1])*length); // preintegrate over [r[2],1] -> [sample3,sample4] val3 = sample(sample3,sample4,(1.0 - r[2])*length); // blend values result.rgba = val0.rgba + vec4(1.0 - val0.a) * (val1.rgba + vec4(1.0 - val1.a) * (val2.rgba + vec4(1.0 - val2.a) * val3.rgba)); } gl_FragColor.rgba = result.rgba; }