Merge pull request #416 from zivy/resampling

Additional resampling example.

Author: zivy

.github/workflows/main.yml

@@ -72,7 +72,7 @@ jobs:
           sudo apt-get update
           sudo apt-get install enchant-2
         elif [ "$RUNNER_OS" == "macOS" ]; then
-          brew update
+          brew update-reset
           brew install enchant
         fi
     - name: Install dependencies

Python/21_Transforms_and_Resampling.ipynb

@@ -22,6 +22,7 @@
    "source": [
     "import SimpleITK as sitk\n",
     "import numpy as np\n",
+    "import itertools\n",
     "\n",
     "%matplotlib inline\n",
     "import gui\n",
@@ -874,7 +875,9 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
     "def resize_and_scale_uint8(image, new_size, outside_pixel_value=0):\n",
@@ -930,7 +933,132 @@
     "# Select the arbitrary new size\n",
     "new_size = [128, 128]\n",
     "resized_images = [resize_and_scale_uint8(image, new_size, 50) for image in images]\n",
-    "gui.multi_image_display2D(resized_images);"
+    "gui.multi_image_display2D(resized_images, figure_size=[8, 4]);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Resampling after registration\n",
+    "\n",
+    "The result of registration is a transformation which maps points from the fixed image coordinate system to the moving image coordinate system. Most often, we then resample the moving image onto the fixed image grid using this transformation.\n",
+    "\n",
+    "This implicitly assumes that the information we are interested in is contained only in the overlapping regions. In some instances, this is not the case. For example, when registration is used to align multiple images to create a composite image which has a larger physical extent. This is common in microscopy where multiple tiles are combined to create the full field of view (the interested reader is referred to the [ITKMontage module](https://github.com/InsightSoftwareConsortium/ITKMontage), not available in SimpleITK).\n",
+    "\n",
+    "In the cells below we illustrate how to combine the result of registration using the common approach, resampling onto the fixed image grid and onto a grid that contains all image data. The only difference between these two is the definition of the resampling grid. Thus, resampling using the same registration result can yield significantly different images."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create an artificial registration dataset and result.\n",
+    "image = sitk.ReadImage(fdata(\"cthead1.png\"))\n",
+    "fixed_image = image[80:200, 100:250]\n",
+    "moving_image = image[20:150, 60:160]\n",
+    "registration_result = sitk.Euler2DTransform([50, 60], -0.785, [70, 80])\n",
+    "# Modify the moving image origin and direction cosine using the\n",
+    "# \"registration\" transformation and then add noise to the transformation\n",
+    "# mimicking registration error.\n",
+    "moving_image.SetOrigin(registration_result.TransformPoint(moving_image.GetOrigin()))\n",
+    "moving_image.SetDirection(registration_result.GetMatrix())\n",
+    "registration_result.SetTranslation(\n",
+    "    [t + np.random.random() * 5 for t in registration_result.GetTranslation()]\n",
+    ")\n",
+    "gui.multi_image_display2D(\n",
+    "    [fixed_image, moving_image], [\"fixed image\", \"moving image\"], figure_size=[5, 4]\n",
+    ");"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Resample the moving image onto the fixed image grid, notice that parts of the original moving image are not included\n",
+    "# because they are outside the physical region occupied by the fixed image.\n",
+    "resampled_moving = sitk.Resample(moving_image, fixed_image, registration_result)\n",
+    "\n",
+    "# Use green/magenta color scheme, image is gray in overlap regions where the registration is aligned.\n",
+    "gui.multi_image_display2D(\n",
+    "    [\n",
+    "        sitk.Cast(\n",
+    "            sitk.Compose(fixed_image, resampled_moving, fixed_image),\n",
+    "            sitk.sitkVectorUInt8,\n",
+    "        )\n",
+    "    ],\n",
+    "    [\"combined image\"],\n",
+    "    figure_size=[4, 4],\n",
+    ");"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Compute the axis aligned bounding box grid for all images. Treat all images the same way, so\n",
+    "# code is concise and applicable to N images of the same dimension, 2D or 3D.\n",
+    "# The choice of spacing is arbitrary and implicitly defines the size of the grid.\n",
+    "images = [fixed_image, moving_image]\n",
+    "transforms = [sitk.Transform(2, sitk.sitkIdentity), registration_result]\n",
+    "dim = images[0].GetDimension()\n",
+    "\n",
+    "boundary_points = []\n",
+    "for image, transform in zip(images, transforms):\n",
+    "    for boundary_index in list(\n",
+    "        itertools.product(*zip([0] * dim, [sz - 1 for sz in image.GetSize()]))\n",
+    "    ):  # Points from the moving image(s) are mapped to the fixed coordinate system using the inverse\n",
+    "        # of the registration_result.\n",
+    "        boundary_points.append(\n",
+    "            transform.GetInverse().TransformPoint(\n",
+    "                image.TransformIndexToPhysicalPoint(boundary_index)\n",
+    "            )\n",
+    "        )\n",
+    "max_coords = np.max(boundary_points, axis=0)\n",
+    "min_coords = np.min(boundary_points, axis=0)\n",
+    "\n",
+    "new_origin = min_coords\n",
+    "# Arbitrarily use the spacing of the first image and its pixel type,\n",
+    "# change these to suite your needs.\n",
+    "new_spacing = images[0].GetSpacing()\n",
+    "new_pixel_type = images[0].GetPixelID()\n",
+    "new_size = (((max_coords - min_coords) / new_spacing).round().astype(int)).tolist()\n",
+    "new_direction = np.identity(dim).ravel()\n",
+    "\n",
+    "# Resample all images onto the common grid.\n",
+    "resampled_images = []\n",
+    "for image, transform in zip(images, transforms):\n",
+    "    resampled_images.append(\n",
+    "        sitk.Resample(\n",
+    "            image,\n",
+    "            new_size,\n",
+    "            transform,\n",
+    "            sitk.sitkLinear,\n",
+    "            new_origin,\n",
+    "            new_spacing,\n",
+    "            new_direction,\n",
+    "            0.0,\n",
+    "            new_pixel_type,\n",
+    "        )\n",
+    "    )\n",
+    "\n",
+    "# Use green/magenta color scheme, image is gray in overlap regions where the registration is aligned\n",
+    "gui.multi_image_display2D(\n",
+    "    [\n",
+    "        sitk.Cast(\n",
+    "            sitk.Compose(resampled_images[0], resampled_images[1], resampled_images[0]),\n",
+    "            sitk.sitkVectorUInt8,\n",
+    "        )\n",
+    "    ],\n",
+    "    [\"combined image\"],\n",
+    "    figure_size=[4, 4],\n",
+    ");"
    ]
   }
  ],

tests/additional_dictionary.txt

@@ -126,6 +126,7 @@ IPython
 IQR
 ITK
 ITK's
+ITKMontage
 ImageFileReader
 ImageFileReader's
 ImageIO