{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## 8. Resampling with ITK" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After image registrations it is often useful to apply the transformation as found by the registration to another image. Maybe you want to apply the transformation to an original (larger) image to gain resolution. Or maybe you need the transformation to apply it to a label image (segmentation). The `itk.resample_image_filter` function will resample an image given a sampling grid, a spatial transformation, and an interpolation function.\n", "\n", "A spatial transformation is defined as a mapping from the fixed image domain to the moving image domain. More information on the precise definition of the transform can be found in the [Elastix Manual](https://elastix.lumc.nl/download/elastix-5.0.1-manual.pdf) or the [ITK Software Guide](https://itk.org/ItkSoftwareGuide.pdf).\n", "\n", "As an alternative, see the [transformix example](ITK_Example08_SimpleTransformix.ipynb)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Elastix" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import itk\n", "from itkwidgets import compare, checkerboard" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Import Images\n", "fixed_image = itk.imread('data/CT_2D_head_fixed.mha', itk.F)\n", "moving_image = itk.imread('data/CT_2D_head_moving.mha', itk.F)\n", "\n", "# Import Default Parameter Maps\n", "parameter_object = itk.ParameterObject.New()\n", "parameter_map_rigid = parameter_object.GetDefaultParameterMap('rigid')\n", "parameter_object.AddParameterMap(parameter_map_rigid)\n", "parameter_map_affine = parameter_object.GetDefaultParameterMap('affine')\n", "parameter_object.AddParameterMap(parameter_map_affine)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Registration with the registration function..." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Call registration function\n", "result_image_elx, result_transform_parameters = itk.elastix_registration_method(\n", " fixed_image, moving_image,\n", " parameter_object=parameter_object)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Convert the Elastix ParameterMap to ITK Transform's." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "euler_transform = itk.Euler2DTransform[itk.D].New()\n", "\n", "pm0 = result_transform_parameters.GetParameterMap(0)\n", "\n", "center = [float(p) for p in pm0['CenterOfRotationPoint']]\n", "fixed_parameters = itk.OptimizerParameters[itk.D](len(center))\n", "for i, p in enumerate(center):\n", " fixed_parameters[i] = p\n", "euler_transform.SetFixedParameters(fixed_parameters)\n", "\n", "elx_parameters = [float(p) for p in pm0['TransformParameters']]\n", "itk_parameters = itk.OptimizerParameters[itk.D](len(elx_parameters))\n", "for i, p in enumerate(elx_parameters):\n", " itk_parameters[i] = p\n", "euler_transform.SetParameters(itk_parameters)\n", "\n", "\n", "dimension = moving_image.GetImageDimension()\n", "affine_transform = itk.AffineTransform[itk.D, dimension].New()\n", "\n", "pm1 = result_transform_parameters.GetParameterMap(1)\n", "\n", "center = [float(p) for p in pm1['CenterOfRotationPoint']]\n", "fixed_parameters = itk.OptimizerParameters[itk.D](len(center))\n", "for i, p in enumerate(center):\n", " fixed_parameters[i] = p\n", "affine_transform.SetFixedParameters(fixed_parameters)\n", "\n", "elx_parameters = [float(p) for p in pm1['TransformParameters']]\n", "itk_parameters = itk.OptimizerParameters[itk.D](len(elx_parameters))\n", "for i, p in enumerate(elx_parameters):\n", " itk_parameters[i] = p\n", "affine_transform.SetParameters(itk_parameters)\n", "\n", "# When creating the composite transform for itk, \n", "# take into account that elastix uses T2(T1(x)) while itk does this the other way around.\n", "# So to get the correct composite transform, add the last transform in elastix first in itk.\n", "composite_transform = itk.CompositeTransform[itk.D, dimension].New()\n", "composite_transform.AddTransform(affine_transform)\n", "composite_transform.AddTransform(euler_transform)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "result_image_itk = itk.resample_image_filter(moving_image,\n", " transform=composite_transform,\n", " use_reference_image=True,\n", " reference_image=fixed_image)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualization\n", "The results of the image transform can be visualized with widgets from the itkwidget library such as the checkerboard and compare widgets." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "3df906dd6b44476690d74150fdfb0429", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(Viewer(annotations=False, interpolation=False, rendered_image=