800+ IT
News
als RSS Feed abonnieren📚 Dealing with DICOM Using ImageIO Python Package
💡 Newskategorie: AI Nachrichten
🔗 Quelle: towardsdatascience.com
Medical Images == DICOM
In this post, we will use the ImageIO Python package to read DICOM files, extract metadata and attributes, and plot image slices using interactive slider widgets using Ipywidgets.
Background Info
DICOM is the bread and butter of medical imaging systems. If you are a Biomedical Engineer, an IT Specialist in the healthcare field, or a Healthcare Data Scientist/Analyst, you’ve probably used or at least heard about DICOM because it’s everywhere related to medical imaging systems.
DICOM, Digital Imaging and Communications in Medicine, is a set of standards that are created to allow communication across multiple modalities, e.g. CT, MRI, Ultrasound, etc., between multiple manufacturers so that all medical machines, that are DICOM-compliant of course, can speak the same language when sending information across a network.
DICOM files usually have information other than just data of pixels. DICOM Standards provide a set of attributes to standardize the data stored in DICOM files. Each attribute has a specific data type (e.g. string, integer, float, etc.) with a fixed value representation. These attributes are called metadata. Examples of DICOM metadata:
- Patient-related attributes: name, age, sex, etc.
- Modality-related attributes: modality (CT, MRI, Ultrasound, etc.), manufacturer, acquisition date, etc.
- Image-related attributes: shape, sampling, aspect ratio, pixel data, etc.
Note that the image-related attributes I mentioned above are very important to understand as they are useful to implement real-world calculations on medical images. We will heavily rely on them down below. It’s preferred that you understand what each represents.
There are dozens of attributes that can characterize DICOM files. So, you cannot just read all of them. You need to focus only on the attributes that you might counter during your work on a DICOM-related thing. For this purpose, I highly recommend this great DICOM Standard Browser which was built by Innolitics, you just need to search for the attribute you want to learn about.
Now let’s get down to our business… We will use python to easily surf DICOM files. We will focus on using the (ImageIO) package to present some of the basics about DICOM files. We will discuss the following:
- Read DICOM files using ImageIO.
- DICOM attributes, Metadata.
- Access a specific DICOM attribute.
- Image representation using Matplotlib package.
- Stack and read multiple slices.
- Understand Pixel Spacing, Shape, Slice Thickness, Aspect Ratio, and the field of view.
- Build an interactive image representation along the axial, coronal, and sagittal planes using Ipywidgets.
Let’s Code
First, let’s import the required packages. We will use (ImageIO) to deal with DICOM files, (NumPy) as the pixel data are read as a NumPy-array, and (matplotlib) to visualize the images. An additional package (Ipywidgets) is also used to build an interactive slider that we can use to scroll among multiple image slices.
I’m using Google Colab for coding in Python
import ipywidgets as widgets
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import imageio
DICOM Dataset
There are alot of DICOM samples on the internet. And you can use your own DICOM files. For this blog, I chose a dataset of 99 slices of chest-CT scan for one patient. You can find it on Kaggle. I saved the dataset on my Drive so I can easily have access through GoogleColab.
Reading DICOM Files
You can easily read a DICOM file and store it in a variable using (.imread) from ImageIO package.
#Reading a DICOM Image
im = imageio.imread('/content/drive/MyDrive/Datasets/Kaggle/DICOM/dicom_lung/000009.dcm')
DICOM Attributes
Using (.meta) outputs a dictionary that contains the attributes, the metadata, of this DICOM file.
#DICOM Metadata
im.meta
Dict([('TransferSyntaxUID', '1.2.840.10008.1.2.1'),
('SOPClassUID', '1.2.840.10008.5.1.4.1.1.2'),
('SOPInstanceUID',
'1.3.6.1.4.1.14519.5.2.1.7085.2626.397047382026069586801778973091'),
('StudyDate', '20100227'),
('SeriesDate', '20100227'),
('AcquisitionDate', '20100227'),
('ContentDate', '20100227'),
('StudyTime', '161937.171'),
('SeriesTime', '162536.14 '),
('AcquisitionTime', '162208.162527 '),
('ContentTime', '162208.162527 '),
('Modality', 'CT'),
('Manufacturer', 'SIEMENS'),
('StudyDescription', 'CT CHEST W IV CONTRAST'),
('SeriesDescription', 'LUNG 3.0 B70f'),
('PatientName', 'C3N-00247'),
('PatientID', 'C3N-00247'),
('PatientBirthDate', ''),
('PatientSex', 'F '),
('PatientAge', '077Y'),
('StudyInstanceUID',
'1.3.6.1.4.1.14519.5.2.1.7085.2626.258626612405225511766549337110'),
('SeriesInstanceUID',
'1.3.6.1.4.1.14519.5.2.1.7085.2626.242193282649561899185427104083'),
('SeriesNumber', 2),
('AcquisitionNumber', 2),
('InstanceNumber', 1),
('ImagePositionPatient', (-143.2177734375, -287.2177734375, 6.0)),
('ImageOrientationPatient', (1.0, 0.0, 0.0, 0.0, 1.0, 0.0)),
('SamplesPerPixel', 1),
('Rows', 512),
('Columns', 512),
('PixelSpacing', (0.564453125, 0.564453125)),
('BitsAllocated', 16),
('BitsStored', 12),
('HighBit', 11),
('PixelRepresentation', 0),
('RescaleIntercept', -1024.0),
('RescaleSlope', 1.0),
('PixelData',
b'Data converted to numpy array, raw data removed to preserve memory'),
('shape', (512, 512)),
('sampling', (0.564453125, 0.564453125))])Since the metadata is stored as a dictionary, you can see the keys, the attributes’ names, of the DICOM file:
# The Attributes of the DICOM File
im.meta.keys()
odict_keys(['TransferSyntaxUID', 'SOPClassUID', 'SOPInstanceUID',
'StudyDate', 'SeriesDate', 'AcquisitionDate', 'ContentDate',
'StudyTime', 'SeriesTime', 'AcquisitionTime', 'ContentTime',
'Modality', 'Manufacturer', 'StudyDescription', 'SeriesDescription',
'PatientName', 'PatientID', 'PatientBirthDate', 'PatientSex',
'PatientAge', 'StudyInstanceUID', 'SeriesInstanceUID',
'SeriesNumber', 'AcquisitionNumber', 'InstanceNumber',
'ImagePositionPatient', 'ImageOrientationPatient',
'SamplesPerPixel', 'Rows', 'Columns', 'PixelSpacing',
'BitsAllocated', 'BitsStored', 'HighBit', 'PixelRepresentation',
'RescaleIntercept', 'RescaleSlope', 'PixelData', 'shape',
'sampling'])
Note: The attributes we see here might not be the whole attributes contained in this DICOM file. This is because ImageIO deals with DICOM images as light as possible. An optional reading is to check this link to see the dictionary of attributes that ImageIO supports.
Access a Specific DICOM Attribute
Accessing a specific attribute, “Modality” for example, can be done like the following:
#The Modality Attribute
im.meta['Modality']
'CT'
Image Representation Using Matplotlib Package
Usually, the “Pixel Data” attribute has the values of pixels. Let’s access this attribute just like before:
# Access the Pixel Data attribute
im.meta['PixelData']
b'Data converted to numpy array, raw data removed to preserve memory'
As it says, the pixel values are stored as a NumPy array. And this is useful because NumPy is suitable for dealing with arrays and their calculation fastly. Now let’s see the pixel values:
#print the image Numpy-array
im
Array([[ -985, -990, -999, ..., -1017, -1008, -971],
[-1016, -984, -963, ..., -1000, -1009, -999],
[-1024, -1008, -996, ..., -979, -1021, -987],
...,
[ -920, -942, -944, ..., -893, -917, -955],
[ -871, -879, -905, ..., -895, -869, -867],
[ -876, -855, -873, ..., -933, -982, -936]], dtype=int16)
We only see a bunch of numbers where each number represents the pixel value in the Hounsfield Unit (HU) form. You can slice the image array using NumPy-array indices:
#Slicing the first 5 rows and first 5 columns
im[0:5, 0:5]
Array([[ -985, -990, -999, -982, -1011],
[-1016, -984, -963, -978, -1005],
[-1024, -1008, -996, -969, -992],
[-1006, -984, -997, -963, -1002],
[ -970, -988, -1006, -992, -999]], dtype=int16)
After understanding what Pixel Data represents, let’s show the image of these pixels.
#Show the image with gray color-map
plt.imshow(im, cmap='gray')
#Don't show tha axes
plt.axis('off')
#Add a title to the plot
plt.title('Axial Slice')
plt.show()
Stack and Read Multiple Slices
A DICOM file can have multiple frames that might be found stacked in the Pixel Data attribute. You can see that in DICOM files which have videos or colored images, RGB channels. Also, sometimes you find multiple DICOM files in one folder in which each DICOM file contains one frame, or slice, of the same patient. In this case, we have to stack these DICOM files ourselves. Fortunately, this is easily done using (.volread) from ImageIO.
#Stacking 99 slices
vol = imageio.volread('/content/drive/MyDrive/Datasets/Kaggle/DICOM/dicom_lung')
Reading DICOM (examining files): 99/99 files (100.0%)
Found 1 correct series.
Reading DICOM (loading data): 99/99 (100.0%)
You can access a specific slice using array indices:
#Access the first slice
vol[0,:,:]
Array([[ -985, -990, -999, ..., -1017, -1008, -971],
[-1016, -984, -963, ..., -1000, -1009, -999],
[-1024, -1008, -996, ..., -979, -1021, -987],
...,
[ -920, -942, -944, ..., -893, -917, -955],
[ -871, -879, -905, ..., -895, -869, -867],
[ -876, -855, -873, ..., -933, -982, -936]], dtype=int16)
Also, image representation for a single specific slice can be done just like before:
#Show the first slice.
plt.imshow(vol[0,:,:], cmap='gray')
#Don't show the axis
plt.axis('off')
#Add a title
plt.title('Axial Slice')
plt.show()
Dealing with stacked slices is very useful for extracting other planes (Sagittal, and Coronal). Also, it is helpful for radiologists to scroll over multiple slices while diagnosing. But to show the three planes and scroll upon slices correctly, we, engineers, have to understand the meaning of the specific attributes. These attributes are Sampling, Shape, and Aspect Ratio.
#Introduce the metadata of the stacked images.
vol.meta
Dict([('TransferSyntaxUID', '1.2.840.10008.1.2.1'),
('SOPClassUID', '1.2.840.10008.5.1.4.1.1.2'),
('SOPInstanceUID',
'1.3.6.1.4.1.14519.5.2.1.7085.2626.397047382026069586801778973091'),
('StudyDate', '20100227'),
('SeriesDate', '20100227'),
('AcquisitionDate', '20100227'),
('ContentDate', '20100227'),
('StudyTime', '161937.171'),
('SeriesTime', '162536.14 '),
('AcquisitionTime', '162208.162527 '),
('ContentTime', '162208.162527 '),
('Modality', 'CT'),
('Manufacturer', 'SIEMENS'),
('StudyDescription', 'CT CHEST W IV CONTRAST'),
('SeriesDescription', 'LUNG 3.0 B70f'),
('PatientName', 'C3N-00247'),
('PatientID', 'C3N-00247'),
('PatientBirthDate', ''),
('PatientSex', 'F '),
('PatientAge', '077Y'),
('StudyInstanceUID',
'1.3.6.1.4.1.14519.5.2.1.7085.2626.258626612405225511766549337110'),
('SeriesInstanceUID',
'1.3.6.1.4.1.14519.5.2.1.7085.2626.242193282649561899185427104083'),
('SeriesNumber', 2),
('AcquisitionNumber', 2),
('InstanceNumber', 1),
('ImagePositionPatient', (-143.2177734375, -287.2177734375, 6.0)),
('ImageOrientationPatient', (1.0, 0.0, 0.0, 0.0, 1.0, 0.0)),
('SamplesPerPixel', 1),
('Rows', 512),
('Columns', 512),
('PixelSpacing', (0.564453125, 0.564453125)),
('BitsAllocated', 16),
('BitsStored', 12),
('HighBit', 11),
('PixelRepresentation', 0),
('RescaleIntercept', -1024.0),
('RescaleSlope', 1.0),
('PixelData', b'Deferred loading of pixel data'),
('shape', (99, 512, 512)),
('sampling', (3.0, 0.564453125, 0.564453125))])Notice the differences in “shape” and “sampling” attributes for the stacked images. This is done by ImageIO reading when using (.volread).
Understand “Shape”, “Sampling” and “PixelAspectRatio” attributes
- Shape: It’s simply the number of rows and columns in each slice. Since we’re dealing with multiple slices, there will be a third dimension equal to the number of slices that are stacked together. In our example, the shape of the stacked images is 99 slices, 512 rows, and 512 columns.
# The shape of the stacked images in each plane
# (Axial, Coronal, and Sagittal, respectively)
n0, n1, n2 = vol.shape
# Print the ouput
print("Number of Slices:\n\t", "Axial=", n0, "Slices\n\t",
"Coronal=", n1, "Slices\n\t",
"Sagittal=", n2, "Slices")
Number of Slices:
Axial= 99 Slices
Coronal= 512 Slices
Sagittal= 512 Slices
- Sampling: If you search for “Sampling” at DICOM Standard, you won’t get a direct answer. Because it’s a special variable built in ImageIO. This Sampling combination of two attributes, SliceThickness, and PixelSpacing. SliceThickness is simply the nominal slice thickness in mm. As for the PixelSpacing attribute, it’s the physical distance, in the patient. It’s specified by a numeric pair: The first value is the row spacing in mm, that is the spacing between the centers of adjacent rows, or vertical spacing. The second value is the column spacing in mm, that is the spacing between the centers of adjacent columns, or horizontal spacing.
# The sampling of the stacked images in each plane
# (Axial, Coronal, and Sagittal, respectively)
d0, d1, d2 = vol.meta['sampling'] # in mm
# Print the output
print("Sampling:\n\t", "Axial=", d0, "mm\n\t",
"Coronal=", d1, "mm\n\t",
"Sagittal=", d2, "mm")
Sampling:
Axial= 3.0 mm
Coronal= 0.564453125 mm
Sagittal= 0.564453125 mm
- Pixel Aspect Ratio: Ratio of the vertical size and horizontal size of the pixels in the image along a specific axis. Notice that “PixelAspectRatio” is not in the metadata of the stacked images. But that’s okay because we can calculate the aspect ratio along each axis by using the “sampling” parameters. The image below represents the PixelAspectRatio attribute.
# The aspect ratio along the axial plane
axial_asp = d1/d2
# The aspect ratio along the sagittal plane
sagittal_asp = d0/d1
# The aspect ratio along the coronal plane
coronal_asp = d0/d2
# Print the output
print("Pixel Aspect Ratio:\n\t", "Axial=", axial_asp, "\n\t",
"Coronal=", coronal_asp, "\n\t",
"Sagittal=", sagittal_asp)
Pixel Aspect Ratio:
Axial= 1.0
Coronal= 5.314878892733564
Sagittal= 5.314878892733564
Field of View
By multiplying the “shape” parameters by the “sampling” parameters, the physical space in mm, along each axis, we get the Field of View along each axis.
print("Field of View:\n\t", "Axial=", n0*d0, "mm\n\t",
"Coronal=", n1*d1, "mm\n\t",
"Sagittal=", n2*d2, "mm")Field of View:
Axial= 297.0 mm
Coronal= 289.0 mm
Sagittal= 289.0 mm
Build an interactive image representation using Ipywidgets
Now, after understanding what sampling and shape mean, let’s use them to do a correct representation of the images with the correct aspect ratio along each axis. We’ll use Ipywidgets to build an integer slider which we can scroll between different slices of the specific planes (Axial, Sagittal, or Coronal).
# Add a slider that starts with 0 and ends at the number of
# slices along the axial plane, n0=99.
@widgets.interact(axial_slice=(0,n0-1))
# Define the function that shows the images of the specified slice number.
# It starts with the 10th slice. And you can scroll over any slice
# using the slider.
def axial_slicer(axial_slice=50):
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
# Show the image of the specified slice number in 'gray' color-map
# and axial aspect ratio
ax.imshow(vol[axial_slice,:,:], cmap='gray', aspect=axial_asp)
# Don't show the axis
ax.axis('off')
Now let’s show the images of three slices each along a specific axis. I intentionally showed the axes to easily understand the shape of the slices along each axis.
# Define a figure with 1 row and 3 columns of plots to show
# the images along the three planes
fig, ax = plt.subplots(1, 3, figsize=(8, 8))
# Axial Plane: show the 10th slice
ax[0].imshow(vol[10,:,:], cmap='gray', aspect= axial_asp)
#ax[0].axis('off')
ax[0].set_title('Axial')
# Coronal Plane: show the slice 100
ax[1].imshow(vol[:,100,:],cmap='gray', aspect= coronal_asp)
#ax[1].axis('off')
ax[1].set_title('Coronal')
# Sagittal Plane: show the slice 100
ax[2].imshow(vol[:,:,100], cmap='gray', aspect= sagittal_asp)
#ax[2].axis('off')
ax[2].set_title('Sagittal')
plt.show()
Let’s try the wrong way… let’s show the images ignoring the aspect ratio of each axis.
# Define a figure with 1 row and 3 columns of plots to show the
# images along the three planes
fig, ax = plt.subplots(1, 3, figsize=(8, 8))
# Axial Plane: show the 10th slice
ax[0].imshow(vol[10,:,:], cmap='gray')
#ax[0].axis('off')
ax[0].set_title('Axial')
# Coronal Plane: show the slice 100
ax[1].imshow(vol[:,100,:],cmap='gray')
#ax[1].axis('off')
ax[1].set_title('Coronal')
# Sagittal Plane: show the slice 100
ax[2].imshow(vol[:,:,100], cmap='gray')
#ax[2].axis('off')
ax[2].set_title('Sagittal')
plt.show()
We can see that only the axial plane is not messy. This is because the argument “aspect” is set to 1 by default, which is the aspect ratio of the axial plane in our case (d1/d2 = 0.564453125 / 0.564453125 = 1).
Our Final Code
Let’s go beyond just one slider with one image plane. We will build a figure of 3 plots that contain the three planes (Axial, Coronal, and Sagittal, respectively). We can specify the slice number using the specific slider for each plane. Don’t forget to add the correct aspect ratio so it won’t get messy just like above.
# Add three sliders that start with 0 and ends at the number of slices
# along each plane.
# Axial: n0=99 slice
# Corornal: n1=512 slice
# Sagittal: n2=512 slice
@widgets.interact(axial_slice=(0,n0-1), coronal_slice=(0,n1-1),\
sagittal_slice=(0,n2-1))
def slicer(axial_slice, coronal_slice, sagittal_slice=100):
fig, ax = plt.subplots(1, 3, figsize=(12, 12))
# Show the specfied slice on the axial plane with 'gray' color-map
# and axial aspect ratio.
ax[0].imshow(vol[axial_slice,:,:], cmap='gray', aspect= axial_asp)
ax[0].axis('off')
ax[0].set_title('Axial')
# Show the specified slice on the coronal plane with 'gray' color-map
# and coronal aspect ratio.
ax[1].imshow(vol[:,coronal_slice,:],cmap='gray', aspect= coronal_asp)
ax[1].axis('off')
ax[1].set_title('Coronal')
# Show the specified slice on the sagittal plane with 'gray' color-map
# and sagittal aspect ratio.
ax[2].imshow(vol[:,:,sagittal_slice], cmap='gray', aspect= sagittal_asp)
ax[2].axis('off')
ax[2].set_title('Sagittal')
Conclusion
- We’ve seen how to use the ImageIO package to read DICOM files. And the power of using (.volread) to stack the images and read them more conveniently.
- We’ve learned that a DICOM file has metadata, a set of attributes, along with Pixel Data. And we’ve introduced some of the basics of specific image-related attributes, such as PixelSapcing, SliceThickness, Shape, and PixelAspectRatio.
- We’ve pointed out the importance of the aspect ratio parameter when plotting the images along the axes and how the image gets messy if it’s not correctly specified.
- We’ve built a figure that shows the specified slice number along the axial, coronal, and sagittal planes using the slider of the Ipywidgets package.
Thanks For Reading…
Recommendation
- For more details about DICOM: refer to the previous blog, What is DICOM?
- DICOM Metadata — A Useful Resource for Big Data Analytics:
This article provides an overview of new ways to represent data by combining patient access and DICOM information, advanced use of medical imaging metadata, analysis of radiation dose and image segmentation, and deep learning for feature engineering to enrich data. - Biomedical Image Analysis in Python Course on DataCamp. This course was my first dive into DICOM images and medical image processing using Python. The information in this blog heavily depends on the first chapter of this course, which is the only free chapter.
- Ipywidget Documentation is a good start to learning about building interactive widgets.
References
[1] DataCamp, First Chapter of Biomedical Image Analysis in Python Course, First Chapter, Exploration.
[2] C. Rossant, IPython Interactive Computing and Visualization Cookbook, (2018), Mastering widgets in the Jupyter Notebook.
[3] ImageIO, GitHub, ImageIO-Plugin-DICOM, (2022), DICOM Attributes, MINIDICT.
[4] Innolitics, DICOM Standard Browser, (2022), Slice Thickness Attribute.
[5] Innolitics, DICOM Standard Browser, (2022), Pixel Spacing Attribute.
[6] Innolitics, DICOM Standard Browser, (2022), Aspect Ratio Attribute.
Dealing with DICOM Using ImageIO Python Package was originally published in Towards Data Science on Medium, where people are continuing the conversation by highlighting and responding to this story.
...