Direct Volume Rendering (DVR): Ray-casting

1. Direct Volume Rendering (DVR): Ray-casting Jian Huang This set of slides references slides used by Prof. Torsten Moeller (Simon Fraser), Prof. Han-Wei Shen (Ohio State).

2. Direct: No conversion to surface geometry Papers • Tuy and Tuy, 1984, IEEE CG & A (one of the earliest volume rendering techniques) • Levoy, 1988 IEEE CG&A, and later improvements • Drebin, Carpenter, Hanrahan, 1988, SIGGRAPH

3. Basic Idea Based on the idea of ray tracing • Trace from eat each pixel • as a ray into object space • Compute color value along • the ray • Assign the value to the • pixel

4. Data Representation • 3D volume data are represented by a finite number of cross sectional slices (hence a 3D raster) • On each volume element (voxel), stores a data value (if it uses only a single bit, then it is a binary data set. Normally, we see a gray value of 8 to 16 bits on each voxel.) N x 2D arraies = 3D array

5. A voxel is a cubic cell, which has a single value cover the entire cubic region A voxel is a data point at a corner of the cubic cell The value of a point inside the cell is determined by interpolation Data Representation (2) What is a Voxel? – Two definitions

6. Viewing • Ray Casting • Where to position the volume and image plane • What is a ‘ray’ • How to march a ray

7. Viewing (1) 1. Position the volume Assuming the volume dimensions is w x w x w We position the center of the volume at the world origin Volume center = [w/2,w/2,w/2] (local space) Translate T(-w/2,-w/2,-w/2) (0,0,0) x (data to world matrix? world to data matrix ) y z

8. (0,0,0) x y z Viewing (2) 2. Position the image plane Assuming the distance between the image plane and the volume center is D, and initially the center of the image plane is (0,0,-D) Image plane

9. Viewing (3) 3. Rotate the image plane A new position of the image plane can be defined in terms of three rotation angle a,b,g with respect to x,y,z axes Assuming the original view vector is [0,0,1], then the new view vector g becomes: cosb 0 -sinb 1 0 0 cosg sing 0 g = [0,0,1]0 1 0 0 cosa sina -sing cosg 0 sinb 0 cosb 0 -sina cosa 0 0 1

10. E0 u0 E v v0 u + S0 S (0,0,0) x y z Viewing (4) B = [0,0,0] S0 = [0,0,-D] u0 = [1,0,0] v0 = [0,1,0] B Now, R: the rotation matrix S = B – D x g U = [1,0,0] x R V = [0,1,0] x R

11. Viewing (5) Image Plane: L x L pixels Then E = S – L/2 x u – L/2 x v So Each pixel (i,j) has coordinates P = E + i x u + j x v S + v u E R: the rotation matrix S = B – D x g U = [1,0,0] x R V = [0,1,0] x R We enumerate the pixels by changing i and j (0..L-1)

12. d p x x x x Q Viewing (6) 4. Cast rays Remember for each pixel on the image plane P = E + i x u + j x v and the view vector g = [0,0,1] x R So the ray has the equation: Q = P + k (d x g) d: the sampling distance at each step K = 0,1,2,…

13. Early Methods • Before 1988 • Did not consider transparency • did not consider sophisticated light transportation theory • were concerned with quick solutions • hence more or less applied to binary data non-binary data - require sophisticated classification/compositing methods!

14. Ray Tracing -> Ray Casting • “another” typical method from traditional graphics • Typically we only deal with primary rays -hence: ray-casting • a natural image-order technique • as opposed to surface graphics - how do we calculate the ray/surface intersection??? • Since we have no surfaces - we need to carefully step through the volume

15. Ray Casting • Stepping through the volume: a ray is cast into the volume, sampling the volume at certain intervals • The sampling intervals are usually equi-distant, but don’t have to be (e.g. importance sampling) • At each sampling location, a sample is interpolated / reconstructed from the grid voxels • popular filters are: nearest neighbor (box), trilinear (tent), Gaussian, cubic spline • Along the ray - what are we looking for?

16. Example: Using the nearest neighbor kernel Q = P + K x V (v=dxg) At each step k, Q is rounded off to the nearest voxel (like the DDA algorithm) Check if the voxel is on the boundary or not (compare against a threshold) If yes, perform shading In tuys’ paper

17. Basic Idea of Ray-casting Pipeline • Data are defined at the corners • of each cell (voxel) • The data value inside the • voxel is determined using • interpolation (e.g. tri-linear) • Composite colors and opacities • along the ray path • Can use other ray-traversal schemes as well c1 c2 c3

18. Ray Traversal Schemes Intensity Max Average Accumulate First Depth

19. Ray Traversal - First Intensity • First: extracts iso-surfaces (again!)done by Tuy&Tuy ’84 First Depth

20. Ray Traversal - Average Intensity • Average: produces basically an X-ray picture Average Depth

21. Ray Traversal - MIP Intensity Max • Max: Maximum Intensity Projectionused for Magnetic Resonance Angiogram Depth

22. Ray Traversal - Accumulate Intensity • Accumulate opacity while compositing colors: make transparent layers visible!Levoy ‘88 Accumulate Depth

23. 1.0 Raycasting volumetric compositing color opacity object (color, opacity)

24. Raycasting Interpolationkernel volumetric compositing color opacity 1.0 object (color, opacity)

25. 1.0 Raycasting Interpolationkernel volumetric compositing color c = c s s(1 - ) + c opacity =  s (1 - ) +  object (color, opacity)

26. Raycasting volumetric compositing color opacity 1.0 object (color, opacity)

27. Raycasting volumetric compositing color opacity 1.0 object (color, opacity)

28. Raycasting volumetric compositing color opacity 1.0 object (color, opacity)

29. Raycasting volumetric compositing color opacity 1.0 object (color, opacity)

30. Raycasting volumetric compositing color opacity object (color, opacity)

31. Volume RenderingPipeline Acquired values Data preparation Prepared values shading classification Voxel colors Voxel opacities Ray-tracing / resampling Ray-tracing / resampling Sample colors Sample opacities compositing Image Pixels

32. DCH DVR Pipeline

33. DCH - Pipeline Original data Material percentage volumes Classification Density volume Opacity volume Color volume Normals Shading Gradient Shaded volume shears Transformed volume Final image compositing

34. Common Components of General Pipeline • Interpolation/reconstruction • Classification or transfer function • Gradient/normal estimation for shading • Question: are normals also interpolated?

35. Levoy - Interpolation

36. Levoy - Interpolation (2) Closest value Weighted average

37. Levoy – Gradient/Normals • Central difference • per voxel Y+1 y-1, z-1 X+1

38. Levoy - Compositing • Image order • back-to-front • using the over operator

39. Levoy - Shading • Phong Shading + Depth Cueing • Cp = color of parallel light source • ka / kd / ks = ambient / diffuse / specular light coefficient • k1, k2 = fall-off constants • d(x) = distance to picture plane • L = normalized vector to light • H = normalized vector for maximum highlight • N(xi) = surface normal at voxel xi

40. Classification Classification: Mapping from data to opacities Region of interest: high opaicity (more opaque) Rest: translucent or transparent The opacity function is typically specified by the user Levoy came up with two formula to compute opacity 1. Isosurface 2. Region boundary (e.g. between bone and fresh)

41. Classification/Transfer Function • Maps raw voxel value into presentable entities: color, intensity, opacity, etc. Raw-data  material (R, G, B, a, Ka, Kd, Ks, ...) • May require probabilistic methods (Drebin). Derive material volume from input. Estimate % of each material in all voxels. Pre-computed. AKA segmentation. • Often use look-up tables (LUT) to store the transfer function that are discovered

42. Levoy - Classification • Usually not only interested in a particular iso-surface but also in regions of “change” • Feature extraction - High value of opacity exists in regions of change • Transfer function (Levoy) - Saliency • Surface “strength”

43. Levoy - Classification • Chemistry Data • only iso-value - loose information of layers • iso-range - could be too narrow or too wide • thickness of region should be constant • hence linear fall off of opacity • wider fall off for larger gradient • Medical Data • assume at most two tissues meet • linear transition between opacities of “neighboring” tissues • reflects linear combination of tissues within one voxel

44. Opacity function (1) • Goal: visualize voxels that have a selected threshold • value fv • - No intermediate geometry is extracted • - The idea is to assign voxels that have value fv the • maximum opacity (say a) • And then create a smooth transition for the surrounding • area from 1 to 0 • Levoy wants to maintain a constant thickness for the • transition area.

45. Opacity function (2) Maintain a constant isosurface thickness Can we assign opacity based on function value instead of distance? (local operation: we don’t know where the isosurface is) Yes – we can based on the value distance f – fv but we need to take into account the local gradient opacity = a opacity = 0

46. Opacity function (3) Assign opacity based on value difference (f-fv) and local gradient gradient: the value fall-off rate grad = Df/Ds Assuming a region has a constant gradient and the isosurface transition has a thickness R F = f(x) Then we interpolate the opacity opacity = a – a * ( fv-f(x))/ (grad * R) opacity = a F = fv opacity = 0 F = fv – grad * R thickness = R

47. Levoy - Classification A

48. Levoy - Classification B

49. DCH - Material Percentage V. • Probabilistic classifier • probability that a voxel has intensity I: • pi - percentage of material • Pi(I) - prob. that material i has value I • Pi(I) given through statistics/physics • pi then given by:

50. DCH - Classification • Like Levoy - assumes only two materials per voxel • that will lead to material percentage volumes • from them we conclude color/opacity: • where Ci=(aiRi, aiGi, aiBi, ai)