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

　　【摘要】传统的图像配准计算强度大，难以满足实时应用的要求．具有较强的可编程性和出色的并行计算能力的GPU，为解决该问题提供了新的途径．根据GPU的自身特点，以薄板样条插值作为变换模型，构建了弹性配准计算平台．对多模态的两组图像进行实验，结果表明，含数据传输时间时，随着标记点个数的增加，GPU的速度远远大于CPU的速度，而且是量级的差距。
　　【关键词】CPU；实验
　　多模医学图像配准是指待配准的两幅图像来源于不同的成像设备。由于扫描设备的原理不同，扫描参数条件各异，所以两种断层图像间不是简单的一一对应关系。例如：CT可以识别密度较大的组织，MR针对软组织，SPECT和PET能够反映人体的功能和代谢信息。因此，在临床应用中，常常需要将CT(或MR)与SPECT(或PET）配准。二者的结合能够较全面提供对象信息，具有临床应用价值。多模医学图像的配准一直是图像配准的重点研究课题，也越来越引起人们的关注。
　　目前，医学图像配准方法主要有矩主轴法和最大互信息配准法两大类．矩主轴法是借用经典力学中物体质量分布的概念,计算两幅图像像素点的质心和主轴,再通过平移和旋转使两幅图像的质心和主轴对齐,从而达到配准的目的。该方法首先将二幅图像标记点间建立对应关系，然后利用插值法求取对应标记点之间的变换，进而将变换作用于整个待配准图像．该方法的缺点是计算时间长、运算量大、应用不方便．
　　图像处理器(graphicsprocessingunits，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
　　智能应用
　　IntelligenceApplication
　　电子制作
　　图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(framebufferobject，FB）的扩展帧缓冲对象直接将结果渲染到纹理之上．
　　2.2计算变换参数
　　在CPU的编程环境中，通常利用循环语句对各数据元素进行操作，而要完成图形的渲染，就要使用多个循环语句．GPU对图形的渲染可以说是手到擒来，其各个顶点数据或像素数据的计算过程也是各自独立的。在渲染的过程中，GPU的片断处理器将同时对原始的每一个纹理元素执行相同的片断函数．因此，GPU的一次渲染过程就能完成CPU中的一轮循环操作．
　　在以纹理形式载人数据后，首先根据Gauss．Jordan定义渲染过程中的纹理运算，其中每个纹理元素对应系数矩阵中的每一项系数．在第k次消元过程中，需要更新数据矩阵中n×(n-k+1)个数据元素，而GPU只需要对该纹理(矩阵)实施一次渲染就可以完成对全部元素的更新。