A few months ago, while clearing up my office, I came across a CD entitled ‘MyBrain’. A colleague of mine was studying for his PhD and had needed volunteers to undergo an MRI scan as part of his research. In return for being one of his guinea pigs I kept a copy of the structural scan element of the data he gathered, and promptly did…er… absolutely nothing with it!

So, having rediscovered the data *years* later, I thought I should finally have a crack at reading & rendering using C# & XNA. I had a number of options to tackle this problem. Amongst the most visually impressing algorithms I have seen, are these demos on Ziggyware. However, I wanted to do more than just produce the best rendering, I wanted to produce a true mesh from the data - for reasons I will get to at the end…

Data Structure

The data came in 2 independent files; a *.img file and associated *.hdr file. The data structures of these files turned out to be pretty straightforward - once I’d found a cool free app, MRIcro, which allowed me to browse the image & header files and provided a lot of useful information! Essentially it was quickly clear I could completely ignore the header (*.hdr) file for the purposes of rendering my particular scan image, since all I needed to know was that the data comprised 64 layers of 256x256 images. The images were each of 16-bit integers representing the density of the tissue at each point in space.

Marching Cubes

The problem I have with creating a watertight, closed, mesh (well formed with properly coincident triangle edges and no gaps) from data of this sort is that there isn’t really a surface over which to create your mesh. One solution to this is the Marching Cubes algorithm. This approach can be effective at creating a closed mesh over at the boundary of a given density value in your data. I started with the algorithm & C++ code as described by Paul Bourke in this extremely good paper. In this algorithm, a grid is created in the 3D space occupied by the density data. Each corner of this grid is then checked against the data, to see if it is greater, or less, than the target density we are fitting the mesh to. Once we’ve done this for each corner of a cube, we can then see if we need to create triangles were we have located a boundary between corners that are inside, and corners that are outside, of our target density.

A *particularly* useful bit of the code provided by Bourke are the look-up tables for all these permutations of triangles – some examples of which are shown above.

In the C++ implementation, the code makes extensive use of pointers and I have to admit that I found this to be the first time I *really* missed them in C#. However, I followed the design pretty closely using lists referenced by index instead. I also wanted to create an indexed triangle mesh for rendering as this saves around (or maybe exactly?) 75% of the required vertices and massively improved performance. I noted from Bourke’s description that there is only one vertex created per cube edge, so by keeping track of the shared edges for each neighbouring cube I was able to reuse vertices fairly easily.

STL File Output

The stereolithography, or STL, file format is an *extremely* simple triangle mesh file format specifically for use in Rapid Prototyping machines. Precise details of the format can be found here, but essentially it has an 80 byte header we can ignore, a 32-bit int for the number of triangles, and then the list of un-indexed triangles comprising 3 vertices (each ending in a 16-bit int, which for all practical purposes can apparently be set to 0). It is probably simplest to look at the code snippet:

public void SaveSTL(string filename, VertexPositionNormalColor[] vertices, int[] indices)
{
    System.IO.FileStream fs = System.IO.File.Create(filename);
    System.IO.BinaryWriter bw = new System.IO.BinaryWriter(fs);

    Byte[] header = new Byte[80];
    bw.Write(header);
    UInt32 facet_count = (UInt32)(indices.Length / 3);
    bw.Write(facet_count);
    UInt16 attribute_byte_count = 0;

    for (int i = 0; i < indices.Length; i += 3)
    {
        bw.Write(0.0f);
        bw.Write(0.0f);
        bw.Write(0.0f);
        bw.Write(vertices[indices[i + 2]].Position.X);
        bw.Write(vertices[indices[i + 2]].Position.Y);
        bw.Write(vertices[indices[i + 2]].Position.Z);
        bw.Write(vertices[indices[i + 1]].Position.X);
        bw.Write(vertices[indices[i + 1]].Position.Y);
        bw.Write(vertices[indices[i + 1]].Position.Z);
        bw.Write(vertices[indices[i]].Position.X);
        bw.Write(vertices[indices[i]].Position.Y);
        bw.Write(vertices[indices[i]].Position.Z);
        bw.Write(attribute_byte_count); 
    }
    bw.Close();
    fs.Close();
}

Ironically (yes I really do it mean that, Tim!), you can see that, having gone to the trouble of producing an indexed triangle mesh for rendering, I now had to undo all the good work to produce a no-nonsense verbose triangle list for the file. :(

Making Use of the Result

So why go to all the bother of creating a mesh at all, let alone an STL file? Well courtesy of this RP service, last week I was able to have a model of my brain built in resin! The physical models are built from the same data used in the XNA rendered in the movie above! This picture show the finished models, which are just over an inch in length.

Potentially you could do this for any mesh model, although obviously they tend to geared for low-poly count in games so I’m not sure how good they’d come out. Nevertheless, I thought this might be interesting for someone…

OK, better get back to some more conventional games programming, now my techno ego trip is complete!!