DICOM: Importowanie danych obrazów medycznych do NumPy za pomocą PyDICOM

 

Czy da się manipulować pojedynczymi warstwami skanów w języku Python? Tak, jest to możliwe! W 2021 roku udało mi się napisać krótki program, który poprzez ręczną zmianę parametrów możliwe jest przeskakiwania pomiędzy poszczególnymi warstwami i o tym będzie dzisiejszy wpis.

Czym są pliki DICOM?

Na samym początku powinniśmy zacząć od tego, czym jest DICOM. DICOM oznacza „Digital Imaging and Communications in Medicine“. Jest to norma międzynarodowa do diagnostyki obrazowej. Plik DICOM składa się z nagłówków oraz zbioru danych obrazu, zawartych w jednym pliku. Informacje dotyczące pacjenta są standaryzowane i identyczne na wszystkich nagłówkach DICOM. Do obszerniejszego wpisu na temat plików DICOM zapraszam Cię tutaj.

Dataset 

Do testów programu użyłam datasetu MyHead, który znajdziesz też na moim githubie: https://github.com/jolapodolszanska/dicom-slicer, ale można je też pobrać stąd. Pomimo faktu, że łatwo jest uzyskać zbiory danych DICOM, większość baz danych zabrania ich redystrybucji przez osoby trzecie. 

Przygotowanie środowiska

Możesz do tego celu użyć dowolnego kompilatora, np. Visual Studio Code, Spyder, PyCharm, Atom itd. Wymagane jest też zainstalowanie Pythona, bo w tym języku będziemy pracować. 

Skoro pracować będziemy na plikach DICOM, to musimy zainstalować bibliotekę, która obsłuży nam tego typu pliki, zatem instalujemy PyDICOM za pomocą pip.

pip install pydicom

Ten pakiet importujemy pod nazwą dicom oraz resztę importów.

import pydicom
import os
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans

Następnie używamy os.path.walk, aby przejść przez katalog MyHead i zebrać wszystkie pliki .dcm na liście o nazwie lstFilesDCM:

PathDicom = "./MyHead dataset/"
lastFilesDCM = []
for dirName, subdirList, fileList in os.walk(PathDicom):
    for filename in fileList:
        if ".dcm" in filename.lower():
            lastFilesDCM.append(os.path.join(dirName,filename))

W PathDicom podajemy ścieżkę do katalogu, gdzie znajdują się pliki z rozszerzeniem .dcm. Jeśli sprawdzisz katalog MyHead, zobaczysz, że pliki .dcm mają nazwy vhf.1501.dcm, vhf.1502.dcm, vhf.1503.dcm itd. Dlatego funkcja przechodząca po kolei przez wszystkie pliki w katalogu zwróci je w odpowiedniej kolejności, ponieważ są one sortowane leksykograficznie przez system operacyjny. Jednak w wielu przypadkach pliki DICOM nie mają tych wszystkich początkowych zer. To spowodowałoby, że lastFilesDCM miałby nazwy plików uporządkowane w taki sposób, jak CT1, CT10, CT100 itp. Ponieważ typowe funkcje sortowania w Pythonie, takie jak sorted i metoda sortowania obiektów listy, są również leksykograficzne, dlatego zawsze jest lepiej użyć pakietu natsort, który można znaleźć na PyPI, instalując za pomocą prostego ).

pip install natsort

Przejdźmy do dalszej części kodu. Godnym uwagi aspektem tego pakietu jest to, że po odczytaniu pliku DICOM tworzy on obiekt dicom.dataset.FileDataset, w którym różne metadane są przypisywane do atrybutów obiektu o tej samej nazwie.

dicomLbl = pydicom.read_file(lastFilesDCM[0])

ConstPixelDims = (int(dicomLbl.Rows), int(dicomLbl.Columns), len(lastFilesDCM))

ConstPixelSpacing = (float(dicomLbl.PixelSpacing[0]), float(dicomLbl.PixelSpacing[1]), float(dicomLbl.SliceThickness))

Na początku ładujemy pierwszy plik DICOM, którego użyjemy jako odniesienia o nazwie dicomLbl, aby wyodrębnić metadane i którego nazwa pliku jest pierwsza na liście lastFilesDCM. Następnie obliczamy całkowite wymiary tablicy 3D NumPy, które są równe (Liczba wierszy pikseli w przekroju) x (Liczba kolumn pikseli w przekroju) x (Liczba przekrojów) wzdłuż osi x i y. Na koniec używamy atrybutów PixelSpacing i SliceThickness do obliczenia odstępów między pikselami w trzech osiach. Wymiary tablicy przechowujemy w ConstPixelDims, a odstępy w ConstPixelSpacing.

x = np.arange(0.0, (ConstPixelDims[0]+1)*ConstPixelSpacing[0], ConstPixelSpacing[0])
y = np.arange(0.0, (ConstPixelDims[1]+1)*ConstPixelSpacing[1], ConstPixelSpacing[1])

Aby otworzyć pliki z rozszerzeniem .dcm wystarczy oprogramowanie microDICOM (tak jak we wpisie, w którym odesłałam Cię do na początku). Jak sprawdzisz metadane, zobaczysz, że Rows, Columns, PixelSpacing i SliceThickness to wszystkie wpisy metadanych. W pydicom po prostu tworzy atrybuty o tych samych nazwach i przypisuje im odpowiednie wartości, dzięki czemu są one łatwo dostępne. Zatem używamy numpy.arange, ConstPixelDims i ConstPixelSpacing do obliczenia osi dla tej tablicy. Następnie pojawia się ostatnia część pydicom:

DicomArray = np.zeros(ConstPixelDims, dtype=dicomLbl.pixel_array.dtype)

print ("Generating image, please wait...")  

for filenameDCM in lastFilesDCM:
    ds = pydicom.read_file(filenameDCM)
    DicomArray[:, :, lastFilesDCM.index(filenameDCM)] = ds.pixel_array  

Na samym końcu używamy stosownych printów do wyświetlania parametrów.

print ("Done!")  
plt.figure(dpi=256)
plt.axes().set_aspect('equal', 'datalim')
plt.pcolormesh(x, y, np.flipud(DicomArray[:, :, 170]))

Zmieniając parametry w 

DicomArray[:, :, 170])

jesteśmy w stanie manipulować widokiem skanów. Przykład poniżej. 

Jak widać, najpierw tworzymy tablicę NumPy o nazwie DicomArray o wymiarach określonych w obliczonym wcześniej ConstPixelDims. Typ dtype tej tablicy jest taki sam jak typ dtype tablicy pixel_array referencyjnego zbioru danych dicomLbl, którego pierwotnie użyliśmy do wyodrębnienia metadanych. Interesujące jest to, że obiekt pixel_array jest czystą tablicą NumPy zawierającą dane pikseli dla konkretnego wycinka/obrazu DICOM. Następnie przechodzimy w pętli przez zebrane nazwy plików DICOM i używamy funkcji dicom.read_file, aby odczytać każdy plik do obiektu dicom.dataset.FileDataset. 

I to wszystko! Mam nadzieję, że to będzie dla Ciebie pierwsze starcie z DICOM. Cały kod jest dostępny na repozytorium: https://github.com/jolapodolszanska/dicom-slicer.