Simple Volumetric Modeling

Setup

Initialization

<< Experimental`

<< Graphics`Animation`

ops  =  {Mesh -> False, Frame -> False};

opsi =  {Mesh -> False, Frame -> False, 
          DisplayFunction -> Identity};

nb  =  
        NotebookCreate[WindowFrame -> "Palette", 
          Visible -> False];

displayToNotebook[nb_NotebookObject, graphic_, 
          opts___] := 
(NotebookWrite[nb, 
            Cell[GraphicsData["PostScript", 
                DisplayString[graphic]], "Graphics", opts]]; 
          graphic);

Read in the dataset

Mathematica 4 includes new and fast algorithms for import of binary data, even unformatted data  like this 450kB file of grayscale values.

\(\(datanew = 
          BinaryImport["C:\desk\40 Demos\head.dat", {"Byte" .. }];\) // Timing\)

Reconstruct the original structure

Because the data was stored in a one-dimensional file, we have to reconstruct the original structure: 57 scans at a resolution of 100x80 pixels. First we partition into 57 sets of  8000 bytes each.

data2  = Partition[datanew, 80*100];

Dimensions[data2]

Next we, I reconstruct the 3D structure and throw away the intermeidate results.

data3  = Table[ Partition[data2[[i]], 80], {i, 1, 57}];  
    Clear[datanew, data2];

We now have the data cube.

Dimensions[data3]

Viewing data

Slicing

To see if everything worked correcly, let's view a slice.

ListDensityPlot[data3[[24, All, All]], ops];

Of course this works for arbitrary directions:

Show[
        GraphicsArray[{ListDensityPlot[data3[[24, All, All]], opsi], 
            ListDensityPlot[data3[[All, 50, All]], opsi], 
            ListDensityPlot[data3[[All, All, 40]], 
              opsi]}], 
DisplayFunction -> $DisplayFunction];

We can also take arbitrary slices:

data4 = Table[data3[[i, All, i]], {i, 57}];

Dimensions[data4]

ListDensityPlot[data4, ops];

Flythroughs

The simplest way to run the animation is in the notebook itself.

Table[
        Show[Graphics[Raster[data3[[i]]/256]], 
          AspectRatio -> 1], {i, 1, 57, 2}];

Of course we can do any number and orientation of flythroughs.

Table[
        Show[GraphicsArray[{Show[
                Graphics[Raster[data3[[Floor[1.14  i]]]/256]], 
                AspectRatio -> 1, DisplayFunction -> Identity], 
              Show[Graphics[
                  Raster[data3[[All, Floor[2  i], All]]/256]], 
                AspectRatio -> 1, DisplayFunction -> Identity], 
              Show[Graphics[
                  Raster[data3[[All, All, Floor[1.6  i]]]/256]], 
                AspectRatio -> 1, DisplayFunction -> Identity]}], 
          DisplayFunction -> $DisplayFunction], {i, 1, 50, 
          2}];

For presentation purposes, we can also display and animate the slices in external windows.

nb  =  
        NotebookCreate[WindowFrame -> "Palette", 
          Visible -> False, WindowSize -> {400, 350}];

Table[displayToNotebook[nb,Show[Graphics[Raster[data3[[i]]/256]], 
                AspectRatio -> 1, DisplayFunction -> Identity],AnimationDisplayTime -> 0.2],{i, 1, 57, 2}];

SelectionMove[nb, All, Notebook]; SetOptions[nb, Visible -> True]; 
    SelectionAnimate[nb]

Interpolation and visualization

For an easy reconstruction like this, we just calculate a three-dimensional linear interpolation of the dataset.

interpol = 
          ListInterpolation[data3, {{0, 1}, {0, 1}, {0, 1}}, 
            InterpolationOrder -> 1]; // Timing

For example

<< Graphics`ContourPlot3D`

ContourPlot3D[interpol[x, y, z], {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, 
      Contours -> {100. }, PlotPoints -> {3, 5}]

Using Mathematic’s 3D capabilities, this reconstruction can be rotated interactively.

<< RealTime3D`

Show[%%, Mesh -> False]