In [ ]:
import matplotlib.pyplot as plt
%matplotlib inline
import SimpleITK as sitk
# Download data to work on
from downloaddata import fetch_data as fdata
SimpleITK has a build in Show method which saves the image to disk and launches a user configurable program ( defaults to ImageJ ), to display the image.
In [ ]:
img1 = sitk.ReadImage(fdata("cthead1.png"))
sitk.Show(img1, title="cthead1")
In [ ]:
img2 = sitk.ReadImage(fdata("VM1111Shrink-RGB.png"))
sitk.Show(img2, title="Visible Human Head")
In [ ]:
nda = sitk.GetArrayFromImage(img1)
plt.imshow(nda)
In [ ]:
nda = sitk.GetArrayFromImage(img2)
ax = plt.imshow(nda)
In [ ]:
def myshow(img):
nda = sitk.GetArrayFromImage(img)
plt.imshow(nda)
In [ ]:
myshow(img2)
In [ ]:
myshow(sitk.Expand(img2, [10]*5))
This image does not appear not bigger.
There are numerous improvements that we can make:
In [ ]:
def myshow(img, title=None, margin=0.05, dpi=80):
nda = sitk.GetArrayFromImage(img)
spacing = img.GetSpacing()
if nda.ndim == 3:
# fastest dim, either component or x
c = nda.shape[-1]
# the the number of components is 3 or 4 consider it an RGB image
if not c in (3,4):
nda = nda[nda.shape[0]//2,:,:]
elif nda.ndim == 4:
c = nda.shape[-1]
if not c in (3,4):
raise Runtime("Unable to show 3D-vector Image")
# take a z-slice
nda = nda[nda.shape[0]//2,:,:,:]
ysize = nda.shape[0]
xsize = nda.shape[1]
# Make a figure big enough to accomodate an axis of xpixels by ypixels
# as well as the ticklabels, etc...
figsize = (1 + margin) * ysize / dpi, (1 + margin) * xsize / dpi
fig = plt.figure(figsize=figsize, dpi=dpi)
# Make the axis the right size...
ax = fig.add_axes([margin, margin, 1 - 2*margin, 1 - 2*margin])
extent = (0, xsize*spacing[1], ysize*spacing[0], 0)
t = ax.imshow(nda,extent=extent,interpolation=None)
if nda.ndim == 2:
t.set_cmap("gray")
if(title):
plt.title(title)
In [ ]:
myshow(sitk.Expand(img2,[2,2]), title="Big Visibile Human Head")
In [ ]:
img1_seg = sitk.ReadImage(fdata("2th_cthead1.png"))
myshow(img1_seg, "Label Image as Grayscale")
In [ ]:
myshow(sitk.LabelToRGB(img1_seg), title="Label Image as RGB")
Most filters which take multiple images as arguments require that the images occupy the same physical space. That is the pixel you are operating must refer to the same location.
In [ ]:
# This fails because the images don't occupy the same physical scale.
myshow(sitk.LabelOverlay(img1, img1_seg), title="Label Overlayed")
Ideally, this meta-data is retained, so images are consistent. However, the situation is not ideal, and image meta-data needs to be fixed.
In [ ]:
img1_seg.CopyInformation(img1)
In [ ]:
myshow(sitk.LabelOverlay(img1, sitk.LabelContour(img1_seg), 1.0))
Now lets move on to visualizing real MRI images with segmentations. The Surgical Planning Laboratory at Brigham and Women's Hospital has a wonderful Multi-modality MRI-based Atlas of the Brain that we can use.
Please note, what is done here is for convenients and is not the correct whay to display images for a radiological work.
In [ ]:
img_T1 = sitk.ReadImage(fdata("nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT1.nrrd"))
img_T2 = sitk.ReadImage(fdata("nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT2.nrrd"))
img_labels = sitk.ReadImage(fdata("nac-hncma-atlas2013-Slicer4Version/Data/hncma-atlas.nrrd"))
In [ ]:
myshow(img_T1)
myshow(img_T2)
myshow(sitk.LabelToRGB(img_labels))
In [ ]:
size = img_T1.GetSize()
myshow(img_T1[:,size[1]//2,:])
In [ ]:
slices =[img_T1[size[0]//2,:,:], img_T1[:,size[1]//2,:], img_T1[:,:,size[2]//2]]
myshow(sitk.Tile(slices, [3,1]), dpi=20)
In [ ]:
nslices = 5
slices = [ img_T1[:,:,s] for s in range(0, size[2], size[0]//(nslices+1))]
myshow(sitk.Tile(slices, [1,0]))
Let's create a version of the show methods which allows the selection of slices to be displayed.
In [ ]:
def myshow3d(img, xslices=[], yslices=[], zslices=[], title=None, margin=0.05, dpi=80):
size = img.GetSize()
img_xslices = [img[s,:,:] for s in xslices]
img_yslices = [img[:,s,:] for s in yslices]
img_zslices = [img[:,:,s] for s in zslices]
maxlen = max(len(img_xslices), len(img_yslices), len(img_zslices))
img_null = sitk.Image([0,0], img.GetPixelIDValue(), img.GetNumberOfComponentsPerPixel())
img_slices = []
d = 0
if len(img_xslices):
img_slices += img_xslices + [img_null]*(maxlen-len(img_xslices))
d += 1
if len(img_yslices):
img_slices += img_yslices + [img_null]*(maxlen-len(img_yslices))
d += 1
if len(img_zslices):
img_slices += img_zslices + [img_null]*(maxlen-len(img_zslices))
d +=1
if maxlen != 0:
if img.GetNumberOfComponentsPerPixel() == 1:
img = sitk.Tile(img_slices, [maxlen,d])
#TODO check in code to get Tile Filter working with VectorImages
else:
img_comps = []
for i in range(0,img.GetNumberOfComponentsPerPixel()):
img_slices_c = [sitk.VectorIndexSelectionCast(s, i) for s in img_slices]
img_comps.append(sitk.Tile(img_slices_c, [maxlen,d]))
img = sitk.Compose(img_comps)
myshow(img, title, margin, dpi)
In [ ]:
myshow3d(img_T1,yslices=range(50,size[1]-50,20), zslices=range(50,size[2]-50,20), dpi=30)
In [ ]:
myshow3d(img_T2,yslices=range(50,size[1]-50,30), zslices=range(50,size[2]-50,20), dpi=30)
In [ ]:
myshow3d(sitk.LabelToRGB(img_labels),yslices=range(50,size[1]-50,20), zslices=range(50,size[2]-50,20), dpi=30)
In [ ]:
# Again the following line fails due to difference in the image location
myshow3d(sitk.LabelOverlay(img_T1,img_labels),yslices=range(50,size[1]-50,20), zslices=range(50,size[2]-50,20), dpi=30)
Why are the following lines needed?
What do they do?
In [ ]:
img_labels.CopyInformation(img_T1)
In [ ]:
img_T1 = sitk.Cast(sitk.RescaleIntensity(img_T1), sitk.sitkUInt8)
This myshow and myshow3d function is really useful. It has been copied into a "myshow.py" file so that it can be imported into other notebooks.