DICOM Processing and Segmentation in Python

DICOM is a pain in the neck.  It also happens to be very helpful.  As clinical radiologists, we expect post-processing, even taking them for granted. However, the magic that occurs behind the scenes is no easy feat, so let’s explore some of that magic.

In this quest, we will be starting from raw DICOM images. We will extract voxel data from DICOM into numpy arrays, and then perform some low-level operations to normalize and resample the data, made possible using information in the DICOM headers.

The remainder of the Quest is dedicated to visualizing the data in 1D (by histogram), 2D, and 3D. Finally, we will create segmentation masks that remove all voxel except for the lungs.

Processing raw DICOM with Python is a little like excavating a dinosaur – you’ll want to have a jackhammer to dig, but also a pickaxe and even a toothbrush for the right situations. Python has all the tools, from pre-packaged imaging process packages handling gigabytes of data at once to byte-level operations on a single voxel.

To follow along, set up your computer using the following Python tutorials: Alternatively, start a free Jupyter notebook from Azure Notebooks.
Howard Chen
(Howard) Po-Hao Chen, MD MBA is a radiology chief resident at Hospital of the University of Pennsylvania. He has an interest in data-driven radiology, quality improvement, and innovation.

Howard will finish training with fellowships in musculoskeletal radiology and nuclear medicine in June 2018 from University of Pennsylvania.

11 Responses to “DICOM Processing and Segmentation in Python

  • Hi Howard,

    Thanks for the tutorial, I have been looking at the images, and I think I understand most of the preprocessing. However, after the resampling (taking into account we have different pixel spacing and slice thickness), we obtain volumes of different dimensions for each patient. How would you recommend resizing in order to get for all patients a volume of dimensions XxYxZ, with spacing of 1x1x1? For me this makes more sense to be able to compare them..
    Thanks in advance,

    Diana

    • Hey Diana!

      If I understand the question correctly, you have a set of DICOM images, each with different real-life size (L * W * H mm), all of which you want to be able to resample to the same pixel dimensions (X * Y * Z) while maintaining 1 x 1 x 1 mm voxel sizes. It turns out it is a natural side effect: resampling isotropically means so all voxels are the same size but each exam will be different sizes (this is a common approach) because patients are different sizes. This is because CT scans are commonly obtained at a constant 512 x 512 matrix. Since each patient is different in size, what changes is the “zoom” (field-of-view), so each voxel represents a different number of mm in real life.

      You can imagine that if we scanned an 85-pound patient at the same “zoom” as a 190-pound patient, you wouldn’t want the scan to occupy only the middle 250 voxels with a wide rim of air – you’d want to zoom in at the time of acquisition so that it makes a full use of the 512 x 512. This means that each CT scan actually represents different dimensions in real life even though they are all 512 x 512 x Z slices.

      This is why when we resample to isotropic 1 mm voxels, they all end up being different sizes.

      Think of each examination as having a fixed millimeter-per-voxel conversion factor which is based on patient size and different from exam to exam. IIf we made all the XxYxZ the same for all exams, then the voxel size can no longer be 1 x 1 x 1 mm, and vice versa.

      • Hi Howard,

        Thanks for your fast response, is much clear now what is happening. My aim is to be able to feed a 3d Neural network with the volumes, and this requires that all of them have the same shape. I understand the trade off you mention in the last paragraph, however, is there a transformation you could suggest to be able to get the images in the shape we want?

        Thanks in advance,

        Diana

        • I see – you may find success in using the resampling method from the blog post and simply hardcoding the new_shape variable (instead of calculating it). This should force all the resampled data to be in the same dimensions.

  • Hi Howard, I’m Minwoo.
    Thank you for this tutorial. I try to do your segmentation tutorial. But I have some problem of your tutorials.
    This problem is that some CT slices don’t make final mask or just one lung mask. Erosion and and dilation process is ok. Then color labels process also is ok. But some CT slices don’t show final mask.
    I’ll waiting for your response.
    Thank you.

    • Hi Minwoo – hope you are making good stride in your DICOM manipulation work!

      The “prop.bbox” code (starting for prop in regions:) in the make_lungmask function is the place you want to take a closer look.

      Essentially the code draws boxes around each of the labeled regions (B = prop.bbox). Then because boxes that represent lungs are more likely to be shaped a certain way than boxes that represent other labels, you can perform the mathematic to determine which label is most likely the lung. In the code, B = (min_row, min_col, max_row, max_col) in absolute pixel locations.

      Therefore, B[2]-B[0] would represent the height of the box that has been drawn.

      Try this reference to understand how bbox works. http://scikit-image.org/docs/dev/api/skimage.measure.html

      The precise numbers are determined empirically, so to get the right masks you may have to try different numbers.

      For instance, if your patients tend to have smaller lungs, then you would adjust the code to get closer to the center of the DICOM image.

      • Thank you Howard.
        I have solved my problem.
        If I have some problem again, See you soon.

  • Eric Voots
    4 weeks ago

    Hi,

    If we loop through all of the images and process them. What is the best way to load them all to be put into a format to analyze the saved .npy files using PCA, neural network etc..? I.e. not sure how to put all .npy files together for analysis.

  • Hey Eric, npy is a good choice for this, and I would go with a numpy.ndarray so you can have a 3D array. Popular packages like sklearn, tensorflow, and Keras all support Numpy variables.

  • Can you support 3D Plotting using vtk?
    Thank so much.

Trackbacks & Pings

Leave a Reply

Your email address will not be published. Required fields are marked *