基于快速凸优化算法配准三维CT与超声肝脏图像的方法

摘要

本发明涉及医学图像后处理领域，旨在提供基于快速凸优化算法配准三维CT与超声肝脏图像的方法。该基于快速凸优化算法配准三维CT与超声肝脏图像的方法包括下述过程：将超声和CT图像分辨率调整到相同；对超声和CT图像的基于刚体变换的粗配准；提取多模态图像配准的统一特征信息；计算当前非刚性形变场u(x)下，数据项中的D(u)和D(u)的梯度场对逐步凸优化方法的每一步进行模型求解，得到形变场最优矫正值h(x)，更新形变场，直到h(x)很小；根据求解的非刚性形变场，对超声图像变换，与CT图像配准。本发明通过建立合理的模型，设计出快速、精确的三维超声-CT肝脏图像配准算法，提高消融手术的精确性、安全性和有效性。

主权项

基于快速凸优化算法配准三维CT与超声肝脏图像的方法，其特征在于，具体包括下述过程：(1)将超声和CT图像分辨率调整到相同；(2)对超声和CT图像的基于刚体变换的粗配准；(3)提取多模态图像配准的统一特征信息；(4)构造快速凸优化算法，计算当前非刚性形变场u(x)下，数据项中的D(u)和D(u)的梯度场(5)对逐步凸优化方法的每一步进行模型求解，得到形变场最优矫正值h(x)，更新形变场，重复(3)，(4)直到h(x)很小；(6)根据求解的非刚性形变场，对超声图像变换，与CT图像配准；所述过程(1)具体包括下述步骤：步骤A：获取所需配准的三维超声图像I^U(x)和三维CT图像I^C(x)，然后将三维超声图像I^U(x)和三维CT图像I^C(x)的图像显示窗位调节至[0，N]，再将三维CT图像I^C(x)的分辨率调节至跟三维超声图像I^U(x)一样；所述N是大于0的整数，所述x表示三维图像中的一个点，图像的定义域为Ω；所述过程(2)具体包括下述步骤：步骤B：分别在三维超声图像I^U(x)和三维CT图像I^C(x)中选取4到6对特征点；步骤C：根据特征点，对三维超声图像I^U(x)和三维CT图像I^C(x)进行刚性配准；所述过程(3)具体包括下述步骤：步骤D：提取独立于不同图像模式的统一特征描述MIND来作为图像配准的相似性度量，对于给定的三维CT图像I^C(x)和三维超声图像I^U(x)，计算图像中每一点x∈Ω的MIND特征，分别记为矢量函数C(x)：＝(c₁(x)，...，c_k(x))^T和R(x)：＝(r₁(x)，...，r_k(x))^T；其中，所述I^C(x)表示步骤A处理过的三维CT图像；所述I^U(x)表示步骤A处理过的三维超声图像；所述x表示三维图像中的一个点，图像的定义域为Ω；所述k是计算MIND信息时设定的局域块的维数，表示MIND特征的维数；所述c_k表示CT图像MIND特征第k维的值；所述r_k表示超声图像MIND特征第k维的值；所述T表示向量(c₁(x)，...，c_k(x))^T，(r₁(x)，...，r_k(x))^T中的向量转置；所述C(x)表示三维CT图像I^C(x)的MIND特征矢量函数；所述R(x)表示三维超声图像I^U(x)的MIND特征矢量函数；所述过程(4)具体包括下述步骤：步骤E：计算在当前形变场为u(x)＝(u₁(x)，u₂(x)，u₃(x))^T时，数据项中的常数$D\left(u\right)\left(x\right)=\underset{i=1}{\overset{k}{\Sigma }}|c_i\left(x\right)-r_i(x+u)|,$梯度场${▿}_{2}D\left(u\right)\left(x\right)=({\partial }_{1}D\left(u\right)\left(x\right),{\partial }_{2}D\left(u\right)\left(x\right),{\partial }_{3}D\left(u\right)\left(x\right));$其中，所述u₁(x)，u₂(x)，u₃(x)分别代表图像中每一点x＝(x₁，x₂，x₃)在x，y，z三个方向的形变量；所述T表示向量转置；所述表示从i＝1...k求和；所述c_i(x)表示过程(3)中计算的图像I^C(x)的MIND特征第i维的值；所述r_i(x+u)表示图像I^U(x)经过形变场u(x)作用后的MIND特征第i维的值；x+u表示对超声图像I^U(x)作用一个形变场u(x)，即每一点x在x，y，z方向加一个形变，得到新的点的位置(x₁(x)+u₁(x)，x₂(x)+u₂(x)，x₃(x)+u₃(x))；所述表示D(u)(x)分别对x，y，z三个方向求偏导，其中符号表示偏导算子；所述过程(5)具体包括下述步骤：步骤F：采取逐步凸优化方法，每一步的具体过程中固定非刚体形变场现值u(x)，优化后的能量泛函是凸的，即得到了一个凸模型为以下形式：${\min }_h{\int }_{\Omega }|D\left(u\right)+{▿}_{u}D·h|\mathrm{dx}+\frac{\alpha }{2}\underset{i=1}{\overset{3}{\Sigma }}{\int }_{\Omega }{|▿(u_i+h_i)|}^2\mathrm{dx}---\left(1\right)$其中，所述h＝(h₁(x)，h₂(x)，h₃(x))^T表示需要求解的最优矫正形变场；凸模型中第一项表示数据项，其中的D(u)和由步骤E已经算得；第二项表示对形变场的光滑正则化；参数α为一个大于0的常数，用于调节数据项和正则项的比重；所述符号表示关于u求梯度；符号表示梯度算子；符号min_h表示关于h求极小值；符号| |表示求绝对值；符号∫_Ω表示在图像区域Ω内求积分；符号dx表示体积元；h₁(x)，h₂(x)，h₃(x)分别表示在图像中每点在x，y，z三个方向的形变量；T表示向量转置；符号·表示向量相乘；符号表示从i＝1，2，3求和；通过原始‑对偶算法将模型(1)的凸模型转化为以下形式：${\min }_h{\max }_{w,q}{L}_c(h,w,q):={\int }_{\Omega }(\mathrm{wD}\left(u\right)+\underset{i=1}{\overset{3}{\Sigma }}{u}_i\mathrm{div}{q}_i)\mathrm{dx}-\frac{1}{2\alpha }\underset{i=1}{\overset{3}{\Sigma }}{\int }_{\Omega }{q_i}^2\mathrm{dx}+\underset{i=1}{\overset{3}{\Sigma }}<h_i,F_i>-\frac{c}{2}\underset{i=1}{\overset{3}{\Sigma }}{\left|\right|F_i\left|\right|}^2$约束于$\left|w\left(x\right)\right|\le 1;F_i\left(x\right):=(w·{\partial }_{u_i}D+{\mathrm{div}q}_i)\left(x\right)---\left(2\right)$$\forall x\in \Omega ,i=\mathrm{1,2,3};$其中，用L_c(h，w，q)表示定义的模型的能量函数，所述w，q＝(q₁，q₂，q₃)表示模型(2)中的对偶变量；所述h表示拉格朗日乘子，即在模型(1)中需要求解的最优矫正形变场；所述α就是模型(1)中用于调节数据项和正则项的参数，是一个大于0的常数；所述c是一个大于0的常数；D(u)，由步骤E已经算得，u＝(u₁，u₂，u₃)是已知的形变场现值；所述符号min_h表示关于h求极小值；max_w，q表示关于w，q求极大值；符号：＝表示“定义为”；符号div表示散度算子；符号∫_Ω表示在图像区域Ω内求积分；符号dx表示体积元；符号表示从i＝1，2，3求和；符号<，>表示求内积；符号表示对u_i求偏导数；符号|| ||²表示求L2‑范数，符号| |表示求绝对值；求解模型(1)等价于求解模型(2)，所以通过优化模型(2)就可以得到模型(1)需要的h(x)；所述过程(6)具体包括下述步骤：步骤G：更新非刚体形变场现值u(x)＝u(x)+h(x)，x∈Ω，对旧的形变场现值加上求得的形变场矫正值，得到新的形变场值；步骤H：根据求解得到的非刚体形变场，对超声图像I^U(x)进行形变I^U(x+u(x))，超声与CT图像即得到配准；I^U(x+u)＝I^U(x₁(x)+u₁(x)，x₂(x)+u₂(x)，x₃(x)+u₃(x))，(x₁(x)，x₂(x)，x₃(x))为像素点x的坐标，(u₁(x)，u₂(x)，u₃(x))为该像素点x在x，y，z方向的形变量。