<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Using matplotlib to display inline images
In 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 SimpleITK images.
Step1:  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.
Step2:  This image does not appear not bigger.
There are numerous improvements that we can make
Step3:  Tips and Tricks for Visualizing Segmentations
Step4:  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.
Step5:  Ideally, this meta-data is retained, so images are consistent. However, the situation is not ideal, and image meta-data needs to be fixed.
Step6:  Tips and Tricks for 3D Image Visualization
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.
Step7:  Let's create a version of the show methods which allows the selection of slices to be displayed.
Step8:  Why are the following lines needed?
What do they do?
<ASSISTANT_TASK:>
Python Code:
import matplotlib.pyplot as plt
%matplotlib inline
import SimpleITK as sitk
# Download data to work on
from downloaddata import fetch_data as fdata
Explanation: Using matplotlib to display inline images
In 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 SimpleITK images.
End of explanation
img1 = sitk.ReadImage(fdata("cthead1.png"))
sitk.Show(img1, title="cthead1")
img2 = sitk.ReadImage(fdata("VM1111Shrink-RGB.png"))
sitk.Show(img2, title="Visible Human Head")
nda = sitk.GetArrayFromImage(img1)
plt.imshow(nda)
nda = sitk.GetArrayFromImage(img2)
ax = plt.imshow(nda)
def myshow(img):
    nda = sitk.GetArrayFromImage(img)
    plt.imshow(nda)
myshow(img2)
myshow(sitk.Expand(img2, [10]*5))
Explanation: 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.
End of explanation
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)
myshow(sitk.Expand(img2,[2,2]), title="Big Visibile Human Head")
Explanation: This image does not appear not bigger.
There are numerous improvements that we can make:
support 3d images
include a title
use physical pixel size for axis labels
show the image as gray values
End of explanation
img1_seg = sitk.ReadImage(fdata("2th_cthead1.png"))
myshow(img1_seg, "Label Image as Grayscale")
myshow(sitk.LabelToRGB(img1_seg), title="Label Image as RGB")
Explanation: Tips and Tricks for Visualizing Segmentations
End of explanation
# This fails because the images don't occupy the same physical scale.
myshow(sitk.LabelOverlay(img1, img1_seg), title="Label Overlayed")
Explanation: 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.
End of explanation
img1_seg.CopyInformation(img1)
myshow(sitk.LabelOverlay(img1, sitk.LabelContour(img1_seg), 1.0))
Explanation: Ideally, this meta-data is retained, so images are consistent. However, the situation is not ideal, and image meta-data needs to be fixed.
End of explanation
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"))
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]))
Explanation: Tips and Tricks for 3D Image Visualization
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.
End of explanation
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)
Explanation: Let's create a version of the show methods which allows the selection of slices to be displayed.
End of explanation
img_labels.CopyInformation(img_T1)
img_T1 = sitk.Cast(sitk.RescaleIntensity(img_T1), sitk.sitkUInt8)
Explanation: Why are the following lines needed?
What do they do?
End of explanation
<END_TASK>