simpleITK - Setup - matplotlib‘s imshow

使用 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 切片的索引是相同的)。

myshowmyshow3d 函数确实很有用。它们已被复制到“myshow.py”文件中,以便可以导入到其他笔记本中。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/42292.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

Github 热点项目 awesome-mcp-servers MCP 服务器合集,3分钟实现AI模型自由操控万物!

【今日推荐】超强AI工具库"awesome-mcp-servers"星数破万! ① 百宝箱式服务模块:AI能直接操作浏览器、读文件、连数据库,比如让AI助手自动整理Excel表格,三分钟搞定全天报表; ② 跨领域实战利器:…

硬件老化测试方案的设计误区

硬件老化测试方案设计中的常见误区主要包括测试周期不足、测试条件过于单一、样品选择不当等方面。其中,测试周期不足尤为突出,容易导致潜在缺陷未被完全暴露。老化测试本质上是通过加速产品老化来模拟长期使用状况,因此测试周期不足会严重削…

CSS学习笔记5——渐变属性+盒子模型阶段案例

目录 通俗易懂的解释 渐变的类型 1、线性渐变 渐变过程 2、径向渐变 如何理解CSS的径向渐变,以及其渐变属性 通俗易懂的解释 渐变属性 1. 形状(Shape) 2. 大小(Size) 3. 颜色停靠点(Color Sto…

Java StringUtils工具类常用方法详解

StringUtils是Apache Commons Lang库中一个极其常用的工具类,它提供了大量处理字符串的静态方法,能够简化我们的日常开发工作,提高代码的可读性和健壮性。下面我将详细介绍StringUtils类中最常用的方法及其使用场景。 一、StringUtils的基本…

设计模式(创建型)- 原型模式

目录 定义 类图 角色 优缺点 优点 缺点 应用场景 案例展示 浅克隆 深克隆 定义 原型模式旨在创建重复的对象,同时确保良好的性能表现。它通过复制现有对象(原型)来创建新对象,而非使用传统的构造函数创建方式。这种设计…

MQ的数据一致性,如何保证?

1 数据一致性问题的原因 这些年在Kafka、RabbitMQ、RocketMQ踩过的坑,总结成四类致命原因: 生产者悲剧:消息成功进Broker,却没写入磁盘就断电。消费者悲剧:消息消费成功,但业务执行失败。轮盘赌局&#x…

Angular由一个bug说起之十五:自定义基于Overlay的Tooltip

背景 工具提示(tooltip)是一个常见的 UI 组件,用于在用户与页面元素交互时提供额外的信息。由于angular/material/tooltip的matTooltip只能显示纯文本,所以我们可以通过自定义Directive来实现一个灵活且功能丰富的tooltip Overlay…

简单介绍一下Unity中的ScriptableObject

ScriptableObject的本质 ScriptableObject是Unity引擎中的一个特殊基类,允许你创建不依附于游戏对象的数据容器,以资产(Asset)形式存储在项目中。这些资产: 可在编辑器中创建和配置 在构建后作为资产打包 可通过Resources或AssetBundle加…

ubuntu24.04.2 NVIDIA GeForce RTX 4060笔记本安装驱动

https://www.nvidia.cn/drivers/details/242281/ 上面是下载地址 sudo chmod x NVIDIA-Linux-x86_64-570.133.07.run # 赋予执行权限把下载的驱动复制到家目录下,基本工具准备,如下 sudo apt update sudo apt install build-essential libglvnd-dev …

LabVIEW 布尔控件回车键触发程序退出

在 LabVIEW 开发过程中,部分用户可能会遇到按下回车键(Enter)后,程序意外退出的问题。该问题主要源于布尔控件的属性设置冲突,包括键分配、数据绑定及 Tab 键行为等。本文将详细分析问题根源,并提供一套完整…

分布式系统面试总结:3、分布式锁(和本地锁的区别、特点、常见实现方案)

仅供自学回顾使用,请支持javaGuide原版书籍。 本篇文章涉及到的分布式锁,在本人其他文章中也有涉及。 《JUC:三、两阶段终止模式、死锁的jconsole检测、乐观锁(版本号机制CAS实现)悲观锁》:https://blog.…

WebWorkers在项目中的使用案例

Worker | 文档 worker 线程的关闭在主线程和 worker 线程都能进行操作,但对 worker 线程的影响略有不同。 // main.js(主线程) const myWorker new Worker(/worker.js); // 创建worker myWorker.terminate(); // 关闭worker 复制代码 // wor…

vue ts+Windi CSS

1、创建vue项目 trae(字节)打开一个空文件夹 npm install -g vue/cli vue create my-project cd my-project vue add typescript npm run serve vue项目创建完成 2、安装windicss vue add windicss vue.config.js配置 npm install vue-router …

【HTML 基础教程】HTML 编辑器

HTML 编辑器推荐 可以使用专业的 HTML 编辑器来编辑 HTML,菜鸟教程为大家推荐几款常用的编辑器: VS Code:Visual Studio Code - Code Editing. RedefinedSublime Text:http://www.sublimetext.com/在线编辑器:HTML/C…

文件上传的小点总结(2)

4.黑名单绕过(.htaccess方法) 源码一打开,遇到这样的黑名单是不是看的头皮发麻,这么多后缀都禁用。 .htaccess可以启用或禁用apache的功能,利用这个特点,我们可以使用该文件来禁用上述黑名单功能,从而上传**文件。 简…

mysql--主从复制--部署

MySQL 主从复制部署教程 一、主节点(Master)配置 1. 创建目录结构 mkdir -p /usr/local/src/mysql_demo/master_replica/{logs,configFile,data}2. 编写主节点的 MySQL 配置文件 my.cnf 路径:/usr/local/src/mysql_demo/master_replica/co…

Qt弹出新窗口并关闭(一个按钮)

参考:Qt基础 练习:弹出新窗口并关闭的两种实现方式(两个按钮、一个按钮)_qt打开一个窗口另一个关闭-CSDN博客 实现: 一个按钮,点击一次,按钮的名字从open window变为close window,…

游戏引擎学习第185天

回顾并计划今天的内容 我们完成了开始整理这些数据的工作,但我们还没有机会真正去查看这些数据的具体内容,因为我们只是刚刚开始了数据整理的基本工作。我们收集了大量的信息,但到目前为止,仍然没有足够的可视化工具来帮助我们理…

《一本书讲透Elasticsearch:原理、进阶与工程实践》读书笔记

1:es的组成部分: Elasticsearch 引擎:核心组件,处理索引和搜索请求 Kibana:es的可视化的数据界面,用于分析和展示数据 Beats(可选)轻量级的日志采集器 2:基本概念 es开…

[React 进阶系列] 组合组件 复合组件

[React 进阶系列] 组合组件 & 复合组件 今天写个人项目练手的时候搜到了一个比价有趣的实现,于是用了一下,发现这个 concept 不是特别的熟,于是上网找了下,返现了一个叫 复合组件(compound components) 的概念。搜索了一下后…