# How to show solid bodies using volumetric rendering?

Note: This question was asked when Mathematica 8 was latest release. Version 9 has built-in support for volume rendering through Image3D.

There's an example on wolfram.com showing some medical data rendered in 3D. Unfortunately there's no code, only an animation.

Note that the data is shown as a solid three-dimensional body, not only 2D slices of it like in the example on the doc page of Texture.

Given a 3D texture, how can such a visualization of a solid 3D body be made using Mathematica? You will find the volumetric data on the doc page of Texture, under Applications -> Volume Rendering.

I notice there's something included with CUDALink. Unfortunately I can't try it (I don't have the hardware), but it looks quite different from the image I linked to above, so I think there must be an independent way to produce the image above.

[With Sjoerd's hardware (it's pretty cool, movements are totally smooth):]

Link to the MathGroup version of this post.

-
 I dropped the output of CUDAVolumetricRender you were talking about in your question. Please remove if you don't like it. It really is an impressive volumetric renderer. – Sjoerd C. de Vries May 30 '12 at 21:56 This is an interesting question for me as well, since a number of the quaternion-based fractals are often rendered through this technique... – J. M.♦ May 31 '12 at 1:06 There is an OpenCLLink set of commands as well, which will work on ATI graphics cards: reference.wolfram.com/mathematica/OpenCLLink/guide/… – Guillochon May 31 '12 at 2:38 @Sjoerd The CUDAVolumetricRender output is actually mouse-rotatable? – Szabolcs May 31 '12 at 7:25 @Szabolcs It is. Moreover, it's amazingly fluent. You really see the power of CUDA acceleration here. – Sjoerd C. de Vries May 31 '12 at 9:13

# Solution 1: Using 3D Texture with Polygons

The idea is to use Polygon with 3D texture supported by Texture, but it requires a bit of undocumented hack to make it smooth.

The original data set is from Stanford Graphics Group website. The dataset that has been used is CThead, 8-bit tiffs (download). Before proceed, make sure that you have a plenty of memory (~500MB would be enough). Also, if you don't have a good graphics card, turning off 3D antialiasing (through "Preference" menu) will be helpful.

### Step 1

filename = "cthead-8bit.tar.gz"; (* Appropriate path to the downloaded file *)
slices = Import[filename, #] & /@ Import[filename];


The images should look like this:

### Step 2

Let's apply some colors and transparency. We apply color function using Colorize, then add alpha channel based on binarized image using ColorCombine. Don't forget to specify the color space "RGB" (otherwise, the last channel will not be interpreted as opacity). Since it is a lot of slices, we will use ParallelMap. I am reversing the result to match the orientation.

colored =
Reverse[ParallelMap[
With[{img = #},
ColorCombine[{
Colorize[img, ColorFunction -> "SunsetColors"], Binarize[img]
}, "RGB"]] &, slices]];


The result will look like this:

We will combine them and pack them so that it will make a nice volumetric color dataset. Again, notice that option DataReversed is used to reverse the orientation.

data = DeveloperToPackedArray[Map[ImageData[#, DataReversed -> True] &, colored]];


### Step 3

Basically, we put a lot of polygons with textures coming from this volumetric color dataset. But to make it work rather smooth, you need to use undocumented option BaseStyle -> {RenderingOptions -> {"DepthPeelingLayers" -> n}} (n being some reasonable number). I won't go deep into how this option works, I will just say that this sets limit on rendering of translucent polygons.

Here is the main code to display it:

(* All are in 0-1 coordinates. *)
(* top is z-value for x-y cutaway. *)
(* right is x value for y-z cutaway. *)
(* step defines how many polygon slices will be inserted. Smaller -> more *)
With[{top = 0.4, right = 0.4, rpad = -.2, step = 0.01},
Graphics3D[{
EdgeForm[],
Opacity[.4], (* Overall transparency of the textured polygons *)
Texture[data], (* Set volumetric texture *)

(* Bottom part up until the variable top *)
With[{pts = Table[{{0, 0, z}, {1, 0, z}, {1, 1, z}, {0, 1, z}}, {z, 0, top, step}]},
Polygon[pts, VertexTextureCoordinates -> pts]
],

(* Top part from top to 1, but from 0 to right in x direction *)
With[{pts = Table[{{x, 0, top}, {x, 1, top}, {x, 1, 1}, {x, 0, 1}}, {x, 0, right, step}]},
Polygon[pts, VertexTextureCoordinates -> pts]
],

(* Cutaway decoration *)
EdgeForm[Directive[Opacity[.6, RGBColor[1., .8, .1]], Thick]],
FaceForm[],
Polygon[{{right, 0, top}, {1 + rpad, 0, top}, {1 + rpad, 1, top},
{right, 1, top}, {right, 1, 1}, {right, 0, 1}}]
},
PlotRange -> {{0, 1 + rpad}, {0, 1}, {0, 1}},
Lighting -> "Neutral", Background -> Black,
RotationAction -> "Clip", SphericalRegion -> True,
BoxStyle -> Directive[Opacity[.2], White],
ImageSize -> 4 {100, 100},
(* Rendering option for smooth volumetric rendering *)
BaseStyle -> {RenderingOptions -> {"DepthPeelingLayers" -> 100}},
BoxRatios -> {1 + rpad, 1, .9}]
]


Here is the result:

Warning: Be very careful of decreasing step or increasing DepthPeelingLayers and ImageSize. It will crash the Front End...

# Solution 2: Using CUDAVolumetricRender

It is a limited solution for the people who has CUDA capable graphics cards. Essentially, it will call CUDAVolumetricRender function which is a part of CUDALink.

### Step 0

<<CUDALink


The following line will make it sure that you get the latest CUDA paclet.

CUDAResourcesInstall[Update -> True]


It takes a long time to download full paclet, so if you are sure that your CUDA distribution is up to date (CUDAResourcesInstall[] returns 8.0.4.2, as of 5/31/2012), then you can skip this step.

### Step 1

filename = "cthead-8bit.tar.gz"; (* Appropriate path to the downloaded file *)
slices = Import[filename, #] & /@ Import[filename];


### Step 2

Create a volume data suitable for the CUDAVolumetricRender. First, the data should be a packed integer array with depth 3, each element should range from 0 to 255 (0 being background). Since original data is already in 0-255 range (8-bit TIFF), you can use ImageData[..., Automatic] to take original values out. Otherwise, you can use ImageData[..., "Byte"] to enforce it to be in byte range.

data = DeveloperToPackedArray[Map[ImageData[#, Automatic] &, slices]];


Here is a slight problem. Apparently there is a bug. To render $d \times h \times w$ voxel data, you need to partition it again, so that the data would have $h \times w \times d$ dimension. We are not talking of transposing it, instead, literally we just need to change dimension, without touching data order. It is hard to explain and no need to understand, so I will just provide a conversion function here:

Clear[prepareCUDAVolumeData];

prepareCUDAVolumeData::arg =
"The argument should be an integer array of depth 3.";

prepareCUDAVolumeData[array_] /; ArrayQ[array, 3, IntegerQ] :=
Module[{x, y, z}, {x, y, z} = Dimensions[array];
DeveloperToPackedArray[
Partition[#, x] & /@ Partition[Flatten[array], x*z]]];

prepareCUDAVolumeData[___] /; (Message[prepareCUDAVolumeData::arg];
False) := Null;


### Step 3

Now, we are ready. Let's try this:

CUDAVolumetricRender[prepareCUDAVolumeData[data]]


Here is the result (with some parameter tweak):

The problem with this approach is that you don't have much flexibility except what is given by the interface.

Note: The volume rendering support in Mathematica 8 is quite limited (without CUDALink). Technically, it is not really volumetric rendering in traditional sense. A proper volumetric rendering requires a transfer function which converts accumulation of density values into certain color per each pixel. Here, we are mimicking the effect by using false coloring and alpha blending. The CUDALink volume rendering is a real deal, but you don't really need CUDA as long as graphics cards supports shader language (which is mostly true nowadays).

Note 2: You can expect much improved (and proper) volumetric support in the future... ;)

-
Could you comment on how one would do this using CUDALink? – rcollyer May 31 '12 at 3:14
@rcollyer , CUDAVolumetricRender is in the documentation (link), although it is a canned function. You don't really have much freedom except what's offered there through the interface... (also, you sometimes get memory allocation error). – Yu-Sung Chang May 31 '12 at 3:24
I've looked at the documentation, but it only gives examples involving data already in the raw format. I'm more interested in how one would take a set of slices, as above, and use CUDAVolumetricRender. – rcollyer May 31 '12 at 3:27
Added the solution using CUDAVolumetricRender. – Yu-Sung Chang May 31 '12 at 5:30
This is really great ;-) – Vitaliy Kaurov May 31 '12 at 5:32

Here is my extremely slow volume rendering solution. It is based on placing a Cuboid[] for each voxel of the data.

I switched from MapThread to ParallelMap as for some reason MapThread cannot be parallelized.

slices = Import["https://dl.dropbox.com/u/3730003/slicedata.wdx"]
data = DeveloperToPackedArray[Map[ImageData, slices]];

(*downsampling slices*)
resized = ImageResize[#, 64] & /@ slices;
data = DeveloperToPackedArray[Map[ImageData, resized]];

(*3D*)
dims = Dimensions@data;
coordswithdata = Table[{data[[x, y, z]], {x, y, z}}, {x, 1, dims[[1]]}, {y, 1, dims[[2]]}, {z, 1, dims[[3]]}];
cubes = {Hue @@ #1, EdgeForm[], Cuboid@#2} &
output = ParallelMap[cubes @@ # &, coordswithdata, {3}];
Graphics3D@output


When I tried to rotate the view Mathematica crashed so here is an ugly screenshot, but it works (sort of)!

(*2D*)
dims = Dimensions@data[[45]];
coordswithdata = Table[{data[[45, x, y]], {x, y, 1}}, {x, 1, dims[[1]]}, {y, 1, dims[[2]]}];
cubes = {Hue @@ #1, EdgeForm[], Cuboid@#2} &
output = ParallelMap[cubes, coordswithdata, {2}];
Graphics3D@output


""

-
This took 7 minutes of kernel time and 1G of kernel memory, then it died in swapping while transferring the data to the front end. Can you post a screenshot just to satisfy my curiosity? – Szabolcs May 31 '12 at 8:09
@Szabolcs, screenshots uploaded. – s0rce May 31 '12 at 15:02
A couple of things — This still isn't using the actual 3D texture, but simply uses Hue to fake colour it. Secondly this won't be smooth and will eat up memory and resources like crazy even for a rough output. Compared to this, Yu-Sung's was fluid. In fact, I have a very similar deleted answer, and the main issue there too was memory and roughness. The missing ingredient is the undocumented DepthPeelingLayers that Yu-Sung showed. – rm -rf May 31 '12 at 15:25
@R.M. You are right, it is very slow and memory intensive. I left it up there in case anyone wanted to see or improve upon it, no doubt, Yu-Sung's answer is much better. However, this doesn't fake the 3D texture using Hue it renders each voxel using the hue value given in the dataset including the opacity so the low-intensity (z) voxels are translucent. – s0rce May 31 '12 at 16:04