一种基于配准的CT图像全心脏自动分割系统

摘要

本发明公开了一种基于配准的CT图像全心脏自动分割系统，该系统利用不同个体间心脏CT的相似性，把配准应用于分割，实现全心脏的自动分割。该系统包括如下模块：输入模块、分割模块、粗配准模块、精配准模块、变换模块和输出模块。该系统能实现全心脏的自动分割并能同时分割出多个心脏腔室。

说明书

技术领域

本发明涉及医学图像分割系统，具体涉及一种基于配准的CT图像全心脏 自动分割系统。

背景技术

医学图像分割就是对感兴趣区域进行提取的一个过程，是进行医学图像三 维可视化前的关键步骤。心脏分割能把心脏各部位分割出来，从而为临床诊断 提供大量的解剖信息，并能为心脏手术的术前计划以及心脏介入手术提供依 据。对分割后的各腔室进行体积计算还能为心功能评价提供量化的信息。

目前，医学图像的分割可分为手动分割、半自动分割和自动分割。手动分 割就是医生利用临床知识在医学图像上勾画出感兴趣区域的轮廓，精度最高。 但效率较低，分割结果不可重复，而且分割结果受到分割者的经验制约。半自 动分割是利用计算机通过人机交互进行分割，分割结果在一定程度仍受到分割 者经验知识的制约。自动分割是由计算机完成医学图像分割的全过程，其分割 结果能够重现。传统的自动分割完全脱离医生的临床经验，因此应用于临床时 有一定的局限性。而本发明把医学专家的临床经验知识与自动分割过程相结 合，能克服传统自动分割的局限性，因此对临床诊断有重要的意义。

发明内容

本发明利用不同个体间心脏CT的相似性，把配准应用于分割中，把医学 专家的临床经验知识与自动分割相结合，实现了全心脏及心脏各腔室的自动分 割。

本发明的全心脏自动分割系统通过以下技术方案实现：

本发明提出一种基于配准的CT图像全心脏自动分割系统，包括：

输入模块100，用于输入心脏CT序列图像模板和待分割的心脏CT图像， 并将心脏CT模板图像传送至分割模块200，将输入心脏CT序列图像模板和 待分割的心脏CT图像传送至粗配准模块300；

分割模块200，用于在CT序列图像模板上分割出心脏的各个部位，为各 部位标记上标记符，制作出Label图像，并传送至变换模块500；该Label图 像能应用于所有分割实例中，无需重复制作；

粗配准模块300，用于基于仿射变换的粗配准，对模板作出平移、缩放、 旋转变换，并将粗配准结果传送给精配准模块400，将粗配准变换参数传送至 变换模块500；

精配准模块400，用于基于B样条变换的精配准，并将精配准变换参数传 送至变换模块500；

变换模块500，根据粗配准和精配准变换参数对Label图像作出变换，并 将变换结果传送至输出模块600，变换结果即自动分割结果；

输出模块600，用于输出自动分割结果，分割结果中不同的标记符即对应 于心脏中不同的部位。

所述粗配准模块300依照步骤(3.1)～(3.8)进行处理：

(3.1)将输入心脏CT序列图像模板作为浮动图像，将待分割的心脏CT 图像作为参考图像，把浮动图像和参考图像分解为若干个不同的分辨率层；

(3.2)在当前分辨率层上对浮动图像进行仿射变换，并记录变换参数；

(3.3)对变换后的浮动图像进行插值；

(3.4)对参考图像与插值后的浮动图像进行相似性度量；

(3.5)利用优化算法搜寻下一次变换参数；

(3.6)重复步骤(3.2)～(3.5)，直至优化算法的迭代次数到达最大值；

(3.7)用当前辨率层的参数结果对下一分辨率层的变换参数进行初始化；

(3.8)在所有分辨率层上重复步骤(3.2)～(3.7)，得到基于仿射变换的 粗配准结果和粗配准变换参数。

所述精配准模块400依照步骤(4.1)～(4.8)进行处理：

(4.1)将基于仿射变换的粗配准的结果分解为若干个不同的分辨率层， 并将其作为新的浮动图像；

(4.2)在当前分辨率层上对新的浮动图像进行B样条变换，并记录变换 参数；

(4.3)对变换后的浮动图像进行插值；

(4.4)对参考图像与插值后的浮动图像进行相似性度量；

(4.5)通过优化算法搜寻下一次变换的参数；

(4.6)重复步骤(4.2)～(4.5)，直至优化算法的迭代次数到达最大值；

(4.7)用当前辨率层的参数结果对下一分辨率层的变换参数进行初始化；

(4.8)在所有分辨率层上重复步骤(4.2)～(4.7)，得到精配准变换参数。

所述粗配准模块300和精配准模块400的多分辨率分解采用金字塔算法； 所述插值方法采用B样条插值法；所述相似性度量采用互信息度量法；所述 优化算法采用自适应随机梯度下降法。

本发明与现有技术相比具有如下优势：

(1)把医学专家的临床经验知识与自动分割过程相结合，在医学专家的 临床经验知识的指导下实现全自动分割。

(2)通过一次自动分割即可同时分割并标记心脏各腔室。

附图说明

图1为基于配准的CT图像全心脏自动分割系统的结构图。

图2为基于配准的CT图像全心脏自动分割系统配准方法的流程图。

图3为几种Sigmoid函数的函数曲线。

图4为Label自动分割结果效果图，图中不同值对应于心脏不同部分。

图5为全心脏的分割结果效果图。

图6为主动脉的分割结果效果图。

图7为左心房的分割结果效果图。

图8为左心室腔的分割结果效果图。

图9为左心室心肌的分割结果效果图。

图10为右心房的分割结果效果图。

图11为右心室的分割结果效果图。

具体实施方式

为了更好地理解本发明的技术方案，以下结合附图和实施例进行详细的描 述，但本发明的实施和保护范围不限于此。

实施例

本实施例的序列模板图像和待分割的心脏图像分辨率为512*512。

如图1所示，本发明提出一种基于配准的CT图像全心脏自动分割系统， 其工作流程如下：

第一步，输入模块100输入心脏CT序列图像模板和待分割的心脏CT图 像。

第二步，输入模块100将心脏CT模板图像传送至分割模块200，在CT 序列图像模板上分割出心脏的各个部位，为各部位标记上标记符，制作出Label 图像，并传送至变换模块500；该Label图像能应用于所有分割实例中，无需 重复制作。

第三步，输入模块100将心脏CT序列图像模板和待分割的心脏CT图像 传送至粗配准模块300，进行粗配准处理，具体步骤如下：

(3.1)将输入心脏CT序列图像模板作为浮动图像，将待分割的心脏CT 图像作为参考图像，把浮动图像和参考图像分解为6不同的分辨率层。

(3.2)在当前分辨率层上对浮动图像进行仿射变换，并记录变换参数；

所述仿射变换可以分解为平移变换、旋转变换以及缩放变换。其数学表达 式为：

T(x)＝RSx+t                (1)

其中x ∈Ω为变换前的点，Ω为图像域，x＝[xyz]^T，R为旋转变换矩阵， S为缩放变换矩阵，t为平移变换矩阵。

平移变换矩阵t的表达式为：

$t=\left(\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right)---\left(2\right)$

其中t_x，t_y，t_z分别为沿x，y，z轴方向平移的距离。

缩放变换矩阵S的表达式为：

$S=\left(\begin{array}{ccc}{S}_x& 0& 0\\ 0& S_y& 0\\ 0& 0& S_z\end{array}\right)---\left(3\right)$

其中S_x，S_y，S_z分别为x，y，z轴方向上的缩放比例。

旋转变换又可分解为绕x，y，z轴旋转，绕各轴旋转的旋转变换矩阵分别 为：

$R_x=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \alpha & \sin \alpha \\ 0& -\sin \alpha & \cos \alpha \end{array}\right)---\left(4\right)$

$R_y=\left(\begin{array}{ccc}\cos \beta & 0& -\sin \beta \\ 0& 1& 0\\ \sin \beta & 0& \cos \beta \end{array}\right)---\left(5\right)$

$R_z=\left(\begin{array}{ccc}\cos \theta & \sin \theta & 0\\ -\sin \theta & \cos \theta & 0\\ 0& 0& 1\end{array}\right)---\left(6\right)$

其中α，β，θ分别为绕x，y，z轴旋转的角度。

整体旋转是绕x，y，z轴旋转的串联，可表示为：

R＝R_xR_yR_z                (7)

把(2)～(7)式代入(1)式可得：

$T\left(x\right)=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right)x+\left(\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right)---\left(8\right)$

由(8)式可知，仿射变换有12个参数，即μ＝[a₁₁，a₁₂，a₁₃，a₂₁，a₂₂，a₂₃，a₃₁，a₃₂，a₃₃， t_x，t_y，t_z]，所述仿射变换根据这12个参数对浮动图像进行变换。

首次仿射变换时，把[t_x，t_y，t_z，S_x，S_y，S_z，α，β，θ]设置为[0，0，0，1，1，1，0，0，0]，并通过 (1)～(8)式求取a₁₁～a₃₃，然后根据(8)式进行仿射变换；此后的各次迭 代中这12个参数由优化过程求取，并根据(8)式进行仿射变换。

(3.3)对变换后的浮动图像进行B样条插值。

在图像中只能直接获取整数点上的值。当一个变换将一个点从一个空间映 射到另一个空间时，一般情况下会被映射到非整数点区域上。此时，需要通过 插值来计算该点的值。

本发明中的插值方法采用B样条插值法。假设图像f(x)由采样点集f_i＝f(x_i) 描述，其中x_i∈Ω是有整数间距的点，Ω为图像域。利用B样条基函数实现 插值，非整数点上的像素值利用下式插值求得：

$f\left(x\right)=\underset{i}{\Sigma }{c}_i{\beta }^{\left(3\right)}(x-x_i)---\left(9\right)$

式中，x为图像域中任意实值像素位置，x＝[xyz]^T；x_i为整数点位置的坐 标矢量，x_i＝[x_iy_iz_i]^T；c_i为B样条的系数，利用递归滤波器计算；β³(x)为可分 离的三次B样条卷积核，β³(x)＝β³(x)β³(y)β³(z)。其中：

${\beta }^3\left(x\right)=\left(\begin{array}{cc}2/3-{\left|x\right|}^2+{\left|x\right|}^3/2,& 0\le \left|x\right|<1\\ {(2-|x\left|\right)}^3/6,& 1\le \left|x\right|<2\\ 0& \left|x\right|\ge 2\end{array}\right)---\left(10\right)$

(3.4)对参考图像与插值后的浮动图像进行相似性度量。

本发明采用互信息度量法，它是基于信息理论的。互信息通常用于描述两 个系统间的相似性，它可通过熵进行计算。熵是对两个系统间的混乱程度的度 量，互信息越大熵越小，两图像间的相似性越大。

假设参考图像和浮动图像的像素的值可看成两个随机变量F和M，则变 量F和M的熵H(F)、H(M)及其联合熵H(F，M)可表示为：

$H\left(F\right)=-\underset{f}{\Sigma }{P}_F\left(f\right)·\log P_F\left(f\right)---\left(11\right)$

$H\left(M\right)=-\underset{m}{\Sigma }{P}_M\left(m\right)·\log P_M\left(m\right)---\left(12\right)$

$H(F,M)=-\underset{f,m}{\Sigma }{P}_{\mathrm{FM}}(f,m)·\log P_{\mathrm{FM}}(f,m)---\left(13\right)$

其中f∈F，m∈M，P_F(f)和P_M(m)分别是F和M完全独立时的概率分布， P_FM(f，m)是F和M的联合概率分布。

利用下式计算参考图像和浮动图像间的互信息：

I(F，M)＝H(F)+H(M)-H(F，M)        (14)

(3.5)利用优化算法搜寻下一次变换参数。

优化过程就是为了寻找一个最佳参数使得代价函数最小，设F为参考图 像，M为浮动图像，T为参数化变换公式，μ为变换参数矢量，配准的任务就 是寻求最佳参数矢量使得代价函数C最小：

$\hat{\mu }=\arg \underset{\mu }{\min }C(F,\mathrm{MoT})---\left(15\right)$

本发明所用优化算法为自适应随机梯度下降法。假设x为图中的坐标，可 表示为x＝[x，y，z]^T；若配准中的参数化变换公式为T(x，μ)，Ω_F为参考图像的图 像域，采用下式作为代价函数：

$C\left(\mu \right)=\Psi \left(\frac{1}{{|\Omega }_F^{\prime }|}\underset{x_i\in {\Omega }_F^{\prime }}{\Sigma }\xi (F\left(x_i\right),M\left(T(x_i,\mu )\right))\right)---\left(16\right)$

其中为参考图像F中抽样用于计算代价函数的点xi的集合，|Ω′F| 为点的数目，此处取2048个抽样点。

其中：

Ψ(u)＝u                                  (17)

ξ(u，v)＝(u-v)²                          (18)

因此得出代价函数的梯度函数：

$g\left(\mu \right)=\frac{\partial C}{\partial \mu }=\frac{1}{{|\Omega }_F^{\prime }|}\underset{x_i\in {\Omega }_F^{\prime }}{\Sigma }\frac{{\partial T}^T}{\partial \mu }\frac{\partial M}{\partial x}\frac{\partial \xi }{\partial v}\frac{\partial \Psi }{\partial u}---\left(19\right)$

然后利用以下迭代优化策略寻找最佳参数：

μ_k+1＝μ_k-γ(t_k)g(μ_k)，k＝0，1，Λ，K   (20)

t_k+1＝[t_k+f(-g(μ_k)^Tg(μ_k-1))]⁺           (21)

γ(t_k)＝a/(t_k+1+A)^α                      (22)

$f\left(x\right)=f_{\mathrm{MIN}}+\frac{f_{\mathrm{MAX}}-f_{\mathrm{MIN}}}{1-(f_{\mathrm{MAX}}/f_{\mathrm{MIN}})e^{-x/\omega }}---\left(23\right)$

式(20)中g(μ_k)为代价函数C在μ_k处的梯度，γ为步长；式(21)中 [x]⁺＝max(x，0)；其中t₀＝0；f为Sigmoid函数，其中几种Sigmoid函数的曲线如 图3所示。式(22)中a＞0，A≥1，0＜α≤1；式(22)中t_k与当前梯度g(μ_k) 和前一梯度g(μ_k-1)的内积有关，如果两个连续迭代步骤的梯度同向，则内积为 正，-g(μ_k)^Tg(μ_k-1)为负，如图3所示，当-g(μ_k)^Tg(μ_k-1)为负时，f(-g(μ_k)^Tg(μ_k-1)) 为负，因此t_k+1相对于t_k减少。从式(22)可见，当t_k+1减少，则步长γ(t_k+1) 增大。自适应随机梯度下降法正是利用这种方法使得迭代步长得到自动调节。

在式(22)中取A＝20，α＝1。利用下式计算a：

$a=a_{\mathrm{MAX}}\frac{E{\left|\right|g\left|\right|}^2}{E{\left|\right|g\left|\right|}^2+E{\left|\right|{ϵ}_k\left|\right|}^2}---\left(24\right)$

其中ε为逼近代价函数梯度g时产生的误差，E为对g作出N次估计时的 期望，a_MAX、E||g||²、E||ε_k||²利用下式计算：

$a_{\mathrm{MAX}}\equiv \frac{\mathrm{A\delta }}{{\sigma }_1}\underset{x_i\in {\Omega }_F^{\prime }}{\min }{[\mathrm{tr}\left(J_i{\mathrm{CJ}}_i^T\right)+2\sqrt{2}{\left|\right|J_i{\mathrm{CJ}}_i^T\left|\right|}_F]}^{-\frac{1}{2}}---\left(25\right)$

$E{\left|\right|g\left|\right|}^2=\mathrm{tr}\left({\sigma }_1^{2}C\right)---\left(26\right)$

$E{\left|\right|{ϵ}_k\left|\right|}^2=\mathrm{tr}\left({\sigma }_2^{2}C\right)---\left(27\right)$

式(25)中δ为参考图像和浮动图像各方向像素间距的平均值，J_i为空间 变换的雅可比矩阵，利用式(28)进行计算：

$J_i=\frac{\partial T}{\partial \mu }(x_i,\mu )---\left(28\right)$

式(25)～(27)中tr(X)为矩阵X的迹，C利用下式进行计算：

$C\equiv \frac{1}{{|\Omega }_F^{\prime }|}\underset{x_i\in {\Omega }_F^{\prime }}{\Sigma }{J}_i^T{J}_i---\left(29\right)$

式(25)～(27)中σ₁、σ₂利用下式计算：

$\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}{\left|\right|g\left({\mu }_n\right)\left|\right|}^2={\sigma }_1^2\mathrm{tr}\left(C\right)---\left(30\right)$

$\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}{\left|\right|{ϵ}_k\left({\mu }_n\right)\left|\right|}^2={\sigma }_2^2\mathrm{tr}\left(C\right)---\left(31\right)$

式(30)(31)中μ_n是从μ_k附近随机抽取的估计参数，其服从正态分布：

${\mu }_n~N({\mu }_k,{\sigma }_3^{2}I)---\left(32\right)$

其中I为单位矩阵，σ₃利用下式计算：

${\sigma }_3^2=\underset{x_i\in {\Omega }_F^{\prime }}{\min }\frac{{\delta }^2}{{\left|\right|J_i\left|\right|}_F^2+2\sqrt{2}{\left|\right|J_i{J}_i^T\left|\right|}_F}---\left(33\right)$

(33)式中首先求取分母的最大值，当分母最大时即能求得σ3。

式(30)(31)中N为对梯度g作出估计的次数，估计的次数N满足下列 各式：

$E+2\sqrt{\mathrm{Var}}<\mathrm{KE}---\left(34\right)$

$E=E\left(\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}g{\left({\mu }_n\right)}^{T}g\left({\mu }_n\right)\right)={\sigma }_1^2\mathrm{tr}\left(C\right)---\left(35\right)$

$\mathrm{Var}=\mathrm{Var}\left(\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}g{\left({\mu }_n\right)}^{T}g\left({\mu }_n\right)\right)=\frac{{2\sigma }_1^4}{N}{\left|\right|C\left|\right|}_F^2---\left(36\right)$

联立(34)～(36)得：

$N>\frac{{8\left|\right|C\left|\right|}_F^2}{{\left[(K-1)\mathrm{tr}\left(C\right)\right]}^2}---\left(37\right)$

若(37)式右边小于2，则N＝2，否则N利用下式计算：

其中K＝1.5。

把各参数代入(25)式中，并求的最大值，则能 求得a_MAX，再通过式(24)求参数a。

式(23)中f_MAX取1，f_MIN及ω利用下式求取：

$f_{\mathrm{MIN}}=\frac{E{\left|\right|g\left|\right|}^2}{E{\left|\right|g\left|\right|}^2+E{\left|\right|{ϵ}_k\left|\right|}^2}-f_{\mathrm{MAX}}---\left(39\right)$

$\omega =\frac{1}{10}\sqrt{{\sigma }_2^4{\left|\right|C\left|\right|}_F^2}---\left(40\right)$

至此，则可利用式(20)～(23)求取各次迭代的变换参数。

(3.6)重复步骤(3.2)～(3.5)，直至优化算法的迭代次数到达250次。

(3.7)用当前辨率层的参数结果对下一分辨率层的变换参数进行初始化。

(3.8)在所有分辨率层上重复步骤(3.2)～(3.7)，得到基于仿射变换的 粗配准结果和粗配准变换参数，粗配准结果传送至精配准模块400，粗配准变 换参数传送至变换模块500。

第四步，粗配准模块300将粗配准结果传送至精配准模块400，进行精配 准处理，具体步骤如下：

(4.1)将基于仿射变换的粗配准的结果分解为6个不同的分辨率层，并 将其作为新的浮动图像。

(4.2)在当前分辨率层上对新的浮动图像进行B样条变换，并记录变换 参数。

本发明利用B样条变换进行精配准，它在形变中采用了B样条基函数。 由于改变某一控制点的位置时只会影响有限的一段B样条曲线，因此改变局 部控制点时不会影响整体的形变，只对局部形变产生影响，所以能利用它进行 局部精配准。

假设有三维体数据V＝{(x，y，z)|0≤x＜X，0≤y＜Y，0≤z＜Z}，Φ是一组网格控制 点φ_i，j，k的集合，其大小为n_x×n_y×n_z，它们在三个方向上的间距分别为δ_x，δ_y， δ_z。利用三次B样条的三维张量积表示形变模型：

$T(x,y,z)=\underset{l=0}{\overset{3}{\Sigma }}\underset{m=0}{\overset{3}{\Sigma }}\underset{n=0}{\overset{3}{\Sigma }}{B}_l\left(u\right)B_m\left(v\right)B_n\left(w\right){\phi }_{i+l,j+m,k+n}---\left(41\right)$

其中：

(41)式中φ_{i+l，j+m，k+n}为网格点(i+l，j+m，k+n)对应的坐标，它就是B样条变 换的参数μ，把初始参数设置为网格点间距为16mm的网格。

B样条的基函数为：

B₀(u)＝(1-u)³/6

B₁(u)＝(3u³-6u²+4)/6                    (44)

B₂(u)＝(1-3u³+3u²+3u)/6

B₃(u)＝u³/6

(4.3)利用如步骤(3.3)所述的方法对变换后的浮动图像进行B样条插 值。

(4.4)利用如步骤(3.4)所述的方法对参考图像与插值后的浮动图像进 行相似性度量。

(4.5)利用如步骤(3.5)所述的方法通过优化算法搜寻下一次迭代的变 换参数。

(4.6)重复步骤(4.2)～(4.5)，直至优化算法的迭代次数到达500次。

(4.7)用当前辨率层的参数结果对下一分辨率层的变换参数进行初始化。

(4.8)在所有分辨率层上重复步骤(4.2)～(4.7)，得到精配准变换参数， 传送至变换模块500。至此完成了整个配准处理过程。

第五步，变换模块500根据上述步骤得到的变换参数对Label图像进行变 换。变换参数记录了配准过程中的浮动图像的变换过程，因此把它应用于Label 图像即能重现浮动图像的变换过程。将变换结果传送至输出模块600，变换结 果即自动分割结果。

第六步，输出模块600输出自动分割结果，分割结果中不同的标记符则对 应于心脏中不同的部位。

至此，即实现整个基于配准的CT图像全心脏自动分割，并得到如图4至 图11所示的分割结果。