Skip to content

Adding advanced image reading code. #488

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Sep 19, 2025
Merged

Adding advanced image reading code. #488

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
194 changes: 193 additions & 1 deletion Python/03_Image_Details.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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": {
Expand All @@ -1237,7 +1429,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.12"
"version": "3.11.10"
}
},
"nbformat": 4,
Expand Down
13 changes: 8 additions & 5 deletions Python/gui.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
(
Expand All @@ -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)

Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand Down Expand Up @@ -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()

Expand Down Expand Up @@ -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()

Expand Down
1 change: 1 addition & 0 deletions tests/additional_dictionary.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -365,6 +365,7 @@ dataframe
dataset
datasets
dciodvfy
dcm
debugOn
defaultPixelValue
deformable
Expand Down
Loading