CN109919987A - 一种基于gpu的三维医学图像配准相似度计算方法 - Google Patents

Abstract

一种基于GPU的三维医学图像配准相似度计算方法，使用具有NVIDIA GeForce系列显卡的操作平台，通过合理分配线程执行核函数，再结合CPU，串并行相结合来计算处理三维图像灰度信息，得出区域互信息，能够对三维医学图像配准进行相似性测度。相比较于常用的互信息方法，本发明可以获得更为精确的评价数据，处理速度通过GPU并行计算得到极大提高，能够为三维医学图像配准过程提供精确、稳定、可重复的方法。

Description

一种基于GPU的三维医学图像配准相似度计算方法

技术领域

本发明涉及信息技术计算机图像学下的医学成像、医学图像配准领域，尤其是一种基于GPU的三维医学图像配准相似度计算方法。

背景技术

随着医学影像技术的发展和计算机设备的革新，医学图像分辨率以及成像精度不断提高，使得医学图像在临床以及医学研究中得到广泛应用。图像配准技术是诸多医学图像应用中不可或缺的重要技术，将多种类型的信息匹配到同一空间下，便于疾病分析和诊断。

三维医学图像配准过程主要由图像预处理、图像变换、灰度插值、目标函数优化和相似性测度计算五个部分组成。图像配准的最终目的就是使待配准两图像间的相似度最大，这通常采用相似性测度值来定量确定。常用的相似性测度方法有误差平方和、均方误差、均方根误差、互相关系数和互信息量等。互信息概念来源于信息论，用于表现一个***中包含另一个***中信息量的多少。在图像配准中，互信息常用于表达一幅图像与另一幅图像的相关程度，当两幅图像在几何上对齐时，图像对应的像素灰度值的统计相关性(即互信息)达到最大。互信息因具有不需对图像进行预处理、配准精度高和良好的鲁棒性等优点而在图像配准技术中受到广泛应用。

互信息是目前常用的相似性测度方法，传统互信息在计算联合熵时，仅考虑两幅图像中单个像素点在相同坐标出现相同灰度的概率，而忽略了该像素在其邻域空间中与其他像素的相对关系，缺少对图像空间信息的把握。在图像有噪声或残缺、分辨率较低的情况下，互信息的计算容易使陷入局部极值。于是为了更好地解决相似性测度问题，出现了区域互信息测度概念来对图像进行配准的方法，将图像领域空间的灰度信息也作为重要衡量标准引入互信息计算，从而提高配准的精确性和鲁棒性。但相对于传统互信息直接读取单个点灰度出现次数，计算统计灰度图和概率密度，区域互信息计算过程更为复杂，存在大量高维矩阵计算，时间消耗远大于传统互信息的计算，付出的时间代价很大，如何更快速地计算高维矩阵是解决问题的关键。

发明内容

为克服上述区域互信息测度中计算高维矩阵耗时长的问题，本发明提供了一种基于GPU的三维医学图像配准相似度计算方法，能够大大缩短运行时间，从而提高配准的精确性和鲁棒性，使得区域互信息测度方法能够更好地应用于配准技术中。

本发明解决其技术问题所采用的技术方案是：

一种基于GPU的三维医学图像配准相似度计算方法，包括如下步骤：

步骤一，读取数据，在CPU内存中开辟内存空间，读入两幅相同分辨率的三维医学图像数据，两幅图像分别记为参考图像R和浮动图像F，对三维医学图像进行标记，三维图像的总体素个数N＝(x-2r)(y-2r)(z-2r)，其中x,y,z为三维图像分辨率，r为邻近点半径；在CPU内存中创建基于邻域空间信息的列向量数组，每个列向量数组包含d＝2(2r+1)³＝54个元素，参考图像R和浮动图像F的体素分别为R_i,j,k和F_i,j,k，将两幅图像每个对应的体素点与其邻域体素点信息连列，获得联合体素V_i,j,k，代表d维空间点P_i，所有P_i组成联合分布，即d×N大小的矩阵P，表示为P＝[P₁,…,P_N]；

步骤二，在显存中开辟空间，传输数据到显存；

在显存中为待运算的数据以及需要存放结果的变量开辟显存空间，CPU将图像数据读入***内存后，获得图像大小，使用cudaMalloc函数在GPU上开辟相同大小的内存空间并初始化，并开辟两个d×N大小的显存空间用于保存联合灰度直方矩阵P和平均值以及一个d×d大小的空间用于保存协方差矩阵C，用cudaMemcpy函数将CPU内存中待运算的图像信息数据拷贝到GPU显存中；

步骤三，计算联合体素的平均值使用GPU中的原子操作，统计求得

步骤四，分配线程，执行核函数，计算联合体素的平均值与联合分布P之差P₀；

每个线程块分配nthreadsPerBlock＝16×16＝256个线程，每个线程格分配个线程块，线程块中的线程索引为threadIndex＝threadIdx.x+threadIdx.y*blockDim.x，线程格中的线程块索引为blockIndex＝blockIdx.x+blockIdx.y*gridDim.x，线程索引为i＝threadIndex+blockIndex*blockDim.x*blockDim.y，每个线程i负责计算：全部线程计算完毕后，K_i组成d×N矩阵P₀；

步骤五，计算邻域每个点的协方差：

通过引入共享内存和矩阵分块概念，将大矩阵分解成若干子矩阵，每个线程块只读取块内共享内存中的子矩阵进行计算，再将共享内存中的矩阵合并得到最终结果，求得协方差C，其中C为d×d矩阵；

步骤六，用cudaMemcpy函数将计算结果传输到CPU内存中，基于CPU计算联合熵；

R^d空间内一组正太分布点的协方差矩阵的信息熵计算如下：

将步骤五求得的协方差C代入上述公式中，求得联合信息熵H(R，F)；

步骤七，求参考图像和浮动图像的信息熵H(R)、H(F)

取C_R为参考图像的协方差矩阵，C_F为浮动图像的协方差矩阵，同样代入公式，求得两幅图像各自的信息熵H(R)、H(F)。其中C_R为协方差矩阵C左上部分大小的子矩阵，C_F为协方差矩阵C的右下部分大小的子矩阵；

步骤八，计算区域互信息RMI

最终区域互信息为RMI(R,F)＝H(R)+H(F)-H(R,F)，通过区域互信息RMI值可以定量分析判断两幅三维医学图像的相似程度，即参考图像和浮动图像的相似度。

由上述步骤可以看出，区域互信息需要使用d×N＝54×(x-2)(y-2)(z-2)个数据进行联合信息熵的计算，需要进行大量的矩阵计算，复杂度高，由于数据处理过程具有独立性，可以使用GPU多线程并行处理，这种基于GPU的三维医学图像配准相似度计算方法可以大大提高配准的效率、精确性以及鲁棒性。

进一步，所述步骤三中，原子操作的核函数线程索引为index＝(threadIdx.x+threadIdx.y×blockDim.x)+(blockIdx.x+blockIdx.y×gridDim.x)×blockDim.x×blockDim.y，其中threadIdx.x和threadIdx.y分别为线程块中x,y维度上的线程索引，blockIdx.x和blockIdx.y分别为线程格中x,y维度上的线程块索引，blockDim.x和blockDim.y分别为线程块x,y维度上的线程数量，gridDim.x和gridDim.y分别为x,y维度上的线程块数量。首先，线程并行读取数据，然后使用原子操作中atomicAdd函数对第d个维度的值进行累加求平均值，此过程循环执行54次，直至求出所有维度的联合体素的平均值

再进一步，所述步骤四中，设置dim3类型线程块，每个线程块选择维度为16×16的大小，即共256个线程，为保证总体素数与分配的线程数是整数倍关系，且分配的线程大于或等于总体素数，以向上取整的策略求出分配线程需要的总线程块数为保证每个体素点都有唯一对应的线程和内核函数参数与计算，每个线程进行独立且相同的操作，根据线程索引i来获取数据和处理数据完成计算。

本发明由于采用上述技术方案，具有以下优点：能够有效提高三维医学图像的配准精确性和鲁棒性；通过GPU并行加速能够提高配准效率；能够为三维医学图像的研究提供稳定、高效、可重复的分析方法。

附图说明

图1是基于GPU的三维医学图像配准相似度计算方法的实现步骤流程图。

图2是区域互信息中图像数据与高维体素及联合体素关系，传统互信息将图像灰度信息以一维体素值分布为基础进行计算，而区域互信息将图像视为N个d维点的联合分布，每个联合体素点包含两幅图像中相同坐标点与其单位邻域的体素灰度分布。

图3是高维体素信息在GPU中的读取及储存，每个线程各自对应一个GPU显存储存区域，用于数据的存储和修改，数据具有独立性，各线程间不发生冲突。

具体实施方式

下面结合附图对本发明作进一步说明。

参照图1～图3，一种基于GPU的三维医学图像配准相似度计算方法，包括如下步骤：

步骤一，读取数据，在CPU内存中开辟内存空间，读入两幅相同分辨率的三维医学图像数据，两幅图像分别记为参考图像R和浮动图像F，对三维医学图像进行标记，由于医学图像中的图像主要信息均位于图像中心，忽略图像边缘体素点对最终信息熵的影响，三维图像的总体素个数N＝(x-2r)(y-2r)(z-2r)，其中x,y,z为三维图像分辨率，r为邻近点半径；

在CPU内存中创建基于邻域空间信息的列向量数组，每个列向量数组包含d＝2(2r+1)³＝54个元素，参考图像R和浮动图像F的体素分别为R_i,j,k和F_i,j,k，将两幅图像每个对应的体素点与其邻域体素点信息连列，获得联合体素V_i,j,k，代表d维空间点P_i，所有P_i组成联合分布，即d×N大小的矩阵P，表示为P＝[P₁,…,P_N]；

步骤二，在显存中开辟空间，传输数据到显存；

在显存中为待运算的数据以及需要存放结果的变量开辟显存空间，CPU将图像数据读入***内存后，获得图像大小，使用cudaMalloc函数在GPU上开辟相同大小的内存空间并初始化，并开辟两个d×N大小的显存空间用于保存联合灰度直方矩阵P和平均值以及一个d×d大小的空间用于保存协方差矩阵C，用cudaMemcpy函数将CPU内存中待运算的图像信息数据拷贝到GPU显存中；

步骤三，计算联合体素的平均值使用GPU中的原子操作，统计求得

步骤四，分配线程，执行核函数，计算联合体素的平均值与联合分布P之差P₀；

每个线程块分配nthreadsPerBlock＝16×16＝256个线程，每个线程格分配个线程块，线程块中的线程索引为threadIndex＝threadIdx.x+threadIdx.y*blockDim.x，线程格中的线程块索引为blockIndex＝blockIdx.x+blockIdx.y*gridDim.x，线程索引为i＝threadIndex+blockIndex*blockDim.x*blockDim.y，每个线程i负责计算：全部线程计算完毕后，K_i组成d×N矩阵P₀；

步骤五，计算邻域每个点的协方差：

协方差计算是矩阵乘法运算，通过引入共享内存和矩阵分块概念，将大矩阵分解成若干子矩阵，每个线程块只读取块内共享内存中的子矩阵进行计算，再将共享内存中的矩阵合并得到最终结果，求得协方差C，其中C为d×d矩阵；

步骤六，用cudaMemcpy函数将计算结果传输到CPU内存中，基于CPU计算联合熵；

R^d空间内一组正太分布点的协方差矩阵的信息熵计算如下：

将步骤五求得的协方差C代入上述公式中，求得联合信息熵H(R，F)；

步骤七，求参考图像和浮动图像的信息熵H(R)、H(F)；

取C_R为参考图像的协方差矩阵，C_F为浮动图像的协方差矩阵，同样代入公式，求得两幅图像各自的信息熵H(R)、H(F)。其中C_R为协方差矩阵C左上部分大小的子矩阵，C_F为协方差矩阵C的右下部分大小的子矩阵；

步骤八，计算区域互信息RMI；

最终区域互信息为RMI(R,F)＝H(R)+H(F)-H(R,F)，通过区域互信息RMI值可以定量分析判断两幅三维医学图像的相似程度，即参考图像和浮动图像的相似度。

由上述步骤可以看出，区域互信息需要使用d×N＝54×(x-2)(y-2)(z-2)个数据进行联合信息熵的计算，需要进行大量的矩阵计算，复杂度高，由于数据处理过程具有独立性，可以使用GPU多线程并行处理，这种基于GPU的三维医学图像配准相似度计算方法可以大大提高配准的效率、精确性以及鲁棒性。

进一步，所述步骤三中，原子操作的核函数线程索引为index＝(threadIdx.x+threadIdx.y×blockDim.x)+(blockIdx.x+blockIdx.y×gridDim.x)×blockDim.x×blockDim.y，其中threadIdx.x和threadIdx.y分别为线程块中x,y维度上的线程索引，blockIdx.x和blockIdx.y分别为线程格中x,y维度上的线程块索引，blockDim.x和blockDim.y分别为线程块x,y维度上的线程数量，gridDim.x和gridDim.y分别为x,y维度上的线程块数量。首先，线程并行读取数据，然后使用原子操作中atomicAdd函数对第d个维度的值进行累加求平均值，此过程循环执行54次，直至求出所有维度的联合体素的平均值

所述步骤四中，设置dim3类型线程块，每个线程块选择维度为16×16的大小，即共256个线程，为保证总体素数与分配的线程数是整数倍关系，且分配的线程大于或等于总体素数，以向上取整的策略求出分配线程需要的总线程块数为保证每个体素点都有唯一对应的线程和内核函数参数与计算，每个线程进行独立且相同的操作，根据线程索引i来获取数据和处理数据完成计算。

Claims (3)

1.一种基于GPU的三维医学图像配准相似度计算方法，其特征在于，所述方法包括以下步骤：

步骤一，读取数据，在CPU内存中开辟内存空间，读入两幅相同分辨率的三维医学图像数据，两幅图像分别记为参考图像R和浮动图像F，对三维医学图像进行标记，由于医学图像中的图像主要信息均位于图像中心，忽略图像边缘体素点对最终信息熵的影响，三维图像的总体素个数N＝(x-2r)(y-2r)(z-2r)，其中x,y,z为三维图像分辨率，r为邻近点半径；

在CPU内存中创建基于邻域空间信息的列向量数组，每个列向量数组包含d＝2(2r+1)³＝54个元素，参考图像R和浮动图像F的体素分别为R_i,j,k和F_i,j,k，将两幅图像每个对应的体素点与其邻域体素点信息连列，获得联合体素V_i,j,k，代表d维空间点P_i，所有P_i组成联合分布，即d×N大小的矩阵P，表示为P＝[P₁,…,P_N]；

步骤二，在显存中开辟空间，传输数据到显存；

在显存中为待运算的数据以及需要存放结果的变量开辟显存空间，CPU将图像数据读入***内存后，获得图像大小，使用cudaMalloc函数在GPU上开辟相同大小的内存空间并初始化，并开辟两个d×N大小的显存空间用于保存联合灰度直方矩阵P和平均值以及一个d×d大小的空间用于保存协方差矩阵C，用cudaMemcpy函数将CPU内存中待运算的图像信息数据拷贝到GPU显存中；

步骤三，计算联合体素的平均值使用GPU中的原子操作，统计求得

步骤四，分配线程，执行核函数，计算联合体素的平均值与联合分布P之差P₀；

每个线程块分配nthreadsPerBlock＝16×16＝256个线程，每个线程格分配个线程块，线程块中的线程索引为threadIndex＝threadIdx.x+threadIdx.y*blockDim.x，线程格中的线程块索引为blockIndex＝blockIdx.x+blockIdx.y*gridDim.x，线程索引为i＝threadIndex+blockIndex*blockDim.x*blockDim.y，每个线程i负责计算：全部线程计算完毕后，K_i组成d×N矩阵P₀；

步骤五，计算邻域每个点的协方差：

协方差计算是矩阵乘法运算，通过引入共享内存和矩阵分块概念，将大矩阵分解成若干子矩阵，每个线程块只读取块内共享内存中的子矩阵进行计算，再将共享内存中的矩阵合并得到最终结果，求得协方差C，其中C为d×d矩阵；

步骤六，用cudaMemcpy函数将计算结果传输到CPU内存中，基于CPU计算联合熵；

R^d空间内一组正太分布点的协方差矩阵的信息熵计算如下：

将步骤五求得的协方差C代入上述公式中，求得联合信息熵H(R，F)；

步骤七，求参考图像和浮动图像的信息熵H(R)、H(F)；

取C_R为参考图像的协方差矩阵，C_F为浮动图像的协方差矩阵，同样代入公式，求得两幅图像各自的信息熵H(R)、H(F)，其中C_R为协方差矩阵C左上部分大小的子矩阵，C_F为协方差矩阵C的右下部分大小的子矩阵；

步骤八，计算区域互信息RMI；

最终区域互信息为RMI(R,F)＝H(R)+H(F)-H(R,F)，通过区域互信息RMI值可以定量分析判断两幅三维医学图像的相似程度，即参考图像和浮动图像的相似度。

2.如权利要求1所述的一种基于GPU的三维医学图像配准相似度计算方法，其特征在于，所述步骤三中，原子操作的核函数线程索引为index＝(threadIdx.x+threadIdx.y×blockDim.x)+(blockIdx.x+blockIdx.y×gridDim.x)×blockDim.x×blockDim.y，其中threadIdx.x和threadIdx.y分别为线程块中x,y维度上的线程索引，blockIdx.x和blockIdx.y分别为线程格中x,y维度上的线程块索引，blockDim.x和blockDim.y分别为线程块x,y维度上的线程数量，gridDim.x和gridDim.y分别为x,y维度上的线程块数量；首先，线程并行读取数据，然后使用原子操作中atomicAdd函数对第d个维度的值进行累加求平均值，此过程循环执行54次，直至求出所有维度的联合体素的平均值

3.如权利要求1或2所述的一种基于GPU的三维医学图像配准相似度计算方法，其特征在于，所述步骤四中，设置dim3类型线程块，每个线程块选择维度为16×16的大小，即共256个线程，为保证总体素数与分配的线程数是整数倍关系，且分配的线程大于或等于总体素数，以向上取整的策略求出分配线程需要的总线程块数为保证每个体素点都有唯一对应的线程和内核函数参数与计算，每个线程进行独立且相同的操作，根据线程索引i来获取数据和处理数据完成计算。