利用GPU 提高多模医学图像配准速度

利用GPU 提高多模医学图像配准速度

宋 斐 宁夏医科大学档案馆 750004

本论文为教育部春晖计划项目（Z2011053），宋斐为项目负责人。

【文章摘要】

传统的图像配准计算强度大，难以满足实时应用的要求．具有较强的可编程性和出色的并行计算能力的GPU，为解决该问题提供了新的途径． 根据GPU 的自身特点，以薄板样条插值作为变换模型，构建了弹性配准计算平台．对多模态的两组图像进行实验，结果表明，含数据传输时间时，随着标记点个数的增加，GPU 的速度远远大于CPU 的速度，而且是量级的差距。

【关键词】

CPU ；实验

多模医学图像配准是指待配准的两幅图像来源于不同的成像设备。由于扫描设备的原理不同，扫描参数条件各异，所以两种断层图像间不是简单的一一对应关系。例如：CT 可以识别密度较大的组织，MR 针对软组织，SPECT 和PET 能够反映人体的功能和代谢信息。因此，在临床应用中， 常常需要将CT( 或MR) 与SPECT( 或PET） 配准。二者的结合能够较全面提供对象信息，具有临床应用价值。多模医学图像的配准一直是图像配准的重点研究课题，也越来越引起人们的关注。

目前，医学图像配准方法主要有矩主轴法和最大互信息配准法两大类．矩主轴法是借用经典力学中物体质量分布的概念 , 计算两幅图像像素点的质心和主轴 , 再通过平移和旋转使两幅图像的质心和主轴对齐 , 从而达到配准的目的。该方法首先将二幅图像标记点间建立对应关系， 然后利用插值法求取对应标记点之间的变换，进而将变换作用于整个待配准图像．该方法的缺点是计算时间长、运算量大、应用不方便．

图像处理器(graphics processing units，GPU) 是显卡的心脏，能够从硬件上支持多边形转换和光源处理（TransformandLighting，T&L）技术。如今的GPU 已经不再局限于3D 图形处理。在浮点运算、并行计算等方面，GPU 可以提供数十倍乃至于上百倍于CPU 的性能。利用GPU 可以完成排序、查找、搜索等基本的数据结构操作，同时在图像分割、图像滤波及色彩变换、数据挖掘、机器学习等诸多领域，同样有着广泛的应用．充分发挥GPU 的优异性能，并在此基础上建立计算模型进行图像配准，将有助于加快处理速度，增强图像处理的实时性．

1 多模医学图像配准方法

基于互信息的配准方法是直接利用图像的灰度值实现两幅图像间的配准。

具体过程如下：

首先根据两幅图像的基本情况预设一个初始参数x0，其中x0(1) 为裁剪旋转x0(3) 角的图像2 行的第一个索引。x0(2) 为裁剪旋转x0(3) 角的图像2 列的第一个索引，x0(3) 为旋转角度，x0(4) 为比例因子。然后按照给定的初始参数对图像2 进行变换，并计算图像1 和图像2 的互信息， 然后利用最优化工具箱中的fminsearch 函数在x0 附近寻找使图像1 和图像2 互信息最大的点，直至搜索到满足精度要求的参数；最后输出配准参数。

1.1 矩主轴法实现全局粗配准

利用特征点自动配对算法，分别找出PET 和CT 图像对应的特征点，然后分别算出两幅图像轮廓线的质心，并把PET 图像的质心移动到与CT 图像质心相同的位置，实现全局粗配准。

对于灰度图像，如果假设背景的灰度值是0，那么图像轮廓的灰度值为非0。对闭合灰度图像边界的轮廓进行跟踪。利用轮廓跟踪法按从上到下、从左到右的顺序查找边缘点，最后得出边界。

基于轮廓的矩主轴配准方法进行的配准，能够得到比较准确的配准结果。这种方法提取图像轮廓特征的方法也较为简单。另外，该方法适用于图像轮廓比较清晰的医学图像的配准。也就是说，如果图像的数据有缺失或图像轮廓不清晰, 运用基于轮廓的矩主轴配准方法的配准效果较差。因此，该方法适用于整个物体完整的出现在两幅图像中的图像配准。

1.2 层次B 样条自适应自由变形法实现细配准

由于不自主的生理运动或患者移动等使其内部的器官和组织的位置、尺寸和形状发生改变，器官脏器的局部位置存在弹性变形。如果直接对图像作弹性变形， 因为对图像任一部分的变形只利用了图像的局部信息，容易产生误配。直观地，可以认为B 样条函数生成一条平滑曲线（或者平滑面） 逼近控制点。考虑到拟合函数的精确性和变形的平滑性，采用出于计算量与配准精度考虑，基于层次B 样条自适应自由变形法对多模态医学图像进行自动细配准，流程如图1 所示。

2 基于GPU 的配准过程

在获得对应标记点之后，我们要建立一个变换模型，然后实现二幅图像的空间对应关系．在多模医学图像的配准过程中，通常要花费大量时间来完成该过程的实现．由于GPU 具有可编程片断处理器的特点，如果以薄板样条插值作为变换模型，就可以实现从浮动图像到参考图像的空间变换．具体过程如图2 所示．

图1 基于B 样条自适应自由变形法实现局部细配准092

智能应用

Intelligence Application

电子制作

图2 基于GPU 的配准流程

给定n 组对应标记点，薄板样条插值为

(1)

以薄板样条插值为变换模型，先将已知标记点数据按式(1) 转换成线性方程组， 然后使用Gauss．Jordan 消元法求解方程组．

2.1 载入标记数据

在CPU 中，数组是常见的数据组织形式，通过下标计算，其寻址方式实质上都是一维线性的；对于GPU，二维数组是天然的数据组织形式，因为它可以以纹理形式存储在GPU 中．根据需要，可采用亮度或RGBA 的形式，在每个纹理元素中存储1 个或4 个数据．因此首先要将已知的对应标记点信息，以纹理形式载入到GPU 中．

假设已知参考图像和浮动图像的n 组对应标记点，步骤如下：

步骤一：根据式(1)，逐一替代n 对已知的对应点，可以获得n 个方程．

步骤二：求解一个方向上的变换参数(x 方向或y 方向)，建立参数个数为n+3 的线性方程组。

步骤三：分别传入边长为的正方形纹理之中．

由于系数矩阵为(n+4)×(n+3)，因此内存空间的大小也要与之对应。GPU 的一大特点是支持正方形形状的纹理，所以最终传入方程组的矩阵是 (n+4)×(n+4) 大小的纹理．为了能够准确控制所访问的纹理数据，需要选择适当的正交投影和观察点。在GPU 实现其良好计算性能的过程中，纹理不仅要作为输入数据，还要作为输出数据．通过使用OpenGL(framebuffer object，FB）的扩展帧缓冲对象直接将结果渲染到纹理之上．

2.2 计算变换参数

在CPU 的编程环境中，通常利用循环语句对各数据元素进行操作，而要完成图形的渲染，就要使用多个循环语句． GPU 对图形的渲染可以说是手到擒来，其各个顶点数据或像素数据的计算过程也是各自独立的。在渲染的过程中，GPU 的片断处理器将同时对原始的每一个纹理元素执行相同的片断函数．因此，GPU 的一次渲染过程就能完成CPU 中的一轮循环操作．

在以纹理形式载人数据后，首先根据Gauss．Jordan 定义渲染过程中的纹理运算，其中每个纹理元素对应系数矩阵中的每一项系数．在第k 次消元过程中，需要更新数据矩阵中n×(n-k+1) 个数据元素， 而GPU 只需要对该纹理( 矩阵) 实施一次渲染就可以完成对全部元素的更新。

Gauss．Jordan 消元是一个多次迭代的过程，即：GPU 需要不断地执行新的渲染过程，而且要将上一次渲染过程所使用的输出纹理转换为下一个渲染过程的输入纹理．

根据前述内容，已知参考图像和浮动图像n 组对应点的情况下，需要建立(n+4)×(n+4) 大小的纹理。迭代了n+4 次后，就可以求出x 方向上的变换参数．同理，求出y 方向上的变换参数．

3 实验对比分析

本实验中，硬件采用的是INTEL(R) Pentium(R)4 2．26 GHz，256 MDDR 400 和Nvidia GeForce5950 Ultrra。

3.1 配准精度

实验选取了同一断层平面，相同模态和不同模态的2 套图像数据分别进行实验，图片的大小均为256×256。

在同模态实验中，以互相关系数作为子区域的相似性测度，实验分别选取了16、25、36 和49 个标记点，其配准结果如图3 所示，并以NMI 作为配准精度的评价标准．配准前后的精度比较如表l 所示． 多模态实验中所选取的二幅图像，整体上没有较大旋转和平移，但是图像若干局部( 例如边缘轮廓) 有明显的差异，为此实验以互信息作为子区域的相似性测度，分别选取了25、64、8l 和100 个标记点，其配准结果如图4 所示，并仍以NMI 作为配准精度的评价标准．配准前后的精度比较如表2 所示．

实验结果分析：

（1）需要选取的标记点数目与图像本身密切相关。

（2）当浮动图像有整体变形时，不需要选择很多标记点。

（3）当浮动图像有明显局部差异时， 只有在局部选取标记点，才能得到较好的配准效果．

3.2 配准速度

对GPU 和CPU 配准速度的比较，有以下几个前提条件：

（1）本实验以多模态的二幅图像为对象，分别选取16、32、64、128、256 和512 个标记点。

（2）计算变换参数的时间近似等于建

图3 同模态配准结果

图4 多模态配准结果

表1 同模态配准精度

表2 多模态配准精度