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
type(dcm)
dcm.pixels
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
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
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
plt.hist(t_bin.numpy(), bins=t_bin, color='c'); plt.plot(t_bin, torch.linspace(0,1,len(t_bin)));
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);
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 brightnessw= window width or range aka contrast
plt.imshow(dcm.windowed(*dicom_windows.brain), cmap=plt.cm.bone);
tensor_ct = TensorCTScan(dcm.pixel_array)
tensor_ct.show();
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.
dcm.show()
dcm.pct_in_window(*dicom_windows.brain)
pct_in_window can be used to check what percentage of the image is composed of meaningful pixels (pixels within the specified window)
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'))
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
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');