IEEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. ?, NO. ?, JAN 2010 1 GPU-Based Ray-Casting of Spherical Functions Applied to High Angular Resolution Diffusion Imaging M. van Almsick, T. H. J. M. Peeters, V. Prˇ ckovska, A. Vilanova, and B. ter Haar Romeny Abstract— Any sufficiently smooth, positive, real-valued function ψ : S 2 → R + on a sphere S 2 can be expanded by a Laplace expansion into a sum of spherical harmonics. Given the Laplace expansion coefficients, we provide a CPU and GPU-based algorithm that renders the radial graph of ψ in a fast and efficient way by ray-casting the glyph of ψ in the fragment shader of a GPU. The proposed rendering algorithm has proven highly useful in the visualization of high angular resolution diffusion imaging (HARDI) data. Our implementation of the rendering algorithm can display simultaneously thousands of glyphs depicting the local diffusivity of water. The rendering is fast enough to allow for interactive manipulation of large HARDI data sets. Index Terms—I.3.3 [Computer Graphics]: Viewing Algorithms—; I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism—; J.3 [Computer Applications]: Life and Medical Sciences— ✦ 1 I NTRODUCTION T HE rendering algorithm introduced in this article implements a new primitive for geometric modeling and visualization. The new primitive can present any smooth, positive, real-valued function ψ that is defined on an ordinary sphere S 2 . ψ : S 2 → R + . (1) We refer to these functions as spherical functions. One can visualize a spherical function by a deformed sphere, where the extent of the radial in- or extrusion is given by ψ. The parameterization of the spherical function domain S 2 is usually given by spherical angles ϑ ∈ [0,π] and ϕ ∈ [0, 2π). Hence, the shape or so-called glyph of ψ is given by the radial graph defined by r = ψ(ϑ,ϕ). Any sufficiently smooth spherical function ψ can be expanded in a set of orthonormal, spherical basis func- tions. The default choice is the orthonormal set of spheri- cal harmonics Y m l . Such an expansion is called a Laplace expansion. We assume, that the function ψ can be given or can be approximated by a finite Laplace expansion. Our goal is to render a spherical function in such an efficient way, that hundreds or even thousands of these glyphs can be displayed at once in real-time. There are basically two techniques to achieve this in a standard graphics pipeline. For each glyph of a spherical function one can either generate a vertex geometry of a sphere via a n th -order tessellation and then in- or extrude the resulting vertices according to ψ. Or, alternatively, • All authors are with the department of Biomedical Technology, Technische Universiteit Eindhoven, MB 5600 Eindhoven, the Netherlands. E-mail: see http://www.bmia.bmt.tue.nl one can apply a ray-casting algorithm to the spherical function and determine the shading of each pixel in a simple polygon which covers the glyphs projection towards the viewport. We opt for the second approach and outline an effi- cient ray-casting algorithm for spherical functions in the following sections. The purpose of this endeavor has been a visualiza- tion task in high angular resolution diffusion imaging (HARDI), a magnetic resonance imaging (MRI) tech- nique, that provides a unique view of the fiber structure of white brain matter in-vivo. HARDI is an extension of the better known diffusion tensor imaging (DTI), a technique introduced in its current form in 1994 by Basser et al. [1]. In DTI, one determines symmetric, positive-definite diffusion tensors (DT) from a minimum of seven diffusion weighted images. This gives a second order approximation of the probability density functions (PDF) which describes the local diffusivity of the water molecules by a Gaussian distribution. One of the most important application of DTI is the prediction of fiber orientation, which is specified by the principal eigen- vector of the DT. The local fiber orientations determine the trajectories of fiber bundles and consequently the neuronal connections between regions in the brain. Due to the crude assumption that a simple 3D Gaus- sian probability density function can capture the Green’s function of a diffusion process [2], DTI is limited. It can only recover structures with at most one direction per voxel. In areas with more complex intra-voxel hetero- geneity (i.e. fiber crossing, kissing, divergence etc.), the second order DT approximation fails, which is a severe limitation particularly when tracking the path of a fiber bundle.