Table of Contents

• Introduction
• Issues
  • Long build time for linux
  • Separate full builds for each python version
  • Saving transformation results
  • There's a bug that prevents composition of transformations
  • SimpleElastix and SimpleITK transformations are incompatible
  • No function to read files with registration evolution (metric values, etc)
  • Colour (multi-channel) images not supported
  • TODO: Masks not working?
• Alternatives

Created by gh-md-toc

Introduction

SimpleElastix is a python project by Kasper Marstal, Erasmus Medical Center, Rotterdam.

It basically extends SimpleITK with an interface to the Elastix image registration toolbox.

Here are some notes on the features, bugs, pros and cons of SimpleElastix regarding Python, as that's what we are interested in in our case. SimpleElastix also provides an interface for C++, Java, R, Ruby, Lua, Tcl and C#.

These notes refer to SimpleElastix v1.0.1, commit abf921e. Later versions may have fixed some of the issues.

My interest in SimpleElastix comes from similarly having developed a light Matlab interface for elastix myself, the ElastixToolbox, that we used to reconstruct a mouse heart from a histology stack in Casero et al., "Transformation diffusion reconstruction of three-dimensional histology volumes from two-dimensional image stacks", MIA, 2017.

Now I'm interested in moving away from Matlab to Python, in part because it forces other users to also buy a Matlab licence, and in part because I'm working more with Machine Learning algorithms, and the field seems to be moving that way too.

That said, however, ElastixToolbox is quite complete (it works with multi-channel images, composition of transformations, etc), so I'm writing this document to point out what I find missing in SimpleElastix, with the idea that myself or somebody else can perhaps use it as a starting point to create a fully functional Python interface for Elastix.

Issues

Long build time for linux

SimpleElastix provides a python wheel with Windows libraries. However, for linux the whole project needs to be built, and this takes ~1 h in a modern computer with 4 cores.

We have found convenient writing a small bash script build_SimpleElastix.sh that downloads the most recent code, installs dependencies in Ubuntu, creates a local conda environment and allows separate builds for different python builds, e.g.

./build_SimpleElastix.sh 2.7
./build_SimpleElastix.sh 3.6

Our build enables OpenMP for speed at run time of elastix, and disables building tests, docs and examples, to speed up the build. It also disables shared libraries to avoid trouble with library linking in local environments.

It is then possible to install those SimpleElastix builds in different local environments running e.g.

source activate my_conda_local_environment
cd ~/Downloads/SimpleElastix/build_3.6/SimpleITK-build/Wrapping/Python/Packaging
python setup.py install

Separate full builds for each python version

I have been unsuccessful at figuring out a way to build the bulk of SimpleElastix once, which is slow, and then quickly generate python bindings for separate python versions.

I opened SimpleElastix issue #142 in case somebody knows how to do it. Kasper also seems to be having trouble with it, as he points out (and links to) "This section in the SimpleITKPythonPackage README file suggests that this should not be the case".

Saving transformation results

Once you have computed a registration, you can extract the transform solution with

affineTransform = affineElastixImageFilter.GetTransformParameterMap()[0]

type(affineTransform)
SimpleITK.SimpleITK.ParameterMap

This type is a SwigPyObject, an interface to a C++ object, that cannot be serialised by pickle, so it cannot be directly saved like this

import pickle
    
with open('/tmp/foo.pickle', 'wb') as f:
    pickle.dump(affineTransform, f, pickle.HIGHEST_PROTOCOL)

TypeError: can't pickle SwigPyObject objects

Instead, you save it as a text parameter file like this

sitk.WriteParameterFile(affineTransform, '/tmp/foo.txt')

You can then reload it like

affineTransform = sitk.ReadParameterFile('/tmp/foo.txt')

There's a bug that prevents composition of transformations

Composition of transformations is important, because repeat interpolation causes cumulative errors.

transform1 -> interpolation -> transform2 -> interpolation -> transform3 -> interpolation

Elastix, using the InitialTransformationFileName parameter, allows to concatenate transformations, so that internally it computes

transform1 -> transform2 -> transform3 -> interpolation

If you want to compose two transformation in Transformix, the first transformation is simply ignored.

transformixFilter = sitk.TransformixImageFilter()
transformixFilter.SetTransformParameterMap(rigidTransform)
transformixFilter.AddTransformParameterMap(affineTransform)
transformixFilter.SetMovingImage(imhisto_gray)
transformixFilter.Execute()

This is reported as bugs in issues #140 and #102.

It is possible though to create a registration object with several levels of transformations (e.g. rigid, affine, B-spline).

SimpleElastix and SimpleITK transformations are incompatible

• SimpleElastix transformations: type SimpleITK.SimpleITK.ParameterMap
• SimpleITK transformations: type e.g. SimpleITK.SimpleITK.AffineTransform

These two classes describe transformations in a different way, and should not be mixed.

No function to read files with registration evolution (metric values, etc)

As far as I can tell, there's no function to read the IterationInfo.0.R0.txt, IterationInfo.0.R1.txt, etc. files produced by Elastix with the evolution of the registration, e.g.

1:ItNr  2:Metric        3a:Time 3b:StepSize     4:||Gradient||  Time[ms]
0       -0.350389       0.000000        0.015852        0.312353        42.0
1       -0.341285       0.000000        0.015852        0.209653        19.0
2       -0.337619       0.000000        0.015852        0.401096        2.8
3       -0.342027       0.000000        0.015852        0.293177        16.8
4       -0.330139       0.000000        0.015852        0.314619        7.9
5       -0.338001       0.000000        0.015852        0.786158        7.9
6       -0.335377       0.000000        0.015852        0.225376        6.8
...

Writing a python function to read them is not difficult, though. For example, in this small function we read the file, read the header to figure out which column contains the metric (it changes depending on the Optimizer used), and then read the numbers in that column

def readIterationInfo(iterationInfoFile):
    with open(iterationInfoFile, 'r') as csvfile:
        dialect = csv.Sniffer().sniff(csvfile.read(1024), delimiters='\t')
        csvfile.seek(0)
        fileReader = csv.reader(csvfile, dialect)
        metric = []
        #header = fileReader.next() 
        for row in fileReader:
            if (fileReader.line_num == 1):
                # find which column contains the metric
                for i, val in enumerate(row):
                    if 'Metric' in val:
                        idxMetric = i
                        break
            else:
                metric.append(float(row[idxMetric]))
    csvfile.close()
    return metric

Colour (multi-channel) images not supported

A colour image has multiple channels (the pixel is a Vector type in ITK). Elastix/Transformix only take greyscale images as inputs.

It is actually possible to register multi-channel images with Elastix following two steps:

• selecting (Registration "MultiMetricMultiResolutionRegistration") in the registration parameters file

• passing each channel as a separate greyscale image to the elastix binary, e.g.

  elastix -f0 channelR.tif -f1 channelG.tif -f2 channelB.tif [...]
  

SimpleElastix does not support multi-channel images. It tries to apply to the input image an sitk.Cast to operate with 32-bit floats, but the cast only works with scalar (1 channel) pixels.

sitk.Cast(im, sitk.sitkFloat32)
RuntimeError: Exception thrown in SimpleITK TransformixImageFilter_Execute: .../SimpleElastix/Code/BasicFilters/src/sitkCastImageFilter.cxx:94:
sitk::ERROR: Filter does not support casting from casting vector of 32-bit float to 32-bit float

when for multi-channel images, the cast should be

sitk.Cast(im, sitk.sitkVectorFloat32)

Apart from that, SimpleElastix would need to be internally updated so that when it detects a multi-channel image, it selects the appropriate registration option, and splits/merges the channels as required.

We implemented this in the Matlab ElastixToolbox interface for elastix

% create text string with the command to run elastix
comm = 'elastix ';
for I = 1:size(fixedfile, 1)
    comm = [...
        comm ...
        ' -f' num2str(I-1) ' "' fixedfile(I, :) '"' ...
        ' -m' num2str(I-1) ' "' movingfile(I, :) '"' ...
        ];
end

and similarly in transformix, we transform each channel separatedly

for I = 1:nchannel
    
    % extract channel from full image
    scich.data = scimat.data(:, :, :, :, I);
    
    % smart estimation of background colour for newly added pixels
    if (opts.AutoDefaultPixelValue)
        
        % we assume as background the typical intensity voxel for this
        % channel
        t.DefaultPixelValue = median(scich.data(:));
        
    end
    
    % save channel to temp file
    scimat_save(tmpfilename, scich);
    
    % apply transform to channel
    sciaux = warp_image(t, tmpfilename, ext, opts);
    
    % aggregate channel
    if (I == 1)
        sciout = sciaux;
    else
        sciout.data(:, :, :, :, I) = sciaux.data;
    end
end

In the meantime, a high-level solution in Python is to split colour images, process each channel separately, combine them again, and cast back to the input pixel type (as images are internally cast to float in transformix)

def transformixColor(transform, im):
    transformixFilter = sitk.TransformixImageFilter()
    transformixFilter.SetTransformParameterMap(transform)

    # grayscale image
    if (im.GetNumberOfComponentsPerPixel()==1):
        transformixFilter.SetMovingImage(im)
        transformixFilter.Execute()
        imout = sitk.Cast(transformixFilter.GetResultImage(), im.GetPixelIDValue())
        return imout
    
    # colour image
    imout = []
    for chan in range(im.GetNumberOfComponentsPerPixel()):
        transformixFilter.SetMovingImage(sitk.VectorIndexSelectionCast(im, chan))
        transformixFilter.Execute()
        imout.append(transformixFilter.GetResultImage())
    
    imout = sitk.Compose(imout)
    imout = sitk.Cast(imout, im.GetPixelIDValue())
    return imout

TODO: Masks not working?

I ran a test where I get a good affine registration. Then I add a mask

affineElastixImageFilter.SetMovingMask(mask)

mask is an Image with pixel type uint8, values {0,1} or {0,255} (I've tried both cases). The registration then seems to produce almost the identity, instead of a similar result.

But I haven't tested this in enough depth to know whether this is a problem with SimpleElastix.

Alternatives

PyElastix is a thin Python wrapper around the Elastix binaries by Almar Klein.

This means that the processing is left to the Elastix/Transformix binaries provided by project Elastix, already compiled for Linux, Windows and MacOS X. PyElastix provides a Python interface that hides the writing and reading of parameter and results files and images.

• The PyElastix approach doesn't require building, and it's arguably easier to upgrade (just download the new binaries for PyElastix vs. make a full rebuild with SimpleElastix).
• However, PyElastix is not a full interface. For example, it's not possible to
  • pass a fixed or moving mask to the registration algorithm
  • compose transformations
  • use colour images
  • work with images in real world coordinates (e.g. Images having spacing, offset and direction values)
• Shortcomings in PyElastix would be relatively easy to overcome adding some code to the project.

This thin interface approach is the same we followed in the [ElastixToolbox] (https://github.com/vigente/gerardus/tree/master/matlab/ElastixToolbox) for Matlab. On the other hand, ElastixToolbox requires buying Matlab, which can be expensive.