Recent Advances in Visualization of Volumetric Data

1
Recent Advances in Visualization of Volumetric
Data

• Ken Brodlie
• Jason Wood
• University of Leeds

2
Applications

• Simulations from computational science
• Medical imaging

3
Applications

• Voxel-Man Anatomical Atlas
• University of Hamburg

4
Scalar Data in 3D

• Rectilinear
• Curvilinear
• Unstructured

5
Reference Model

• Haber-McNabb model describes visualization as a
  pipeline

Draw the geometry
Interpolate the samples
Apply a technique
6
Structure of the Talk
7
Data Enrichment Interpolation

• Rectilinear data
• piecewise trilinear
• 7 linear interpolations
• F(x,y,z) a bx cy dz eyz fzx gxy
  hxyz
• piecewise tricubic
• F(x,y,z) 64 terms!
• .. recently used by Cohen-Or et al (2000) for
  smoothing boundary surfaces

f111
f011
f001
f101
f010
f110
f000
f100
..but deceptively complex, eg surface F
constant is conic surface
8
Data Enrichment Interpolation

• Medical imaging cubic interpolation in
  z-direction only could make sense
• Unstructured data
• piecewise linear within each tetrahedron
• F(x,y,z) a bx cy dz
• other approaches - radial basis functions,
  Shepard methods, (see Nielson review papers)

.. surface F constant really is simple
9
Mapping

• Three approaches
• isosurface
• F(x,y,z) k
• section through dependent variable
• slice
• F(x,y,z) such that (x,y,z) e P
• section through independent variable
• volume rendering
• F(x,y,z) mapped to volume of translucent gel

10
Rendering

• Geometry rendered using computer graphics
  techniques
• polygon rendering and ray casting
• Note interaction with mapping step
• isosurface geometry approximation disguised by
  clever rendering ..

.. but care needed or else..
11
Structure of Talk
12
Isosurfacing - Classical Approach

• Marching Cubes
• Lorensen and Cline (1987)
• Rectilinear grid of cells
• Estimate F(x,y,z) by trilinear interpolant within
  each cell
• Calculate intersections of isosurface with cell
  edges
• Approximate with simple triangulation in interior
• 15 canonical configurations
• Marching Tetrahedra
• Unstructured data

13
Isosurfacing - Improving Marching Cubes

• Surface representation
• naïve interior representation can leave holes
• increase robustness and accuracy
• preserve shapes
• Performance
• few cells contribute to surface
• how can we find them quickly?

14
Structure of Talk
15
Surface representation Robustness

• Getting rid of the dreaded holes!
• Full 256 case table allows a robust algorithm

16
Surface representation Topological correctness

• Ambiguous faces
• asymptotic decider of Nielson and Hamann (1992) -
  resolves ambiguity by matching the bilinear face
  behaviour
• decided by saddle point value
• subcases for each of the 6 ambiguous cases of the
  original 15 -making 27 cases in all
• interior triangulation still simple

17
Surface representation Topological correctness

• Ambiguous interior behaviour
• MC case 5
• increase the data values until surfaces touch and
  then a tunnel appears
• Natarajan (1994) resolved ambiguity by body
  saddle value
• body saddle 3D equivalent of 2D saddle, where
  all first derivatives zero

18
Surface representationTopological correctness

• Definitive work by Chernyaev (1995)
• 15 basic canonical configurations
• Ambiguous faces divide 6 into 18 subcases -
  totalling 27 cases
• Ambiguous interiors (tunnel or no tunnel) divide
  6 into 12 cases - totalling 33 cases
• Marching Cubes 33

19
Surface representation Accuracy

• Can we represent the interior surface within a
  cell more accurately?
• One approach is to approximate the trilinear
  surface with a curved surface
• Hamann et al (1997) used triangular
  rational-quadratic Bezier patches

20
Surface representation Accuracy

• Another approach is to refine the polygonal
  approximation
• Start from contouring
• For bilinear interpolant, contour is hyperbolic
  arc (PRQ)
• Usually approximated by single line (PQ)
• Suppose we use two lines - how do we best do
  this?
• Shoulder point (R) is point on hyperbola furthest
  from chord

21
Surface representation Accuracy

• Lopes (1999) carried this over to isosurfacing
• Greater accuracy from
• adding shoulder points in faces
• adding inflection points in interior
• inflection points points on surface which are
  saddle points on planes parallel to cell faces

22
Surface representation Accuracy

• Inflection points mark significant changes in
  topology
• Also characterise the nature of tunnels

23
Surface representation Accuracy

• An application where triangle count is
  unimportant is manufacturing
• Bailey (2000) describes refining triangles by
  comparing value at mid-points of sides and
  centroid with isosurface value

Similar recent work by Cignoni et al (2000) -
comparison based on distance to isosurface
24
Surface representation Rendering

• Gouraud or Phong shading will improve appearance
  of triangular mesh
• Better results by calculating gradient vector of
  function
• rectilinear grid central differences at
  gridpoints ( normal to best fit plane)
• unstructured grid best fit plane to data - with
  distance weighting

25
Surface representation rendering

• At this conference
• Neumann et al describe new method which fits
  linear function at each grid point using distance
  weighting
• F(x,y,z) a bx cy dz
• choosing a, b, c, d to best fit at each grid
  point
• in context of volume rendering, but could be
  applied to isosurface shading

26
Surface representation Direct Surface Rendering

• Jones and Chen (1995) suggested direct ray
  tracing of surface
• High quality images, with shadows
• No intermediate geometry, view dependent...but is
  this a problem?
• Parker et al (1999) follow up with interactive
  ray tracing
• exploits parallelism on 64 processor SGI Reality
  Monster

27
Isosurfaces - Shape-based

• Early approach was to contour 2D slices -and
  stitch them together into 3D surface- but
  problems in handling branches
• Jones and others have suggested expressing 2D
  slices as distance fields (distance of point from
  contour line)
• Stacking the distance field slices into a volume
  and isosurfacing solves branching problems

distance field characterises the shape of the
contour
28
Isosurfaces - Shape-based

• Now substantial interest in shape-based
  interpolation which uses the distance field idea
• Udupa, University of Pennsylvania

29
Isosurfaces - Shape-based Example

• Visualizing Circle of Willis is a challenging
  problem
• Isosurfacing the raw data not successful
• Segmentation can identify the structure, but on
  closer inspection boundary surface not smooth

30
Isosurfaces - Shape-based

• Improved results from
• extract contours
• smooth
• create distance field
• stack distance fields
• interpolate between slices to increase resolution
• isosurface the distance field data

31
Structure of Talk
32
Surface Extraction - Performance

• Performance issues with the Marching Cubes
  algorithm.
• Visits and tests every cell in the data set for
  surface intersections.
• Produces large numbers triangles, some of them
  extremely small.

33
Performance

• Structured data
• Implicit connectivity between cells.
• Data stored in simple 3D array.
• Unstructured data
• Connectivity explicitly described.
• Data not necessarily stored in a consecutive,
  particularly true of adaptive meshes.

34
Structured data - Octree

• Can take advantage of implicit cell ordering when
  constructing search structures.
• Octree - recursive sub-division of geometric
  space.
• Trees are only optimal when data size is 2nx2nx2n

35
Structured data - Branch on need octree

• Branch On Need Octree (BONO)
• Wilhelms and Van Gelder 1992
• reduces tree overhead by minimising the number of
  subdivisions, particularly at the lowest levels.

36
Structured Data - Multi-resolution

• Multi-Resolution Approach
• Reduce the number of cells to visit by combining
  cells with similar values into larger cells.

• Gives a mesh of varying resolution across the
  data set.

• Leaves cracks in surface where areas of different
  resolution join.

37
Multi-resolution

• Tetrahedral Framework approach -
  Zhou,Chen,Kaufman 1997
• Check data is ((2n1)x(2n1)x(2n1)), pad with
  blank data if required.
• Take centre point and form 6 pyramids with the
  faces as the base of each.
• Recursively sub-divide each pyramid into
  tetrahedra until at original cell level

38
Multi-resolution

• Tetrahedral Framework approach - contd.
• Use some method to allow re-combination of low
  level tetrahedra using an error estimate and
  ensuring no hanging nodes.
• Use Marching Tetrahedra to construct iso-surface.
• Many others working in this area such as Cignoni
  1997, Ohlberger 1998, and others.

39
Different approach for unstructured data

• Structured data approaches use a geometric
  approach to simplification that is not easy to
  apply to unstructured grids.
• Many of the approaches for unstructured data
  construct value based search structures.

40
Performance - Extrema Graphs

• Itoh Koyamada 1995.
• Find maxima and minima within data set
• reduce to single cells
• Join maxima and minima points by arcs of
  connected cells
• Direct arc where possible, otherwise use
  polygonal arcs.

41
Extrema Graphs

• Extrema Graphs contd.
• if direct connection fails due to crossing
  boundary of data set pick different point to
  connect to.

If no cells on direct path then must use
polygonal arc, expensive.
42
Extrema Graph

• Extrema Graph contd.
• to create isosurface, search through boundary
  lists and connection lists looking for seed
  points.
• Use a surface propagation algorithm to grow
  surface from these seed points.

43
Performance - Kd-Trees

• Livnat, Shen Johnson 1996
• previous approaches create structures ordered by
  either min or max, or have a structure for each.
• Combines min max lists into single tree data
  structure.
• Kd-Tree is multi-dimensional binary tree.

44
Kd-Trees

• Representation using Span Space
• allows us to geometrically understand range based
  methods.

• plot each range as a single point in span space.

• for a given threshold, T, can easily show which
  cells are active.

45
Kd-Trees

• Constructing the Kd-tree
• First find min max range of each cell - can be
  used for any cell type.

• Partially sort cells by min value to find median
  cell.
• For each sub-tree, partially sort by max value to
  find median.
• Repeat . . .

46
Kd-Trees

• Searching Kd-Tree
• compare threshold (T) with min of root cell,
• If greater, then test cell at this node, then
  test both children.

• if less need only progress down one half of tree.

47
Performance - Lattice Sub-Division

• Shen, Hansen, Livnat Johnson 1996.
• Uses the idea that each cell range can be
  represented as a single point in the span space.

• Subdivide span space into a lattice and place
  each data cell into a lattice element.

48
Lattice Sub-Division

• Searching throws up 5 possible cases
• case 1 - no intersection
• case 2 - definite intersection
• case 3 - test max only
• case 4 - test min only
• case 5 - test min and max

49
Lattice Sub-Division

• In practice, lattice elements of equal size do
  not contain equal numbers of cells, so elements
  are allowed to be of unequal size.
• This sub-division method suitable for use on
  parallel machines, elements are assigned to a
  different processor using a round-robin method.

50
Performance - Interval Trees

• Cignoni, Marino, Montani, Puppo and Scopigno
  1997.
• Applies the optimally efficient Interval Trees
  data structure of Edelsbrunner to searching for
  active cells.
• Find range for each cell and create a set of
  ranges.
• find the middle value of the overall range of the
  data set, dr, and use it as the discriminating
  value for the root node.

51
Interval Trees

• Using dr, at the root of the tree place all the
  intervals that contain dr, into two lists, one,
  AL, ordered by min, the other list, DR, ordered
  by max.

• For each of the 2 children, of the root create
  similar lists with respect to their
  discriminating value.
• Repeat . . .

52
Interval Trees

• To search the interval tree for a value T
• if T lt dr, scan list AL until list value gt T, use
  all these cells, then search left subtree only.
• if T gt dr, scan list DR until list value lt T, use
  these cells, then search right subtree only.
• If T dr, just use cells in AL.
• This method has been described in use with
  Structured data.

53
Alternatives for visualization over the network

• Reduce size of surface to be sent to viewer by
  reducing number of triangles.
• Can be done using mesh simplification methods.
• Other alternatives exist . . .

54
Web Vis

• Engel, Westermann Ertl,
• Reduce triangle count by creating triangle
  strips.
• Use a multi-resolution approach, but rather than
  reducing the data by means of an error based
  model, take the users focus of interest and a
  given radius for hi quality mesh, the rest at
  lower quality.

55
View Dependent Isosurface

• Livnat and Hansen 1998
• Only create isosurface from cells that will be
  visible from users viewpoint, 3 step algorithm.
• Step 1 - find active cells by using Octree
  representation plus front to back traversal.

• Step 2 - coarse software visibility tests to
  further reduce cells used.
• Step 3 - send triangles to hardware.

56
Surface Simplification

• Discretized Marching Cubes - Montani et al 1994
• Uses own lookup table
• doesnt interpolate along edges.

• combines co-planar triangles into larger polygons

57
Slicing

• Geometric intersection of a plane with the cells
  of the data set.
• Simple approach is to test each cell against the
  plane to look for intersections.

• Most cells are not required - hence performance
  improvements can be found by reducing number of
  cell tests.

58
Slicing

• Performance improvements easy for structured data
  when slice is axially aligned.
• More difficult when cell has arbitrary
  orientation.

59
Slicing

• More difficult for unstructured data.
• Possible to use range based methods from
  isosurfacing since problem can be reduced to a
  single range.
• For axially aligned planes, construct tree using
  appropriate x/y/z coordinate.
• For arbitrary orientation, apply suitable
  rotation then construct tree using appropriate
  x/y/z coordinate.

60
Structure of Talk
61
Modelling the Data as Translucent Gel

• Basic concept is to model data as a translucent
  gel
• Classification step
• maps data values to opacity via opacity transfer
  function
• .. and to colour via colour transfer function

62
Classical Approach - Volume Rendering Integral
C(s)light reflected at point s

• Cast rays through image plane into volume, and
  measure light received
• Kajiya and von Hertzen (1984)
• Max (1995)

m(s) density at point s
s
L
I ?L0C(s)m(s) exp -?s0 m(t)dt ds
light
density
attenuation
63
Simplifying the Integral
I ?L0C(s)m(s) exp -?s0 m(t)dt ds

• Approximate using Riemann sums (n number of
  steps)
• Approximate exponential by Taylor series and
  introduce opacity, a, and unit spacing
• Calculate recursively front-to-back as...

I Sni0 C(iDs)m(iDs)Ds Pi-1j0 exp -m
(jDs)Ds
I Sni0 C(i)a(i) Pi-1j0 (1 - a(j))
Cout Cin (1-ain)aiCi aout ain ai(1 - ain)
stop when a 1
Compositing associative but not commutative ie
can group but cannot re-order
64
The Two Approaches

• Image order
• from image to volume
• classical ray casting method of Levoy (1988)
• Recent attention
• integration
• interpolation
• different meshes
• fast traversal

• Object order
• from volume to image
• classical splatting method of Westover (1989)
• Recent attention
• better splatting
• shear warp rendering
• different meshes
• texture mapping
• hardware advances

65
Structure of Talk
66
Integration

• Riemann sums can be replaced by more accurate
  techniques (eg Simpsons rule)
• Recent work by Jung et al (1998) has found
  semi-analytical solution
• C(s), m(s) expressed as 3rd degree polynomials
  via trilinear interpolation
• numerical approx to exponential term

I ?L0C(s)m(s) exp -?s0 m(t)dt ds
67
Maximum Intensity Projection

• When performance rather than accuracy is the
  goal, we can avoid compositing altogether and
  approximate I by maximum intensity along ray
• MIP Maximum Intensity Projection
• Often used in angiography...

68
Maximum Intensity Projection

• Performance is major issue
• lack of shading in image drives need for
  real-time rotation
• fast identification of maximum becomes important

- analytical maximum in each cell along ray -
maximum of samples along ray - skip cells below
maximum
.. but need to achieve high quality too See Mroz
et al (EG2000) for fast, high quality MIP
69
Interpolation

• Sample points occur within cells, not at grid
  points, so we need to interpolate
• Do we
• classify at grid points, then interpolate
• interpolate, then classify
• ?

• Classify - interpolate
• classification done as pre-processing
• smoothing effect can obscure detail
• Interpolate - classify
• classification now within the inner loop of the
  ray cast (sample points are view dependent)
• in return, fine detail can be picked out

See Gasparakis (1999)
70
Classify - Interpolate

• Wittenbrink et al (1998) point out a danger in
  interpolation after classification
• Naïve colour interpolation would assign 3 parts
  yellow, 1 part blue to centre point

but if opacity of bottom left is zero?
Correct approach is to weight according to
opacity, so colour at centre is yellow!
71
Different Meshes

• Fundamental problem is volume rendering
  compositing is non-commutative
• Hence order matters
• Lesser problem for rectilinear grid where cells
  are naturally ordered
• Big complexity problem for curvilinear and
  unstructured grids

72
Different Meshes - Curvilinear

• Fruhauf (1994) algorithm transforms to
  rectilinear grid to ray cast, then back
• Hong and Kaufman (1998) ray cast into curved
  volume
• find first cell, entry and exit point
• accumulate colour and opacity
• find next cell

• Key is to find the sequence of cells efficiently
• cell faces projected to image plane
• bucket sort to get depth ordering
• example of 3D complexity reduced to 2D problem

73
Different Meshes - Unstructured

• Giertsen (1992) introduced idea of sweep plane
• Sweep plane contains all rays for 1 scan line
• Find cells intersecting plane and order (reduced
  to 2D problem)
• Keep set of active cells to exploit coherence
• Silva and Mitchell (1997) lazy sweep ray casting
  algorithm

image plane
intersection with tetrahedron
74
Fast Traversal

• Template-based ray traversal (Yagel et al, 1992)
• Pre-process volume to identify regions of
  significance
• octree decomposition (Parker et al, 1999)
• boundary method

• Boundary method
• Wan et al (1999)
• project boundary cells to image plane
• create nearest and furthest buffers
• only process rays which intersect, and only
  process from nearest to farthest

75
Structure of Talk
76
Better Splatting

• Original splatting does shading then interpolate
  - causing smoothing effect
• Recent work at Ohio has re-ordered to allow
  shading after interpolation, getting better
  detail
• Mueller et al (IEEE Vis99)

77
Shear Warp Rendering

• To get fast traversal, shear volume by
  translating each slice then can resample as
  shown
• Project front-to-back to get intermediate image
• Then warp image
• Note parallelised versions

78
Unstructured Grids

• Cells are projected onto image plane
• For all pixels covered by a cell, compositing
  operation applied
• Ordering of cells is challenging computational
  geometry problem
• Williams, Max and Stein (1998) describe high
  quality algorithm

79
Texture Mapping

• Exploits texture mapping and blending provided in
  OpenGL environments
• Volume is sliced parallel to viewing plane
• texture painted on to rectangle slices
• textured rectangles are composited

80
Texture Mapping

• SGI OpenGL Volumizer is software which exploits
  3D texture mapping

81
Texture Mapping

• With 2D texture hardware, approach is still
  possible
• Generate views parallel to co-ordinate planes
• Choose closest to viewing direction
• Example using VRML for medical volume
  visualization

Hendin, John, Schochet (1998)
82
Hardware Advances

• Holy grail real-time volume rendering
• Main searcher has been Kaufman through Cube
  architectures
• Also major European effort
• VIZARD (Tubingen)
• VolumePro System now commercially available from
  Mitsubishis RealTime Visualization

83
Conclusions

• Volume Visualization
• 1980s basic algorithms
• 1990s enhancements in terms of robustness,
  accuracy and performance
• 2000 hardware solutions
• Isosurface, slice or volume render?
• the winner is the user