diff --git a/Python/03_Image_Details.ipynb b/Python/03_Image_Details.ipynb index 602edcaf..271d44ea 100644 --- a/Python/03_Image_Details.ipynb +++ b/Python/03_Image_Details.ipynb @@ -68,6 +68,9 @@ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", + "import random\n", + "import tempfile\n", + "import shutil\n", "\n", "from ipywidgets import interact, fixed\n", "import os\n", @@ -1218,6 +1221,195 @@ "print(image_opt1.GetMetaDataKeys())\n", "print(image_opt2.GetMetaDataKeys())" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Reading non-uniformly spaced images\n", + "\n", + "ITK/SimpleITK images have uniform pixel spacing. When reading an image as a set of slices in DICOM format there are situations where the spacing between slices is non uniform. In such cases the image is read and a warning, \"Non uniform sampling or missing slices detected...\" is issued. The resulting image's spatial extent is incorrect, but the per slice intensity information is correct.\n", + "\n", + "Reading such images into a spatially correct image requires interpolation and sometimes extrapolation. The following code enables reading of such non-uniformly spaced images. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def read_dcm_volume(sorted_file_names, dcmIO=\"GDCMImageIO\"):\n", + " \"\"\"\n", + " Read a volume from a list of sorted DICOM file names, usually obtained via a call to\n", + " sitk.ImageSeriesReader_GetGDCMSeriesFileNames. In ITK/SimpleITK image spacing is required\n", + " to be uniform. When spacing is non-uniform ITK issues a warning 'Non uniform sampling or\n", + " missing slices detected...' but the image is read and the user needs to deal with this.\n", + " This function resolves the problem for non-uniformly spaced z slices, in which case the volume\n", + " is resampled to have a uniform z spacing equal to the minimal z spacing in the original dataset.\n", + " If spacing is uniform, the volume is read and returned as is.\n", + "\n", + " The approach taken here works well when the interpolated regions are not too\n", + "\n", + " The function assumes orthonormal axes and z scan direction.\n", + " \"\"\"\n", + " reader = sitk.ImageFileReader()\n", + " reader.SetImageIO(dcmIO)\n", + "\n", + " # Get the z coordinate for all slices\n", + " original_z = []\n", + " for file_name in sorted_file_names:\n", + " reader.SetFileName(file_name)\n", + " reader.ReadImageInformation()\n", + " original_z.append(reader.GetOrigin()[2])\n", + " original_z_set = set(original_z)\n", + "\n", + " z_spacings = np.array(original_z[1:]) - np.array(original_z[0:-1])\n", + "\n", + " # If a uniformly spaced dataset, return as is.\n", + " # Comparing spacings is safe because the original information from the DICOM\n", + " # Image Position (Patient) tag, 0020|0032, is stored as a string with limited precision.\n", + " if len(set(z_spacings)) == 1:\n", + " return sitk.ReadImage(sorted_file_names, imageIO=dcmIO)\n", + "\n", + " # Minimal z spacing is the new uniform spacing\n", + " new_z_spacing = np.min(z_spacings)\n", + "\n", + " # Get the z coordinates for the new volume by identifying the overlap with the largest\n", + " # number of original slices. The resulting volume may occupy a slightly\n", + " # different spatial region compared to the original. This is due to the use of\n", + " # uniform spacing.\n", + " max_intersection_size = 0\n", + " for current_z in original_z:\n", + " if not max_intersection_size == len(original_z):\n", + " z_min = (\n", + " current_z\n", + " - int((current_z - original_z[0]) / new_z_spacing + 0.5) * new_z_spacing\n", + " )\n", + " z_max = (\n", + " current_z\n", + " + int((original_z[-1] - current_z) / new_z_spacing + 0.5)\n", + " * new_z_spacing\n", + " )\n", + " slices_z = np.linspace(\n", + " z_min, z_max, int((z_max - z_min) / new_z_spacing) + 1\n", + " )\n", + " intersection_with_orig = original_z_set.intersection(set(slices_z))\n", + " current_intersection_size = len(intersection_with_orig)\n", + " if current_intersection_size > max_intersection_size:\n", + " slices_z_best = slices_z\n", + " max_intersection_size = current_intersection_size\n", + "\n", + " # Configure the resample filter using the fixed information\n", + " image_slice = sitk.ReadImage(sorted_file_names[0], imageIO=dcmIO)\n", + " resample_size = image_slice.GetSize()\n", + " resample_origin = image_slice.GetOrigin()\n", + " resample_direction = image_slice.GetDirection()\n", + " resample_spacing = image_slice.GetSpacing()\n", + "\n", + " resample_filter = sitk.ResampleImageFilter()\n", + " resample_filter.SetInterpolator(sitk.sitkLinear)\n", + " resample_filter.SetTransform(sitk.Transform())\n", + " resample_filter.UseNearestNeighborExtrapolatorOn()\n", + " resample_filter.SetSize(resample_size)\n", + " resample_filter.SetOutputSpacing(resample_spacing)\n", + " resample_filter.SetOutputDirection(resample_direction)\n", + "\n", + " # Create the output volume with appropriate image metadata\n", + " # original size, origin and spacing in x and y, and the modified\n", + " # information in the z direction.\n", + " result = sitk.Image(resample_size[0:2] + (len(slices_z_best),), reader.GetPixelID())\n", + " result.SetOrigin(resample_origin[0:2] + (slices_z_best[0],))\n", + " result.SetSpacing(resample_spacing[0:2] + (new_z_spacing,))\n", + " result.SetDirection(resample_direction)\n", + "\n", + " # Read consecutive pairs of images and use this mini-volume to linearly\n", + " # interpolate slices if they are not present in the original dataset,\n", + " # otherwise read original slices.\n", + " i = 0\n", + " curr_vol = None\n", + " for j, curr_z in enumerate(slices_z_best):\n", + " if i < len(original_z) - 1 and curr_z > original_z[i + 1]:\n", + " i += 1\n", + " curr_vol = None # current sub-volume no longer relevant\n", + " try: # if there is an original slice at this z location use it\n", + " z_file_name = sorted_file_names[original_z.index(curr_z)]\n", + " result[..., j] = sitk.ReadImage(z_file_name, imageIO=dcmIO)[..., 0]\n", + " except (\n", + " ValueError\n", + " ): # an exception is raised when the current z is not in the original list\n", + " if (\n", + " curr_vol == None\n", + " ): # resampling is needed, read the new sub-volume only once\n", + " curr_vol = sitk.ReadImage(sorted_file_names[i : i + 2], imageIO=dcmIO)\n", + " resample_origin = resample_origin[0:2] + (curr_z,)\n", + " resample_filter.SetOutputOrigin(resample_origin)\n", + " result[..., j] = resample_filter.Execute(curr_vol)\n", + " return result" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Read a DICOM volume and compare to a non-uniform version of that volume." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Read the first DICOM series found in the data_directory\n", + "data_directory = os.path.dirname(fdata(\"CIRS057A_MR_CT_DICOM/readme.txt\"))\n", + "file_list = sitk.ImageSeriesReader_GetGDCMSeriesFileNames(data_directory)\n", + "uniform_image = sitk.ReadImage(file_list)\n", + "\n", + "# Copy a random subset of the images to a temporary directory, creating a non-uniform\n", + "# volume and read it using the read_dcm function.\n", + "percentage_to_remove = 0.3\n", + "num_files = len(file_list)\n", + "image_indexes_to_remove = random.sample(\n", + " range(num_files), int(num_files * percentage_to_remove + 0.5)\n", + ")\n", + "\n", + "# Uncomment the following to see an example where this approach generates a uniform volume\n", + "# from the non-uniformly spaced one, but the contents are not accurate due to the interpolation\n", + "# based approach (removing 50 slices from the middle of the volume is too much for interpolation\n", + "# to work).\n", + "# Scroll to slice 75 in the original image and click on any point. You will see that the slice\n", + "# obtained via resampling is very different from the original, not surprising given that the\n", + "# interpolated slice is created using slices that are far away from it.\n", + "# image_indexes_to_remove = range(50,100)\n", + "\n", + "with tempfile.TemporaryDirectory() as tmpdirname:\n", + " for i in range(num_files):\n", + " if i not in image_indexes_to_remove:\n", + " shutil.copy(file_list[i], tmpdirname)\n", + " nonuniform_image = read_dcm_volume(\n", + " sitk.ImageSeriesReader_GetGDCMSeriesFileNames(tmpdirname)\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib widget\n", + "import gui\n", + "\n", + "gui.RegistrationPointDataAquisition(\n", + " uniform_image,\n", + " nonuniform_image,\n", + " figure_size=(8, 4),\n", + " known_transformation=sitk.Transform(),\n", + " image1_title=\"original image\",\n", + " image2_title=\"image from non-uniform slices\",\n", + ");" + ] } ], "metadata": { @@ -1237,7 +1429,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.12" + "version": "3.11.10" } }, "nbformat": 4, diff --git a/Python/gui.py b/Python/gui.py index 42ba8a4f..83341761 100644 --- a/Python/gui.py +++ b/Python/gui.py @@ -41,6 +41,8 @@ def __init__( moving_window_level=None, figure_size=(10, 8), known_transformation=None, + image1_title="fixed image", + image2_title="moving image", ): self.fixed_image = fixed_image ( @@ -61,7 +63,8 @@ def __init__( ) # Keep a history of user point localizations, enabling undo of last localization. self.known_transformation = known_transformation # If the transformation is valid (not None) then corresponding points are automatically added. self.text_and_marker_color = "red" - + self.image1_title = image1_title + self.image2_title = image2_title ui = self.create_ui() display(ui) @@ -118,7 +121,7 @@ def create_ui(self): self.fixed_slider = self.moving_slider = None if self.fixed_npa.ndim == 3: self.fixed_slider = widgets.IntSlider( - description="fixed image z slice:", + description=f"{self.image1_title} z slice:", min=0, max=self.fixed_npa.shape[0] - 1, step=1, @@ -128,7 +131,7 @@ def create_ui(self): self.fixed_slider.observe(self.on_slice_slider_value_change, names="value") self.moving_slider = widgets.IntSlider( - description="moving image z slice:", + description=f"{self.image2_title} z slice:", min=0, max=self.moving_npa.shape[0] - 1, step=1, @@ -225,7 +228,7 @@ def update_display(self): color=self.text_and_marker_color, ) self.fixed_axes.set_title( - f"fixed image - localized {len(self.fixed_point_indexes)} points" + f"{self.image1_title} - {len(self.fixed_point_indexes)} points" ) self.fixed_axes.set_axis_off() @@ -264,7 +267,7 @@ def update_display(self): color=self.text_and_marker_color, ) self.moving_axes.set_title( - f"moving image - localized {len(self.moving_point_indexes)} points" + f"{self.image2_title} - {len(self.moving_point_indexes)} points" ) self.moving_axes.set_axis_off() diff --git a/tests/additional_dictionary.txt b/tests/additional_dictionary.txt index 28df9aed..cac924f5 100644 --- a/tests/additional_dictionary.txt +++ b/tests/additional_dictionary.txt @@ -365,6 +365,7 @@ dataframe dataset datasets dciodvfy +dcm debugOn defaultPixelValue deformable