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.

42 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
    2 years 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.

    • Hi khiem, as you mentioned, VTK does support 3D plotting, and does a very good job at it.

      However, there is no easy way for me to show it on this blog because Jupyter does not directly support VTK, making it difficult to share the outputs.

      • Hi Howard Chen,
        Thank for reply, can you share code in here i will run it.

        • Hey khiem – it’s a little out of scope for this blog post, so maybe in the future I’ll write up something in full about VTK. It’s on my list of things to explore for web-based outputs.

        • You can use 3D Slicer, it support VTK and Python, 3D visualization is more simple than using 3D plotting.

  • Wogayehu
    1 year ago

    HI
    thank you for the tutorials since am new for deep neural networks i got clear idea about preprocessing and segmentation but does it mean we can feed the result of segmented region of lung to any deep Learning model(deepCNN, DBN, SVM) so on or only for deep CNN?

    • Hey Wogayehu,

      It should be useful for any number of the learning algorithms. The concept of doing segmentation and preprocessing in radiology is to standardize and focus on only the portion of the images you really care about (not always the right thing to do, but it often is). Any learning algorithm can potentially benefit from doing this.

      If you have a background in other learning algorithms like SVM and have used it for statistical learning with standard datasets, you may recall that data preprocessing, normalization, and filtering is often a good thing to do beforehand. Conceptually this may be though of as the imaging equivalent of that.

  • Wogayehu
    1 year ago

    thank you very much Howard for reply
    do you give me python pseudo-code resources sites for DBNs in image recognition specially for medical images?

  • Kassahun G
    1 year ago

    Dear Howard,
    Greetings!! Thank you for sharing such a nice tutorial. I would be happy if I can get the PDF version.
    Kassahun

    • Hi,

      I’m glad you’ve found it helpful! We don’t typically publish PDF versions of blog posts. Maybe you can try printing the page from your web browser to a PDF file.

      Best,
      Howard

  • Youssef Emad
    12 months ago

    Hey could u please where is the cancer in slide 97 and 112 .. i can’t see it

    • Hi Youssef,

      First let’s take at look at the right-sided lung (that’s actually the patient’s LEFT lung, but it’s just the way CT is displayed in America by convention).

      The cancer is not just on slice 97 and 112, it’s on slices from 97 through 112 (all the slices in between). Remember lung cancer is a 3D object so you should expect to see it on multiple slices. It is best seen on slice 100 as a cloud-looking round thing in the lung. Hope this helps!

  • Hi Chen,

    This is indeed a very useful tutorial. Could you please let me know whether, the images ready to hit the CNN are data saved by

    np.save(output_path + “maskedimages_%d.npy” % (id), imgs)

  • Hi Howard and thank you for this great tutorial,

    I’m working on an automated segmentation and 3D surface reconstruction script for
    cancellous and cortical bone. My problem is that as this is for potential total hip replacement
    cases, it usually involves patients with severe osteoarthritis. Segmentation using thresholding
    results in incomplete bone contours, which in turn result in holes on my reconstructed 3D model. Do you have any suggestions on how I should go about to tackle this issue?

    Many thanks!

  • captcoma
    11 months ago

    Hi Howard,

    Thank you very much for your great tutorial. I am new to Python and I do not know how to implement the “Import Packages” section.
    Do I have to run the code as a script?

    Many thanks for your help

    • import is Python syntax that includes packages. In short, it’s similar to import statements in Java and other languages. When you’re in Jupyter, the notebook will automatically execute your Python code without your having to save it separately as a script. Hope this helps.

  • You can use 3D Slicer, it support VTK and Python, 3D visualization is more simple than using 3D plotting.

  • Amazing tutorial! Thank you for this!

  • Interesting

  • Hi Howard,

    I‘m studying segmentation of MR brain-tissue(including segment White Matter、Gray Matter、Cerebrospinal fluid from brain-tissue),and I want to use support vertor machine to segment, I have got the feature vector from pixel, but i don’t konw how to get the labels, because SVM need the labels to complete, In other words, I don’t konw how to extract the labels of each tissue from DICOM files.Many thanks for your help!!!

    • Hi Jeffery – do you mean the ground truth labels for white, gray, and CSF for each pixel? (FYI – it’s technically called a voxel because the pixels have 3 dimensions).

      Unfortunately, radiology tissue type labels don’t come pre-labeled in the DICOM normally, which is why segmentation AI is valuable. Establishing ground truths typically require a human expert who hand-draw regions on each slice. This information is then compiled into a data format the script can read. Then, with the ground truth labels for every image, you can separate your input files into training, validating, and testing sets.

      As you can imagine, creating annotated ground truth is laborious, making annotated datasets very valuable. If you are open to using publically available datasets, take a look at MRBrainS by searching for it on Google.

      Best of luck,
      Howard

  • Hi Sir,

    Great tutorial Helped a lot, can you please also help how to use convolution neural network to classify stages of lung cancer and increase accuracy….
    Please mail similar kind of tutorial to train the data and classify stages..

    kazisajid023@gmail.com

  • Hi,
    I am trying to run this code using python 3.6.4 (installed using anaconda) but there is no compatible version of pydicom for 3.6.4.
    I have also tried it with Python 2.7 but then I run into errors while installing sci-kit.

    What should I do?

    • Hi Areeb,
      Anaconda allows you to install a different version of Python. See here. In short, you would just create a new conda environment, like this:

      conda create --name raddq python=3.5

      Then depending on your operating system you can activate it accordingly via either “activate” (Win) or “source activate” (Linux/Mac) commands.

  • Hunar A.Ahmed
    4 months ago

    I have a folder contains 5 Dicom images:
    Shape before resampling (5, 512, 512)
    Shape after resampling (175, 340, 340)
    I know that before resampling number of images is 5 and each image 512×512 (height x width), but after resampling, it showing me 175, is this mean the number of images is now 175 and each image is 340×340 (height x width)? or anything else, I am confusing with that?

    • Yes you probably have 175 resampled slices. If you did isometric resampling, it probably means the distance between your first and last slice was pretty big, so the algorithm tried to fill in the distance in between. Unfortunately since it only had 5 source slices my guess is your resampled images might have some quality issues.

      • Hunar A.Ahmed
        4 months ago

        thank you for your replay Mr.Howard, in your replay,
        you said ‘distance between your first and last slice was pretty big’, is the slice thickness distance between first and last slice or distance between two slices? also, you said ‘algorithm tried to fill in the distance in between’ what you mean by this?
        I use these 5 images folder for test only because I have a low computing power Pc, I have the complete folder image with 133 slices (from LIDC-IDIR) when displaying slice thickness with 5 folder images it shows 30 mm but when I use the 133 folder image it shows 2.5 mm, please can you explain these for me? my resampling code is same as your code.

        • Hunar A.Ahmed
          4 months ago

          sorry, when displaying slice thickness with 5 folder images it shows 35 mm not 30 mm.

  • Bhakti Raul Palkar
    3 months ago

    Why my 3d graph shows everything in white?

  • Thanks for the detailed tutorial. I can no longer access the image set from Kaggle website. Is there an alternate location from where I can download the data?

Trackbacks & Pings

Leave a Reply

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