Medical Imaging

Helpers for working with DICOM files

Patching

---

source

get_dicom_files

def get_dicom_files(
    path, recurse:bool=True, folders:NoneType=None
):

Get dicom files in path recursively, only in folders, if specified.

---

source

Path.dcmread

def dcmread(
    fn:Path, force:bool=False
):

Open a DICOM file

fastai.medical.imaging uses pydicom.dcmread to read a DICOM file. To view the header of a DICOM, specify the path of a test file and call dcmread.

TEST_DCM = Path('images/sample.dcm')
dcm = TEST_DCM.dcmread()
dcm

Dataset.file_meta -------------------------------
(0002, 0000) File Meta Information Group Length  UL: 176
(0002, 0001) File Meta Information Version       OB: b'\x00\x01'
(0002, 0002) Media Storage SOP Class UID         UI: CT Image Storage
(0002, 0003) Media Storage SOP Instance UID      UI: 9999.180975792154576730321054399332994563536
(0002, 0010) Transfer Syntax UID                 UI: Explicit VR Little Endian
(0002, 0012) Implementation Class UID            UI: 1.2.40.0.13.1.1.1
(0002, 0013) Implementation Version Name         SH: 'dcm4che-1.4.38'
-------------------------------------------------
(0008, 0018) SOP Instance UID                    UI: ID_e0cc6a4b5
(0008, 0060) Modality                            CS: 'CT'
(0010, 0020) Patient ID                          LO: 'ID_a107dd7f'
(0020, 000d) Study Instance UID                  UI: ID_6468bdd34a
(0020, 000e) Series Instance UID                 UI: ID_4be303ae64
(0020, 0010) Study ID                            SH: ''
(0020, 0032) Image Position (Patient)            DS: [-125.000, -122.268, 115.936]
(0020, 0037) Image Orientation (Patient)         DS: [1.000000, 0.000000, 0.000000, 0.000000, 0.978148, -0.207912]
(0028, 0002) Samples per Pixel                   US: 1
(0028, 0004) Photometric Interpretation          CS: 'MONOCHROME2'
(0028, 0010) Rows                                US: 256
(0028, 0011) Columns                             US: 256
(0028, 0030) Pixel Spacing                       DS: [0.488281, 0.488281]
(0028, 0100) Bits Allocated                      US: 16
(0028, 0101) Bits Stored                         US: 16
(0028, 0102) High Bit                            US: 15
(0028, 0103) Pixel Representation                US: 1
(0028, 1050) Window Center                       DS: "40.0"
(0028, 1051) Window Width                        DS: "100.0"
(0028, 1052) Rescale Intercept                   DS: "-1024.0"
(0028, 1053) Rescale Slope                       DS: "1.0"
(7fe0, 0010) Pixel Data                          OW: Array of 131072 elements

type(dcm)

pydicom.dataset.FileDataset

---

source

TensorDicom

def TensorDicom(
    args:VAR_POSITIONAL, kwargs:VAR_KEYWORD
):

Inherits from TensorImage and converts the pixel_array into a TensorDicom

---

source

PILDicom

def PILDicom(
    
)->None:

Base class for a Pillow Image that can show itself and convert to a Tensor

---

source

Path.png16read

def png16read(
    
):

---

Dataset.pixels

def pixels(
    
):

pixel_array as a tensor

dcm.pixels

tensor([[-1024., -1024., -1024.,  ..., -1024., -1024., -1024.],
        [-1024., -1024., -1024.,  ..., -1024., -1024., -1024.],
        [-1024., -1024., -1024.,  ..., -1024., -1024., -1024.],
        ...,
        [-1024., -1024., -1024.,  ..., -1024., -1024., -1024.],
        [-1024., -1024., -1024.,  ..., -1024., -1024., -1024.],
        [-1024., -1024., -1024.,  ..., -1024., -1024., -1024.]])

---

Dataset.scaled_px

def scaled_px(
    
):

pixels scaled by RescaleSlope and RescaleIntercept

scaled_px uses RescaleSlope and RescaleIntercept values to correctly scale the image so that they represent the correct tissue densities. You can observe what scaled_px does by viewing the the pixel distribution of a dicom image. The histogram below displays the current pixel distribution which shows a pixel range between -1133 and 2545.

plt.hist(dcm.pixels.flatten().numpy());

As shown in the header of the test image the RescaleIntercept has a value of -1024.0 and a RescaleSlope value of 1.0. scaled_px will scale the pixels by these values.

plt.hist(dcm.scaled_px.flatten().numpy());

The pixel distibution is now between -2157 and 1521

---

source

array_freqhist_bins

def array_freqhist_bins(
    n_bins:int=100
):

A numpy based function to split the range of pixel values into groups, such that each group has around the same number of pixels

---

source

Tensor.freqhist_bins

def freqhist_bins(
    n_bins:int=100
):

A function to split the range of pixel values into groups, such that each group has around the same number of pixels

For example with n_bins set to 1 this means the bins will be split into 3 distinct bins (the beginning, the end and the number of bins specified by n_bins.

t_bin = dcm.pixels.freqhist_bins(n_bins=1)
t_bin

tensor([-1076.,    40.,  2375.])

plt.hist(t_bin.numpy(), bins=t_bin, color='c')
plt.plot(t_bin, torch.linspace(0,1,len(t_bin)));

with n_bins at 100

t_bin = dcm.pixels.freqhist_bins(n_bins=100)
t_bin

tensor([-1076., -1026., -1024., -1021.,    28.,    30.,    31.,    32.,    33.,
           34.,    35.,    36.,    37.,    38.,    39.,    40.,    41.,    42.,
           44.,    48.,    52.,    58.,    66.,    72.,    76.,    80.,    85.,
           91.,    94.,    98.,   103.,   111.,   123.,   161.,   219.,   478.,
          829.,   999.,  1027.,  1038.,  1044.,  1047.,  1049.,  1050.,  1051.,
         1052.,  1053.,  1054.,  1055.,  1056.,  1057.,  1058.,  1059.,  1060.,
         1062.,  1066.,  1108.,  1265.,  1453.,  1616.,  1741.,  1838.,  1943.,
         2051.,  2220.,  2375.])

plt.hist(t_bin.numpy(), bins=t_bin, color='c'); plt.plot(t_bin, torch.linspace(0,1,len(t_bin)));

---

source

Tensor.hist_scaled_pt

def hist_scaled_pt(
    brks:NoneType=None
):

---

source

Tensor.hist_scaled

def hist_scaled(
    brks:NoneType=None
):

Scales a tensor using freqhist_bins to values between 0 and 1

The test image has pixel values that range between -1000 and 2500

plt.hist(dcm.pixels.flatten().numpy(), bins=100);

hist_scaled provides a way of scaling the input pixel values to between 0 and 1

tensor_hists = dcm.pixels.hist_scaled()
plt.hist(tensor_hists.flatten().numpy(), bins=100);

---

Dataset.hist_scaled

def hist_scaled(
    brks:NoneType=None, min_px:NoneType=None, max_px:NoneType=None
):

Pixels scaled to a min_px and max_px value

data_scaled = dcm.hist_scaled()
plt.imshow(data_scaled, cmap=plt.cm.bone);

data_scaled = dcm.hist_scaled(min_px=100, max_px=1000)
plt.imshow(data_scaled, cmap=plt.cm.bone);

Dicom images can contain a high amount of voxel values and windowing can be thought of as a means of manipulating these values in order to change the apperance of the image so particular structures are highlighted. A window has 2 values:

• l = window level or center aka brightness
• w = window width or range aka contrast

---

source

Tensor.windowed

def windowed(
    w, l
):

Scale pixel intensity by window width and window level

---

Dataset.windowed

def windowed(
    w, l
):

plt.imshow(dcm.windowed(*dicom_windows.brain), cmap=plt.cm.bone);

---

source

TensorCTScan

def TensorCTScan(
    args:VAR_POSITIONAL, kwargs:VAR_KEYWORD
):

Inherits from TensorImageBW and converts the pixel_array into a TensorCTScan

tensor_ct = TensorCTScan(dcm.pixel_array)
tensor_ct.show();

---

source

PILCTScan

def PILCTScan(
    
)->None:

Base class for a Pillow Image that can show itself and convert to a Tensor

---

Dataset.show

def show(
    scale:bool=True,
    cmap:LinearSegmentedColormap=<matplotlib.colors.LinearSegmentedColormap object at 0x7f363f3f2cc0>,
    min_px:int=-1100, max_px:NoneType=None, ax:NoneType=None, figsize:NoneType=None, title:NoneType=None,
    ctx:NoneType=None, norm:NoneType=None, aspect:NoneType=None, interpolation:NoneType=None, alpha:NoneType=None,
    vmin:NoneType=None, vmax:NoneType=None, colorizer:NoneType=None, origin:NoneType=None, extent:NoneType=None,
    interpolation_stage:NoneType=None, filternorm:bool=True, filterrad:float=4.0, resample:NoneType=None,
    url:NoneType=None, data:NoneType=None
):

Display a normalized dicom image by default

scales = False, True, dicom_windows.brain, dicom_windows.subdural
titles = 'raw','normalized','brain windowed','subdural windowed'
for s,a,t in zip(scales, subplots(2,2,imsize=4)[1].flat, titles):
    dcm.show(scale=s, ax=a, title=t)

dcm.show(cmap=plt.cm.gist_ncar, figsize=(6,6))

Some dicom datasets such as the The Thyroid Segmentation in Ultrasonography Dataset is a dataset where each image has multiple frames per file (hundreds in this case). By default the show function will display 1 frame but if the dataset has multiple frames you can specify the number of frames to view.

---

Dataset.show

def show(
    frames:int=1, scale:bool=True,
    cmap:LinearSegmentedColormap=<matplotlib.colors.LinearSegmentedColormap object at 0x7f363f3f2cc0>,
    min_px:int=-1100, max_px:NoneType=None, kwargs:VAR_KEYWORD
):

Adds functionality to view dicom images where each file may have more than 1 frame

dcm.show()

---

Dataset.pct_in_window

def pct_in_window(
    dcm:DcmDataset, w, l
):

% of pixels in the window (w,l)

dcm.pct_in_window(*dicom_windows.brain)

0.19049072265625

pct_in_window can be used to check what percentage of the image is composed of meaningful pixels (pixels within the specified window)

---

source

uniform_blur2d

def uniform_blur2d(
    x, s
):

Uniformly apply blurring

ims = dcm.hist_scaled(), uniform_blur2d(dcm.hist_scaled(), 20), uniform_blur2d(dcm.hist_scaled(), 50)
show_images(ims, titles=('original', 'blurred 20', 'blurred 50'))

---

source

gauss_blur2d

def gauss_blur2d(
    x, s
):

Apply gaussian_blur2d kornia filter

ims = dcm.hist_scaled(), gauss_blur2d(dcm.hist_scaled(), 20), gauss_blur2d(dcm.hist_scaled(), 50)
show_images(ims, titles=('original', 'gauss_blur 20', 'gauss_blur 50'))

Images are often affected by random variations in intensity values, called noise. Gaussian noise contains variatons in intensity that are drawn from a Gaussian or normal distribution. A Guassian filter is usually used to blur edges and remove smaller or thinner areas in order to preserve the most important information

---

source

Tensor.mask_from_blur

def mask_from_blur(
    x:Tensor, window, sigma:float=0.3, thresh:float=0.05, remove_max:bool=True
):

Create a mask from the blurred image

---

Dataset.mask_from_blur

def mask_from_blur(
    x:DcmDataset, window, sigma:float=0.3, thresh:float=0.05, remove_max:bool=True
):

Create a mask from the blurred image

mask = dcm.mask_from_blur(dicom_windows.brain, sigma=0.9, thresh=0.1, remove_max=True)
wind = dcm.windowed(*dicom_windows.brain)

_,ax = subplots(1,3)
show_image(wind, ax=ax[0], title='window')
show_image(mask, alpha=0.5, cmap=plt.cm.Reds, ax=ax[1], title='mask')
show_image(wind, ax=ax[2])
show_image(mask, alpha=0.5, cmap=plt.cm.Reds, ax=ax[2], title='window and mask');

---

source

mask2bbox

def mask2bbox(
    mask
):

bbs = mask2bbox(mask)
lo,hi = bbs
show_image(wind[lo[0]:hi[0],lo[1]:hi[1]]);

---

source

crop_resize

def crop_resize(
    x, crops, new_sz
):

px256 = crop_resize(to_device(wind[None]), bbs[...,None], 128)[0]
show_image(px256)
px256.shape

torch.Size([1, 128, 128])

Comparing the original image with the image from using the mask and crop_resize function

_,axs = subplots(1,2)
dcm.show(ax=axs[0])
show_image(px256, ax=axs[1]);

---

source

Tensor.to_nchan

def to_nchan(
    x:Tensor, wins, bins:NoneType=None
):

---

Dataset.to_nchan

def to_nchan(
    x:DcmDataset, wins, bins:NoneType=None
):

to_nchan takes a tensor or a dicom as the input and returns multiple one channel images (the first depending on the choosen windows and a normalized image). Setting bins to 0 only returns the windowed image.

show_images(dcm.to_nchan([dicom_windows.brain], bins=0))

show_images(dcm.to_nchan([dicom_windows.brain], bins=None))

---

source

Tensor.to_3chan

def to_3chan(
    x:Tensor, win1, win2, bins:NoneType=None
):

---

Dataset.to_3chan

def to_3chan(
    x:DcmDataset, win1, win2, bins:NoneType=None
):

show_images(dcm.to_nchan([dicom_windows.brain,dicom_windows.subdural,dicom_windows.abdomen_soft]))

---

Dataset.save_jpg

def save_jpg(
    x:Tensor | DcmDataset, path, wins, bins:NoneType=None, quality:int=90
):

Save tensor or dicom image into jpg format

---

source

Tensor.save_jpg

def save_jpg(
    x:Tensor | DcmDataset, path, wins, bins:NoneType=None, quality:int=90
):

Save tensor or dicom image into jpg format

---

Dataset.to_uint16

def to_uint16(
    x:Tensor | DcmDataset, bins:NoneType=None
):

Convert into a unit16 array

---

source

Tensor.to_uint16

def to_uint16(
    x:Tensor | DcmDataset, bins:NoneType=None
):

Convert into a unit16 array

---

Dataset.save_tif16

def save_tif16(
    x:Tensor | DcmDataset, path, bins:NoneType=None, compress:bool=True
):

Save tensor or dicom image into tiff format

---

source

Tensor.save_tif16

def save_tif16(
    x:Tensor | DcmDataset, path, bins:NoneType=None, compress:bool=True
):

Save tensor or dicom image into tiff format

_,axs=subplots(1,2)
with tempfile.TemporaryDirectory() as f:
    f = Path(f)
    dcm.save_jpg(f/'test.jpg', [dicom_windows.brain,dicom_windows.subdural])
    show_image(Image.open(f/'test.jpg'), ax=axs[0])
    dcm.save_tif16(f/'test.tif')
    show_image(Image.open(str(f/'test.tif')), ax=axs[1]);

---

Dataset.set_pixels

def set_pixels(
    px
):

---

Dataset.zoom

def zoom(
    ratio
):

Zoom image by specified ratio

Check to see the current size of the dicom image

dcm.pixel_array.shape

(256, 256)

dcm.zoom(7.0)
dcm.show(); dcm.pixel_array.shape

---

Dataset.zoom_to

def zoom_to(
    sz
):

Change image size to specified pixel size

dcm.zoom_to(200); dcm.pixel_array.shape

(200, 200)

---

Dataset.shape

def shape(
    
):

Returns the shape of a dicom image as rows and columns

dcm2 = TEST_DCM.dcmread()
dcm2.zoom_to(90)
test_eq(dcm2.shape, (90,90))

dcm2 = TEST_DCM.dcmread()
dcm2.zoom(0.25)
dcm2.show()

---

Dataset.as_dict

def as_dict(
    px_summ:bool=True, window:tuple=(80, 40)
):

Convert the header of a dicom into a dictionary

as_dict takes in 2 parameters: px_summ which by default is set to True and this returns additional stats such as minimal pixel value, maximum pixel value, the mean pixel value and the image standard deviation. The window parameter calculates the pct_in_window value depending on the window that is specified.

dcm.as_dict(px_summ=True, window=dicom_windows.brain);

Creating a dataframe of the values within the header of the dicom

pneumothorax_source = untar_data(URLs.SIIM_SMALL)
items = get_dicom_files(pneumothorax_source, recurse=True,  folders='train')
dicom_dataframe = pd.DataFrame.from_dicoms(items, window=dicom_windows.brain, px_summ=True)
dicom_dataframe.head(2).T.tail(5)

	0	1
img_min	0	0
img_max	254	250
img_mean	160.398	114.525
img_std	53.8549	70.7523
img_pct_window	0.0870295	0.326269

---

source

DicomSegmentationDataLoaders

def DicomSegmentationDataLoaders(
    loaders:VAR_POSITIONAL, # [`DataLoader`](https://docs.fast.ai/data.load.html#dataloader) objects to wrap
    path:str | Path='.', # Path to store export objects
    device:NoneType=None, # Device to put [`DataLoaders`](https://docs.fast.ai/data.core.html#dataloaders)
):

Basic wrapper around DICOM DataLoaders with factory methods for segmentation problems

path = untar_data(URLs.TCGA_SMALL)
codes = np.loadtxt(path/'codes.txt', dtype=str)
fnames = get_dicom_files(path/'dicoms')
label_func = lambda o: path/'labels'/f'{o.stem}.png'

dls = DicomSegmentationDataLoaders.from_label_func(path, fnames, label_func, codes=codes, bs=4)
dls.show_batch()