CN111080681A - 一种基于LoG算子的3D/2D医学图像配准方法 - Google Patents

Abstract

本发明公开了一种基于LoG算子的3D/2D医学图像配准方法，先获取2D和3D医学图像并进行去噪预处理，再根据X射线成像计算模型和刚体变换模型建立投影坐标系；然后利用GPU的并行计算完成DRR图像的快速生成过程，并基于LoG算子建立相似性测度函数，最后利用混合优化算法对相似性测度函数进行迭代寻优，找出满足精度配准时的刚体变换参数，则迭代停止，从而完成配准过程。

Description

一种基于LoG算子的3D/2D医学图像配准方法

技术领域

本发明属于图像配准技术领域，更为具体地讲，涉及一种基于LoG算子的3D/2D医学图像配准方法。

背景技术

随着医学成像技术的发展，临床需要不同维度的医学图像为手术提供导航，术中拍摄的2D图像能够给医生提供实时信息，却丢失了器官组织的三维空间信息，这些空间信息对于手术影像导航至关重要，但术中难以获取三维图像。为了能在术中有效利用术前的3D图像信息，需要确定术前3D图像和术中2D图像的位置信息和方向信息，并通过空间变换，调整3D图像位姿，使术前3D图像投影产生的图像与术中2D图像一致，这里就用到了2D图像与3D图像的配准技术。配准后的3D图像与2D图像，能够为影像导航手术实时提供空间信息，提高影像导航手术的安全性与可靠性。

目前的配准方法可分为基于灰度的方法、基于特征的方法和基于深度学习的方法。基于特征的方法又可分为基于内部特征和外部特征，基于内部特征的方法通过提取对图像进行分割并提取特定的特征信息，这就导致损失了一定的图像信息，且配准结果取决于分割精度与特征准确度；基于外部特征的方法通过在患者体内植入标记物，根据标记物的距离测度来实现配准，但植入物可能对人体造成损害，总的来说基于特征的配准方法存在较大的误差。基于灰度的方法利用相似性测度，通过迭代算法逐个比较像素之间的差距，仍然没有挖掘和应用图像中像素之间的内在联系，计算量大且存在精度误差。基于深度学习的配准方法通过神经网络实现端到端的参数预测，难以保证较高的配准精度，且需要大量的配对的真实医学图像进行训练，这在实际中是难以获取的。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于LoG算子的3D/2D医学图像配准方法，实现了术前3D图像和术中2D图像实时配准，从而为影像导航手术提供高精度的实时空间信息，确保手术的安全性与可靠性。

为实现上述发明目的，本发明一种基于LoG算子的3D/2D医学图像配准方法，其特征在于，包括以下步骤：

(1)、获取2D和3D医学图像；

获取术中X射线图像作为配准过程的2D参考图像，获取术前CT序列作为配准过程的待配准3D浮动图像，然后分别对2D参考图像和3D浮动图像进行去噪处理；

(2)、初始化混合优化算法；

初始迭代步数k＝0，设置最大迭代步数为Max；初始刚体变换参数为T₀；设置迭代阈值为ε，迭代终止温度为Temp；设置迭代算法启用次数count＝0，迭代算法最大启用次数为Time，初始化迭代标记flag＝0；设置相似性测度阈值F；

(3)、建立投影坐标系；

以CT中心为原点，根据X射线成像计算模型和刚体变换模型建立投影坐标系；其中，投影坐标系中，用T＝t_x,t_y,t_z,r_x,r_y,r_z表示六维的刚体变换参数；其中，t_x表示在X轴上的平移参数，t_y表示在Y轴上的平移参数，t_z表示在Z轴上的平移参数，r_x表示沿X轴的旋转参数，r_y表示沿Y轴的旋转参数，r_z表示沿Z轴的旋转参数；

(4)、通过迭代更新的方式实现3D/2D医学图像配准；

(4.1)、判断迭代标记flag，若flag＝0，则进入步骤(4.2)；若flag＝1，则进入步骤(4.3)；

(4.2)、利用Powell优化算法对相似性测度函数进行迭代寻优；

(4.2.1)、当迭代步数为k时，从步骤(1)中随机选取一张3D浮动图像，并输入至投影坐标系，根据初始刚体变换参数为T_k，投影坐标系投影产生DRR_k图像，再将DRR_k图像和2D参考图像的像素大小统一到m×n，像素值归一化到同一灰度范围，处理完成后，将2D参考图像标记为I_f，将DRR_k图像标记为I_k；

(4.2.2)、基于LoG算子，构建图像I_f和I_k的相似性测度函数S(I_f,I_k)；

S(I_f,I_k)＝-[N₁(I_f,I_k)+N₂(L_f,L_k)]

L_f(i,j)＝LoG(i,j)*I_f(i,j)，L_k(i,j)＝LoG(i,j)*I_k(i,j)

其中，符号*表示卷积，I_f(i,j)为参考图像在像素(i,j)处的像素值，I_k(i,j)为DRR_k图像在像素(i,j)处的像素值，

为参考图像的像素平均值，

为DRR_k图像的像素平均值，L_f为I_f经过LoG算子卷积后的图像，L_k为I_k经过LoG算子卷积后的图像，

为图像L_f的像素平均值，

为图像L_k的像素平均值，N₁(I_f,I_k)表示参考图像与DRR图像的NCC值，N₂(L_f,L_k)表示图像L_f和图像L_k的NCC值，LoG算子中的参数σ为高斯分布的标准差；

(4.2.3)、利用Powell优化算法对相似性测度函数进行寻优，得到本轮迭代后的相似性测度函数值以及下一轮迭代时的刚体变换参数为T_k+1；

(4.2.4)、若相邻两次迭代的相似性测度函数值之差的绝对值小于阈值ε，或k＞Max，则停止迭代，输出刚体变换参数T_k，同时更新T₀＝T_k，flag＝1，k＝0，count＝count+1，并跳转至步骤(4.4)；否则，令k＝k+1，再返回至步骤(4.2.1)进行下一轮迭代寻优；

(4.3)、利用模拟退火算法对相似性测度函数进行迭代寻优；

(4.3.1)、按照步骤(4.2.1)-(4.2.2)所述方法，构建图像I_f和I_k的相似性测度函数S(I_f,I_k)；

(4.3.2)、利用模拟退火算法对相似性测度函数进行寻优，得到本轮迭代后的相似性测度函数值以及下一轮迭代时的刚体变换参数为T_k+1；

(4.3.3)、若模拟退火算法的温度小于Temp，或k＞Max，则停止迭代，输出刚体变换参数T_k，同时更新T₀＝T_k，flag＝0，k＝0，count＝count+1，并跳转至步骤(4.4)；否则，令k＝k+1，再返回至步骤(4.3.1)进行下一轮迭代寻优；

(4.4)、若相似性测度函数值满足|S(I_f,I_k)-F|＜ε或count＞Time，则结束配准，将刚体变换参数T_k作为最终的输出；否则，令F＝S(I_f,I_k)，再返回到步骤(4.1)；

(4.5)、根据输出的变换参数T_k去调整3D浮动图像，从而到达3D、2D图像的配准。

本发明的发明目的是这样实现的：

本发明一种基于LoG算子的3D/2D医学图像配准方法，先获取2D和3D医学图像并进行去噪预处理，再根据X射线成像计算模型和刚体变换模型建立投影坐标系；然后利用GPU的并行计算完成DRR图像的快速生成过程，并基于LoG算子建立相似性测度函数，最后利用混合优化算法对相似性测度函数进行迭代寻优，找出满足精度配准时的刚体变换参数，则迭代停止，从而完成配准过程。

同时，本发明一种基于LoG算子的3D/2D医学图像配准方法还具有以下有益效果：

(1)、相似性测度函数中，在基于像素值信息的归一化互相关的基础上，引入基于LoG算子的归一化互相关，即引入了图像的方向信息和边缘信息，能够有效利用图像的灰度信息和特征信息，提高配准的鲁棒性，将配准精度提高了一个量级，而且目标函数的收敛速度更快，所需迭代步数更少。

(2)、本发明采用Powell算法和模拟退火算法的混合优化策略，结合Powell算法局部收敛快和模拟退火算法全局寻优能力强的优点，当Powell算法收敛域局部最优时，模拟退火算法能够跳出局部最优去寻找全局最优，限制局部最优值，提高全局最优能力，有效解决容易陷入局部最优值的问题。

(3)、利用GPU的并行计算，能够有效解决3D图形计算带来的时间效率问题。

(4)、本发明实现了术前3D图像和术中2D图像实时配准，从而为影像导航手术提供高精度的实时空间信息，确保手术的安全性与可靠性。

附图说明

图1是本发明一种基于LoG算子的3D/2D医学图像配准方法流程图；

图2是通过混合优化算法迭代更新来实现3D/2D图像的配准流程图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

为了方便描述，先对具体实施方式中出现的相关专业术语进行说明：

GPU(Graphic Processing Unit)：图形处理器；

DRR(Digitally Reconstructed Radiographs)：数字重建放射影像；

NCC(Normalized Cross Correlation)：归一化互相关；

LoG(Laplacian of Gaussian)：高斯拉普拉斯。

图1是本发明一种基于LoG算子的3D/2D医学图像配准方法流程图。

在本实施例中，如图1所示，本发明一种基于LoG算子的3D/2D医学图像配准方法，包括以下步骤：

S1、获取2D和3D医学图像；

获取术中X射线图像作为配准过程的2D参考图像，获取术前CT序列作为配准过程的待配准3D浮动图像，然后分别对2D参考图像和3D浮动图像进行去噪处理；

S2、初始化混合优化算法；

初始迭代步数k＝0，设置最大迭代步数为Max；初始刚体变换参数为T₀；设置迭代阈值为ε，迭代终止温度为Temp；设置迭代算法启用次数count＝0，迭代算法最大启用次数为Time，初始化迭代标记flag＝0；设置相似性测度阈值F；

在本实施例中，初始变换参数为T₀＝(0,0,0,0,0,0)，设最大迭代步数为Max＝100，迭代阈值为ε＝0.001，迭代终止温度为Temp＝10^-8，迭代算法最大启用次数为Time＝10，初始F＝10。

S3、建立投影坐标系；

以CT中心为原点，根据X射线成像计算模型和刚体变换模型建立投影坐标系；其中，投影坐标系中，用T＝t_x,t_y,t_z,r_x,r_y,r_z表示六维的刚体变换参数；其中，t_x表示在X轴上的平移参数，t_y表示在Y轴上的平移参数，t_z表示在Z轴上的平移参数，r_x表示沿X轴的旋转参数，r_y表示沿Y轴的旋转参数，r_z表示沿Z轴的旋转参数；

在本实施例中，投影坐标系中，光源点、CT中心、投影面板中心位于同一条直线上；

S4、如图2所示，通过迭代更新的方式实现3D/2D医学图像配准；

S4.1、判断迭代标记flag，若flag＝0，则进入步骤S4.2；若flag＝1，则进入步骤S4.3；

S4.2、利用Powell优化算法对相似性测度函数进行迭代寻优；

S4.2.1、当迭代步数为k时，从步骤S1中随机选取一张3D浮动图像，并输入至投影坐标系，根据初始刚体变换参数为T_k，利用GPU并行计算，在投影坐标系中对3D浮动图像投影产生DRR_k图像，再将DRR_k图像和2D参考图像的像素大小统一到m×n，像素值归一化到同一灰度范围，在本实施例中，m＝512，n＝512，像素值归一化到灰度范围[0,255]；处理完成后，将2D参考图像标记为I_f，将DRR_k图像标记为I_k；

S4.2.2、在本实施例中，采用结合像素值和LoG算子的归一化互相关NCC作为相似性测度函数，因此，构建图像I_f和I_k的相似性测度函数S(I_f,I_k)；

L_f(i,j)＝LoG(i,j)*I_f(i,j)，L_k(i,j)＝LoG(i,j)*I_k(i,j)

其中，符号*表示卷积，I_f(i,j)为参考图像在像素(i,j)处的像素值，I_k(i,j)为DRR_k图像在像素(i,j)处的像素值，

为参考图像的像素平均值，

为DRR_k图像的像素平均值，L_f为I_f经过LoG算子卷积后的图像，L_k为I_k经过LoG算子卷积后的图像，

为图像L_f的像素平均值，

为图像L_k的像素平均值，N₁(I_f,I_k)表示参考图像与DRR图像的NCC值，N₂(L_f,L_k)表示图像L_f和图像L_k的NCC值，LoG算子中的参数σ为高斯分布的标准差，参数σ的取值需要根据不同医学图像器官组织来选取不同的值；在本实施例中，在对头骨进行配准的实验中发现，当LoG算子的σ＝1时，滤波后的图像质量较好；

S4.2.3、利用Powell优化算法对相似性测度函数进行寻优，得到本轮迭代后的相似性测度函数值以及下一轮迭代时的刚体变换参数为T_k+1；

S4.2.4、若相邻两次迭代的相似性测度函数值之差的绝对值小于阈值ε，或k＞Max，则停止迭代，输出刚体变换参数T_k，同时更新T₀＝T_k，flag＝1，k＝0，count＝count+1，并跳转至步骤S4.4；否则，令k＝k+1，再返回至步骤S4.2.1进行下一轮迭代寻优；

S4.3、利用模拟退火算法对相似性测度函数进行迭代寻优；

S4.3.1、按照步骤S4.2.1-S4.2.2所述方法，构建图像I_f和I_k的相似性测度函数S(I_f,I_k)；

S4.3.2、利用模拟退火算法对相似性测度函数进行寻优，得到本轮迭代后的相似性测度函数值以及下一轮迭代时的刚体变换参数为T_k+1；

S4.3.3、若模拟退火算法的温度小于Temp，或k＞Max，则停止迭代，输出刚体变换参数T_k，同时更新T₀＝T_k，flag＝0，k＝0，count＝count+1，并跳转至步骤S4.4；否则，令k＝k+1，再返回至步骤S4.3.1进行下一轮迭代寻优；

S4.4、若相似性测度函数值满足|S(I_f,I_k)-F|＜ε或count＞Time，则结束配准，将刚体变换参数T_k作为最终的输出；否则，令F＝S(I_f,I_k)，再返回到步骤S4.1；

S4.5、根据最终输出的变换参数T_k去调整3D浮动图像，从而到达3D、2D图像的配准。

综上所述，本发明所涉及的基于LoG算子的3D/2D医学图像配准方法，其理论基础为LoG算子和混合优化策略。通过在NCC相似性测度中引入LoG算子，结合图像的像素信息、边缘信息和方向信息，目标函数收敛速度更快，解决配准的精度问题；Powell算法和模拟退火算法的混合优化策略，有效解决优化陷入局部最优的问题；GPU的并行计算能力，有效解决3D图像计算带来的时间效率问题。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (1)

1.一种基于LoG算子的3D/2D医学图像配准方法，其特征在于，包括以下步骤：

(1)、获取2D和3D医学图像；

获取术中X射线图像作为配准过程的2D参考图像，获取术前CT序列作为配准过程的待配准3D浮动图像，然后分别对2D参考图像和3D浮动图像进行去噪处理；

(2)、初始化混合优化算法；

初始迭代步数k＝0，设置最大迭代步数为Max；初始刚体变换参数为T₀；设置迭代阈值为ε，迭代终止温度为Temp；设置迭代算法启用次数count＝0，迭代算法最大启用次数为Time，初始化迭代标记flag＝0；设置相似性测度阈值F；

(3)、建立建立投影坐标系；

以CT中心为原点，根据X射线成像计算模型和刚体变换模型建立投影坐标系；其中，投影坐标系中，用T＝t_x,t_y,t_z,r_x,r_y,r_z表示六维的刚体变换参数；其中，t_x表示在X轴上的平移参数，t_y表示在Y轴上的平移参数，t_z表示在Z轴上的平移参数，r_x表示沿X轴的旋转参数，r_y表示沿Y轴的旋转参数，r_z表示沿Z轴的旋转参数；

(4)、通过迭代更新的方式实现3D/2D医学图像配准；

(4.1)、判断迭代标记flag，若flag＝0，则进入步骤(4.2)；若flag＝1，则进入步骤(4.3)；

(4.2)、利用Powell优化算法对相似性测度函数进行迭代寻优；

(4.2.1)、当迭代步数为k时，从步骤(1)中随机选取一张3D浮动图像，并输入至投影坐标系，根据初始刚体变换参数为T_k，投影坐标系投影产生DRR_k图像，再将DRR_k图像和2D参考图像的像素大小统一到m×n，像素值归一化到同一灰度范围，处理完成后，将2D参考图像标记为I_f，将DRR_k图像标记为I_k；

(4.2.2)、基于LoG算子，构建图像I_f和I_k的相似性测度函数S(I_f,I_k)；

S(I_f,I_k)＝-[N₁(I_f,I_k)+N₂(L_f,L_k)]

L_f(i,j)＝LoG(i,j)*I_f(i,j)，L_k(i,j)＝LoG(i,j)*I_k(i,j)

其中，符号*表示卷积，I_f(i,j)为参考图像在像素(i,j)处的像素值，I_k(i,j)为DRR_k图像在像素(i,j)处的像素值，

为参考图像的像素平均值，

为DRR_k图像的像素平均值，L_f为I_f经过LoG算子卷积后的图像，L_k为I_k经过LoG算子卷积后的图像，

为图像L_f的像素平均值，

为图像L_k的像素平均值，N₁(I_f,I_k)表示参考图像与DRR图像的NCC值，N₂(L_f,L_k)表示图像L_f和图像L_k的NCC值，LoG算子中的参数σ为高斯分布的标准差；

(4.2.3)、利用Powell优化算法对相似性测度函数进行寻优，得到本轮迭代后的相似性测度函数值以及下一轮迭代时的刚体变换参数为T_k+1；

(4.2.4)、若相邻两次迭代的相似性测度函数值之差的绝对值小于阈值ε，或k＞Max，则停止迭代，输出刚体变换参数T_k，同时更新T₀＝T_k，flag＝1，k＝0，count＝count+1，并跳转至步骤(4.4)；否则，令k＝k+1，再返回至步骤(4.2.1)进行下一轮迭代寻优；

(4.3)、利用模拟退火算法对相似性测度函数进行迭代寻优；

(4.3.1)、按照步骤(4.2.1)-(4.2.2)所述方法，构建图像I_f和I_k的相似性测度函数S(I_f,I_k)；

(4.3.2)、利用模拟退火算法对相似性测度函数进行寻优，得到本轮迭代后的相似性测度函数值以及下一轮迭代时的刚体变换参数为T_k+1；

(4.3.3)、若模拟退火算法的温度小于Temp，或k＞Max，则停止迭代，输出刚体变换参数T_k，同时更新T₀＝T_k，flag＝0，k＝0，count＝count+1，并跳转至步骤(4.4)；否则，令k＝k+1，再返回至步骤(4.3.1)进行下一轮迭代寻优；

(4.4)、若相似性测度函数值满足|S(I_f,I_k)-F|＜ε或count＞Time，则结束配准，将刚体变换参数T_k作为最终的输出；否则，令F＝S(I_f,I_k)，再返回到步骤(4.1)；

(4.5)、根据输出的变换参数T_k去调整3浮动图像，从而到达3D、2D图像的配准。

Images