Skip to content

Commit cb10928

Browse files
zivyzivy
authored
Merge pull request #488 from zivy/readDICOMNonuniformSpacing
Adding advanced image reading code.
2 parents 54750e0 + 8accc1d commit cb10928

File tree

3 files changed

+202
-6
lines changed

3 files changed

+202
-6
lines changed

‎Python/03_Image_Details.ipynb‎

Lines changed: 193 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -68,6 +68,9 @@
6868
"%matplotlib inline\n",
6969
"import matplotlib.pyplot as plt\n",
7070
"import numpy as np\n",
71+
"import random\n",
72+
"import tempfile\n",
73+
"import shutil\n",
7174
"\n",
7275
"from ipywidgets import interact, fixed\n",
7376
"import os\n",
@@ -1218,6 +1221,195 @@
12181221
"print(image_opt1.GetMetaDataKeys())\n",
12191222
"print(image_opt2.GetMetaDataKeys())"
12201223
]
1224+
},
1225+
{
1226+
"cell_type": "markdown",
1227+
"metadata": {},
1228+
"source": [
1229+
"#### Reading non-uniformly spaced images\n",
1230+
"\n",
1231+
"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",
1232+
"\n",
1233+
"Reading such images into a spatially correct image requires interpolation and sometimes extrapolation. The following code enables reading of such non-uniformly spaced images. "
1234+
]
1235+
},
1236+
{
1237+
"cell_type": "code",
1238+
"execution_count": null,
1239+
"metadata": {},
1240+
"outputs": [],
1241+
"source": [
1242+
"def read_dcm_volume(sorted_file_names, dcmIO=\"GDCMImageIO\"):\n",
1243+
" \"\"\"\n",
1244+
" Read a volume from a list of sorted DICOM file names, usually obtained via a call to\n",
1245+
" sitk.ImageSeriesReader_GetGDCMSeriesFileNames. In ITK/SimpleITK image spacing is required\n",
1246+
" to be uniform. When spacing is non-uniform ITK issues a warning 'Non uniform sampling or\n",
1247+
" missing slices detected...' but the image is read and the user needs to deal with this.\n",
1248+
" This function resolves the problem for non-uniformly spaced z slices, in which case the volume\n",
1249+
" is resampled to have a uniform z spacing equal to the minimal z spacing in the original dataset.\n",
1250+
" If spacing is uniform, the volume is read and returned as is.\n",
1251+
"\n",
1252+
" The approach taken here works well when the interpolated regions are not too\n",
1253+
"\n",
1254+
" The function assumes orthonormal axes and z scan direction.\n",
1255+
" \"\"\"\n",
1256+
" reader = sitk.ImageFileReader()\n",
1257+
" reader.SetImageIO(dcmIO)\n",
1258+
"\n",
1259+
" # Get the z coordinate for all slices\n",
1260+
" original_z = []\n",
1261+
" for file_name in sorted_file_names:\n",
1262+
" reader.SetFileName(file_name)\n",
1263+
" reader.ReadImageInformation()\n",
1264+
" original_z.append(reader.GetOrigin()[2])\n",
1265+
" original_z_set = set(original_z)\n",
1266+
"\n",
1267+
" z_spacings = np.array(original_z[1:]) - np.array(original_z[0:-1])\n",
1268+
"\n",
1269+
" # If a uniformly spaced dataset, return as is.\n",
1270+
" # Comparing spacings is safe because the original information from the DICOM\n",
1271+
" # Image Position (Patient) tag, 0020|0032, is stored as a string with limited precision.\n",
1272+
" if len(set(z_spacings)) == 1:\n",
1273+
" return sitk.ReadImage(sorted_file_names, imageIO=dcmIO)\n",
1274+
"\n",
1275+
" # Minimal z spacing is the new uniform spacing\n",
1276+
" new_z_spacing = np.min(z_spacings)\n",
1277+
"\n",
1278+
" # Get the z coordinates for the new volume by identifying the overlap with the largest\n",
1279+
" # number of original slices. The resulting volume may occupy a slightly\n",
1280+
" # different spatial region compared to the original. This is due to the use of\n",
1281+
" # uniform spacing.\n",
1282+
" max_intersection_size = 0\n",
1283+
" for current_z in original_z:\n",
1284+
" if not max_intersection_size == len(original_z):\n",
1285+
" z_min = (\n",
1286+
" current_z\n",
1287+
" - int((current_z - original_z[0]) / new_z_spacing + 0.5) * new_z_spacing\n",
1288+
" )\n",
1289+
" z_max = (\n",
1290+
" current_z\n",
1291+
" + int((original_z[-1] - current_z) / new_z_spacing + 0.5)\n",
1292+
" * new_z_spacing\n",
1293+
" )\n",
1294+
" slices_z = np.linspace(\n",
1295+
" z_min, z_max, int((z_max - z_min) / new_z_spacing) + 1\n",
1296+
" )\n",
1297+
" intersection_with_orig = original_z_set.intersection(set(slices_z))\n",
1298+
" current_intersection_size = len(intersection_with_orig)\n",
1299+
" if current_intersection_size > max_intersection_size:\n",
1300+
" slices_z_best = slices_z\n",
1301+
" max_intersection_size = current_intersection_size\n",
1302+
"\n",
1303+
" # Configure the resample filter using the fixed information\n",
1304+
" image_slice = sitk.ReadImage(sorted_file_names[0], imageIO=dcmIO)\n",
1305+
" resample_size = image_slice.GetSize()\n",
1306+
" resample_origin = image_slice.GetOrigin()\n",
1307+
" resample_direction = image_slice.GetDirection()\n",
1308+
" resample_spacing = image_slice.GetSpacing()\n",
1309+
"\n",
1310+
" resample_filter = sitk.ResampleImageFilter()\n",
1311+
" resample_filter.SetInterpolator(sitk.sitkLinear)\n",
1312+
" resample_filter.SetTransform(sitk.Transform())\n",
1313+
" resample_filter.UseNearestNeighborExtrapolatorOn()\n",
1314+
" resample_filter.SetSize(resample_size)\n",
1315+
" resample_filter.SetOutputSpacing(resample_spacing)\n",
1316+
" resample_filter.SetOutputDirection(resample_direction)\n",
1317+
"\n",
1318+
" # Create the output volume with appropriate image metadata\n",
1319+
" # original size, origin and spacing in x and y, and the modified\n",
1320+
" # information in the z direction.\n",
1321+
" result = sitk.Image(resample_size[0:2] + (len(slices_z_best),), reader.GetPixelID())\n",
1322+
" result.SetOrigin(resample_origin[0:2] + (slices_z_best[0],))\n",
1323+
" result.SetSpacing(resample_spacing[0:2] + (new_z_spacing,))\n",
1324+
" result.SetDirection(resample_direction)\n",
1325+
"\n",
1326+
" # Read consecutive pairs of images and use this mini-volume to linearly\n",
1327+
" # interpolate slices if they are not present in the original dataset,\n",
1328+
" # otherwise read original slices.\n",
1329+
" i = 0\n",
1330+
" curr_vol = None\n",
1331+
" for j, curr_z in enumerate(slices_z_best):\n",
1332+
" if i < len(original_z) - 1 and curr_z > original_z[i + 1]:\n",
1333+
" i += 1\n",
1334+
" curr_vol = None # current sub-volume no longer relevant\n",
1335+
" try: # if there is an original slice at this z location use it\n",
1336+
" z_file_name = sorted_file_names[original_z.index(curr_z)]\n",
1337+
" result[..., j] = sitk.ReadImage(z_file_name, imageIO=dcmIO)[..., 0]\n",
1338+
" except (\n",
1339+
" ValueError\n",
1340+
" ): # an exception is raised when the current z is not in the original list\n",
1341+
" if (\n",
1342+
" curr_vol == None\n",
1343+
" ): # resampling is needed, read the new sub-volume only once\n",
1344+
" curr_vol = sitk.ReadImage(sorted_file_names[i : i + 2], imageIO=dcmIO)\n",
1345+
" resample_origin = resample_origin[0:2] + (curr_z,)\n",
1346+
" resample_filter.SetOutputOrigin(resample_origin)\n",
1347+
" result[..., j] = resample_filter.Execute(curr_vol)\n",
1348+
" return result"
1349+
]
1350+
},
1351+
{
1352+
"cell_type": "markdown",
1353+
"metadata": {},
1354+
"source": [
1355+
"Read a DICOM volume and compare to a non-uniform version of that volume."
1356+
]
1357+
},
1358+
{
1359+
"cell_type": "code",
1360+
"execution_count": null,
1361+
"metadata": {},
1362+
"outputs": [],
1363+
"source": [
1364+
"# Read the first DICOM series found in the data_directory\n",
1365+
"data_directory = os.path.dirname(fdata(\"CIRS057A_MR_CT_DICOM/readme.txt\"))\n",
1366+
"file_list = sitk.ImageSeriesReader_GetGDCMSeriesFileNames(data_directory)\n",
1367+
"uniform_image = sitk.ReadImage(file_list)\n",
1368+
"\n",
1369+
"# Copy a random subset of the images to a temporary directory, creating a non-uniform\n",
1370+
"# volume and read it using the read_dcm function.\n",
1371+
"percentage_to_remove = 0.3\n",
1372+
"num_files = len(file_list)\n",
1373+
"image_indexes_to_remove = random.sample(\n",
1374+
" range(num_files), int(num_files * percentage_to_remove + 0.5)\n",
1375+
")\n",
1376+
"\n",
1377+
"# Uncomment the following to see an example where this approach generates a uniform volume\n",
1378+
"# from the non-uniformly spaced one, but the contents are not accurate due to the interpolation\n",
1379+
"# based approach (removing 50 slices from the middle of the volume is too much for interpolation\n",
1380+
"# to work).\n",
1381+
"# Scroll to slice 75 in the original image and click on any point. You will see that the slice\n",
1382+
"# obtained via resampling is very different from the original, not surprising given that the\n",
1383+
"# interpolated slice is created using slices that are far away from it.\n",
1384+
"# image_indexes_to_remove = range(50,100)\n",
1385+
"\n",
1386+
"with tempfile.TemporaryDirectory() as tmpdirname:\n",
1387+
" for i in range(num_files):\n",
1388+
" if i not in image_indexes_to_remove:\n",
1389+
" shutil.copy(file_list[i], tmpdirname)\n",
1390+
" nonuniform_image = read_dcm_volume(\n",
1391+
" sitk.ImageSeriesReader_GetGDCMSeriesFileNames(tmpdirname)\n",
1392+
" )"
1393+
]
1394+
},
1395+
{
1396+
"cell_type": "code",
1397+
"execution_count": null,
1398+
"metadata": {},
1399+
"outputs": [],
1400+
"source": [
1401+
"%matplotlib widget\n",
1402+
"import gui\n",
1403+
"\n",
1404+
"gui.RegistrationPointDataAquisition(\n",
1405+
" uniform_image,\n",
1406+
" nonuniform_image,\n",
1407+
" figure_size=(8, 4),\n",
1408+
" known_transformation=sitk.Transform(),\n",
1409+
" image1_title=\"original image\",\n",
1410+
" image2_title=\"image from non-uniform slices\",\n",
1411+
");"
1412+
]
12211413
}
12221414
],
12231415
"metadata": {
@@ -1237,7 +1429,7 @@
12371429
"name": "python",
12381430
"nbconvert_exporter": "python",
12391431
"pygments_lexer": "ipython3",
1240-
"version": "3.8.12"
1432+
"version": "3.11.10"
12411433
}
12421434
},
12431435
"nbformat": 4,

‎Python/gui.py‎

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,8 @@ def __init__(
4141
moving_window_level=None,
4242
figure_size=(10, 8),
4343
known_transformation=None,
44+
image1_title="fixed image",
45+
image2_title="moving image",
4446
):
4547
self.fixed_image = fixed_image
4648
(
@@ -61,7 +63,8 @@ def __init__(
6163
) # Keep a history of user point localizations, enabling undo of last localization.
6264
self.known_transformation = known_transformation # If the transformation is valid (not None) then corresponding points are automatically added.
6365
self.text_and_marker_color = "red"
64-
66+
self.image1_title = image1_title
67+
self.image2_title = image2_title
6568
ui = self.create_ui()
6669
display(ui)
6770

@@ -118,7 +121,7 @@ def create_ui(self):
118121
self.fixed_slider = self.moving_slider = None
119122
if self.fixed_npa.ndim == 3:
120123
self.fixed_slider = widgets.IntSlider(
121-
description="fixed image z slice:",
124+
description=f"{self.image1_title} z slice:",
122125
min=0,
123126
max=self.fixed_npa.shape[0] - 1,
124127
step=1,
@@ -128,7 +131,7 @@ def create_ui(self):
128131
self.fixed_slider.observe(self.on_slice_slider_value_change, names="value")
129132

130133
self.moving_slider = widgets.IntSlider(
131-
description="moving image z slice:",
134+
description=f"{self.image2_title} z slice:",
132135
min=0,
133136
max=self.moving_npa.shape[0] - 1,
134137
step=1,
@@ -225,7 +228,7 @@ def update_display(self):
225228
color=self.text_and_marker_color,
226229
)
227230
self.fixed_axes.set_title(
228-
f"fixed image - localized {len(self.fixed_point_indexes)} points"
231+
f"{self.image1_title} - {len(self.fixed_point_indexes)} points"
229232
)
230233
self.fixed_axes.set_axis_off()
231234

@@ -264,7 +267,7 @@ def update_display(self):
264267
color=self.text_and_marker_color,
265268
)
266269
self.moving_axes.set_title(
267-
f"moving image - localized {len(self.moving_point_indexes)} points"
270+
f"{self.image2_title} - {len(self.moving_point_indexes)} points"
268271
)
269272
self.moving_axes.set_axis_off()
270273

‎tests/additional_dictionary.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -365,6 +365,7 @@ dataframe
365365
dataset
366366
datasets
367367
dciodvfy
368+
dcm
368369
debugOn
369370
defaultPixelValue
370371
deformable

0 commit comments

Comments
 (0)