Opacity Function Generation

Mathematical Boundary Analysis

In order to develop a method for opacity function generation that uses our boundary model and the information stored in the histogram volume, it is helpful to look at the equation we have used to describe the ideal boundary data value as a function of position:

$$v = f(x) = v_{min}+(v_{max}-v_{min})\frac{1+\operatorname{erf}(\frac{x}{\sigma\sqrt{2}})}{2}$$ (6)

Just as we have taken ``position'' to always be along an axis pointing in the gradient direction, we define zero to always be the position of the inflection point in the boundary. $v_{min}$ and $v_{max}$ are the data values of the materials on either side of the boundary. As $\operatorname{erf}()$ ranges from $-1$ to $1$, $v$ ranges from $v_{min}$ to $v_{max}$. The parameter controlling the amount of boundary blurring is $\sigma$. The first and second derivatives of $f$ are as follows:

$$f'(x) = \frac{v_{max}-v_{min}}{\sigma\sqrt{2\pi}}\exp(-\frac{x^2}{2\sigma^2})$$ (7)

$$f''(x) = -\frac{x(v_{max}-v_{min})}{\sigma^3\sqrt{2\pi}}\exp(-\frac{x^2}{2\sigma^2})$$ (8)

Our choice of boundary parameterization means that $f'(x)$ is a Gaussian, with $\sigma$ being the usual standard deviation. Since the Gaussian has inflection points at $\pm\sigma$, this is where $f''(x)$ attains its extrema. The same positions can serve as artificial delimiters for the extent of the boundary -- we define the ``thickness'' of the boundary to be $2\sigma$. Note that the thickness of a boundary can be recovered if the maximum values of $f'$ and $f''$ are known:

$$\frac{f'(0)}{f''(-\sigma)} = \sigma \sqrt{e}$$ (9)

More importantly, once $\sigma$ is known, we can recover the position $x$ knowing only the values of $f'$ and $f''$:

$$\frac{f''(x)}{f'(x)} = -\frac{x}{\sigma^2}$$ (10)

Opacity functions of data value

Before using Eqn. 10 as the basis for opacity function generation, we define some important functions of data value: $g(v)$ is the average first directional derivative of $f$ over all the positions $\mathbf{x}$ at which $f(\mathbf{x})=v$, and $h(v)$ is likewise the average second directional derivative at value $v$. These two functions can be obtained from the histogram volume by slicing it at data value $v$, and finding the centroid of the scatterplot of $f'$ and $f''$ at that value. The $f'$ axis coordinate of the centroid is $g(v)$, and the $f''$ axis coordinate is $h(v)$.

Knowing $g(v)$ and $h(v)$ for all $v$, one can find the ratio of their maxima to recover $\sigma$ with Eqn. 9, assuming that $g$ attains its maxima at $f(0)$, and that $h$ attains its maxima at $f(-\sigma)$. With this information, we define a mapping $p(v)$ from data value to an approximate position along a boundary:

$$p(v) = \frac{-\sigma^2 h(v)}{g(v)}$$ (11)

$$\approx \frac{-\sigma^2 f''(f^{-1}(v))}{f'(f^{-1}(v))} = x$$

Roughly speaking, $p(v)$ tells us on which side of the nearest boundary a data value $v$ tends to fall. For values closer to $v_{min}$, $p(v)$ will be negative; for values closer to $v_{max}$, $p(v)$ will be positive. In practice, we have found it useful to modify Eqn. 11 to account for the fact that the gradient magnitude at the interior of materials is rarely exactly zero. Knowing how it differs from zero is again a matter of experience, but assuming one can find a $g_{thresh}$ which approximates the ambient gradient magnitude, Eqn. 11 is re-formulated, with a slight loss of mathematical accuracy, as

$$p(v) = \frac{-\sigma^2 h(v)}{\max(g(v) - g_{thresh}, 0)}$$ (12)

The user supplies the last piece of information needed: a function $b(x)$ we term the boundary emphasis function, which maps from position along a boundary to opacity. As the intent is to make only boundaries visible in the rendering, $b(x)$ should be non-zero only near zero. For this reason, we have not been especially careful to prevent $p(v)$ from attaining infinite values due to a low $g(v)$; such a data value $v$ should not contribute to the final image. With $b(x)$, the user can directly control whether rendered boundaries will appear thick or thin, sharp or fuzzy, and the proximity of the rendered boundary to the object interior. The final opacity function $\alpha (v)$ is then defined as

$$\alpha(v) = b(p(v))$$ (13)

Instead of exploring the parameter space of all possible opacity functions, the user explores the parameter space of $b(x)$ and lets the information from the histogram volume, embodied in $p(v)$, constrain the search to those opacity functions which display object boundaries. Defining opacity as a function of position within a boundary then becomes a more intuitive task than defining opacity as a function of data value, as there is a more predictable relationship between changes made to the boundary emphasis function and the corresponding change in rendered results. As there is a potential for inaccuracy in the calculation of $\sigma$ from Eqn. 9, the user may need to experiment with different scalings in the domain of $b(x)$. Fig. 11 shows how the choice of boundary emphasis function affects the opacity function and the rendered image, for a synthetically created dataset containing two concentric spheres at distinct data values. It should be stressed that the user does not set the location of the peaks in $\alpha (v)$, since this is determined by the information in $p(v)$, but the user can influence the location, as well as the width, height, and shape of the peaks.   This is the main benefit of this method: if the histogram volume has successfully captured information about the boundaries in the dataset, the user enjoys high-level control over the character of the rendered boundaries without having to worry about the exact specification of $\alpha (v)$. Yet, the $\alpha (v)$ so generated is sensible enough that it could be edited by hand if desired. For example, since this technique will attempt to make all boundaries opaque, a useful supplement to the interface would be a feature which allows supression of the peaks in $\alpha (v)$ for one or more boundaries.

Figure 11: Relationship between $b(x)$, $\alpha (v)$, and the rendered result.

Even though we have made some strong assumptions about the boundary characteristics in the volume dataset, the technique described here typically works well even if the material boundaries are not ``ideal''. Essentially, by taking the ratio of the second and first derivatives, and by having $b(x)$ assign opacity to positions around zero, we are more apt to make opaque those data values associated with both low second derivatives and high first derivatives. Or, even if $p(v)$ is not a perfect indicator of ``position relative to boundary'', the sign change in $f''$ around its zero-crossing affords us some control over whether we want to emphasize regions closer to or further from the object's interior. Fig. 12 (colorplate) on the accompanying colorplate shows a rendering of an MRI dataset which does not have ideal boundaries but for which this technique still works.

Opacity functions of data value and gradient magnitude

So far the opacity functions under consideration have assigned opacity based on data value alone. Higher quality renderings can sometimes be obtained, however, if the opacity is assigned as a function of both data value and gradient magnitude. Defining these two-dimensional opacity functions by hand is especially challenging because there are even more degrees of freedom than in one-dimensional, value-based opacity functions. Fortunately, the ideas presented so far easily generalize to allow semi-automatic generation of two-dimensional opacity functions.

Analogous to the definition of $h(v)$, we define $h(v,g)$ to be average second derivative over all locations where the data value is $v$ and the gradient magnitude is $g$; this is also easily calculated from the histogram volume. We similarly define a new position function, and from that an opacity function:

$$p(v,g) = \frac{-\sigma^2 h(v,g)}{\max(g - g_{thresh}, 0)}$$ (14)

$$\alpha(v,g) = b(p(v,g))$$ (15)

$\sigma$ is calculated as before; the ratio of the extremum of the average first and second derivatives. The benefit of this kind of opacity function is that it can distinguish between boundaries that have overlapping ranges of values. For instance, the nested cylinders volume in Fig. 9 and the engine block volume in Fig. 10 each have one boundary which overlaps the two other boundaries in data value, spanning from the higher of the two material values to the background value. Selectively rendering this single boundary is impossible with a value-based opacity function, but because the boundary has a distinct curve in the plot of data value versus first derivative, it is possible to create an opacity function which selects only the voxels comprising this boundary. As it did in the case of value-based opacity functions, the technique presented here will generate two-dimensional opacity functions which make all detected boundaries opaque; a simple ``lasso'' tool could then be used to select different regions in the two-dimensional opacity function to render one boundary at a time. In Fig. 13 (colorplate), the feet of the female Visible Human CT dataset [17] are rendered with four different two-dimensional opacity functions. Using a modification of an automatically generated opacity function, one rendering shows almost exclusively the registration cord laced around the body prior to scanning.

The lower right rendering in Fig. 13 demonstrates another advantage of two-dimensional opacity functions -- the ability to accurately render the surface of a material which attains a wide range of data values, as is the case for the bone tissue in this same CT scan. Different parts of the bone surface are more radio-opaque than others, leading to a wide range of data values associated with bone, which in turn causes a wide range of gradient magnitudes within the boundary region between bone and soft tissue. Knowing the average second derivative for each location in $(v,g)$ space, we can make opaque only those voxels near the middle of the boundary (near the zero-crossing in $f''$), regardless of the bone data value. As is visible in the opacity function generated with Eqn. 14, this implies that as gradient magnitude increases, there is an upward shift in the data values which should be made most opaque. This kind of careful opacity assignment is not possible with a simple value-based opacity function, though it is reminiscent of the two-dimensional opacity functions described by Levoy [11]. Although space does not permit a detailed comparison between our approach and Levoy's, the main difference is that (ideally) the measured first and second derivative information serves to constrain the opacity function generation so as to only show boundaries, while in Levoy's method the user still has to experiment to find the right parameter settings.