ITK (Insight Toolkit)

convert_elastix_to_itk_transform.py

from pathlib import Path

import itk
import SimpleITK as sitk

config = {
    'transform' : {
        'affine': 'AffineTransform'
        , 'bspline': 'BSplineTransform'
    }
}

def convert_elastix_itk(elastix_transform_folder_url):
    if Path(elastix_transform_folder_url).exists():
        
        # Step 1 - Get Elastix Parameters
        parameter_object = itk.ParameterObject.New()
        for filepath in Path(elastix_transform_folder_url).glob('*.txt'):
            parameter_object.AddParameterFile(str(filepath))
        
        
        # Step 2 - Convert to ITK Transform
        affine = None
        bspline = None
        for i in range(parameter_object.GetNumberOfParameterMaps()):
            transform_obj = parameter_object.GetParameterMap(i)
            if transform_obj['Transform'][0] == config['transform']['affine']:
                transform_vals = [float(each) for each in transform_obj['TransformParameters']]
                affine = sitk.AffineTransform(3)
                affine.SetTranslation(transform_vals[-3:])
                affine.SetMatrix(transform_vals[:-3])
                affine.SetCenter([float(each) for each in transform_obj['CenterOfRotationPoint']])
                affine_path = Path(elastix_transform_folder_url).joinpath('affine.tfm')
                sitk.WriteTransform(affine, str(affine_path))
                
            elif transform_obj['Transform'][0] == config['transform']['bspline']:
                # Ref: http://insightsoftwareconsortium.github.io/SimpleITK-Notebooks/Python_html/22_Transforms.html#BSpline
                # Ref: https://www.slicer.org/wiki/Documentation/4.3/FAQ/Registration#What_is_the_Meaning_of_.27Fixed_Parameters.27_in_the_transform_file_.28.tfm.29_of_a_BSpline_registration_.3F
                bspline = sitk.BSplineTransform(3,1)
                bspline.SetTransformDomainOrigin([float(each) for each in transform_obj['Origin']]) # from fixed image
                bspline.SetTransformDomainPhysicalDimensions([int(each) for each in transform_obj['Size']]) # from fixed image
                bspline.SetTransformDomainDirection([float(each) for each in transform_obj['Direction']]) # from fixed image
                # bspline.SetTransformDomainMeshSize([int(each)-1 for each in transform_obj['GridSize']])
                
                fixedParams = [int(each) for each in transform_obj['GridSize']]
                fixedParams += [float(each) for each in transform_obj['GridOrigin']]
                fixedParams += [float(each) for each in transform_obj['GridSpacing']]
                fixedParams += [float(each) for each in transform_obj['GridDirection']]
                bspline.SetFixedParameters(fixedParams)
                bspline.SetParameters([float(each) for each in transform_obj['TransformParameters']])
                
                bspline_path = Path(elastix_transform_folder_url).joinpath('bspline.tfm')
                sitk.WriteTransform(bspline, str(bspline_path))
        
        # Step 3 - Save as .tfm
        if affine is not None and bspline is not None:
            composite = sitk.CompositeTransform([affine, bspline])
            composite_path = Path(elastix_transform_folder_url).joinpath('composite.tfm')
            print (composite)
            sitk.WriteTransform(composite, str(composite_path))
            
    else:
        print (' - Error with path')


elastix_transform_folder_url = './' # should contain some elastix TransformParameter{i}.txt files 
convert_elastix_itk(elastix_transform_folder_url)

elastix.py

"""
 - Ref
  - https://github.com/InsightSoftwareConsortium/ITKElastix
"""

import itk

fixed_image_url = 'Accession Nu.0/img.mha'
moving_image_url = 'Accession Nu.3/img.mha'
fixed_image = itk.imread(fixed_image_url) # itk.ImageFileReader.GetTypes()
moving_image = itk.imread(moving_image_url)

# itk.itkImagePython.itkImageSS3 = Signed Short (for 3D data)
fixed_image_new = itk.CastImageFilter[itk.itkImagePython.itkImageSS3, itk.itkImagePython.itkImageD3].New(fixed_image) # itk.CastImageFilter.GetTypes()
moving_image_new = itk.CastImageFilter[itk.itkImagePython.itkImageSS3, itk.itkImagePython.itkImageD3].New(moving_image)

registered_image, params = itk.elastix_registration_method(fixed_image_new, moving_image_new) # itk.ElastixRegistrationMethod.GetTypes()



# reader = itk.ImageFileReader.New(FileName=fixed_image_url)
# reader.Update()
# fixed_image = reader.GetOutput()

elastix_jupyterlab.py

import itk
import time
from itkwidgets import view, compare, checkerboard

def register(fixed_url, moving_url, viz=False):
    
    fixed, moving = read(fixed_url, moving_url)
    
    t0 = time.time()
    registered_moving, transform_parameters = itk.elastix_registration_method(fixed, moving, log_to_console=False)
    print (' - Registration Time taken: ', time.time() - t0,'s')
    print (' - registered_moving: ', registered_moving.shape)
    
    if viz:
        viz_obj = compare(fixed, registered_moving, mode='z', ui_collapsed=True)
    else:
        viz_obj = None
        
    return viz_obj, registered_moving, transform_parameters

def visualize_checkerboard(fixed_url, moving_url):
    fixed, moving = read(fixed_url, moving_url)
    return checkerboard(fixed, moving, mode='z', ui_collapsed=False)

def visualize_compare(fixed_url, moving_url):
    fixed, moving = read(fixed_url, moving_url)
    return compare(fixed, moving, mode='z', ui_collapsed=True)
    
def read(fixed_url, moving_url):
    fixed = itk.imread(fixed_url, itk.F)
    moving = itk.imread(moving_url, itk.F)
    print (' - shape: ',  fixed.shape, moving.shape)
    print (' - spacing: ', fixed.GetSpacing(), moving.GetSpacing())
    print (' - origin: ', fixed.GetOrigin(), moving.GetOrigin())
    return fixed, moving

#################################################
fixed_url = 'data/CT_3D_lung_fixed.mha'
moving_url = 'data/CT_3D_lung_moving.mha'

# visualize_compare(fixed_url, moving_url)
viz_obj, registered_moving, transform_parameters = register(fixed_url, moving_url, viz=True)
viz_obj