A Multi-platform DICOM Toolbox for Academic Radiologists¶
The Digital Imaging and Communications in Medicine (DICOM) standard is used to store and transmit medical images. Radiology workstations attached to a picture archiving and communication systems (PACS) are used to display these images for diagnostic purposes.
However, diagnostic workstations are sometimes not optimal for academic research purposes. It is often necessary to anonymize DICOM images for research studies or for publication.
This Jupyter notebook contains code to accomplish several useful tasks with DICOM images.
- How do I open and view a DICOM image?
- How do I view and edit DICOM data elements?
- How do I anonymize a DICOM image?
- How do I anonymize a whole folder of DICOM images?
- How do I save a DICOM image as some other format (.JPG, .PNG, .TIFF)?
- How do I window/level and otherwise optimize a DICOM image?
- How do I invert a gray-scale image?
- How do I flip an image about its horizontal and vertical axes?
- How do I open a color DICOM image?
- How do I redact protected health information that is embedded within an image (i.e. US or radionuclide images)?
This toolbox includes several HIPAA-compliant DICOM images for your delectation. We have done our best to replace every bit of the original protected health information (PHI) with fictitious entries, using names and locations of characters from the works of Mark Twain.
Installing pydicom¶
This notebook assumes that you have already installed Jupyter and scientific Python. If not, we highly recommend the Anaconda Python distribution. The Anaconda distribution makes it easy to conveniently install Python, the Jupyter Notebook, and many other commonly used packages for scientific computing and data science.
While Anaconda contains many wonderful packages for science and mathematics, it does not know how to deal with DICOM images. However, there is a great and free add-on package called pydicom, and located at the Anaconda website.
You can install pydicom on your system very simply via this Jupyter notebook. Simply select the next cell in this notebook and hit shift-Return.
Please note that this only has to be done once.¶
!conda install -c conda-forge pydicom=0.9.9 --yes
Having installed pydicom, it is now time to work with some images.
Loading a few useful libraries¶
First, we will ask Python to import a few useful libraries that will do the heavy lifting for us:
- pydicom --- each time we start up our notebook, we need to reload this library. Oddly enough, loading this requires that we type "import dicom", rather than "import pydicom".
- matplotlilb --- we will use this rather powerful library to do most of the actual plotting.
- numpy --- an amazing library of mathematical routines
- scikit-image library --- a number of useful image-processing routines
- os --- this module to perform certain operating system functions, such as reading or writing a file
- glob --- this module finds all the pathnames matching a specified pattern, which is useful in performing operations on every file found in a given folder
- datetime --- this module allows one to manipulate dates and times in both simple and complex ways
# this following line tells Jupyter to display images here in the browser,
# rather than in separate window.
%matplotlib inline
# import pydicom library
import dicom
# import matplotlib and numpy
import matplotlib.pyplot as plt
import matplotlib.image as mpimage
import numpy as np
# load some handy functions from the scikit-image library
from skimage import exposure
import skimage.morphology as morp
from skimage.filters import rank
# import operating system and glob libraries
import os, glob
# import some useful date functions
from datetime import datetime
1. Reading and viewing a DICOM image¶
We have included several sample DICOM images in this tutorial. This example uses an image called: "Knee_MR.dcm". We will now show you how to load and view this image using Jupyter.
ds1 = dicom.read_file("Knee_MR.dcm") # read the DICOM image into memory
# tell matplotlib to display our images as a 6 x 6 inch image, with resolution of 100 dpi
plt.figure(figsize = (6,6), dpi=100)
# tell matplotlib to display our image, using a gray-scale lookup table.
plt.imshow(ds1.pixel_array, cmap=plt.cm.gray)
Voilà ! A coronal MR image is now visible for our delectation.
Now, let's look at the DICOM header to learn a bit more about this patient.
2. Examining data elements in the DICOM image header¶
The DICOM file format is extremely large and complex. It contains not only the image pixel data, but also a large image header with many DICOM_data_elements. These data elemenets contain patient identifiers and detailed information about the imaging device and how the image was obtained.
Thankfully, most DICOM files that we use will not include every possible data element from the DICOM standard. However, there will be a bunch of them. To see what is lurking under the hood on this image, select the code field below, and hit shift-return. Then, take a moment and scroll through some of the many data fields that were included in this image file.
ds1
Most of the time, we are only interested in a few of these fields. Rather than dumping the whole list, we can pick and choose the ones we want.
Since one may not always remember the exact name of a data element, the pydicom library provides a handy dir() method, useful during interactive sessions at the python prompt. Using this command, we can quickly find the names of all data elements that have the word "patient" in them thusly:
ds1.dir("patient")
We can now examine any of these one by one, as shown below:
ds1.PatientName # show patient name
ds1.PatientSex # show patient gender
ds1.AcquisitionDate # show date of study acquisition
ds1.PatientAge # show patient's age at time of scan
Changing Data Elements in a DICOM file¶
Sometimes, data (such as age) is sometimes entered incorrectly by a data entry minion. One way to check this is to see if the age entered is congruent with the date of birth and the date of study. We can do this with the following Python code:
DOB = datetime.strptime(ds1.PatientBirthDate, "%Y%m%d")
StudyDate = datetime.strptime(ds1.StudyDate, "%Y%m%d")
Age = (StudyDate - DOB).days/365
print ("Age = ", Age)
Sometimes a lot of time has passed since an image was done, and we might like to know how old our patient is currently. We can calculate this by subtracting the date of birth from the current date. In this case, the DICOM header suggests that this patient was born on January 1, 1870. The following Python code will grab today's date and do the calculation:
DOB = datetime.strptime(ds1.PatientBirthDate, "%Y%m%d")
now = datetime.now()
Age = (now - DOB).days/365
print ("Current age = ", Age, "years")
For an timeless literary character such as Huckleberry Finn, 147+ years does not seem excessive. For most actual patients, a result like this suggests a data entry error. If a quick audit were to show that the actual birthdate were July 4, 1970, we could fix the incorrect birthdate with the following code:
PatientBDActual = '19700704'
ds1.PatientBirthDate = PatientBDActual # change birthdate
We can change any other DICOM data field in a similar manner.
Once we have made all of the changes we desire, our next step is to save all of our changes in a new DICOM image file, which we will call "Knee_MR_corrected.dcm".
ds1.save_as("Knee_MR_corrected.dcm") # save as new file with corrected birthdate
3. Anonymizing a DICOM image file¶
In order to comply with the Health Insurance Portability and Accountability Act (HIPAA), it is occasionally necessary to anonymize the protected health information (PHI) contained in a DICOM image. Guidance as to which data elements constitute PHI can be found at the U. S. Health & Human Services website.
One can anonymize a DICOM image file by changing the PHI fields to blank or neutral values, as shown in the section above.
However, one can also completely delete the PHI fields from the DICOM image file using the "del" command. The lines of code shown below demonstrate how one might delete the fields for the following PHI elements:
- Patient name
- Accession number
- Medical record number
- Date of Birth
# delete the following DICOM data elements
del ds1.PatientName
del ds1.AccessionNumber
del ds1.PatientID
del ds1.PatientBirthDate
ds1.save_as("Knee_MR_anonymized.dcm") # save as new file with deleted PHI
4. Anonymizing a folder of DICOM image files¶
Anonymizing a single DICOM image file by hand can be tedious enough. Anonymizing a whole folder of DICOM images is cruel and unusual punishment. The following code will allow you to anonymize an entire folder of DICOM images without breaking a sweat. You're welcome.
The following code will look in a folder called "files_to_anonymize", which is at the same location as this DICOM toolbox notebook.
When you run the following code, it will do the following:
- create a list of every file in the "files_to_anonymize" folder
- sequentially open each DICOM file
- delete four PHI fields from the DICOM header
- save the newly anonymized DICOM file, overwriting the original
- print out the name of each DICOM file as it is processed
Caveat: If this program doesn't find any DICOM images in the folder called "files_to_anonymize", it will give an error message. It will also give an error message if a DICOM image lackss any of the four PHI data elements for name, accession number, patient ID or birthdate (i.e. if they have already been anonymized).
You are welcome to edit the list of data elements in the code below to delete any data elements you desire.
# grab the name of every file in our folder
for filename in glob.iglob('files_to_anonymize/*.dcm', recursive=True):
# read in a DICOM file into memory
ds = dicom.read_file(filename, force=True)
# get the name of the folder the dicomfiles are stored in
foldername=os.path.basename(os.path.dirname(os.path.dirname(filename)))
#delete patient name, accesion number, patient ID and patient birthdate
del ds.PatientName
del ds.AccessionNumber
del ds.PatientID
del ds.PatientBirthDate
# overwrite the original file
dicom.write_file(filename,ds)
# output progress to the screen
print(filename)
5. Saving a DICOM image in another format¶
There are a lot of good reasons for converting a DICOM image into another format, such as the .JPG, .PNG or .TIFF formats. The .JPG and .PNG formats are perfect for use with a website. When submitting an image for publication, journals often request that the image file be submitted in .TIFF format.
The following code will take the image data from the DICOM pixel array and save a copy in the .JPG, .PNG and .TIFF formats:
plt.imsave('knee_MR.jpg', ds1.pixel_array, cmap=plt.cm.gray)
plt.imsave('knee_MR.png', ds1.pixel_array, cmap=plt.cm.gray)
plt.imsave('knee_MR.tif', ds1.pixel_array, cmap=plt.cm.gray)
6. Optimizing a DICOM image¶
The Scientific Python (SciPy) library contains a number of handy image processing routines as part of the scikit-image library. These routines will come pre-installed as part of the Anaconda or Enthought Python distributions. In addition to scikit-image, we will also use some of the routines from the Numerical Python (numpy) library.
Let us look again at our original knee MR, and see how we might optimze it.
ds1 = dicom.read_file("Knee_MR.dcm") # read the knee MR image into memory
# tell matplotlib to display our images as 6 x 6 inch image, with resolution of 100 dpi
plt.figure(figsize = (6,6), dpi=100)
# tell matplotlib to display our image, using a gray-scale lookup table.
plt.imshow(ds1.pixel_array, cmap=plt.cm.gray)
Image histograms can be useful¶
This knee MR is nice-looking, but it is a bit dark. Let's have a look at the image histogram. Histograms can often give us some helpful clues as the best way to optimize an image.
plt.hist(ds1.pixel_array.flatten(), bins=64) # calculate a histogram of our image
plt.show() # display that histogram
print ("pixel array = ", ds1.pixel_array.shape)
print("minimum value = ", np.amin(ds1.pixel_array)) # find minimum pixel value in the image array
print("maximum value = ", np.amax(ds1.pixel_array)) # find maximum pixel value in the image array
Our MR image has 512 * 512 = 262,144 pixels. The pixels with the lowest value (0 in this case) are displayed as black. The pixels with the highest value (513 in this case) are displayed as white. All other values are assigned a value somewhere on a linear scale between black and white.
This histogram of the our image shows that about 67,000 of the pixels are black (value = 0), as shown by the tallest bar over to the far left. Although we know that our highest pixel value is 513, this histogram actually doesn't show much visible image information to the right of 260. The histogram also shows us that pixel values of 260 will be mapped to a shade of gray halfway between black and white.
Window and level settings (aka contrast stretching)¶
Besides useful anatomical details, medical images often contain other things that are less relevant (the air outside the patient, right/left markers, etc.). Depending on how the image was obtained, a number of gray shades may be devoted to these less interesting structures. We would therefore like to adjust our image so that it is optimized for the stuctures we consider useful.
As a rule of thumb, gray-scale images tend to look better when the left-hand edge of the useful pixel values is mapped to black and the right-hand edge of the useful pixel values is mapped to white. When one does this on a workstation, it is called "windowing" and "leveling". One might also call this "contrast stretching".
We will use the exposure.rescale_intensity function from the sci-kit image library to do the heavy lifting for us. We will set the lower and upper edges of our window to various values and see how that affects the image. We could use actual pixel values to set the edges of our window. However, rather than having to worry about the exact values, I find it easier to set the window edges using percentiles of the pixel data.
In the example below, we will leave the lower edge of the window at the 0th percentile, but drop the upper edge of the window to the 99.5th percentile.
# Contrast stretching try 1
plt.figure(figsize = (6,6), dpi=100)
p_lo, p_hi = np.percentile(ds1.pixel_array, (0, 99.5))
img_rescale_1 = exposure.rescale_intensity(ds1.pixel_array, in_range=(p_lo, p_hi))
figure = plt.imshow(img_rescale_1, cmap=plt.cm.gray)
That rescaled image looks a lot snappier to my eye.
Let's see what the histogram looks like for this re-windowed image:
plt.hist(img_rescale_1.flatten(), bins=64)
plt.show()
print("p_lo = ", p_lo)
print("p_hi = ",p_hi)
From the histogram, we can see that the 99.5th percentile is equivalent to a pixel value of 253. All pixels with values greater than or equal to 253 are mapped to white. The histogram now shows that the window edges are indeed now very close to the edges of the visible image data.
Next, let's try dropping the upper edge of the window down to the 95th percentile.
# Contrast stretching try 2
plt.figure(figsize = (6,6), dpi=100)
p_lo, p_hi = np.percentile(ds1.pixel_array, (0, 95))
img_rescale_2 = exposure.rescale_intensity(ds1.pixel_array, in_range=(p_lo, p_hi))
figure = plt.imshow(img_rescale_2, cmap=plt.cm.gray)