git clone https://gitweb.michael.orlitzky.com/spline3.git
Spline3 is an implementation of the 3D interpolation scheme described in paper, Local quasi-interpolation by cubic C^1 splines on type-6 tetrahedral partitions by Sorokina and Zeilfelder.
It takes volumetric data as input, and allows you “zoom in” on it, producing higher-resolution data as a result. We “fill in the gaps” through interpolation.
The program is written in Haskell and is novel because the main algorithm is purely functional. Nowhere in the main algorithm are any state or global variables manipulated. This has a unique benefit: the program can parallelize itself. Because the algorithm is a pure function, the compiler knows that it's safe to partition the computation across all of the available processors.
In fact, our results show close-to-perfect gains. In other words, running on two processors is essentially twice is fast as running on one. The program will run on all available cores.
Spline3 was featured in the paper, Guiding Parallel Array Fusion with Indexed Types by Ben Lippmeier, Manuel M. T. Chakravarty, Gabriele Keller, and Simon Peyton Jones.
Ben Lippmeier presented a talk based on the paper at Haskell Symposium 2012. A video of the talk is available on YouTube.
The input data is volumetric. Basically, it's just pixels in space. The data used came from the Stanford Volume Data Archive. Still, this data needs to be preprocessed a little before spline3 will accept it. Note: The spline3 package already includes a preprocessed copy of the MRI data in data/mri.bin.
First, we download the tarball from the website:
user $ wget -q https://graphics.stanford.edu/data/voldata/MRbrain.tar.gz
Then, extract it and remove the tarball.
user $ tar -xf MRbrain.tar.gz
user $ rm MRbrain.tar.gz
Now, we're left with 109 data files. We want to concatenate all of them together. Fortunately, they're named sequentially—but not in alphabetical order. We can use a little shell magic to concatenate them in the right order:
user $ rm -f mri.bin
user $ for x in $(seq 1 109); do cat MRbrain.$x >> mri.bin; done;
The result will be a file named mri.bin containing all 109 layers. Other data from the website can be combined similarly.
In all cases, you will need to supply a height, width, and depth to the program so that is knows the dimensions of its data. For the MRI data, this can be found on the website (although the program's defaults already assume you're using the MRI data):
109 slices of 256 x 256 pixels
So, the correct program invocation would be,
user $ spline3 --depth=109 --height=256 --width=256 <input> <output>
The names for the dimensions are somewhat arbitrary. We deviate a little from the traditional x-y-z plane terminology; instead, the data can be thought of as a stack of 256x256 images. When you're looking at one of the images, the depth (third coordinate) would naturally be towards or away from you.
The scale factor (default: 2) specifies how far you want to zoom in on the data. Higher values produce larger images, but take longer to compute. Only integral values (2, 3, 4, etc.) are supported. For example,
user $ spline3 --slice=50 --scale=8 data/mri.bin out.bmp
would produce an image 8x larger than the source slice.
There are two modes that the program can run in. The first is two-dimensional. In two dimensions, the algorithm is not very impressive: similar results can be achieved with Photoshop. However, it's useful for testing the program.
Since the input data is three-dimensional, we choose one “slice” of it to work on in 2D. If you pass the --slice flag to the program, it will cut that slice out of the input data and operate on it in 2D. For example,
user $ ./spline3 --slice=50 data/mri.bin output.bmp
In two dimensions, the program will output a bitmap as a result. This can be viewed in any image viewer.
In 3D, things are a little trickier. The output format is the same as the input format, which means that you'll need to jump through some hoops to view it.
First of all, to interpolate in 3D, just don't pass the --slice argument. An example, zooming in to the default of 2x:
user $ ./spline3 data/mri.bin output.bin
To view the volumetric data, you'll need to use Mayavi. A python script is provided for the MRI data: util/view-mri-data.py. To view the input data, for example, you would run,
user $ util/view-mri-data.py data/mri.bin
To view an output file (created previously),
user $ util/view-mri-data.py output.bin
Beware that this will consume a ton of memory!
First, a two-dimensional result, from slice #50 of the MRI data. We have applied a heat map (i.e. color ramp) to the result, making the visualization more attractive. Each image is linked to its full-resolution (2800x2800) version:
Now, the three-dimensional results. Each of the images below is linked to its full-resolution (1000x1000) version. Higher magnifications would require more computing power.
There are a number of properties specified in the paper that we can test automatically in the code.
For example, Lemma 5.1 states that the resulting spline should reproduce trilinears. In our tests, we construct 500 arbitrary trilinears (at random), apply the method to them, and check that the result is equal to the input.
There are a large number of these properties, all of which you can see in the test output below. The tests are integrated with Cabal, but you can still run them manually (and the output is much prettier):
user $ runghc -isrc test/TestSuite.hs
Cardinal Tests: c-tilde_2100 rotation correct: [OK] FunctionValues Tests: test directions: [OK] Grid Tests: trilinear c0 t0: coefficients: c0030 is correct: [OK] c0003 is correct: [OK] c0021 is correct: [OK] c0012 is correct: [OK] c0120 is correct: [OK] c0102 is correct: [OK] c0111 is correct: [OK] c0210 is correct: [OK] c0201 is correct: [OK] c0300 is correct: [OK] c1020 is correct: [OK] c1002 is correct: [OK] c1011 is correct: [OK] c1110 is correct: [OK] c1101 is correct: [OK] c1200 is correct: [OK] c2010 is correct: [OK] c2001 is correct: [OK] c2100 is correct: [OK] c3000 is correct: [OK] face0 vertices: v0 is correct: [OK] v1 is correct: [OK] v2 is correct: [OK] v3 is correct: [OK] p. 80, Section (2.9) Properties: c0120 identity: [OK, passed 500 tests] c0111 identity: [OK, passed 500 tests] c0201 identity: [OK, passed 500 tests] c0102 identity: [OK, passed 500 tests] c0210 identity: [OK, passed 500 tests] c0300 identity: [OK, passed 500 tests] cube indices within bounds: [OK, passed 500 tests] Misc Tests: flatten (1): [OK] Tetrahedron Tests: tetrahedron1 geometry: volume1: [OK] tetrahedron2 geometry: volume1: [OK] Cube Properties: p. 79, Section (2.6) Properties: c0120 identity1: [OK, passed 500 tests] c0120 identity2: [OK, passed 500 tests] c0120 identity3: [OK, passed 500 tests] c0120 identity4: [OK, passed 500 tests] c0120 identity5: [OK, passed 500 tests] c0120 identity6: [OK, passed 500 tests] c0120 identity7: [OK, passed 500 tests] c0210 identity1: [OK, passed 500 tests] c0300 identity1: [OK, passed 500 tests] c1110 identity: [OK, passed 500 tests] c1200 identity1: [OK, passed 500 tests] c2100 identity1: [OK, passed 500 tests] p. 79, Section (2.7) Properties: c0102 identity1: [OK, passed 500 tests] c0201 identity1: [OK, passed 500 tests] c0300 identity2: [OK, passed 500 tests] c1101 identity: [OK, passed 500 tests] c1200 identity2: [OK, passed 500 tests] c2100 identity2: [OK, passed 500 tests] p. 79, Section (2.8) Properties: c3000 identity: [OK, passed 500 tests] c2010 identity: [OK, passed 500 tests] c2001 identity: [OK, passed 500 tests] c1020 identity: [OK, passed 500 tests] c1002 identity: [OK, passed 500 tests] c1011 identity: [OK, passed 500 tests] Edge Incidence Tests: t0 shares edge with t6: [OK, passed 500 tests] t0 shares edge with t1: [OK, passed 500 tests] t0 shares edge with t3: [OK, passed 500 tests] t1 shares edge with t2: [OK, passed 500 tests] t1 shares edge with t19: [OK, passed 500 tests] t2 shares edge with t3: [OK, passed 500 tests] t2 shares edge with t12: [OK, passed 500 tests] t3 shares edge with t21: [OK, passed 500 tests] t4 shares edge with t5: [OK, passed 500 tests] t4 shares edge with t7: [OK, passed 500 tests] t4 shares edge with t10: [OK, passed 500 tests] t5 shares edge with t6: [OK, passed 500 tests] t5 shares edge with t16: [OK, passed 500 tests] t6 shares edge with t7: [OK, passed 500 tests] t7 shares edge with t20: [OK, passed 500 tests] opposite octant tetrahedra are disjoint (1): [OK, passed 500 tests] opposite octant tetrahedra are disjoint (2): [OK, passed 500 tests] opposite octant tetrahedra are disjoint (3): [OK, passed 500 tests] opposite octant tetrahedra are disjoint (4): [OK, passed 500 tests] opposite octant tetrahedra are disjoint (5): [OK, passed 500 tests] opposite octant tetrahedra are disjoint (6): [OK, passed 500 tests] all volumes positive: [OK, passed 500 tests] all volumes exact: [OK, passed 500 tests] v0 all equal: [OK, passed 500 tests] interior values all identical: [OK, passed 500 tests] c-tilde_2100 rotation correct: [OK, passed 500 tests] c-tilde_2100 correct: [OK, passed 500 tests] Tetrahedron Properties: p. 78, Section (2.4) Properties: c3000 identity: [OK, passed 500 tests] c2100 identity: [OK, passed 500 tests] c1110 identity: [OK, passed 500 tests] b0_v0_always_unity: [OK, passed 500 tests] b0_v1_always_zero: [OK, passed 500 tests] b0_v2_always_zero: [OK, passed 500 tests] b0_v3_always_zero: [OK, passed 500 tests] b1_v1_always_unity: [OK, passed 500 tests] b1_v0_always_zero: [OK, passed 500 tests] b1_v2_always_zero: [OK, passed 500 tests] b1_v3_always_zero: [OK, passed 500 tests] b2_v2_always_unity: [OK, passed 500 tests] b2_v0_always_zero: [OK, passed 500 tests] b2_v1_always_zero: [OK, passed 500 tests] b2_v3_always_zero: [OK, passed 500 tests] b3_v3_always_unity: [OK, passed 500 tests] b3_v0_always_zero: [OK, passed 500 tests] b3_v1_always_zero: [OK, passed 500 tests] b3_v2_always_zero: [OK, passed 500 tests] swapping_vertices_doesnt_affect_coefficients1: [OK, passed 500 tests] swapping_vertices_doesnt_affect_coefficients2: [OK, passed 500 tests] swapping_vertices_doesnt_affect_coefficients3: [OK, passed 500 tests] swapping_vertices_doesnt_affect_coefficients4: [OK, passed 500 tests] Misc Properties: factorial greater: [OK, passed 500 tests] Cardinal Properties: ccwx rotation changes direction: [OK, passed 500 tests] cwx rotation changes direction: [OK, passed 500 tests] ccwy rotation changes direction: [OK, passed 500 tests] cwy rotation changes direction: [OK, passed 500 tests] ccwz rotation changes direction: [OK, passed 500 tests] cwz rotation changes direction: [OK, passed 500 tests] ccwx rotation result unique: [OK, passed 500 tests] cwx rotation result unique: [OK, passed 500 tests] ccwy rotation result unique: [OK, passed 500 tests] cwy rotation result unique: [OK, passed 500 tests] ccwz rotation result unique: [OK, passed 500 tests] cwz rotation result unique: [OK, passed 500 tests] four cwx is identity: [OK, passed 500 tests] four ccwx is identity: [OK, passed 500 tests] four cwy is identity: [OK, passed 500 tests] four ccwy is identity: [OK, passed 500 tests] four cwz is identity: [OK, passed 500 tests] four ccwz is identity: [OK, passed 500 tests] Slow Tests: trilinear reproduced: [OK] trilinear9x9x9 reproduced: [OK] zeros reproduced: [OK] Properties Test Cases Total Passed 100 32 132 Failed 0 0 0 Total 100 32 132
Spline3 is a Haskell program using the Cabal build system. The file spline3.cabal lists all dependencies and it is recommended that you use Cabal to build the project. For convenience, a makefile is provided to build the project (in the Git repo, at least).
Cabal handles the build, so do whatever you normally do to install cabal packages. If you just want to install it for your user,
user $ runghc Setup.hs configure --user
user $ runghc Setup.hs build
user $ runghc Setup.hs install
should do it.
Email them to me at michael@orlitzky.com.
