这个例子展示了如何使用基于强度的配准来自动对齐两个体积图像。
这个示例使用Imregister.
,imregtform
和imwarp
自动对齐两个容积数据集:CT图像和从同一患者收集的T1加权MR图像在不同时间。不像其他一些技术,Imregister.
和imregtform
不要寻找功能或使用控制点。基于强度的注册通常非常适合于医学和遥感图像。
本例中使用的三维CT和MRI数据集由迈克尔·菲茨帕特里克博士作为的一部分回顾性图像配准评估(RIRE)数据集.
本例使用了同一患者头部的两张3d图像。在配准问题中,我们将一幅图像视为固定图像,另一幅图像视为运动图像。配准的目的是使运动图像与固定图像对齐。在本例中,固定图像是T1加权MRI图像。我们要配准的运动图像是CT图像。数据以回顾性图像配准评估(RIRE)项目使用的文件格式存储。使用multibandread
读取包含图像数据的二进制文件。使用helperReadHeaderRIRE
函数获取与每个图像关联的元数据。你可使用以下连结查阅更多有关RIRE档案格式的资料:RIRE数据格式
固定主人= HelperReadheaderre(“rirePatient007MRT1.header”);movingHeader = helperReadHeaderRIRE (“rirePatient007CT.header”);fixedVolume = multibandread (“rirePatient007MRT1.bin”,...[fixedHeader。行,fixedHeader。列,fixedHeader。片),...“int16 = >单”0,“bsq”,“ieee-be”);movingVolume = multibandread (“rirePatient007CT.bin”,...[movingHeader。行,movingHeader。列,movingHeader。片),...“int16 = >单”0,“bsq”,“ieee-be”);
的helperVolumeRegistration
函数是一个辅助函数,用来帮助判断三维配准结果的质量。您可以交互地旋转视图,两个轴将保持同步。
helperVolumeRegistration (fixedVolume movingVolume);
你也可以用imshowpair
从固定和移动的体块来看单个平面,以获得体块整体对齐的感觉。在重叠的图像中imshowpair
在美国,灰色区域对应的是具有相似强度的区域,而品红和绿色区域显示的是其中一个图像比另一个图像更亮的地方。使用imshowpair
观察图像体积沿轴向切片通过每个体积的中心的错误配准。
centerFixed =大小(fixedVolume) / 2;centerMoving =大小(movingVolume) / 2;图imshowpair (movingVolume (:,:, centerMoving (3)), fixedVolume (:,:, centerFixed (3)));标题(“未登记的轴向片”)
的imregconfig
函数可以轻松地选择正确的优化器和度量配置Imregister.
.这两张图像来自两种不同的模式,MRI和CT,所以“多模式”选择是合适的。
[优化器,度量标准] = imregconfig(“多通道”);
所使用的算法Imregister.
当指定了关于输入图像的分辨率和/或位置的空间参考信息时,将更快地收敛到更好的结果。在本例中,CT和MRI数据集的分辨率在图像元数据中定义。使用此元数据来构造imref3d
我们将作为注册输入参数传递的空间引用对象。
rfixed = imref3d(大小(固定volume),fixedheader.pixelsize(2),fixedheader.pixelsize(1),固定主机.Slicethickness);rmoving = imref3d(size(thopingvolume),motherheader.pixelsize(2),motherheader.pixelsize(1),搬家主机.Slicethickness);
空间引用对象的属性定义了相关图像体积在世界坐标系中的位置以及每个维度的像素范围。Rmoving的XWorldLimits属性定义了移动体在X维中的位置。PixelExtentInWorld属性定义了X维度(沿列)中每个像素的世界单位的大小。移动体积在世界X坐标系中从0.3269扩展到334.97,每个像素的范围为0.6536mm。单位以毫米为单位,因为用于构建空间引用的标题信息以毫米为单位。ImageExtentInWorldX属性以世界单位确定移动图像体积的全部范围。
Rmoving。XWorldLimits
ans =1×20.3268 - 334.9674
rmoving.pixelextentinworldx.
ans = 0.6536
Rmoving。ImageExtentInWorldX
ans = 334.6406.
两卷之间的错位包括翻译,缩放和旋转。使用相似性转换来注册图像。
首先使用Imregister.
获得一个配准输出图像量,可以直接查看和观察配准结果的质量。
为优化器的InitialRadius属性指定一个非默认设置,以实现更好的注册结果收敛。
优化器。InitialRadius = 0.004;movingRegisteredVolume = imregister(movingVolume,Rmoving, fixedVolume,Rfixed,'死板的'优化器,指标);
使用imshowpair
再次重复检查轴向切片的对齐过程,通过注册体的中心来了解注册的成功程度。
图imshowpair (movingRegisteredVolume (:,:, centerFixed (3)), fixedVolume (:,:, centerFixed (3)));标题('注册量的轴向切片')
从上面的轴向切片来看,似乎配准成功了。使用helperVolumeRegistration
再次查看注册量,继续判断注册是否成功。
helperVolumeRegistration (fixedVolume movingRegisteredVolume);
的imregtform
函数可以用来当你对几何变换估计感兴趣的时候Imregister.
形成已注册的输出图像。imregtform
使用相同的算法Imregister.
接受和。相同的输入参数Imregister.
.由于目视检查所产生的体积Imregister.
表示注册成功,可以拨打imregtform
使用相同的输入参数来获得与该配准结果相关联的几何变换。
geomtform = imregtform(movingVolume,Rmoving, fixedVolume,Rfixed,)'死板的'、优化、指标)
geomtform = affine3d,属性:T: [4x4 double] dimension: 3
的结果imregtform
是一个几何变换对象。这个对象包含一个属性T,它定义了三维仿射变换矩阵。
geomtform。T
ans =4×40.9704 -0.0143 -0.2410 0 0.0228 0.9992 0.0324 0 0.2404 -0.0369 0.9700 0 -15.8648 -17.5692 29.1830 1.0000
的transformPointsForward
函数可用于确定点[u,v,w]在配准后在运动图像映射中的位置。因为空间引用的输入被指定为imregtform
,几何变换映射在世界坐标系中从移动到固定的点。的transformPointsForward
函数用于确定移动图像的中心在世界坐标系中的变换位置。
centerXWorld =意味着(Rmoving.XWorldLimits);centerYWorld =意味着(Rmoving.YWorldLimits);centerZWorld =意味着(Rmoving.ZWorldLimits);[xWorld, yWorld zWorld] = transformPointsForward (geomtform、centerXWorld centerYWorld, centerZWorld);
你可以使用worldToSubscript
功能要确定与移动卷中心对齐的固定卷的元素。
(r、c、p) = worldToSubscript (Rfixed、xWorld yWorld, zWorld)
r = 116
c = 132
p = 13
的imwarp
函数可用于应用从imregtform到三维体的几何变换估计。'OutputView' name-value参数用于定义一个空间引用参数,该参数确定重新采样的输出图像的世界限制和分辨率。你可以得到相同的结果Imregister.
通过使用与固定图像相关联的空间引用对象。这产生了输出量,其中固定和移动图像的世界限制和分辨率是相同的。一旦世界的限制和分辨率都是相同的,就会有像素对应于移动和固定卷的每个样本之间的像素对应关系。
vishrogisteredvolume = imwarp(移动volume,rmoving,geomtform,“双三次的”,“OutputView”, Rfixed);
使用imshowpair
再次查看轴向切片通过中心的注册体积产生imwarp
.
图imshowpair (movingRegisteredVolume (:,:, centerFixed (3)), fixedVolume (:,:, centerFixed (3)));标题('注册量的轴向切片')