使用 matplotlib
显示内联图像
在此笔记本中,我们将探索使用 matplotlib
显示笔记本中的图像,并致力于开发可重复使用的函数来显示 SimpleITK 图像的 2D、3D、颜色和标签叠加层。
我们还将研究使用需要输入图像重叠的图像过滤器的微妙之处。
%matplotlib inline
import matplotlib.pyplot as plt
import SimpleITK as sitk# Download data to work on
%run update_path_to_download_script
from downloaddata import fetch_data as fdata
SimpleITK 有一个内置的“Show”方法,可以将图像保存到磁盘并启动用户可配置的程序(默认为 ImageJ)来显示图像。
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.GetArrayViewFromImage(img1)
plt.imshow(nda)
nda = sitk.GetArrayViewFromImage(img2)
ax = plt.imshow(nda)
def myshow(img):nda = sitk.GetArrayViewFromImage(img)plt.imshow(nda)
myshow(img2)
myshow(sitk.Expand(img2, [10] * 5))
此图像看起来并不大。
我们可以进行许多改进:
- 支持 3D 图像
- 包含标题
- 使用物理像素大小作为轴标签
- 将图像显示为灰度值
def myshow(img, title=None, margin=0.05, dpi=80):nda = sitk.GetArrayViewFromImage(img)spacing = img.GetSpacing()if nda.ndim == 3:# fastest dim, either component or xc = nda.shape[-1]# the the number of components is 3 or 4 consider it an RGB imageif 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-slicenda = nda[nda.shape[0] // 2, :, :, :]ysize = nda.shape[0]xsize = nda.shape[1]# Make a figure big enough to accommodate an axis of xpixels by ypixels# as well as the ticklabels, etc...figsize = (1 + margin) * ysize / dpi, (1 + margin) * xsize / dpifig = 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")
可视化分割的技巧和窍门
我们首先加载分割图像。由于分割只是具有积分数据的图像,因此我们可以像显示任何其他图像一样显示标签。
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")
大多数以多幅图像作为参数的过滤器都要求图像占据相同的物理空间。也就是说,您正在操作的像素必须指向相同的位置。幸运的是,我们的图像和标签确实占据相同的物理空间,这使我们能够将分割叠加到原始图像上。
myshow(sitk.LabelOverlay(img1, img1_seg), title="Label Overlayed")
我们还可以将标签叠加为轮廓。
myshow(sitk.LabelOverlay(img1, sitk.LabelContour(img1_seg), 1.0))
3D 图像可视化的技巧和窍门
现在让我们继续使用分割来可视化真实的 MRI 图像。布莱根妇女医院的外科规划实验室提供了一套出色的基于 MRI 的多模态大脑图谱,可供我们使用。
请注意,此处所做的操作是为了方便,并不是放射学工作中显示图像的常用方式。
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")
)
Fetching nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT1.nrrd
Fetching nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT2.nrrd
Fetching 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]))
让我们创建一个显示方法的版本,允许显示切片的选择。
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.GetPixelID(), img.GetNumberOfComponentsPerPixel())img_slices = []d = 0if len(img_xslices):img_slices += img_xslices + [img_null] * (maxlen - len(img_xslices))d += 1if len(img_yslices):img_slices += img_yslices + [img_null] * (maxlen - len(img_yslices))d += 1if len(img_zslices):img_slices += img_zslices + [img_null] * (maxlen - len(img_zslices))d += 1if maxlen != 0:if img.GetNumberOfComponentsPerPixel() == 1:img = sitk.Tile(img_slices, [maxlen, d])# TO DO check in code to get Tile Filter working with vector imageselse: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,
)
接下来,我们将 T1 图像与标签叠加进行可视化。
# Why doesn't this work? The images do overlap in physical space.
myshow3d(sitk.LabelOverlay(img_T1, img_labels),yslices=range(50, size[1] - 50, 20),zslices=range(50, size[2] - 50, 20),dpi=30,
)
解决我们问题的两种方法:(1) 将标签重新采样到图像网格上 (2) 将图像重新采样到标签网格上。从计算的角度来看,两者之间的差异取决于网格大小以及用于估计非网格位置值的插值器。
请注意,使用可以生成非标签值的插值器对标签图像进行插值是有问题的,因为您最终可能会得到比原始图像具有更多类别/标签的图像。这就是为什么我们在处理标签图像时只使用最近邻插值器的原因。
# Option 1: Resample the label image using the identity transformation
resampled_img_labels = sitk.Resample(img_labels,img_T1,sitk.Transform(),sitk.sitkNearestNeighbor,0.0,img_labels.GetPixelID(),
)
# Overlay onto the T1 image, requires us to rescale the intensity of the T1 image to [0,255] and cast it so that it can
# be combined with the color overlay (we use an alpha blending of 0.5).
myshow3d(sitk.LabelOverlay(sitk.Cast(sitk.RescaleIntensity(img_T1), sitk.sitkUInt8),resampled_img_labels,0.5,),yslices=range(50, size[1] - 50, 20),zslices=range(50, size[2] - 50, 20),dpi=30,
)
# Option 2: Resample the T1 image using the identity transformation
resampled_T1 = sitk.Resample(img_T1, img_labels, sitk.Transform(), sitk.sitkLinear, 0.0, img_T1.GetPixelID()
)
# Overlay onto the T1 image, requires us to rescale the intensity of the T1 image to [0,255] and cast it so that it can
# be combined with the color overlay (we use an alpha blending of 0.5).
myshow3d(sitk.LabelOverlay(sitk.Cast(sitk.RescaleIntensity(resampled_T1), sitk.sitkUInt8), img_labels, 0.5),yslices=range(50, size[1] - 50, 20),zslices=range(50, size[2] - 50, 20),dpi=30,
)
为什么上面的两个显示不同?(提示:在对“myshow3d”函数的调用中,y 和 z 切片的索引是相同的)。
myshow
和 myshow3d
函数确实很有用。它们已被复制到“myshow.py”文件中,以便可以导入到其他笔记本中。