ICP Transform to SimpleITK images in Python

Support
1

Hello !

I am new to VTK and I am trying to use VTK ICP algorithm in a rigid way to get a transform to use in SimpleITK then. I do that in Python.
Here is how I proceed:

  1. Import nifty file and apply Discrete Marching Cube Algorithm to obtain surfaces to register
  2. Feed source and target surface to vtkIterativeClosestPointTransform() and compute
  3. Apply the transform to surfaces through vtkTransformPolyDataFilter() -> Everything looks great, my surface are aligned in the viewer
  4. Retrieve the icp transform matrix from vtkIterativeClosestPointTransform() and convert it into physical coordinates (at least the translation)
  5. Apply this transform to my itk images through a sitk.Euler3DTransform() -> Images are not aligned at all
    I have also tried with the inverse transform just in case but results are bad too. I have tried to do everything in VTK too but same results.

Here is my code (I skip the part where I read DICOM files in SITK to get fixed_image and moving_image):

def extract_surface_from_nifty(mask, image, option_visu=0):
    source_reader=vtk.vtkNIFTIImageReader()
    source_reader.SetFileName(mask)
    origin=image.GetOrigin()
    spacing=image.GetSpacing()
    source_reader.SetDataSpacing(spacing)
    source_reader.SetDataOrigin(origin)
    
    print("Start Discrete Marching Cube algo")
    source_dmc = vtk.vtkDiscreteMarchingCubes()
    source_dmc.SetInputConnection(source_reader.GetOutputPort())
    source_dmc.SetValue(0, 255)
    source_dmc.Update()
    return source_dmc

def icp_transform(fixed_mask, moving_mask, fixed_image, moving_image):
    target=extract_surface_from_nifty(fixed_mask, fixed_image, option_visu=0)
    source=extract_surface_from_nifty(moving_mask, moving_image, option_visu=0)
    icp = vtk.vtkIterativeClosestPointTransform()
    icp.SetSource(source.GetOutput())
    icp.SetTarget(target.GetOutput())
    icp.GetLandmarkTransform().SetModeToRigidBody()
    icp.DebugOn()
    icp.SetMaximumNumberOfIterations(20)
    icp.StartByMatchingCentroidsOn()
    icp.Modified()
    icp.Update()
    icpTransformFilter = vtk.vtkTransformPolyDataFilter()
    icpTransformFilter.SetInputData(source.GetOutput())
    icpTransformFilter.SetTransform(icp)
    icpTransformFilter.Update()
    transformedSource = icpTransformFilter.GetOutput()
    return target, source, transformedSource, icp, icpTransformFilter

target, source, transformedSource, icp, icpTx = icp_transform(path_fixed_mask, path_moving_mask, fixed_image, moving_image)

translation=[icp.GetMatrix().GetElement(0,3),icp.GetMatrix().GetElement(1,3),icp.GetMatrix().GetElement(2,3)]
w_translation = np.zeros([1,3])
spacing=moving_image.GetSpacing()
origin=moving_image.GetOrigin()
w_translation = np.array(origin)+np.array(spacing)*np.array(translation)
outTx_itk=sitk.Euler3DTransform()
outTx_itk.SetTranslation(w_translation)
moving_resampled = sitk.Resample(moving_image, fixed_image, outTx_itk, sitk.sitkLinear, 0.0, moving_image.GetPixelID())

Do you have an idea why it isn’t working ? I suppose it has something to do with coordinate systems but I run through solutions… I have also tried to use Slicer with the parameters of the transform but it provides same results as SimpleITK… :roll_eyes:

2

You have read image origin and spacing, but missed the image direction matrix. Until recently, VTK only supported axis-aligned images (LPS axes had to be the same direction as IJK). Can you copy-paste here the IJK to LPS transform of your input image to check if your problems are indeed due to having a non-identity direction matrix?

If that’s the problem then you can pre-transform the volume to be axis-aligned and then transform the end result back to the original space. If you want to get started quickly, then you can use 3D Slicer for this: figure out a good registration workflow using the GUI, then automate by scripting it in Python. You can use VTK, SimpleITK, etc. the same way in Slicer’s Python environment, there are just additional functions for bridging gaps like this one.

3

Hi Andras,

Thanks a lot for your answer. I will try to get this transform but for the moment I can’t find how… Is it possible that this transform is something else than a diagonal matrix ?

4

In medical images quite often image axes are rotated or permuted.

GetQFormMarrix and GetQFormMarrix methods of the reader should provide you the direction matrices you need.