Fast 3D matrix re-slicing
Question: I have a 3D matrix of size _(X, Y, Z)_ which is stored in a data structure as
_Z_ matrices, each _X x Y_ in size. I would like to re-slice these matrices to
obtain _X_ slices, each _Y x Z_ in size. In other words, I want to reslice a
3D matrix stored as XY slices in the YZ plane. The use case is to reslice
axial CT images into sagittal images. I am working inside a **browser
environment.**

Here's an example of what I am trying to achieve:

![http://ieeexplore.ieee.org/ieee_pilot/articles/06/ttg2009061407/assets/img/article_1/fig_2/large.jpg\]\(series\)!\[enter
image description here](http://i.stack.imgur.com/prSnW.jpg)

I have implemented the naive (iterative) solution in Python, which takes O(Y *
Z) per slice. I haven't even bothered writing out the corresponding JavaScript
implementation, because this approach is too slow by several orders of
magnitude.

    
    
    import glob
    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.misc import imread
    
    height, width, depth = 512, 512, 100
    
    volume = np.zeros((height, width, depth))
    
    s = 0
    
    for filename in glob.iglob('./*.jpg'):
      volume[:,:,s] =  imread(filename)[...,0]/255.0
      s += 1
    
    reslice = np.zeros((depth, height, width))
    
    for s in xrange(0, width):
      current = np.zeros((depth, height))
      for i in xrange(0, height):
        for j in xrange(0, depth):
          current[j,i] = volume[i,s,j]
      reslice[:,:,s] = current
    

This algorithm seems to be amenable to parallelization. For example, in CUDA,
one could load the 3D data into global memory, create one thread per pixel,
then iterate for every slice in the new direction, and on each iteration ask
the right pixels to fire in order to fill out the current slice. This would be
a trivial kernel to write, and would be approximately O(1) per slice. However,
I don't have access to CUDA in the browser.

Mapping from CUDA to WebCL is relatively straightforward, but WebCL is out of
question given inexistent vendor support ATM. Therefore, I'm thinking WebGL is
the ideal solution.

I'm not too sure how this would be done in the "WebGL" paradigm, but I'm sure
it can be done, and I suspect it is fairly trivial as well. I can't seem to
find where to start, however, and resources on doing general-purpose
computations with OpenGL are extremely scarce. How would I go about using
OpenGL to speed up reslicing of a 3D matrix inside the browser?


Answer: You don't have necessarily to use webGL to be fast enough.

If you use a 3D array, JavaScript might be too slow but by using a **_flat_**
array the time to reslice is in fact similar to the time it takes to create
the array!

Another trick is to use a **_typed array_** to reduce memory usage and improve
performances (Uint8Array).

I created a small class to handle such a flat array and to slice it.

I think the most relevant thing you want is in fact to get a **_view_** ,
either on (x, y) axes or (y, z) axes.

Since Array creation is very slow, you want to build the view on place within
a fixed buffer. And since you want also a sliced view, you have to create a
buffer and method also for the sliced view. It's fast: creating a view for
your 512X512x100 example take **_less than 5 ms!_** (So in fact, the
putImageData you'll have to do afterward will quite take more time! )

Fiddle is here: <http://jsfiddle.net/n38mwh95/1/>

Here's the class handling the data, you'll have to change the constructor so
it accepts the real raw data:

    
    
    function Array3D(xSize, ySize, zSize) {
        this.xSize = xSize;
        this.ySize = ySize;
        this.zSize = zSize;
        var xyMultiplier = xSize * ySize;
    
        this.array = new Uint8Array(xSize * ySize * zSize);
    
        this.view = new Uint8Array(xSize * ySize);
    
        this.slicedView = new Uint8Array(ySize * zSize);
    
        this.valueAt = function (x, y, z) {
            return this.array[x + xSize * (y + z * ySize)];
        };
    
        this.setValueAt = function (x, y, z, val) {
            return this.array[x + xSize * (y + z * ySize)] = val;
        };
    
        this.buildView = function (z) {
            var src = this.array;
            var view = this.view;
            for (var x = 0; x < xSize; x++) {
                for (var y = 0; y < ySize; y++) {
                    view[x + xSize * y] = src[x + xSize * (y + z * ySize)];
                }
            }
            return view;
        };
    
        this.buildSlicedView = function (x) {
            var src = this.array;
            var sView = this.slicedView;
            for (var y = 0; y < ySize; y++) {
                for (var z = 0; z < zSize; z++) {
                    sView[y + ySize * z] = src[x + xSize * (y + z * ySize)];
                }
            }
            return sView;
        };
    }
    

In use:

    
    
    var xSize = 512;
    var ySize = 512;
    var zSize = 100;
    
    var t1, t2;
    
    t1 = performance.now();
    var testArray = new Array3D(xSize, ySize, zSize);
    t2 = performance.now();
    console.log('created in :' + (t2 - t1));
    
    t1 = performance.now();
    var resliced = testArray.buildView(10);
    t2 = performance.now();
    
    console.log('building view in :' + (t2 - t1));
    
    var x = 80;
    t1 = performance.now();
    var resliced = testArray.buildSlicedView(x);
    t2 = performance.now();
    
    console.log('building sliced view in :' + (t2 - t1));
    

Results:

    
    
    created in :33.92199998779688 (index):73
    building view in :2.7559999871300533 (index):79
    building sliced view in :5.726000003051013 
    

At the end of the code I also added some code to render the view.

Don't forget to cache the canvas imageData: create it only once then re-use it
for best performance.

You could easily have a real-time rendering in fact.

