SimpleITK  
DicomSeriesFromArray/DicomSeriesFromArray.py
1# =========================================================================
2#
3# Copyright NumFOCUS
4#
5# Licensed under the Apache License, Version 2.0 (the "License");
6# you may not use this file except in compliance with the License.
7# You may obtain a copy of the License at
8#
9# http://www.apache.org/licenses/LICENSE-2.0.txt
10#
11# Unless required by applicable law or agreed to in writing, software
12# distributed under the License is distributed on an "AS IS" BASIS,
13# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14# See the License for the specific language governing permissions and
15# limitations under the License.
16#
17# =========================================================================
18
19""" A SimpleITK example demonstrating how to write a DICOM series. """
20
21import sys
22import time
23import os
24import numpy as np
25import SimpleITK as sitk
26
27pixel_dtypes = {"int16": np.int16, "float64": np.float64}
28
29
30def writeSlices(series_tag, in_image, out_dir, i):
31 """ Write slices to output directory """
32 image_slice = in_image[:, :, i]
33
34 # Tags shared by the series.
35 list(
36 map(
37 lambda tag_value: image_slice.SetMetaData(tag_value[0], tag_value[1]),
38 series_tag,
39 )
40 )
41
42 # Slice specific tags.
43 # Instance Creation Date
44 image_slice.SetMetaData("0008|0012", time.strftime("%Y%m%d"))
45 # Instance Creation Time
46 image_slice.SetMetaData("0008|0013", time.strftime("%H%M%S"))
47
48 # Setting the type to CT so that the slice location is preserved and
49 # the thickness is carried over.
50 image_slice.SetMetaData("0008|0060", "CT")
51
52 # (0020, 0032) image position patient determines the 3D spacing between
53 # slices.
54 # Image Position (Patient)
55 image_slice.SetMetaData(
56 "0020|0032",
57 "\\".join(map(str, in_image.TransformIndexToPhysicalPoint((0, 0, i)))),
58 )
59 # Instance Number
60 image_slice.SetMetaData("0020|0013", str(i))
61
62 # Write to the output directory and add the extension dcm, to force
63 # writing in DICOM format.
64 writer.SetFileName(os.path.join(out_dir, str(i) + ".dcm"))
65 writer.Execute(image_slice)
66
67
68if len(sys.argv) < 3:
69 print(
70 "Usage: python "
71 + __file__
72 + " <output_directory> ["
73 + ", ".join(pixel_dtypes)
74 + "]"
75 )
76 sys.exit(1)
77
78# Create a new series from a numpy array
79try:
80 pixel_dtype = pixel_dtypes[sys.argv[2]]
81except KeyError:
82 pixel_dtype = pixel_dtypes["int16"]
83
84new_arr = np.random.uniform(-10, 10, size=(3, 4, 5)).astype(pixel_dtype)
85new_img = sitk.GetImageFromArray(new_arr)
86new_img.SetSpacing([2.5, 3.5, 4.5])
87
88# Write the 3D image as a series
89# IMPORTANT: There are many DICOM tags that need to be updated when you modify
90# an original image. This is a delicate operation and requires
91# knowledge of the DICOM standard. This example only modifies some.
92# For a more complete list of tags that need to be modified see:
93# http://gdcm.sourceforge.net/wiki/index.php/Writing_DICOM
94# If it is critical for your work to generate valid DICOM files,
95# It is recommended to use David Clunie's Dicom3tools to validate
96# the files:
97# http://www.dclunie.com/dicom3tools.html
98
99writer = sitk.ImageFileWriter()
100# Use the study/series/frame of reference information given in the meta-data
101# dictionary and not the automatically generated information from the file IO
102writer.KeepOriginalImageUIDOn()
103
104modification_time = time.strftime("%H%M%S")
105modification_date = time.strftime("%Y%m%d")
106
107# Copy some of the tags and add the relevant tags indicating the change.
108# For the series instance UID (0020|000e), each of the components is a number,
109# cannot start with zero, and separated by a '.' We create a unique series ID
110# using the date and time. Tags of interest:
111direction = new_img.GetDirection()
112series_tag_values = [
113 ("0008|0031", modification_time), # Series Time
114 ("0008|0021", modification_date), # Series Date
115 ("0008|0008", "DERIVED\\SECONDARY"), # Image Type
116 (
117 "0020|000e",
118 "1.2.826.0.1.3680043.2.1125." + modification_date + ".1" + modification_time,
119 ), # Series Instance UID
120 (
121 "0020|0037",
122 "\\".join(
123 map(
124 str,
125 (
126 direction[0],
127 direction[3],
128 direction[6],
129 direction[1],
130 direction[4],
131 direction[7],
132 ),
133 )
134 ),
135 ), # Image Orientation
136 # (Patient)
137 ("0008|103e", "Created-SimpleITK"), # Series Description
138]
139
140if pixel_dtype == np.float64:
141 # If we want to write floating point values, we need to use the rescale
142 # slope, "0028|1053", to select the number of digits we want to keep. We
143 # also need to specify additional pixel storage and representation
144 # information.
145 rescale_slope = 0.001 # keep three digits after the decimal point
146 series_tag_values = series_tag_values + [
147 ("0028|1053", str(rescale_slope)), # rescale slope
148 ("0028|1052", "0"), # rescale intercept
149 ("0028|0100", "16"), # bits allocated
150 ("0028|0101", "16"), # bits stored
151 ("0028|0102", "15"), # high bit
152 ("0028|0103", "1"),
153 ] # pixel representation
154
155# Write slices to output directory
156list(
157 map(
158 lambda i: writeSlices(series_tag_values, new_img, sys.argv[1], i),
159 range(new_img.GetDepth()),
160 )
161)
162
163# Re-read the series
164# Read the original series. First obtain the series file names using the
165# image series reader.
166data_directory = sys.argv[1]
167series_IDs = sitk.ImageSeriesReader.GetGDCMSeriesIDs(data_directory)
168if not series_IDs:
169 print(
170 'ERROR: given directory "'
171 + data_directory
172 + '" does not contain a DICOM series.'
173 )
174 sys.exit(1)
175series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
176 data_directory, series_IDs[0]
177)
178
179series_reader = sitk.ImageSeriesReader()
180series_reader.SetFileNames(series_file_names)
181
182# Configure the reader to load all of the DICOM tags (public+private):
183# By default tags are not loaded (saves time).
184# By default if tags are loaded, the private tags are not loaded.
185# We explicitly configure the reader to load tags, including the
186# private ones.
187series_reader.LoadPrivateTagsOn()
188image3D = series_reader.Execute()
189print(image3D.GetSpacing(), "vs", new_img.GetSpacing())
190sys.exit(0)
Write out a SimpleITK image to the specified file location.
Read series of image files into a SimpleITK image.