Reaction-Diffusion Textures

Introduction

Our goal in this section is to use reaction-diffusion textures as a means of visualizing three-dimensional diffusion tensor data. We start by describing a simple model of reaction-diffusion texture that works in two and three dimensions, and then discuss how to modify its calculation to make the texture reflect measured diffusion tensor data. Then, we describe how to render the three-dimensional textures as a stand-alone method for diffusion tensor visualization, as well as how to integrate them into the rendering methods described in previous sections. The use of reaction-diffusion textures for this purpose is closely related to previous work which tuned spot noise to portray local characteristics of scalar and vector fields [36], or which used three-dimensional line integral convolution of spot noise to perform flow visualization of volumetric vector data [14,27,6]. In our case, instead of tuning noise, we are tuning what emerges as a well-organized pattern. The pattern closely corresponds to a field of ellipsoids, the traditional means of diffusion tensor visualization.

The origin of reaction-diffusion textures is a paper by Alan Turing [30] that sought to mathematically model the formation of the immense variety of growth and pigmentation patterns found in the animal kingdom. Turing's paper describes a pair of non-linear partial differential equations modeling the reactions between two chemicals (called ``morphogens''), which diffuse at different rates and interact according to certain rules of activation and inhibition. Reaction-diffusion became a popular method in computer graphics for generating textures with the development of methods for normalizing the density of texture on parameterized surface and polygonal models, and for generating a rich variety of texture patterns [31,38]. Reaction-diffusion textures have been applied to a wide variety of other contexts as well [22].

The reaction-diffusion equations Turing proposed are quite simple. The concentrations of the two morphogens are represented by $a$ and $b$, and the differential equations tell how to increment $a$ and $b$ as function of their reaction and diffusion. The initial condition at $t = 0$ is that $a = b = 4$ everywhere.

$$\frac{\partial a}{\partial t} = s (k (16 - ab) + d_a \nabla^{2} a)$$ (12)

$$\frac{\partial b}{\partial t} = s (k (ab - b - 12 + \beta) + d_b \nabla^{2} b)$$

The scaling factor $k$ controls the size of the reaction terms of the equations relative to the diffusion terms, determining the size of the emergent patterns. Larger values produce patterns that stabilize more quickly, but have smaller characteristic size. The diffusion rates $d_a$ and $d_b$ control how fast the two chemicals can spread in the medium. The overall speed of the pattern's emergence is controlled by $s$: higher values make the system run faster, but values too large can lead the system into divergent instability. The remain ingredient is $\beta$, a pattern of uniformly distributed random values in a small interval centered around 0. It is this pattern that pushes the system away from the unstable equilibrium of the initial conditions, and towards the final texture. In practice, these reaction-diffusion systems are simulated on a regular two-dimensional discrete grid, in which case the Laplacian $\nabla^{2}$ of a chemical $c$ can be measured with a discrete convolution mask $L$:

$$\nabla^{2} c = L * c = \left[ \begin{array}{ccc} 0 & 1.0 & 0 \\ 1.0 & -4.0 & 1.0 \\ 0 & 1.0 & 0 \end{array} \right] * c$$ (13)

Figure 11 shows a simple two-dimensional texture that was generated using the equation above and the indicated parameter settings.

Figure 11: Amount of chemical $a$ in solution to Equation 12 on a $100 \times 100$ grid with indicated parameter settings; black is about $3.0$, white is about $7.5$.

An important property of Equation 12 is that they are general with respect to dimension. Specifically, they work equally well to create a volumetric solid texture [24]. Such a texture can be calculated once and then mapped onto any surface placed within the volume. Or, in the context of volume rendering, the texture volume value can modulate the material color calculated at each sample point. The only implementation change is that the Laplacian is now measured with a three-dimensional convolution mask, shown in Figure 12.

Figure 12: Three-dimensional convolution mask $L$ for measuring Laplacian $\nabla ^2 c$.

Figure 13: Amount of chemical $a$ in solution to Equation 12 on a $100 \times 100 \times 100$ grid with indicated parameter settings, volume rendered with colored lights and depth cueing.

The result of running the three-dimensional reaction-diffusion equation is shown in Figure 13. Note that none of the simulation parameters have changed, and the resulting texture is a close three-dimensional analog to what was seen in Figure 11 (though some of the spots have joined together into more distended blobs).

Tuning the Texture with Tensor Data

Suppose we have a single chemical $c$ in an isotropic medium with diffusivity $d$ (a scalar). The rate of change in $c$ due to non-steady-state diffusion is governed by Fick's second law [21]:

$$\frac{\partial c}{\partial t} = \nabla \cdot (d \nabla c)$$ (14)

$$= d \nabla^2 c$$

Equation 14 says that the amount of chemical $c$ changes according to the divergence of its concentration gradient, and the diffusivity $d$. Since the diffusion is isotropic, the scalar $d$ can be brought outside the divergence, and we get the $\nabla ^2 c$ factor which appears in Equation 12. In an anisotropic medium, however, $d$ is replaced by the diffusion tensor matrix $D$, which transforms how the gradient of concentration $\nabla c$ determines flux. $D$ is the matrix calculated from the measured diffusion-weighted MRI images; it comprises the diffusion tensor field we wish to visualize by volume rendering. Generalizing Equation 14 by replacing $d$ with $D$, we derive a convolution mask $M$ that makes calculating the reaction-diffusion simulation simple.

$$\frac{\partial c}{\partial t} = \nabla \cdot (D \nabla c)$$

$$= \nabla \cdot \left( \left[ \begin{array}{ccc} D_{xx} & D_{xy} & D... ...al c}{\partial y} \\ \frac{\partial c}{\partial z} \end{array} \right] \right)$$

$$= \nabla \cdot \left( \begin{array}{c} D_{xx}\frac{\partial c}{\par... ...artial c}{\partial y} + D_{zz}\frac{\partial c}{\partial z} \end{array} \right)$$

$$= D_{xx}\frac{\partial ^2 c}{\partial x^2} + D_{yy}\frac{\partial ^2 c}{\partial y^2} + D_{zz}\frac{\partial ^2 c}{\partial z^2} +$$ (15)

$$2 D_{xy}\frac{\partial ^2 c}{\partial x \partial y} + 2 D_{xz}\fr... ...c}{\partial x \partial z} + 2 D_{yz}\frac{\partial ^2 c}{\partial y \partial z}$$

$$= M * c$$

Equation 15 should be viewed as the linear combination of six different derivatives of the concentration $c$, each of which can be implemented with a convolution mask defined on the same $3 \times 3 \times 3$ grid seen in Figure 12. Using a combination of first and second central differences to evaluate the derivatives, we arrive at the mask $M$ shown in Figure 14. This mask is different at each location in the field because it is built directly from the components of the diffusion tensor matrix.

Figure 14: Three-dimensional convolution mask $M$ for measuring $\nabla \cdot (D \nabla c)$ to determine amount of diffusion. $tr(D)$ is the trace of $D$, $D_{xx} + D_{yy} + D_{zz}$.

Simply substituting the position-independent $L$ of Figure 12 in Equation 12 with the position-dependent $M$ shown in Figure 14 creates a reaction-diffusion texture that visualizes anisotropic diffusion tensor data. Depending on the magnitudes of values in $D$, one may have to adjust the $s$ and $k$ parameters in Equation 12 in order to produce a simulation that converges on spots of an appropriate size. Instead of having approximately spherical spots, the spots will be stretched into ellipsoids that reflect the diffusion tensor data in their local neighborhood.

This simple substitution, however, is not exactly what we want. The random-walk nature of diffusion dictates that if we drop some ink into an anisotropic medium at a location where the diffusion tensor matrix has eigenvalues $\lambda_1$, $\lambda_2$, and $\lambda_3$, then the shape of the ink spot as it grows over time will approximate an ellipsoid whose axes are proportional to $\sqrt{\lambda_1}$, $\sqrt{\lambda_2}$, and $\sqrt{\lambda_3}$ [21]. That is, the ellipsoid discussed in previous sections- the image of the unit sphere under the tensor matrix- is not the same as the ellipsoid produced by running a diffusion simulation on an initial point of dye in a uniformly anisotropic medium. Because of this, we must run the reaction-diffusion simulation on a different tensor field, one in which the diffusion tensor matrix with eigenvalues $\lambda_1$, $\lambda_2$, $\lambda_3$ is replaced by a matrix with eigenvalues $\lambda_1^2, \lambda_2^2$, $\lambda_3^2$ (the eigenvectors are unchanged). This way the ellipsoids appearing in the results of the reaction-diffusion simulation will have the correct aspect ratio.

Figure 15: Comparison of regular array of ellipses (a) with anisotropic reaction-diffusion textures (b) for both a synthetic dataset (top row) and a portion of measured diffusion tensor data in the human brain (bottom row). Only $2 \times 2$ sub-matrices of the $3 \times 3$ diffusion tensor matrices were used to steer the textures in column (b).

The above discussion also holds for the simpler case of two-dimensional diffusion tensor data. Instead of the whole convolution mask shown in Figure 14, we use its slice at $z = 0$. Figure 15 demonstrates how well two-dimensional reaction diffusion textures can represent diffusion tensor data. There are two major advantages to the use of reaction diffusion textures over a regular grid of ellipses. The first advantage is that the texture spots are packed together according to their size. Unlike with the ellipse arrays, there are no large gaps in which the tensor data is not visualized, nor do the texture spots ever overlap. While it is certainly the case that an algorithm could be developed for intelligently placing ellipses on the field, the benefit of using these textures is that the packing of texture spots arises automatically from the reaction-diffusion simulation. This benefit applies equally well to three-dimensional reaction-diffusion textures.

The second advantage of the reaction-diffusion textures is that the ellipses created in them are placed stochastically in a way that better allows the natural structure of the data to be seen. The reaction-diffusion texture facilitates visually tracking curved, twisting features more easily than on a regular grid. This issue is also addressed by research into how to position streamlines for the most effective flow visualization [32]. Because the ellipses are elongated, if they are placed in a way such that they approximately line up end to end (as happens in the circular synthetic dataset in Figure 15, column (a)), a non-existent linear structure spanning multiple ellipses is perceived, which accentuates the placement scheme of the ellipses, rather than the orientation of the underlying data.

Visualizing the Reaction-Diffusion Texture

Visualizing a three-dimensional texture generated from diffusion tensor data would be most useful if we had a way of removing the texture spots that occur in isotropic regions so that they do not obscure our view of the more interesting features. Fortunately, removing isotropic spots is a relatively easy operation to perform. Because of the uniform boundary and brightness of all the texture spots, it is trivial to choose a threshold for making the texture into a binary image in which all the spots are separated from each other. Next, we perform connected component analysis, defining adjacency by face neighbors. Since the spots are generally convex and thicker than a single voxel, each spot is correctly detected as one connected component. Then, we can determine the average value of a barycentric opacity map inside each spot, since we have defined the texture pattern as overlaying the tensor data. Finally, spots with average opacity below $0.5$ are removed from the texture.

Figure 16: Sample reaction-diffusion for three-dimensional tensor data. The texture is shown before (a) and after (b) segmenting. (c) shows a rendering using the same barycentric opacity map which was used for spot segmentation. The image in (b) shows both the structures of (c) and the structures' anisotropy orientation.

Figure 16 shows a small texture volume before and after segmenting out the isotropic spots, as well as structures that were selected by the barycentric opacity map. Renderings of segmented textures are useful visualizations in their own right, since as a whole, they can show the shape of the anisotropic structures in a flexible way, while the individual spots in the structure represent the local properties of the diffusion tensor field.

Figure 17: Texture-mapping with reaction-diffusion texture.

Figure 17(a) shows a segmented texture for half of the brain dataset shown in previous figures, and Figure 17(b) shows the (unsegmented) texture applied to a volume rendering of the tensor data. For this volume rendering, applying the texture was a simple matter of modulating the calculated material color at each point by the corresponding value in the texture volume; a more sophisticated technique like bump-mapping is not as straight-forward to accomplish. The benefit of texture-mapping the reaction-diffusion pattern onto the surface is that now the direction of anisotropy is indicated on what would be an otherwise isotropic surface.