基于B样条和水平集方法的医学图像联合分割与配准

吕凯; 吴军

doi:doi:10.3788/LOP57.101007

激光与光电子学进展, 2020, 57 (10): 101007, 网络出版: 2020-05-08

基于B样条和水平集方法的医学图像联合分割与配准下载： 808次

Joint Segmentation and Registration of Medical Image Based on B-Spline and Level Set Method

论文大纲

吕凯 ^*吴军

作者单位

江西理工大学信息工程学院, 江西赣州 341000

分层B样条水平集图像配准图像分割 hierarchical B-spline level set image registration image segmentation

AI 词云图 AI一句话精读 AI短摘要

注：本部分内容由 AI 自动生成，请您知悉。

摘要

针对分层B样条非刚性配准存在局部极值以及水平集分割方法不适用于噪声图像分割的问题,提出了一种基于局部更新分层B样条双向变换和水平集方法的医学图像联合分割与配准方法。该方法在分割算法中加入配准变换,在配准中融入图像分割的结构信息。使用B样条水平集函数对变换和分割的图像进行平滑表示,并在配准中引入双向变换以提高配准的精度和平滑性。在水平集方法的基础上,引入双向分层B样条变换构造分割与配准的联合能量泛函,并结合梯度下降法将能量泛函最小化以优化目标函数。实验结果表明:本方法与单独的图像分割方法相比,Dice度量均在99%以上;与单独的图像配准方法相比,均方误差下降了30%,能够提高图像的配准精度,且在分割噪声图像时有较好的鲁棒性。

Abstract

Aim

ing at the problems that non-rigid registration of layered B-spline exists local extremum and the level set segmentation method is not suitable for noisy image segmentation, a joint segmentation and registration of medical image based on locally updated hierarchical B-spline bidirectional transformation and level set method is proposed. The proposed method adds a registration transformation to the segmentation algorithm, and the structure information of image segmentation is incorporated into the registration. The B-spline level set function is used to smooth the transformed and segmented image, and a two-way transformation is introduced in the registration to improve the accuracy and smoothness of the registration. Based on the level set method, the bi-directional layered B-spline transform is introduced to construct the joint energy functional of segmentation and registration, and the gradient descent method is used to minimize the energy functional to optimize the objective function. Experimental results show that the Dice metric is always above 99% and the mean square error is reduced by 30% compared with the single image segmentation method. The proposed method improves the registration accuracy and has better robustness in noise image segmentation.

1 引言

图像配准和图像分割作为医学图像处理领域的关键技术,对临床诊断和手术定位有重要意义。传统方法都是将分割与配准独立进行,但随着医学图像处理的要求日益严格,独立使用分割^[1]和配准方法^[2]的局限性逐渐突出。

医学图像的配准方法可分为刚性配准和非刚性配准两大类,对未发生形变的图像使用刚性配准,相反图像间存在较大形变时就要使用非刚性配准。非刚性配准方法主要有基于样条和基于物理模型的方法^[3],常用的样条方法有薄板样条法、B样条法等。B样条法因具有良好的局部变换能力,成为国内外应用最广泛的样条方法。Gu等^[4]使用B样条法对形变较大的图像进行配准;Xie等^[5]将B样条控制网格逐层细化提出了分层B样条(HBS)配准方法;Pawar等^[6]为去除HBS中重叠的部分提出了截断分层B样条(THBS)配准方法;刘晨等^[7]提出了一种将尺度不变特征转换(SIFT)和分层B样条相结合的配准方法。虽然分层B样条能有效提高图像的配准精度^[8],但配准形变复杂的医学图像时,局部极值问题很突出。图像分割方法主要有基于阈值的分割方法^[9]、基于区域的活动轮廓分割方法^[10]、基于边缘^[11]的活动轮廓分割方法,在活动轮廓模型中,水平集方法应用较为广泛。Li等^[12]提出了一种基于区域的水平集图像分割方法用于医学图像分割; Chan等^[11]提出了一种基于Mumford-Shah函数和水平集的活动轮廓分割模型,对图像中的复杂结构进行分割;韩哲等^[13]提出了一种基于改进遗传模糊聚类和水平集方法的图像分割算法,提高了对目标边缘的分割精度,但对有噪声图像的分割效果欠佳。

为解决分层B样条配准的局部极值和水平集分割算法鲁棒性差的问题,引入了联合分割与配准方法,目前对于联合分割与配准方法的研究较少。Droske^[14]和Vemuri^[15]等提出了联合分割与配准的理论,Swierczynski等^[16]提出了一种将水平集方法和密集位移场相结合的联合分割与配准方法,提高了肺部电子计算机断层扫描(CT)图像的配准精度,证明了联合方法的可行性。本文结合了基于分层B样条的图像配准方法和基于水平集的分割方法,并在局部更新的分层B样条配准中引入双向变换,提出了一种基于局部更新分层B样条双向变换和水平集方法的医学图像联合分割与配准方法。先从待配准图像中分割出其水平集轮廓,然后使用局部更新的分层B样条变换模型将其变形并与参考图像进行匹配,从而实现对参考图像的分割;使用局部更新分层B样条变换进行配准的同时,利用分割得到的目标边界信息对图像特征进行精确变形,并引入双向变换提高配准的精度和平滑性;在水平集方法的基础上引入了双向分层B样条变换构造了分割与配准的联合能量泛函,实现联合分割与配准;结合梯度下降法将能量泛函最小化以优化目标函数,得到最终的分割与配准结果。

2 分割与配准模型

2.1 B样条变换

设二维图像中控制点的位置为φ_i_,_j,控制网格的间距为δ_x×δ_y,则图像上任意一点(x,y)的B样条变换可表示为

\begin{matrix} T (x, y) = \overset{3}{\sum_{k = 0}} \overset{3}{\sum_{n = 0}} B_{k} (u) B_{n} (v) φ_{i + k, j + n}, (1) \end{matrix}

式中,φ_i₊_k_,_j₊_n为邻近控制点的坐标位置, $\begin{matrix} i = [\frac{x}{δ_{x}}] - 1, j = [\frac{y}{δ_{y}}] - 1 分别为控制点的位置索引, [] 为取整函数, u = [\frac{x}{δ_{x}}] - \frac{x}{δ_{x}}, v = [\frac{y}{δ_{y}}] - \frac{y}{δ_{y}} 分别为 (x, y) 对应单元控制网格的位置, B_{k} (u) 为对应位置 (x 方向) 的三阶 B 样条基函数, B_{n} (v) 为对应位置 (y 方向) 的三阶 B 样条基函数。 \end{matrix}$ 以B_k(u)为例

\begin{matrix} \{\begin{matrix} B_{0} (u) = \frac{{(1 - u)}^{2}}{6} \\ B_{1} (u) = \frac{3 u^{3} - 6 u^{2} + 4}{6} \\ B_{2} (u) = \frac{- 3 u^{3} - 6 u^{2} + 4}{6} \\ B_{3} (u) = \frac{u^{3}}{6} \end{matrix} 。 (2) \end{matrix}

2.2 局部更新的分层B样条变换

分层B样条配准是通过逐层增加控制网格的密度直到某一网格下的配准误差小于设定的阈值,可以提高配准精度,但计算量大,因此将控制网格局部细化对提高配准效率有重要作用。Bornemann等^[17]提出了一种针对分层B样条控制网格局部细化的更新策略,其原理如图1所示,图中Ω₀~Ω₃分别表示不同细化级别下的控制网格。

图 1. HBS局部细化示意图

Fig. 1. Local refinement diagram of HBS

下载图片查看所有图片

局部更新分层B样条不是在每一个细化级别均匀的细化网格,而是只细化图像中变化较大的区域。首先在待配准图像上构造一个4×4的控制网格,然后通过调整控制点的位置实现图像的变形,对于图像变化复杂的区域,其控制网格会在局部进行细化和调整,从而完善图像局部的变形。局部更新分层B样条中,细化级别l上的B样条基函数都可以由细化级别l+1上的子B样条基函数的线性组合表示为

\begin{matrix} B_{m}^{l, 3} (u) = \overset{Nc}{\sum_{i = 1}} c_{m, i} B_{i}^{(l + 1), 3} (u), (3) \end{matrix}

式中,c_m_,_i为通过节点插入算法确定的控制网格细分系数, $\begin{matrix} B_{m}^{l, 3} \end{matrix}$ (u)为第m个控制网格细化级别l上的B样条基函数, $\begin{matrix} B_{i}^{(l + 1), 3} \end{matrix}$ (u)为 $\begin{matrix} B_{m}^{l, 3} \end{matrix}$ (u)在细化级别l+1中的子基函数,N_c为子B样条基函数的数量。

2.3 水平集方法

水平集方法的基本原理是通过零水平集在三维连续函数曲面中嵌入二维平面闭合曲线,根据二维平面闭合曲线的方程演化出水平集函数方程,以此对图像目标轮廓进行分割。设二维平面的曲线方程α=β(r),三维连续函数的曲面方程μ(r,s)=α-β(r),则μ(r,s)=0就隐式表达了二维平面闭合曲线,即用曲面方程中与平面曲线具有相同函数值的曲线隐式表示平面曲线,通过不断更新水平集函数,演化隐含在水平集函数中的闭合曲线,即使闭合曲线的拓扑结构发生了改变,水平集函数依然有效。

3 本文方法

3.1 配准中引入双向变换

使用的联合方法在局部更新的分层B样条配准过程中引入双向变换以增加图像配准的精确性和对称性。首先给定待配准图像F(X)和参考图像R(X),F(X)和R(X)之间的空间变换函数为T(X),通过迭代最小化能量泛函再求出最优T(X)。由当前迭代的T_n(X)和下一次迭代的变形场V_n₊₁可得到下一迭代的T_n₊₁(X)

\begin{matrix} T_{n + 1} (X) = T_{n} (X) \times (X + V_{n + 1}) 。 (4) \end{matrix}

前向配准中,将F(X)通过变换函数f(X)变形为R(X),后向配准中,将R(X)通过变换函数b(X)变形为F(X)。在双向变换过程中,利用更新的控制点位置可计算出更新的前向变形场V_f(X)和后向变形场V_b(X)

\begin{matrix} \begin{matrix} V_{f} (X) = \overset{Nb}{\sum_{m = 1}} C_{m}^{f} B_{m} (X) - X, (5) \\ V_{b} (X) = \overset{Nb}{\sum_{m = 1}} C_{m}^{b} B_{m} (X) - X, (6) \end{matrix} \end{matrix}

式中, $\begin{matrix} C_{m}^{f} \end{matrix}$ 和 $\begin{matrix} C_{m}^{b} \end{matrix}$ 分别为与基函数B_m(X)相关的一组控制点集,N_b为基函数的总数。

3.2 图像的平滑

采用B样条水平集函数对变换和分割后的图像进行C²平滑表示,平滑表示不仅能减少噪声的影响,还能提高数值计算的精度和收敛性。B样条水平集函数S(X)可表示为

\begin{matrix} S (X) = \overset{n 2}{\sum_{j = 1}} \overset{n 1}{\sum_{i = 1}} a_{i, j} N_{i, j}^{P} (x, y), (7) \end{matrix}

式中n₁,n₂分别为每个方向上单变量基函数的数量, $\begin{matrix} N_{i, j}^{P} \end{matrix}$ (x,y)为X在(x,y)处的P阶B样条基函数,a_i_,_j为对应位置的水平集系数,可表示为

\begin{matrix} a_{i, j} = \frac{\int N_{i, j}^{P} (x, y) g (x, y) dΩ}{\int N_{i, j}^{P} (x, y) dΩ}, (8) \end{matrix}

式中,g(x,y)为图像的强度值,Ω为图像域。

通过滤波的方法将图像进行平滑表示,将(2)式中的B样条基函数离散化得到离散B样条基函数B'(x)

\begin{matrix} B' (x) = \{\begin{matrix} \frac{{(2 + x)}^{3}}{6}, - 2 \leq x < - 1 \\ \frac{2}{3} - x^{2} - \frac{x^{3}}{2}, - 1 \leq x < 0 \\ \frac{2}{3} - x^{2} + \frac{x^{3}}{2}, 0 \leq x < 1 \\ \frac{{(2 - x)}^{3}}{6}, 1 \leq x < 2 \\ 0, ot h erwise \end{matrix} 。 (9) \end{matrix}

将B'(x)作为数字滤波器,通过卷积运算对图像进行平滑表示。则在图像X的坐标(x,y)处的平滑表示为

\begin{matrix} S (X) = \overset{m 1 + 3}{\sum_{i = m 1}} \overset{m 2 + 3}{\sum_{j = m 2}} a_{i, j} B' (x - i) B' (y - j), (10) \end{matrix}

式中,m₁=[x-2],m₂=[y-2]。因此(4)式可改写为

\begin{matrix} \begin{matrix} T_{n + 1} (X) = T_{n} (X) \times (X + V_{n + 1}) = \\ \overset{m 1 + 3}{\sum_{i = m 1}} \overset{m 2 + 3}{\sum_{j = m 2}} A_{i, j}^{n} N_{i, j}^{P} (X + V_{n + 1}), (11) \end{matrix} \end{matrix}

式中, $\begin{matrix} A_{i, j}^{n} \end{matrix}$ =[ $\begin{matrix} a_{i, j}^{x} \end{matrix}$ , $\begin{matrix} a_{i, j}^{y} \end{matrix}$ ]为当前迭代中更新的B样条水平集系数, $\begin{matrix} N_{i, j}^{P} (X + V_{n + 1}) 为 (X + V_{n + 1}) \end{matrix}$ 的P阶B样条基函数。

3.3 联合分割与配准

在水平集方法的基础上引入双向分层B样条变换构造分割与配准的联合能量泛函,从而实现联合分割与配准。通过L²梯度下降法最小化能量泛函以优化目标函数,从而更新控制点位置、变形场和变换等。联合方法的能量泛函由两部分组成,可表示为

\begin{matrix} E [f (X), b (X)] = E_{fid} [f (X), b (X)] + E_{reg} [V_{f} (X), V_{b} (X)] 。 (12) \end{matrix}

用来平衡配准和分割的部分E_fid $\begin{matrix} [f (X), b (X)] \end{matrix}$ 可表示为

\begin{matrix} \begin{matrix} E_{fid} [f (X), b (X)] = θ_{1} \int \{{[R (X) - c_{in}]}^{2} H_{ε} [\tilde{φ} (X)] + {[R (X) - c_{out}]}^{2} \{1 - H_{ε} [\tilde{ϕ} (X)]\}\} dΩ + \\ θ_{2} \int \{g^{f} (X) {[\tilde{F} (X) - R (X)]}^{2} + g^{b} (X) {[\tilde{R} (X) - F (X)]}^{2}\} dΩ, (13) \end{matrix} \end{matrix}

式中,θ₁和θ₂为加权参数,用于平衡分割和配准, $\begin{matrix} \tilde{F} \end{matrix}$ (X)=f $\begin{matrix} [F (X)] \end{matrix}$ 为变换后的待配准图像, $\begin{matrix} \tilde{R} \end{matrix}$ (X)=b $\begin{matrix} [R (X)] \end{matrix}$ 为变换后的参考图像。c_in和c_out分别为分割后参考图像的内部和外部区域图像的平均强度,可表示为

\begin{matrix} \begin{matrix} c_{in} = \frac{\int R (X) H_{ε} [\tilde{ϕ} (X)] dΩ}{\int H_{ε} [\tilde{ϕ} (X)] dΩ}, (14) \\ c_{out} = \frac{\int R (X) \{1 - H_{ε} [(\tilde{ϕ} (X)]\} dΩ}{\int \{1 - H_{ε} [\tilde{ϕ} (X)]\} dΩ} 。 (15) \end{matrix} \end{matrix}

与θ₁相乘的部分表示待配准图像的水平集轮廓变形后对参考图像进行水平集轮廓分割。 $\begin{matrix} \tilde{ϕ} \end{matrix}$ (X)=f $\begin{matrix} [ϕ (X)] \end{matrix}$ 为变换后的水平集轮廓,H_ε(X)为正则化的Heaviside函数,ε为正则化参数。

\begin{matrix} H_{ε} (X) = \frac{1}{2} [1 + \arctan (\frac{X}{ε})] 。 (16) \end{matrix}

与θ₂相乘的部分为最小化的两幅图像灰度值强度的误差平方和,g^f(X)和g^b(X)分别用于减少图像边界的配准,可表示为

\begin{matrix} \begin{matrix} g^{f} (X) = \frac{1}{\sqrt[]{γ + {(\frac{\partial \tilde{F} (X)}{\partial u})}^{2} + {(\frac{\partial \tilde{F} (X)}{\partial v})}^{2}}}, (17) \\ g^{b} (X) = \frac{1}{\sqrt[]{γ + {(\frac{\partial \tilde{R} (X)}{\partial u})}^{2} + {(\frac{\partial \tilde{R} (X)}{\partial v})}^{2}}}, (18) \end{matrix} \end{matrix}

式中,γ是为了防止分母变为0而加入的自定义参数,此处设为10^-11~10^-12。

用L²梯度下降法将E_fid $\begin{matrix} [f (X), b (X)] \end{matrix}$ 部分最小化后得到分层B样条变换的控制点集 $\begin{matrix} C_{m}^{f} \end{matrix}$ 和 $\begin{matrix} C_{m}^{b} \end{matrix}$ 的更新表达式为

\begin{matrix} \begin{matrix} \frac{C_{m}^{f} - C_{m, 0}^{f}}{ε_{t}} = - δ E_{m, fid}^{f} [f (X), b (X)], (19) \\ \frac{C_{m}^{b} - C_{m, 0}^{b}}{ε_{t}} = - δ E_{m, fid}^{b} [f (X), b (X)], (20) \end{matrix} \end{matrix}

式中, $\begin{matrix} C_{m}^{f} \end{matrix}$ 为当前迭代的前向配准控制点集, $\begin{matrix} C_{m}^{b} \end{matrix}$ 为当前迭代的后向配准控制点集, $\begin{matrix} C_{m, 0}^{f} \end{matrix}$ 和 $\begin{matrix} C_{m, 0}^{b} \end{matrix}$ 分别为初始转换时定义的控制点集,ε_t为t时刻的迭代步长。 $\begin{matrix} δ E_{m, fid}^{f} [f (X), b (X)] 和 δ E_{m, fid}^{b} [f (X), b (X)] 分别为 E_{fid} [f (X), b (X)] 相对控制点 C_{m}^{f} 和 C_{m}^{b} 的一阶变化, \end{matrix}$ 可表示为

\begin{matrix} \begin{matrix} δ E_{m, fid}^{f} [f (X), b (X)] = θ_{1} \int \{{[R (X) - c_{in}]}^{2} - \\ {[R (X) - c_{out}]}^{2}\} H'_{ε} [\tilde{ϕ} (X)] \nabla \tilde{ϕ} (X) B_{m} (X) dΩ + \\ 2 θ_{2} \int g^{f} (X) [\tilde{F} (X) - I_{t} (X)] \nabla \tilde{F} (X) B_{m} (X) dΩ, (21) \\ δ E_{m, f id}^{b} [f (X), b (X)] = 2 θ_{2} \int g^{b} (X) [\tilde{R} (X) - F (X)] \nabla \tilde{R} (X) B_{m} (X) dΩ, (22) \end{matrix} \end{matrix}

式中,Ñ为梯度。能量泛函中第二部分E_reg $\begin{matrix} [V_{f} (X), V_{b} (X)] \end{matrix}$ 用来规范双向变形场,可表示为

\begin{matrix} \begin{matrix} E_{reg} [V_{f} (X), V_{b} (X)] = λ_{1} \int ‖ X + V_{f, u} (X) ‖_{2}^{2} + ‖ X + V_{f, v} (X) ‖_{2}^{2}) dΩ + λ_{2} \int ‖ X + V_{b, u} (X) ‖_{2}^{2} + \\ ‖ X + V_{b, v} (X) ‖_{2}^{2}) dΩ + λ_{3} \int \{{[l_{f} (X) \times l_{b} (X) - X]}^{2} + {[l_{b} (X) \times l_{f} (X) - X]}^{2}\} dΩ, (23) \end{matrix} \end{matrix}

式中, $\begin{matrix} l_{f} (X) = X + V_{f} (X), l_{b} (X) = X + V_{b} (X), \end{matrix} \begin{matrix} ‖\dots ‖_{2}^{2} \end{matrix}$ 为求2范数的平方,λ₁、λ₂、λ₃分别为优化过程中的3个正则化参数。与λ₁和λ₂对应的正则项使图像在每个参数方向上保持平滑变形,令λ₁=λ₂以保持优化过程中的变换的对称性。与λ₃对应的正则项为双向变换的组合,用于将前向变换和后向变换约束为同一性,以此保证图像平滑逼真地变形,从而有效减少医学图像配准的偏差。

使用L²梯度下降法将E_reg $\begin{matrix} [V_{f} (X), V_{b} (X)] \end{matrix}$ 最小化,用以更新控制点集和变形场,控制点集 $\begin{matrix} P_{m}^{f} \end{matrix}$ 和 $\begin{matrix} P_{m}^{b} \end{matrix}$ 的更新表达式为

\begin{matrix} \begin{matrix} \frac{P_{m}^{f} - C_{m}^{f}}{ε_{t}} = - δ E_{m, reg}^{f} [V_{f} (X), V_{b} (X)], (24) \\ \frac{P_{m}^{b} - C_{m}^{b}}{ε_{t}} = - δ E_{m, reg}^{b} [V_{f} (X), V_{b} (X)], (25) \end{matrix} \end{matrix}

式中, $\begin{matrix} C_{m}^{f} 、 C_{m}^{b} 和 P_{m}^{f} 、 P_{m}^{b} 分别为能量泛函第一部分和第二部分最小化得出的控制点集。能量泛函 E_{reg} [V_{f} (X), V_{b} (X)] \end{matrix}$ 的一阶变化可表示为

\begin{matrix} \begin{matrix} δ E_{m, reg}^{f} [V_{f} (X), V_{b} (X)] = 2 λ_{1} \int [V_{f, u} (X) B_{m, u} (X) + V_{f, v} (X) B_{m, v} (X)] dΩ + \\ 2 λ_{3} \int [l_{f} (X) \times l_{b} (X) - X] B_{m} [l_{b} (X)] dΩ + 2 λ_{3} \int [l_{b} (X) \times l_{f} (X) - X] D \{l_{b} [l_{f} (X)]\} B_{m} (X) dΩ, (26) \\ δ E_{m, reg}^{b} [V_{f} (X), V_{b} (X)] = 2 λ_{2} \int [V_{b, u} (X) B_{m, u} (X) + V_{b, v} (X) B_{m, v} (X)] dΩ + \\ 2 λ_{3} \int [l_{f} (X) \times l_{b} (X) - X] D \{l_{f} [l_{b} (X)]\} B_{m} (X) dΩ + 2 λ_{3} \int [l_{b} (X) \times l_{f} (X) - X] B_{m} [l_{f} (X)] dΩ, (27) \end{matrix} \end{matrix}

式中, $\begin{matrix} D \{l_{b} [l_{f} (X)]\} 和 D \{l_{f} [l_{b} (X)]\} 分别为关于 l_{f} (X) 和 l_{b} (X) 的雅可比矩阵, B_{m, u} (X) 和 B_{m, v} (X) 分别为 B 样条基函数 B_{m} (X) 关于 u 、 v 的偏导数。 \end{matrix}$ 将新的控制点代入(5)式和(6)式计算出新的变形场,再代入(11)式可得到最新迭代的T(X)。

3.4 算法流程

基于局部更新分层B样条双向变换和水平集方法的医学图像联合分割与配准算法的流程如图2所示。

图 2. 算法流程图

Fig. 2. Flow chart of algorithm

下载图片查看所有图片

主要步骤:

1) 输入分辨率相同的两幅医学图像,分别记为待配准图像F(X)和参考图像R(X)。

2) 初始化分层B样条模型参数、迭代次数以及最小误差等参数。

3) 分割出待配准图像的水平集轮廓。

4) 用B样条水平集函数平滑表示初步变换图像和分割图像。

5) 构造分割与配准的联合能量泛函。

6) 计算分割精度的Dice相似度(DS)和待配准图像之间的相似性测度均方误差(MSE)。

7) 判断参数是否达到最优,若是则退出迭代,否则返回步骤5)。通过梯度下降法最小化能量泛函更新参数,首先使用L²梯度下降法最小化E_fid $\begin{matrix} [f (X), b (X)] \end{matrix}$ 更新控制点集 $\begin{matrix} C_{m}^{f} \end{matrix}$ 和 $\begin{matrix} C_{m}^{b} \end{matrix}$ 的位置,再使用L²梯度下降法最小化E_reg $\begin{matrix} [V_{f} (X), V_{b} (X)] \end{matrix}$ 更新控制点集 $\begin{matrix} P_{m}^{f} \end{matrix}$ 和 $\begin{matrix} P_{m}^{b} \end{matrix}$ ,并代入(5)式和(6)式计算出新的变形场V_f(X)和V_b(X),最后按(11)式更新双向变换函数后继续迭代,达到优化器最大迭代次数后退出优化。

8) 输出配准结果和分割结果。

4 实验结果与分析

4.1 联合分割与配准结果

实验环境:CPU为Intel Core i7-7700,内存为8 GB,操作系统为Windows10版本,软件平台为Matlab2016a。由于实验使用的联合分割与配准方法计算量较大,因此通过Matlab中的并行计算工具箱提高计算速度,计算量特别大的部分将MEX函数转换为C++函数,从而提高了计算效率。

对一组脑部核磁共振图像(MRI)和一组肺部CT图像使用本方法进行分割和配准的结果如图3所示。图3(a)为待配准图像,图3(b)为待配准图像对应的参考图像,图3(c)为两组图像的初始图像差,图3(d)为待配准图像的水平集轮廓,从图3(e)~图3(g)中控制网格的局部细化过程可以看出,对形变复杂的部分本方法能有效地局部细化控制网格,提高配准效率。图3(h)、图3(i)分别为前向和后向的变形网格,从图3(j)可以看出,本方法能准确分割脑部MRI中的脑白质凹凸区和薄区以及肺部CT图像中的肺边界,原因是局部细化的控制网格捕获了更精细的特征。图3(k)~图3(l)分别为两组图像的配准结果和配准后的图像差。

图 3. 脑部MRI和肺部CT图像的分割与配准结果。(a)待配准图像;(b)参考图像;(c)初始图像差;(d)分割的轮廓;(e)~(g)控制网格;(h)前向变形网格;(i)后向变形网格;(j)分割结果;(k)配准结果;(l)配准后图像差

Fig. 3. Segmentation and registration results of brain MRI and lung CT images. (a) Image to be registered; (b) image of target; (c) initial image difference; (d) segmented contours; (e)-(g) control grid; (h) forward deformed mesh; (i) backward deformed mesh; (j) result of segmentation; (k) result of registration; (l) image difference after registration

下载图片查看所有图片

4.2 分割和配准的评价标准

为客观评价分割和配准的效果,用MSE评估配准性能,用DS评估分割结果。MSE可表示为

\begin{matrix} X_{MSE} = \frac{\overset{M}{\sum_{x = 0}} \overset{N}{\sum_{y = 0}} [R (x, y) - F {(x, y)]}^{2}}{M \times N}, (28) \end{matrix}

式中,R(x,y)和F(x,y)分别为图像R和F在(x,y)处的像素灰度,M×N为图像的分辨率,MSE值越小表明配准效果越好。

DS可表示为

\begin{matrix} X_{DS} = \frac{2 N_{C}}{N_{R} + N_{F}} \times 100 %, (29) \end{matrix}

式中,N_C为待配准图像和参考图像共同区域的特定标签的像素数量,N_R和N_F分别为参考图像和待配准图像的像素数量。DS值越接近100%,表明分割的准确度越高。

4.3 配准和分割结果分析

两组不同形状图像的分割结果对比图如图4所示,图4(a)为三角形和C形的参考图像,图4(b)为圆形的待配准图像,图4(c)为待配准图像的分割结果,图4(d)~图4(f)分别为对参考图像分别添加强度为0.02、0.05、0.1的椒盐噪声得到的参考图像。为验证本方法在分割部分的鲁棒性,分别用本方法和水平集分割方法对三组图像进行分割,图4(g)~图4(i)为水平集分割方法的分割结果,图4(j)~图4(l)为本方法的分割结果。

对比图4(g)~图4(i)和图4(j)~图4(l)可以看出,与单独的水平集分割方法相比,本方法抗噪声干扰性强,具有较好的分割边界。两种方法的分割结果的DS如表1所示,表中数据均为5次对比实验的平均值。可以看出,水平集分割方法在噪声强度达到0.05时两组图像的DS就低于99.00%,而本方法分割的两组图像的DS始终高于水平集分割方法得到的结果,且都高于99.00%,这表明本方法分割含有噪声的图像时鲁棒性较高。

图 4. 分割结果。(a)无噪声参考图像;(b)待配准图像;(c)待配准图像的分割结果;(d)~(f)噪声强度为0.02、0.05、0.1的参考图像;(g)~(i)水平集方法对图(d)~图(f)的分割结果;(j)~(l)本方法对图(d)~图(f)的分割结果

Fig. 4. Results of segmentation. (a) Target images without noise-free; (b) images to be registered; (c) segmentation result of images to be registered; (d)-(f) target images with noise intensities of 0.02, 0.05, and 0.1; (g)-(i) segmentation results of Fig. (d)-Fig. (f) by level set method; (j)-(l) segmentation results of Fig. (d)-Fig. (f) by proposed method

下载图片查看所有图片

为了验证双向变换对配准的影响,选取三组医学图像将分层B样条方法和加入了双向变换的分层B样条方法(DHBS)进行对比实验,两种方法的配准结果如图5所示。图5(a)为三组待配准图像,图5(b)为三组参考图像。由图5(c)、图5(e)的配准结果和图5(d)、图5(f)中配准后的图像差可以看出,加入了双向变换方法的配准结果在图像轮廓处的配准效果较优,原因是双向变换的约束减少了变换过程中的偏差。两种方法的MSE和配准时间如表2所示,可以看出,加入双向变换方法的MSE值更小,三组图像的MSE值平均下降了约5%,这证明了在配准中加入双向变换可有效提高配准精度。

表 1. 两种分割方法的DS值

Table 1. DS value of two segmentation methods

Noise intensity	Level set method	Proposed method
Noise intensity	Triangle /%		C-shape /%	Triangle /%	C-shape /%
0.02	99.32	99.21	99.33	99.36
0.05	96.76	95.99	99.27	99.23
0.10	91.83	90.06	99.12	99.08

查看所有表

图 5. 两种方法配准结果。(a)待配准图像;(b)参考图像;(c)分层B样条方法配准结果;(d)分层B样条方法配准后图像差;(e)分层B样条双向变换方法配准结果;(f)分层B样条双向变换方法配准后图像差

Fig. 5. Registration results of two methods. (a) Images to be registered; (b) images of target; (c) registration results of HBS; (d) final image differences of HBS; (e) registration results of DHBS; (f) final image differences of DHBS

下载图片查看所有图片

表 2. 两种配准方法下MSE和配准时间

Table 2. MSE and registration time of two registration methods

Image	MSE	Registration time /s
Image	HBS		DHBS	HBS	DHBS
First group	35.98	33.86	762.44	799.68
Second group	45.52	42.37	892.54	927.59
Third group	37.56	36.15	810.26	831.63

查看所有表

为了验证本方法在图像配准方面的精确性,选取四组不同的医学图像将本方法、分层B样条方法以及文献[ 16]中的联合分割与配准方法进行对比实验,配准结果如图6所示。图6(a)为四组待配准图像,图6(b)为四组参考图像。对比图6(c)、图6(e)的配准结果,以及图6(f)、图6(h)中的配准后图像差可以看出,与单独的配准方法相比本方法对图像轮廓的配准效果较好,原因是本方法在配准中引入了分割轮廓的结构信息,为配准提供了良好的初始变换,并在配准中引入了双向变换以增加配准的对称性,有效减少了配准变换过程中的偏差。本方法、分层B样条和文献[ 16]方法的配准结果MSE和配准时间如表3所示,表中数据均为5次对比实验的平均值。

表 3. 三种配准结果的MSE和配准时间

Table 3. MSE and registration time of three registration methods

Image	MSE	Registration time /s
Image	HBS		Joint method	Proposed method	HBS	Joint method	Proposed method
First group	35.98	26.82	23.97	762.44	813.56	831.26
Second group	45.52	35.17	31.38	892.54	939.17	960.78
Third group	37.56	29.33	26.89	810.26	890.35	899.62
Fourth group	72.31	52.16	51.37	997.68	1105.29	1123.60

查看所有表

从表3中可以看出,虽然本方法的配准时间比分层B样条方法和文献[ 16]中的方法长,但本方法的MSE更小。相比分层B样条方法,本方法的配准结果MSE平均下降了约30%,这表明虽然加入了双向变换和分割的步骤导致配准时间有所延长,但是分割引入的边界轮廓信息以及配准中双向变换的约束提升了配准精度;相比文献[ 16]中的方法,本方法的配准结果MSE平均下降了8%,虽然配准肺部图像时两种方法差距不大,但配准脑部图像时本方法的优势更加突出,综合表明本方法的配准精度优于文献[ 16]中使用的方法。

为了检验本方法在配准时的抗噪性,将加入不同强度椒盐噪声的脑部图像作为待配准图像,分别使用本方法和文献[ 16]中的方法进行配准,结果如图7所示。图7(a)为三组不同噪声强度的待配准图像,图7(b)为无噪声的参考图像。不同噪声强度下图像配准结果的MSE和配准时间如表4所示,表中数据均为5次对比实验的平均值。

图 6. 配准结果。(a)待配准图像;(b)参考图像;(c)分层B样条方法配准结果;(d)联合方法配准结果;(e)本方法配准结果;(f)分层B样条方法配准后图像差;(g)联合方法配准后图像差;(h)本方法配准后图像差

Fig. 6. Results of registration. (a) Images to be registered; (b) images of target; (c) registration results of HBS; (d) registration results of joint method; (e) registration results of proposed method; (f) final image differences of HBS; (g) final image differences of joint method; (h) final image differences of proposed method

下载图片查看所有图片

图 7. 噪声实验对比图。(a)待配准图像;(b)参考图像;(c)联合方法配准结果;(d)本方法配准结果;(e)联合方法配准后图像差;(f)本方法配准后图像差

Fig. 7. Noise experiment comparison chart. (a) Images to be registered; (b) images of target; (c) registration results of joint method; (d) registration results of proposed method; (e) final image differences of joint method; (f) final image differences of proposed method

下载图片查看所有图片

表 4. 不同噪声强度下MSE和配准时间

Table 4. MSE and registration time under different noise intensities

Noise intensity	MSE	Registration time /s
Noise intensity	Joint method		Proposed method	Joint method	Proposed method
0	26.82	23.97	813.56	831.26
0.02	29.69	25.63	861.78	895.21
0.05	33.56	28.29	901.33	942.45

查看所有表

由表4可知,相比无噪声的图像,当噪声强度为0.02时,文献[ 16]使用的方法MSE值提高了10%,而本方法只提高了7%;当噪声强度为0.05时,文献[ 16]使用的方法MSE值提高了25%,而本方法只提高了18%。这表明相比于文献[ 16]中的方法,本方法在图像配准时的抗噪性更好。

5 结论

提出了一种基于局部更新分层B样条双向变换和水平集方法的医学图像联合分割与配准的方法,减少了噪声对图像分割的影响和局部极值对配准精度的影响。实验结果表明,本方法能够提高图像配准的精度,并在分割噪声图像时具有较好的鲁棒性,但由于该方法将分割和配准联合起来,所以相比单独的图像分割或图像配准方法所需时间更长,因此提高联合方法的配准和分割效率是接下来的研究重点。

参考文献

[1] 姚红兵, 卞锦文, 丛嘉伟, 等. 基于局部稀疏形状表示的医学图像分割模型[J]. 激光与光电子学进展, 2018, 55(5): 051011.

Yao H B, Bian J W, Cong J W, et al. Medical image segmentation model based on local sparse shape representation[J]. Laser & Optoelectronics Progress, 2018, 55(5): 051011.

[2] Cazoulat G, Owen D, Matuszak M M, et al. Biomechanical deformable image registration of longitudinal lung CT images using vessel information[J]. Physics in Medicine & Biology, 2016, 61(13): 4826-4839.

[3] 刘加林, 王慧琴, 王可, 等. 基于快速稳健特征最大子矩阵的光谱图像配准方法[J]. 激光与光电子学进展, 2019, 56(6): 063002.

Liu J L, Wang H Q, Wang K, et al. Spectral image registration method based on SURF maximum submatrix[J]. Laser & Optoelectronics Progress, 2019, 56(6): 063002.

[4] Gu S C, Meng X, Sciurba F C, et al. Bidirectional elastic image registration using B-spline affine transformation[J]. Computerized Medical Imaging and Graphics, 2014, 38(4): 306-314.

[5] Xie Z Y, Farin G E. Image registration using hierarchical b-splines[J]. IEEE Transactions on Visualization and Computer Graphics, 2004, 10(1): 85-94.

[6] Pawar A, Zhang Y J, Jia Y, et al. Adaptive FEM-based nonrigid image registration using truncated hierarchical B-splines[J]. Computers & Mathematics with Applications, 2016, 72(8): 2028-2040.

[7] 刘晨, 张龙波, 王雷, 等. 基于超像素重建的多尺度B样条医学图像配准[J]. 智能计算机与应用, 2019, 9(1): 24-27.

Liu C, Zhang L B, Wang L, et al. Multi-scale B-spline medical image registration based on pixel reconstruction[J]. Intelligent Computer and Applications, 2019, 9(1): 24-27.

[8] Pawar A, Zhang Y J, Anitescu C, et al. DTHB3D_Reg: dynamic truncated hierarchical B-spline based 3D nonrigid image registration[J]. Communications in Computational Physics, 2018, 23(3): 877-898.

[9] 周莉莉, 姜枫. 图像分割方法综述研究[J]. 计算机应用研究, 2017, 34(7): 1921-1928.

Zhou L L, Jiang F. Survey on image segmentation methods[J]. Application Research of Computers, 2017, 34(7): 1921-1928.

[10] 谭光鸿, 侯进, 韩雁鹏, 等. 基于卷积神经网络的低参数量实时图像分割算法[J]. 激光与光电子学进展, 2019, 56(9): 091003.

Tan G H, Hou J, Han Y P, et al. Low-parameter real-time image segmentation algorithm based on convolutional neural network[J]. Laser & Optoelectronics Progress, 2019, 56(9): 091003.

[11] Chan T F, Sandberg B, Vese L A. Active contours without edges for vector-valued images[J]. Journal of Visual Communication and Image Representation, 2000, 11(2): 130-141.

[12] Li C M, Huang R, Ding Z H, et al. A level set method for image segmentation in the presence of intensity inhomogeneities with application to MRI[J]. IEEE Transactions on Image Processing, 2011, 20(7): 2007-2016.

[13] 韩哲, 李灯熬, 赵菊敏, 等. 基于改进遗传模糊聚类和水平集的图像分割算法[J]. 计算机工程与设计, 2019, 40(5): 1390-1393, 1412.

Han Z, Li D A, Zhao J M, et al. Image segmentation algorithm based on improved genetic fuzzy clustering and level set[J]. Computer Engineering and Design, 2019, 40(5): 1390-1393, 1412.

[14] Droske M, Rumpf M. Multiscale joint segmentation and registration of image morphology[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2007, 29(12): 2181-2194.

[15] VemuriB, Chen YM. Joint image registration and segmentation[M] ∥Geometric Level Set Methods in Imaging, Vision, and Graphics. New York: Springer-Verlag, 2003: 251- 269.

[16] Swierczynski P. Papie z· B W, Schnabel J A, et al. A level-set approach to joint image segmentation and registration with application to CT lung imaging[J]. Computerized Medical Imaging and Graphics, 2018, 65: 58-68.

[17] Bornemann P B, Cirak F. A subdivision-based implementation of the hierarchical b-spline finite element method[J]. Computer Methods in Applied Mechanics and Engineering, 2013, 253: 584-598.

吕凯, 吴军. 基于B样条和水平集方法的医学图像联合分割与配准[J]. 激光与光电子学进展, 2020, 57(10): 101007. Kai Lü, Jun Wu. Joint Segmentation and Registration of Medical Image Based on B-Spline and Level Set Method[J]. Laser & Optoelectronics Progress, 2020, 57(10): 101007.

基于B样条和水平集方法的医学图像联合分割与配准下载： 808次

1 引言

2 分割与配准模型

2.1 B样条变换

2.2 局部更新的分层B样条变换

图 1. HBS局部细化示意图

Fig. 1. Local refinement diagram of HBS

2.3 水平集方法

3 本文方法

3.1 配准中引入双向变换

3.2 图像的平滑

3.3 联合分割与配准

3.4 算法流程

图 2. 算法流程图

Fig. 2. Flow chart of algorithm

4 实验结果与分析

4.1 联合分割与配准结果

图 3. 脑部MRI和肺部CT图像的分割与配准结果。(a)待配准图像;(b)参考图像;(c)初始图像差;(d)分割的轮廓;(e)~(g)控制网格;(h)前向变形网格;(i)后向变形网格;(j)分割结果;(k)配准结果;(l)配准后图像差

4.2 分割和配准的评价标准

4.3 配准和分割结果分析

图 4. 分割结果。(a)无噪声参考图像;(b)待配准图像;(c)待配准图像的分割结果;(d)~(f)噪声强度为0.02、0.05、0.1的参考图像;(g)~(i)水平集方法对图(d)~图(f)的分割结果;(j)~(l)本方法对图(d)~图(f)的分割结果

表 1. 两种分割方法的DS值

Table 1. DS value of two segmentation methods

图 5. 两种方法配准结果。(a)待配准图像;(b)参考图像;(c)分层B样条方法配准结果;(d)分层B样条方法配准后图像差;(e)分层B样条双向变换方法配准结果;(f)分层B样条双向变换方法配准后图像差

Fig. 5. Registration results of two methods. (a) Images to be registered; (b) images of target; (c) registration results of HBS; (d) final image differences of HBS; (e) registration results of DHBS; (f) final image differences of DHBS

表 2. 两种配准方法下MSE和配准时间

Table 2. MSE and registration time of two registration methods

表 3. 三种配准结果的MSE和配准时间

Table 3. MSE and registration time of three registration methods

图 6. 配准结果。(a)待配准图像;(b)参考图像;(c)分层B样条方法配准结果;(d)联合方法配准结果;(e)本方法配准结果;(f)分层B样条方法配准后图像差;(g)联合方法配准后图像差;(h)本方法配准后图像差

图 7. 噪声实验对比图。(a)待配准图像;(b)参考图像;(c)联合方法配准结果;(d)本方法配准结果;(e)联合方法配准后图像差;(f)本方法配准后图像差

Fig. 7. Noise experiment comparison chart. (a) Images to be registered; (b) images of target; (c) registration results of joint method; (d) registration results of proposed method; (e) final image differences of joint method; (f) final image differences of proposed method

表 4. 不同噪声强度下MSE和配准时间

Table 4. MSE and registration time under different noise intensities

5 结论

Article Outline

关于本站 Cookie 的使用提示

全站搜索

基于B样条和水平集方法的医学图像联合分割与配准 下载： 808次

1 引言

2 分割与配准模型

2.1 B样条变换

2.2 局部更新的分层B样条变换

图 1. HBS局部细化示意图

Fig. 1. Local refinement diagram of HBS

2.3 水平集方法

3 本文方法

3.1 配准中引入双向变换

3.2 图像的平滑

3.3 联合分割与配准

3.4 算法流程

图 2. 算法流程图

Fig. 2. Flow chart of algorithm

4 实验结果与分析

4.1 联合分割与配准结果

图 3. 脑部MRI和肺部CT图像的分割与配准结果。(a)待配准图像;(b)参考图像;(c)初始图像差;(d)分割的轮廓;(e)~(g)控制网格;(h)前向变形网格;(i)后向变形网格;(j)分割结果;(k)配准结果;(l)配准后图像差

4.2 分割和配准的评价标准

4.3 配准和分割结果分析

图 4. 分割结果。(a)无噪声参考图像;(b)待配准图像;(c)待配准图像的分割结果;(d)~(f)噪声强度为0.02、0.05、0.1的参考图像;(g)~(i)水平集方法对图(d)~图(f)的分割结果;(j)~(l)本方法对图(d)~图(f)的分割结果

表 1. 两种分割方法的DS值

Table 1. DS value of two segmentation methods

图 5. 两种方法配准结果。(a)待配准图像;(b)参考图像;(c)分层B样条方法配准结果;(d)分层B样条方法配准后图像差;(e)分层B样条双向变换方法配准结果;(f)分层B样条双向变换方法配准后图像差

Fig. 5. Registration results of two methods. (a) Images to be registered; (b) images of target; (c) registration results of HBS; (d) final image differences of HBS; (e) registration results of DHBS; (f) final image differences of DHBS

表 2. 两种配准方法下MSE和配准时间

Table 2. MSE and registration time of two registration methods

表 3. 三种配准结果的MSE和配准时间

Table 3. MSE and registration time of three registration methods

图 6. 配准结果。(a)待配准图像;(b)参考图像;(c)分层B样条方法配准结果;(d)联合方法配准结果;(e)本方法配准结果;(f)分层B样条方法配准后图像差;(g)联合方法配准后图像差;(h)本方法配准后图像差

图 7. 噪声实验对比图。(a)待配准图像;(b)参考图像;(c)联合方法配准结果;(d)本方法配准结果;(e)联合方法配准后图像差;(f)本方法配准后图像差

Fig. 7. Noise experiment comparison chart. (a) Images to be registered; (b) images of target; (c) registration results of joint method; (d) registration results of proposed method; (e) final image differences of joint method; (f) final image differences of proposed method

表 4. 不同噪声强度下MSE和配准时间

Table 4. MSE and registration time under different noise intensities

5 结论

Article Outline

相关论文

相关资讯

关于本站 Cookie 的使用提示

全站搜索

基于B样条和水平集方法的医学图像联合分割与配准下载： 808次