In [4]:
nda = sitk.GetArrayFromImage(img1)
imshow(nda)
Out[4]:
matplotlib to display inline imagesIn this notebook we will explore using matplotlib to display images in our notebooks, and work towards developing a reusable function to display 2D,3D, color, and label overlays for our SimpleITK images.
import SimpleITK as sitk
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.
img1 = sitk.ReadImage("Data/cthead1.png")
sitk.Show(img1, title="cthead1")
img2 = sitk.ReadImage("Data/VM1111Shrink-RGB.png")
sitk.Show(img2, title="Visible Human Head")
nda = sitk.GetArrayFromImage(img1)
imshow(nda)
nda = sitk.GetArrayFromImage(img2)
ax = imshow(nda)
def myshow(img):
nda = sitk.GetArrayFromImage(img)
imshow(nda)
myshow(img2)
myshow(sitk.Expand(img2, [10]*5))
This image is not bigger as expected.
There are numerous improvements that we can make:
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 = 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)
myshow(sitk.Expand(img2,[2,2]), title="Big Visibile Human Head")
img1_seg = sitk.ReadImage("Data/2th_cthead1.png")
myshow(img1_seg, "Label Image as Grayscale")
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.
# This fails because the images don't occupy the same physical scale.
myshow(sitk.LabelOverlay(img1, img1_seg), title="Label Overlayed")
Anyone responsible doing image processing will maintain this metadata, so images are consistent. However, sometime things are not ideal and need to be fixed.
img1_seg.CopyInformation(img1)
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 wonder 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.
img_T1 = sitk.ReadImage("Data/nac-brain-atlas-1.0/volumes/A1_grayT1.nrrd")
img_T2 = sitk.ReadImage("Data/nac-brain-atlas-1.0/volumes/A1_grayT2.nrrd")
img_labels = sitk.ReadImage("Data/nac-brain-atlas-1.0/volumes/A1_labels.nrrd")
myshow(img_T1)
myshow(img_T2)
myshow(sitk.LabelToRGB(img_labels))
size = img_T1.GetSize()
myshow(img_T1[:,size[1]//2,:])
slices =[img_T1[size[0]//2,:,:], img_T1[:,size[1]//2,:], img_T1[:,:,size[2]//2]]
myshow(sitk.Tile(slices, [3,1]), dpi=20)
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.
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)
myshow3d(img_T1,yslices=range(50,size[1]-50,20), zslices=range(50,size[2]-50,20), dpi=30)
myshow3d(img_T2,yslices=range(50,size[1]-50,30), zslices=range(50,size[2]-50,20), dpi=30)
myshow3d(sitk.LabelToRGB(img_labels),yslices=range(50,size[1]-50,20), zslices=range(50,size[2]-50,20), dpi=30)
# 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?
img_labels.CopyInformation(img_T1)
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.