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 = Developer`ToPackedArray[Map[ImageData, slices]]; (*downsampling slices*) resized = ImageResize[#, 64] & /@ slices; data = Developer`ToPackedArray[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)! ![screenshot][1] (*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 "![Mathematica graphic](https://i.sstatic.net/UpIk1.png)" [1]: https://i.sstatic.net/IPhag.jpg