注册多模态3-D医学图像
这个例子展示了如何使用基于强度的配准自动对齐两个体积图像。
这个例子使用了imregister
,imregtform
而且imwarp
自动对齐两个容量数据集:同一患者在不同时间收集的CT图像和T1加权MR图像。不像其他技术,imregister
而且imregtform
不要寻找特性或使用控制点。基于强度的登记通常非常适合于医疗和遥感图像。
本例中使用的三维CT和MRI数据集由Michael Fitzpatrick博士作为一部分回顾性图像配准评估(RIRE)数据集.
步骤1:加载图像
这个例子使用了同一病人头部的两张3d图像。在配准问题中,我们认为一个图像是固定图像,另一个图像是运动图像。配准的目标是使运动图像与固定图像对齐。在本例中,固定图像为T1加权MRI图像。我们想要配准的运动图像是CT图像。数据以回顾性图像配准评估(RIRE)项目使用的文件格式存储。使用multibandread
读取包含图像数据的二进制文件。使用helperReadHeaderRIRE
函数获取与每个图像关联的元数据。
fixedHeader = helperReadHeaderRIRE(“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 = size(fixedVolume)/2;centerMoving = size(movingVolume)/2;图imshowpair(movingVolume(:,:,centerMoving(3)), fixedVolume(:,:,centerFixed(3)));标题(“未注册轴向切片”)
步骤2:设置初始注册
的imregconfig
函数可以很容易地选择正确的优化器和指标配置来使用imregister
.这两张图片来自两种不同的方式,核磁共振成像和CT,所以“多通道”
选项是合适的。
[optimizer,metric] = imregconfig(“多通道”);
所使用的算法imregister
当指定了关于输入图像的分辨率和/或位置的空间引用信息时,将更快地收敛到更好的结果。在这种情况下,CT和MRI数据集的分辨率在图像元数据中定义。使用此元数据来构造imref3d
空间引用对象以作为注册的输入参数传递。
Rfixed = imref3d(size(fixedVolume),fixedHeader.PixelSize(2),fixedHeader.PixelSize(1),fixedHeader.SliceThickness);Rmoving = imref3d(size(movingVolume),movingHeader.PixelSize(2),movingHeader.PixelSize(1),movingHeader.SliceThickness);
空间引用对象的属性定义了相关图像体积在世界坐标系中的位置以及每个维度的像素范围。的XWorldLimits
的属性Rmoving
定义移动体在X维中的位置。的PixelExtentInWorld
属性定义了X维世界单位中每个像素的大小(沿着列)。移动体积在世界X坐标系中从0.3269到334.97 mm,每个像素的范围为0.6536 mm。单位以毫米为单位,因为用于构建空间引用的标题信息是以毫米为单位的。的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
再次重复检查通过注册体中心拍摄的轴向切片的对齐过程,以了解注册的成功程度。
figure imshowpair(movingRegisteredVolume(:,:,centerFixed(3)), fixedVolume(:,:,centerFixed(3)));标题(“登记体积轴向切片”)
从上面的轴向切片中,您可以看到卷已经成功对齐。使用helperVolumeRegistration
再次查看注册量,继续判断注册是否成功。
helperVolumeRegistration (fixedVolume movingRegisteredVolume);
步骤3:获得移动与固定对齐的3- d几何变换。
的imregtform
函数,当您对所使用的几何变换估计感兴趣时,可以使用imregister
形成注册后的输出图像。imregtform
使用相同的算法imregister
并接受相同的输入参数imregister
.自目视检查所产生的体积imregister
表示注册成功,可以拨打电话imregtform
使用相同的输入参数来获得与该配准结果相关的几何变换。
geomtform = imregtform(movingVolume,Rmoving, fixedVolume,Rfixed,“刚性”,优化器,度量)
geomtform = affine3d属性:T: [4x4 double]维度: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 = mean(remove . xworldlimits);centerYWorld = mean(删除。yworldlimits);centerZWorld = mean(remove . 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
变成一个3-D体积。的“OutputView”
名称-值参数用于定义空间引用参数,该参数确定输出重采样图像的世界限制和分辨率。你可以得到相同的结果imregister
通过使用与固定图像相关联的空间引用对象。这将创建一个输出量,其中固定和移动图像的世界限制和分辨率是相同的。一旦两个体积的世界极限和分辨率相同,移动和固定体积的每个样本之间就有像素到像素的对应关系。
movingRegisteredVolume = imwarp(movingVolume,Rmoving,geomtform,“双三次的”,“OutputView”, Rfixed);
使用imshowpair
再次查看轴向切片通过注册体的中心产生的imwarp
.
figure imshowpair(movingRegisteredVolume(:,:,centerFixed(3)), fixedVolume(:,:,centerFixed(3)));标题(“登记体积轴向切片”)