MINC/Tutorials/NumericPyminc

Numeric Pyminc, or how python made me love fortran, or eat my shorts, Matlab!
The previous tutorials taught some pyminc basics, but left the basic question of why use python and pyminc at all open. If C and perl were good enough for Peter Neelin, why care? The main reason is, in my opinion, that you can use python as a programming language that's a pleasure to use (unlike, gag, Matlab), is freely distributable (Matlab? Heh.), yet can give you speeds close to native C. And there's more than one way to do it. Read on.

Cortical thickness, or an ode to Pierre-Simon, the Marquis de Laplace
One of the ways we have of measuring cortical thickness is with Laplace's equation to create streamlines between the inside and the outside of the cortex. So, for the numerical examples below, I will take a core bit of that problem and solve it in all the different numerical ways accessible to python. In short, given defined boundaries of inside and outside (see labels.mnc in the link above) it creates isopotential lines between them. This is done through a simple iterative six-point averaging style. As a lookahead, the sample labels and the code can all be found here:

http://www.bic.mni.mcgill.ca/users/jason/pyminc-tutorial/

Standard python
The simple solution given the boundaries is the following:

Iterate over every voxel, sum up all neighbours and divide by six. That's it. On the example grid, this takes about 7 seconds on my computer, and there be the rub. Python is an interpreted language, and simply not very efficient when iterating over voxels as described above.

Vectorized python
As any Matlab nerd will tell you, one solution is to vectorize the problem. That can be done in python too, resulting in the following code:

Much faster - about 0.4 seconds. Uglier to read, however, and not every problem is amenable to these kinds of tricks.

Enter weave
So you could rewrite the whole shebang in C or C++. But then you'd have to deal with the pain of coding in either. So why not just take the inner loop and rewrite that? Simple in python - enter a nice extension known as weave. Here a C++ string can be embedded in python. The first time it sees it, it compiles it on the fly, and unless the code changes will from then on reuse the compiled output.

Again takes about 0.35 seconds. Not bad.

do 300, v0=1, nv0 - rocking it old school
Along with that C++ interface in weave, there's a fantastic Fortran interface to python in f2py. Wouldn't you know it, python makes Fortran loveable again. While any Fortran dialect can be used, I stuck to F77 for the masochism of it all. Take the following Fortran code:

Then, on the command line, compile it with 'f2py -c laplace.f -m laplace', and call it like so in python:

And done! 0.315 seconds now. Note that if you look at the full source code I had to add a transpose due to the joys of row versus column major arrays.

Purdy cython
But what if you wanted to write clean python code and not delve in Fortran? Yet another extension to the rescue - cython. It's a dialect of python which lets you add type information in order to get near native speeds. The code is the following:

Then this bit is compiled - notice how the code looks almost identical to the simple and slow python loop except for the addition of types - and called from the main python code as follows:

Takes about 0.32 seconds.

Conclusions
Above are all the reasons to use pyminc - many options for how to write efficient code in a language that's a pleasure to work with. The full code used is below - have fun: