Skip to content

Commit afaed76

Browse files
jonasteuwenjonasteuwen
committed
First examples commit
1 parent d134c76 commit afaed76

File tree

6 files changed

+468
-0
lines changed

6 files changed

+468
-0
lines changed

‎.gitignore

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,6 @@
1+
# emacs
2+
.#*
3+
14
# Byte-compiled / optimized / DLL files
25
__pycache__/
36
*.py[cod]

‎README.md

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
SimpleITK examples in Python
2+
============================
3+
4+
In this repository you will find a couple of examples on how to use SimpleITK with Python. If you want any specific example, please open an issue.
5+
An open source viewer (and more) for medical images is [Slicer 3D](http://slicer.org).
6+
7+
8+
Current:
9+
- *dcm_to_nrrd.py*: reads a complete dicom series from a folder and converts this to a nrrd file. Reads the window and level tags, and crops the image to this.
10+
- *resample_isotropically.py*: An example to read in a image file, and resample the image to a new grid.
11+
- *resample_tests.py*: Several ways to downsample an image.
12+
- *apply_lut.py*: in some DICOM files there is a tag VOILUTFunction. This is an example on how to apply this.
13+
14+
TODO
15+
----
16+
- Juypter notebooks
17+
- More examples

‎examples/apply_lut.py

Lines changed: 97 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,97 @@
1+
# encoding: utf-8
2+
"""Some dicom images require a VOILUTFunction to be applied before display.
3+
This example script reads in the image, apply the function and writes out
4+
back to dicom.
5+
6+
Reads:
7+
- http://dicom.nema.org/MEDICAL/dicom/2014c/output/chtml/part03/sect_C.11.2.html
8+
- http://dicom.nema.org/medical/dicom/2014c/output/chtml/part03/sect_C.7.6.3.html#sect_C.7.6.3.1.2
9+
"""
10+
import dicom
11+
import numpy as np
12+
import SimpleITK as sitk
13+
14+
15+
def display_metadata(dcm_file, expl):
16+
"""Reads a dicom file, and takes an explanation variable
17+
18+
Parameters
19+
----------
20+
dcm_file : str
21+
path to dicom file
22+
expl : str
23+
name of corresponding Window Center Width Explanation.
24+
25+
Returns
26+
-------
27+
center, width, VOILUTFunction, total_bits
28+
"""
29+
dcm = dicom.read_file(dcm_file)
30+
explanation = dcm.WindowCenterWidthExplanation
31+
idx = explanation.index(expl)
32+
center = dcm.WindowCenter[idx]
33+
width = dcm.WindowWidth[idx]
34+
35+
if hasattr(dcm, 'VOILUTFunction'):
36+
lut_func = dcm.VOILUTFunction.strip()
37+
else:
38+
lut_func = ''
39+
40+
if dcm.SamplesPerPixel == 1:
41+
if dcm.PhotometricInterpretation == 'MONOCHROME1':
42+
raise NotImplementedError('The image needs to be inverted. Not implemented.')
43+
44+
return center, width, lut_func, dcm.BitsStored
45+
46+
47+
def apply_window_level(dcm_file):
48+
image = sitk.ReadImage(dcm_file)
49+
center, width, lut_func, bits_stored = display_metadata(dcm_file, 'SOFTER')
50+
51+
if lut_func == 'SIGMOID':
52+
print('Applying `SIGMOID`.')
53+
image = sigmoid_lut(image, bits_stored, center, width)
54+
elif lut_func == '':
55+
print('Applying `LINEAR`.')
56+
image = linear_lut(image, bits_stored, center, width)
57+
else:
58+
raise NotImplementedError('`VOILUTFunction` can only be `SIGMOID`, `LINEAR` or empty.')
59+
60+
return image
61+
62+
63+
def sigmoid_lut(image, out_bits, center, width):
64+
"""If VOILUTFunction in the dicom header is equal to SIGMOID,
65+
apply this function for visualisation.
66+
"""
67+
array = sitk.GetArrayFromImage(image).astype(np.float)
68+
array = np.round((2**out_bits - 1)/(1 + np.exp(-4*(array - center)/width))).astype(np.int)
69+
ret_image = sitk.GetImageFromArray(array)
70+
ret_image.SetSpacing(image.GetSpacing())
71+
ret_image.SetOrigin(image.GetOrigin())
72+
ret_image.SetDirection(image.GetDirection())
73+
74+
return ret_image
75+
76+
77+
def linear_lut(image, out_bits, center, width):
78+
"""If VOILUTFunction in the dicom header is equal to LINEAR_EXACT,
79+
apply this function for visualisation.
80+
"""
81+
lower_bound = center - (width - 1)/2
82+
upper_bound = center + (width - 1)/2
83+
84+
min_max = sitk.MinimumMaximumImageFilter()
85+
min_max.Execute(image)
86+
87+
image = sitk.IntensityWindowing(image,
88+
lower_bound, upper_bound,
89+
min_max.GetMinimum(),
90+
min_max.GetMaximum())
91+
92+
return image
93+
94+
95+
def normalize_dicom(dcm_file, dcm_out):
96+
image = apply_window_level(dcm_file)
97+
sitk.WriteImage(image, dcm_out)

‎examples/dcm_to_nrrd.py

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
# encoding: utf-8
2+
import SimpleITK as sitk
3+
import dicom
4+
from tqdm import tqdm
5+
6+
7+
def dcm_to_nrrd(folder, to_path, intensity_windowing=True, compression=False):
8+
"""Read a folder with DICOM files and convert to a nrrd file.
9+
Assumes that there is only one DICOM series in the folder.
10+
11+
Parameters
12+
----------
13+
folder : string
14+
Full path to folder with dicom files.
15+
to_path : string
16+
Full path to output file (with .nrrd extension). As the file is
17+
outputted through SimpleITK, any supported format can be selected.
18+
intensity_windowing: bool
19+
If True, the dicom tags 'WindowCenter' and 'WindowWidth' are used
20+
to clip the image, and the resulting image will be rescaled to [0,255]
21+
and cast as uint8.
22+
compression : bool
23+
If True, the output will be compressed.
24+
"""
25+
reader = sitk.ImageSeriesReader()
26+
series_ids = reader.GetGDCMSeriesIDs(folder)
27+
28+
assert len(series_ids) == 1, 'Assuming only one series per folder.'
29+
30+
filenames = reader.GetGDCMSeriesFileNames(folder, series_ids[0])
31+
reader.SetFileNames(filenames)
32+
image = reader.Execute()
33+
34+
if intensity_windowing:
35+
dcm = dicom.read_file(filenames[0])
36+
assert hasattr(dcm, 'WindowCenter') and hasattr(dcm, 'WindowWidth'),\
37+
'when `intensity_windowing=True`, dicom needs to have the `WindowCenter` and `WindowWidth` tags.'
38+
center = dcm.WindowCenter
39+
width = dcm.WindowWidth
40+
41+
lower_bound = center - (width - 1)/2
42+
upper_bound = center + (width - 1)/2
43+
44+
image = sitk.IntensityWindowing(image,
45+
lower_bound, upper_bound, 0, 255)
46+
image = sitk.Cast(image, sitk.sitkUInt8) # after intensity windowing, not necessarily uint8.
47+
48+
writer = sitk.ImageFileWriter()
49+
if compression:
50+
writer.UseCompressionOn()
51+
52+
writer.SetFileName(to_path)
53+
writer.Execute(image)
54+
55+
56+
def main():
57+
folders = ['/data/folder1/', '/data/folder2']
58+
for folder in tqdm(folders):
59+
dcm_to_nrrd(folder, '~', intensity_windowing=True)
60+
61+
if __name__ == '__main__':
62+
main()

‎examples/resample_isotropically.py

Lines changed: 123 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,123 @@
1+
# encoding: utf-8
2+
"""A an example to read in a image file, and resample the image to a new grid.
3+
"""
4+
5+
import SimpleITK as sitk
6+
import os
7+
from glob import glob
8+
from tqdm import tqdm
9+
import numpy as np
10+
11+
12+
# https://github.com/SimpleITK/SlicerSimpleFilters/blob/master/SimpleFilters/SimpleFilters.py
13+
_SITK_INTERPOLATOR_DICT = {
14+
'nearest': sitk.sitkNearestNeighbor,
15+
'linear': sitk.sitkLinear,
16+
'gaussian': sitk.sitkGaussian,
17+
'label_gaussian': sitk.sitkLabelGaussian,
18+
'bspline': sitk.sitkBSpline,
19+
'hamming_sinc': sitk.sitkHammingWindowedSinc,
20+
'cosine_windowed_sinc': sitk.sitkCosineWindowedSinc,
21+
'welch_windowed_sinc': sitk.sitkWelchWindowedSinc,
22+
'lanczos_windowed_sinc': sitk.sitkLanczosWindowedSinc
23+
}
24+
25+
26+
27+
def resample_sitk_image(sitk_image, spacing=None, interpolator=None,
28+
fill_value=0):
29+
"""Resamples an ITK image to a new grid. If no spacing is given,
30+
the resampling is done isotropically to the smallest value in the current
31+
spacing. This is usually the in-plane resolution. If not given, the
32+
interpolation is derived from the input data type. Binary input
33+
(e.g., masks) are resampled with nearest neighbors, otherwise linear
34+
interpolation is chosen.
35+
36+
Parameters
37+
----------
38+
sitk_image : SimpleITK image or str
39+
Either a SimpleITK image or a path to a SimpleITK readable file.
40+
spacing : tuple
41+
Tuple of integers
42+
interpolator : str
43+
Either `nearest`, `linear` or None.
44+
fill_value : int
45+
46+
Returns
47+
-------
48+
SimpleITK image.
49+
"""
50+
51+
if isinstance(sitk_image, str):
52+
sitk_image = sitk.ReadImage(sitk_image)
53+
num_dim = sitk_image.GetDimension()
54+
55+
if not interpolator:
56+
interpolator = 'linear'
57+
pixelid = sitk_image.GetPixelIDValue()
58+
59+
if pixelid not in [1, 2, 4]:
60+
raise NotImplementedError(
61+
'Set `interpolator` manually, '
62+
'can only infer for 8-bit unsigned or 16, 32-bit signed integers')
63+
if pixelid == 1: # 8-bit unsigned int
64+
interpolator = 'nearest'
65+
66+
orig_pixelid = sitk_image.GetPixelIDValue()
67+
orig_origin = sitk_image.GetOrigin()
68+
orig_direction = sitk_image.GetDirection()
69+
orig_spacing = np.array(sitk_image.GetSpacing())
70+
orig_size = np.array(sitk_image.GetSize(), dtype=np.int)
71+
72+
if not spacing:
73+
min_spacing = orig_spacing.min()
74+
new_spacing = [min_spacing]*num_dim
75+
else:
76+
new_spacing = [float(s) for s in spacing]
77+
78+
assert interpolator in _SITK_INTERPOLATOR_DICT.keys(),\
79+
'`interpolator` should be one of {}'.format(_SITK_INTERPOLATOR_DICT.keys())
80+
81+
sitk_interpolator = _SITK_INTERPOLATOR_DICT[interpolator]
82+
83+
new_size = orig_size*(orig_spacing/new_spacing)
84+
new_size = np.ceil(new_size).astype(np.int) # Image dimensions are in integers
85+
new_size = [int(s) for s in new_size] # SimpleITK expects lists, not ndarrays
86+
87+
resample_filter = sitk.ResampleImageFilter()
88+
89+
resampled_sitk_image = resample_filter.Execute(sitk_image,
90+
new_size,
91+
sitk.Transform(),
92+
sitk_interpolator,
93+
orig_origin,
94+
new_spacing,
95+
orig_direction,
96+
fill_value,
97+
orig_pixelid)
98+
99+
return resampled_sitk_image
100+
101+
def main(REGEX_TO_IMAGES, WRITE_TO):
102+
for image in tqdm(glob(REGEX_TO_IMAGES)):
103+
tqdm.write('Resampling {}'.format(image))
104+
resampled_image = resample_sitk_image(
105+
image, spacing=[1, 1, 1],
106+
interpolator=None, fill_value=-1000
107+
)
108+
base_name = 'resampled_' + os.path.basename(image)
109+
write_to = os.path.join(WRITE_TO, base_name)
110+
tqdm.write('Writing resampled image to {}'.format(write_to))
111+
sitk.WriteImage(resampled_image, write_to)
112+
113+
114+
if __name__ == '__main__':
115+
base_path = './datasets/'
116+
write_path = './resampled/'
117+
118+
subpaths = ['path1/mhd/*.mhd', 'path2/*.nii.gz']
119+
write_paths = ['path1/', 'path2/']
120+
121+
for i, sub in enumerate(subpaths):
122+
main(base_path + sub, write_path + write_paths[i])
123+

0 commit comments

Comments
 (0)