基于增强弹性形变的CARTO电解剖图与CT图像配准方法和装置

摘要

本发明公开了一种基于增强弹性形变的CARTO电解剖图与CT图像的配准方法和装置，包括：获取CARTO电解剖图与CT图像；采用刚体变换模型，以基于主轴的方法，粗配准CARTO电解剖图与CT图像；在粗配准的基础上，以距离平方和为相似性度量函数，采用增强自由形变模型和金字塔分层方法，精配准CARTO电解剖图与CT图像；将精配准形变后的CT图像叠加在CARTO电解剖图上，得到CARTO电解剖图与CT图像的最终配准图像。本发明的配准方法和装置为房颤消融手术中精确快速导航消融导管提供了保障。

说明书

技术领域

本发明涉及手术导航中图像配准技术领域，尤其涉及心房颤动(房颤)射 频消融介入手术导航中的CARTO电解剖图(Electro-anatomicalMap,EAM) 与CT图像的配准方法和装置。

背景技术

CARTO系统指导下的房颤消融是当前治疗房颤最常用也最为有效的手 段。手术时，医生将导管插入心房，依X线透视在心腔内壁依次标测近百个 点，模拟构建心腔内表面三维解剖模型，即CARTO电解剖图(electro-anatomic map，EAM)；医生在CARTO电解剖图上导航导管，设计并定位消融线；导管 头释放电能，逐点消融，把导致房颤的纤维组织烧掉。由于标测点有限，CARTO 电解剖图分辨率低，而心房解剖结构复杂，医生很难在CARTO电解剖图上准 确定位并有效隔离肺静脉，因而极大影响消融疗效。若术前采集清晰的心房 CT图像，从中分割出作为消融靶区的心腔内表面(简称CT曲面)，然后将 CARTO系统采集的近百个标测点叠加到CT曲面对应位置上，实现两者在空 间位置上最佳叠合，即配准CARTO电解剖图与CT曲面，便可帮助医生从多 个角度仔细观察消融靶区的解剖结构，从而增加消融位点的精确性，减少手术 并发症，缩短手术和X线曝光时间。

然而，现有的CARTO电解剖图与CT图像配准方法都不能很好地模拟实 际CARTO电解剖图与CT曲面之间的形变差异。我们知道，CT图像采集于 术前，CARTO电解剖图采集于术中，病人的躺卧姿势、呼吸深浅、心跳状况、 心脏状态等，均存在差异，使得最终生成的CARTO电解剖图与CT曲面形态 差异很大。同时，消融手术中，导管接触心腔内壁的力度、角度随着医生手术 时操作的特点不同而有所不同，心脏的形态也会随之产生不同程度的扭曲。因 此，现有方法常用的刚体变换和仿射变换，都不能很好地模拟实际CARTO电 解剖图与CT曲面之间的形变差异，即使偶有使用弹性形变模型，其形变量也 不足，致使配准精度难以满足临床实际应用需求。

发明内容

有鉴于此，本发明提供了一种基于增强弹性形变的CARTO电解剖图与 CT图像配准方法和装置，以更为精确细致地描述CARTO电解剖图与CT图 像之间的形变差异，从而提高CARTO电解剖图与CT图像的配准精度。

为了达到上述发明目的，本发明采用了如下技术方案：

一种基于增强弹性形变的CARTO电解剖图与CT图像的配准方法，包括：

获取CARTO电解剖图与CT图像，所述CT图像是对CT原始图像进行 重建分割之后获得的心房内壁的曲面图像；

采用刚体变换模型，以基于主轴的方法，粗配准CARTO电解剖图与CT 图像；粗配准过程中，CT图像为基准图像，CARTO电解剖图为浮动图像；

在粗配准的基础上，以距离平方和为相似性度量函数，采用增强自由形变 模型和金字塔分层方法，精配准CARTO电解剖图与CT图像；精配准过程中， CARTO电解剖图为基准图像，CT图像为浮动图像；

将精配准形变后的CT图像叠加在CARTO电解剖图上，得到CARTO电 解剖图与CT图像的最终配准图像。

基于上述所述的配准方法，本发明还提供了基于增强弹性形变的CARTO 电解剖图与CT图像的配准装置，其包括：

获取单元，用于获取CARTO电解剖图与CT图像，所述CT图像是对CT 原始图像进行重建分割之后获得的心房内壁的曲面图像；

粗配准单元，用于采用刚体变换模型，以基于主轴的方法，粗配准CARTO 电解剖图与CT图像；粗配准过程中，CT图像为基准图像，CARTO电解剖图 为浮动图像；

精配准单元，用于在粗配准的基础上，以距离平方和为相似性度量函数， 采用增强自由形变模型和金字塔分层方法，精配准CARTO电解剖图与CT图 像；精配准过程中，CARTO电解剖图为基准图像，CT图像为浮动图像；

叠加单元，用于将精配准形变后的CT图像叠加在CARTO电解剖图上， 得到CARTO电解剖图与CT图像的最终配准图像。

相较于现有技术，本发明具有以下有益效果：

本发明提供的基于增强弹性形变的CARTO电解剖图与CT图像的配准方 法，采用增强自由形变模型和金字塔分层技术，可以对或大或小的弹性形变建 模，构建的形变场平滑、连续，能够保持拓扑结构和一一对应关系，能更为精 确细致地描述CARTO电解剖图与CT图像之间的形变差异，从而大幅提高 CARTO电解剖图与CT图像的配准精度；同时，金字塔分层技术能减少因弹 性形变产生的大的计算量，在保证配准精度的同时，也能获得较快的配准速度。 通过本发明的配准方法和装置得到的配准图像，能够为房颤消融手术中精确快 速导航消融导管提供保障。

附图说明

为了清楚地理解本发明的技术方案，下面对描述本发明具体实施方式时用 到的附图进行简要说明。显而易见，这些附图仅是本发明的一部分附图，本领 域普通技术人员在不付出创造性劳动的前提下，还可以获得其它附图。

图1是本发明实施例提供的基于增强弹性形变的CARTO电解剖图与CT 图像的配准方法流程示意图；

图2是本发明实施例提供的CARTO电解剖图与CT图像的粗配准方法流 程示意图；

图3是本发明实施例提供的CARTO电解剖图与CT图像的精配准方法流 程示意图；

图4是本发明实施例提供的获取每个控制点形变量的方法流程示意图；

图5是本发明实施例提供的由形变后的控制点位置，确定覆盖其下的CT 图像中各点位置的方法流程示意图；

图6是本发明实施例提供的基于增强弹性形变的CARTO电解剖图与CT 图像的配准装置结构示意图。

具体实施方式

为使本发明的目的、技术手段和有益效果更加清楚完整，下面结合附图对 本发明的具体实施方式进行描述。

图1是本发明实施例一提供的基于增强弹性形变的CARTO电解剖图与CT 图像的配准方法流程示意图。如图1所示，该配准方法包括以下步骤：

S11、获取CARTO电解剖图与CT图像，所述CT图像是对CT原始图像 进行重建分割之后获得的心房内壁的曲面图像。

S12、采用刚体变换模型，以基于主轴的方法，粗配准CARTO电解剖图 与CT图像；粗配准过程中，CT图像为基准图像，CARTO电解剖图为浮动图 像：

图2示出了步骤S12的具体实现过程。如图2所示，本发明实施例所述的粗 配准过程具体包括以下步骤：

S121、计算CT图像的质心C_CT、惯性矩阵I_CT和特征矩阵E_CT；计算CARTO 电解剖图EAM的质心C_EAM、惯性矩阵I_EAM和特征矩阵E_EAM。

S122、计算与CT图像具有主轴对应关系的CARTO电解剖图EAM’：

惯性矩阵中的每一个特征向量对应一个主轴，坐标系中的任何坐标轴旋转 180度，并不会影响图像的惯性矩阵。因此，设定CT图像的特征矩阵E_CT中的 特征向量排列固定不变，与之对应，构成CARTO电解剖图的特征矩阵E_EAM的 特征向量排列方式有种，即CARTO电解剖图的特征矩阵E_EAM共 有48种。对每一个E_EAM,分别计算对应的CARTO电解剖图EAM’，计算公式 如下：

${\mathrm{EAM}}^{\prime }=E_{CT}*E_{EAM}^T*(I_{EAM}-C_{EAM})+C_{CT}---\left(1\right)$

其中，是特征矩阵E_EAM的转置矩阵。

S123、计算并比较48个CARTO电解剖图EAM’到所述CT图像的距离， 与CT图像距离最小的CARTO电解剖图EAM’即是粗配准后的CARTO电解 剖图EAM*。

在本发明实施例中，CARTO电解剖图EAM’到所述CT图像的距离可以 定义为所有CARTO标测点到CT图像最近点欧式距离的均方根。因为仅需遍 历48种主轴的对应关系，寻找其中使得两图距离最小的一种作为最终的配准 结果，所以粗配准算法中不涉及任何优化搜索算法。

S13、在粗配准的基础上，以距离平方和为相似性度量函数，采用增强自 由形变模型和金字塔分层方法，精配准CARTO电解剖图与CT图像，得到形 变后的CT图像；精配准过程中，CARTO电解剖图为基准图像，CT图像为浮 动图像：

需要说明的是，以自由形变模型(FreeFormDeformation,FFD)配准 CARTO电解剖图与CT图像的基本思路是：将由M×N×O个控制点P_m,n,o构成的 控制网格覆盖在CT图像所在空间，调整形变参数，操控由控制点组成的控制 网格，实现对CT图像形状的改变。

然而，FFD将每个控制点的最大位移约束在两个控制点间隔的一半以内， 所以其形变自由度较低，能描述的形变量有限。为了实现更为精确细致的 CARTO电解剖图与CT图像配准，本发明采用增强自由形变模型(Incremental FreeFormDeformation，IFFD)，增强自由形变模型为每个控制点P_m,n,o增加了 形变量δP_m,n,o；进一步，本发明采用了基于B样条的局部增强自由形变模型。

如此，每个控制点的空间位置满足以下关系：

P_m,n,o＝P_m,n,o⁰+δP_m,n,o(2)

其中，P_m,n,o⁰是控制点P_m,n,o未发生形变时的空间位置， (m,n,o)∈[1,M]×[1,N]×[1,O]；δP_m,n,o是控制点的形变量， δP_m,n,o＝{(δP^x_m,n,o,δP^y_m,n,o,δP^z_m,n,o)}。

依据上述定义，给定控制点从P_m,n,o⁰形变到P_m,n,o，形变量为δP_m,n,o，控制网 格覆盖下的CT图像各点Y＝(x,y,z)亦随之发生形变，形变量为 点Y的最终空间位置为L(Y)＝Y+δL(Y)。

这里，

表示向下取整，B_k(u)表示B样条的第k个基函数， δP_i+k,j+l,q+r((k,l,r)∈[0,3]×[0,3]×[0,3])是覆盖在CT图像点Y＝(x,y,z)上的相邻64个 控制点的形变量。

从上所述，CT图像中任意点Y的位置，完全由覆盖其上的相邻64个控制点 的形变量δP_m,n,o决定，改变δP_m,n,o即可实现CT图像和CARTO电解剖图之间的配 准。因此，基于IFFD的CARTO电解剖图与CT图像配准，即等同于寻找控制网 格的形变量δP_m,n,o，使得形变后的CT图像与CARTO电解剖图空间位置达到一 致。

本发明实施例提供了步骤S13的具体实现方式。图3是本发明实施例中精配 准CARTO电解剖图与CT图像的方法流程示意图，如图3所示，该精配准方法 包括以下步骤：

S131、设定控制网格分辨率M×N×O和金字塔层数s的初始值，令 M＝N＝O＝3，s＝1：

为了减少弹性形变产生的大计算量，本发明实施例采用了金字塔分层技 术，随着金字塔层数的增加，其控制网格的分辨率逐渐增加。

S132、对CARTO电解剖图每一个点X，运用迭代最近点方法寻找CT图 像上与之对应的点L(δP_m,n,o；Χ)；其中，CARTO电解剖图位于空间Ω中， Ω＝{X}＝{(x,y,z)|1≤x≤X,1≤y≤Y,1≤z≤Z}；

S133、构建由M×N×O个控制点均匀分布构成的控制网格P，将控制网格 P覆盖在CT图像空间Γ之上：

所述控制网格P具体描述如下：

P＝{P_m,n,o}＝{(P^x_m,n,o,P^y_m,n,o,P^z_m,n,o)}((m,n,o)∈[1,M]×[1,N]×[1,O])

本发明实施例采用增强自由形变模型，每个控制点的空间位置满足以下关 系：

P_m,n,o＝P_m,n,o⁰+δP_m,n,o(4)

P_m,n,o为控制网格的控制点发生形变后的空间位置；

P_m,n,o⁰是控制点P_m,n,o未发生形变时的空间位置；

δP_m,n,o是控制点P_m,n,o的形变量，δP_m,n,o＝{(δP^x_m,n,o,δP^y_m,n,o,δP^z_m,n,o)}；令控制点 的初始形变量δP_m,n,o＝0；

S134、以梯度下降法优化目标函数，获得每个控制点P_m,n,o的形变量δP_m,n,o， 然后确定形变后控制点的位置：

本发明的目标函数的表达式如下：

E(δP_m,n,o)＝E_similarity(δP_m,n,o)+αE_smoothness(δP_m,n,o)(5)

E_similarity(δP_m,n,o)为相似性度量函数，本发明以距离平方和度量相似性时，该 相似性度量函数的表达式为：

$F_{similarity}\left({\mathrm{\delta P}}_{m,n,o}\right)=\underset{\Omega }{\int \int \int }{(X-L({\mathrm{\delta P}}_{m,n,o};X\left)\right)}^{2}dX---\left(6\right)$

X表示CARTO电解剖图中的点，L(δP_m,n,o；Χ)表示CT图像中与X相对应 的点。

E_smoothness(δP_m,n,o)为平滑约束项，用于避免出现形状扭曲，保持配准网格的 规整性。E_smoothness(δP_m,n,o)的计算公式为：

$E_{smoothness}=\underset{\Omega }{\int \int \int }\left(\right|\left|\frac{\partial \delta L({\mathrm{\delta P}}_{m,n,o};X)}{\partial x}\right|{|}^2+\left|\right|\frac{\partial \delta L({\mathrm{\delta P}}_{m,n,o};X)}{\partial y}|{|}^2+|\left|\frac{\partial \delta L({\mathrm{\delta P}}_{m,n,o};X)}{\partial z}\right|{|}^2)dX---\left(7\right)$

α是平衡E_simialrity和E_smoothness两项的一个常量，α越小，收敛越快；α越大， 则形变区域越平滑。

图4进一步显示了S134的具体过程如下：

S1341、对控制网格上的每一个控制点的形变量δP_m,n,o，根据公式(8)计 算目标函数的梯度

$▿E=\frac{\partial E\left({\mathrm{\delta P}}_{m,n,o}\right)}{\partial {\mathrm{\delta P}}_{m,n,o}}---\left(8\right)$

其中，

$\begin{array}{c}\frac{\partial E\left({\mathrm{\delta P}}_{m,n,o}\right)}{\partial {\mathrm{\delta P}}_{m,n,o}}=-2\underset{\Omega }{\int \int \int }(X-L({\mathrm{\delta P}}_{m,n,o};X\left)\right)(▿L({\mathrm{\delta P}}_{m,n,o};X)·\frac{\partial \delta L({\mathrm{\delta P}}_{m,n,o};X)}{\partial {\mathrm{\delta P}}_{m,n,o}})dX+\\ 2\alpha \underset{\Omega }{\int \int \int }(\frac{\partial \delta L({\mathrm{\delta P}}_{m,n,o};X)}{\partial x}·\frac{\partial }{\partial {\mathrm{\delta P}}_{m,n,o}}(\frac{\partial \delta L({\mathrm{\delta P}}_{m,n,o};X)}{\partial x})+\frac{\partial \delta L({\mathrm{\delta P}}_{m,n,o};X)}{\partial y}·\frac{\partial }{\partial {\mathrm{\delta P}}_{m,n,o}}(\frac{\partial \delta L({\mathrm{\delta P}}_{m,n,o};X)}{\partial y}\\ \frac{\partial \delta L({\mathrm{\delta P}}_{m,n,o};X)}{\partial z}·\frac{\partial }{\partial {\mathrm{\delta P}}_{m,n,o}}\left(\frac{\partial \delta L({\mathrm{\delta P}}_{m,n,o};X)}{\partial z}\right))dX\end{array})+---\left(9\right)$

S1342、判断是否成立，如果否，按照公式(10)更新该控制点形 变量δP_m,n,o，并返回执行步骤S1341；如果是，则δP_m,n,o为该控制点P_m,n,o的形变 量δP_m,n,o，输出δP_m,n,o；

${\mathrm{\delta P}}_{m,n,o}={\mathrm{\delta P}}_{m,n,o}+\mu \frac{▿E}{\left|\right|▿E\left|\right|}---\left(10\right)$

其中，μ表示每一步沿梯度向量迭代的步幅，ε表示阈值。

S135、由形变后的控制点位置，确定覆盖其下的CT图像中各点的位置。

图5显示S135的具体过程如下：

S1351、根据所有控制点P_m,n,o的形变量δP_m,n,o，由B样条的自由形变模型计 算得出CT图像上每一个点Y＝(x,y,z)的形变量δL(Y)，计算公式如下

$\delta L\left(Y\right)=\underset{k=0}{\overset{3}{\Sigma }}\underset{l=0}{\overset{3}{\Sigma }}\underset{r=0}{\overset{3}{\Sigma }}{B}_k\left(u\right)B_l\left(v\right)B_r\left(w\right){\mathrm{\delta P}}_{i+k,j+l,q+r}---\left(11\right)$

B_k(u)表示B样条的第k个基函数，δP_i+k,j+l,q+r((k,l,r)∈[0,3]×[0,3]×[0,3])是覆盖 在CT图像点Y＝(x,y,z)上的相邻64个控制点的形变量；

S1352、根据CT图像上每一个点Y＝(x,y,z)的形变量δL(Y)，确定形变 后CT图像各点Y的最终空间位置L(Y)，其中，L(Y)＝Y+δL(Y)。

S136、更新控制网格的分辨率以及金字塔层数，令M＝M+1，N＝N+1， O＝O+1，s＝s+1。

S137、循环执行步骤S132至步骤S136，直至s＝s_阈值：

需要说明的是，在本发明实施例中，s_阈值可以根据实际需要选取。选择s_阈值需要考虑计算量和配准精度之间的平衡，一般情况下，s_阈值选择4或5即可。

S14、将精配准形变后的CT图像叠加在CARTO电解剖图上，得到CARTO 电解剖图与CT图像的最终配准图像。

以上即为本发明实施例提供的基于增强弹性形变的CARTO电解剖图与 CT图像的配准方法的具体实施方式。在该配准方法中，采用增强自由形变模 型和金字塔分层技术，可以对或大或小的弹性形变建模，构建的形变场平滑、 连续，能够保持拓扑结构和一一对应关系，能更为精确细致地描述CARTO电 解剖图与CT图像之间的形变差异，从而大幅提高CARTO电解剖图与CT图 像的配准精度；同时，金字塔分层技术能减少因弹性形变产生的大的计算量， 在保证配准精度的同时，也能获得较快的配准速度。通过本发明提供的配准方 法得到的配准图像，能够为房颤消融手术中精确快速导航消融导管提供保障。

基于上述实施例提供的基于增强弹性形变的CARTO电解剖图与CT图像的 配准方法，本发明实施例还提供了基于增强弹性形变的CARTO电解剖图与CT 图像的配准装置，具体参见以下实施例。

图6为本发明实施例提供的基于增强弹性形变的CARTO电解剖图与CT图 像的配准装置的结构示意图。如图6所示，该配准装置包括以下单元：

获取单元61，用于获取CARTO电解剖图与CT图像，所述CT图像是对 CT原始图像进行重建分割之后获得的心房内壁的曲面图像；

粗配准单元62，用于采用刚体变换模型，以基于主轴的方法，粗配准 CARTO电解剖图与CT图像；粗配准过程中，CT图像为基准图像，CARTO 电解剖图为浮动图像；

精配准单元63，用于在粗配准的基础上，以距离平方和为相似性度量函 数，采用增强自由形变模型和金字塔分层方法，精配准CARTO电解剖图与 CT图像；精配准过程中，CARTO电解剖图为基准图像，CT图像为浮动图像；

叠加单元64，用于将精配准形变后的CT图像叠加在CARTO电解剖图上， 得到CARTO电解剖图与CT图像的最终配准图像。

本发明提供的基于增强弹性形变的CARTO电解剖图与CT图像的配准装 置，采用增强自由形变模型和金字塔分层技术，可以对或大或小的弹性形变建 模，构建的形变场平滑、连续，能够保持拓扑结构和一一对应关系，能更为精 确细致地描述CARTO电解剖图与CT图像之间的形变差异，从而大幅提高 CARTO电解剖图与CT图像的配准精度；同时，金字塔分层技术能减少因弹 性形变产生的大的计算量，在保证配准精度的同时，也能获得较快的配准速度。 通过本发明提供的配准装置得到的配准图像，能够为房颤消融手术中精确快速 导航消融导管提供保障。

作为本发明的一个具体实施例，所述增强自由形变模型为基于B样条的 增强自由形变模型。

作为本发明的一个具体实施例，所述精配准单元63包括：

初始值设定子单元631，用于设定控制网格分辨率M×N×O和金字塔层 数s的初始值，令M＝N＝O＝3，s＝1；

寻找同名点对子单元632，用于对CARTO电解剖图每一个标测点X，运 用迭代最近点方法确定CT图像上与之对应的点L(δP_m,n,o；Χ)；CARTO电解剖图 位于空间Ω中，Ω＝{X}＝{(x,y,z)|1≤x≤X,1≤y≤Y,1≤z≤Z}；

构建控制网格子单元633，用于将由M×N×O个控制点P_m,n,o均匀分布构 成的控制网格P覆盖在CT图像空间Γ之上；

所述控制网格P描述如下：

P＝{P_m,n,o}＝{(P^x_m,n,o,P^y_m,n,o,P^z_m,n,o)}((m,n,o)∈[1,M]×[1,N]×[1,O])

其中，P_m,n,o＝P_m,n,o⁰+δP_m,n,o；

P_m,n,o⁰是控制点P_m,n,o未发生形变时的空间位置；

δP_m,n,o是控制点P_m,n,o的形变量，δP_m,n,o＝{(δP^x_m,n,o,δP^y_m,n,o,δP^z_m,n,o)}；令控制点 的初始形变量δP_m,n,o＝0；

P_m,n,o是控制点发生形变后的空间位置；

计算控制点形变量子单元634，用于以梯度下降法优化目标函数，获得每 个控制点P_m,n,o的形变量δP_m,n,o，然后确定形变后控制点的位置；

确定CT图像形变后各点位置子单元635，用于由形变后的控制点位置， 确定覆盖其下的CT图像中各点的位置；

更新金字塔层子单元636，用于更新控制网格的分辨率及金字塔分层的层 数，令M＝M+1，N＝N+1，O＝O+1，s＝s+1；并返回寻找同名点对子单元，依 次执行，直至s＝s_阈值。

作为本发明的一个具体实施例，需要说明的是，本发明的目标函数的表 达式如下：

E(δP_m,n,o)＝E_similarity(δP_m,n,o)+αE_smoothness(δP_m,n,o)(5)

E_similarity(δP_m,n,o)为相似性度量函数，采用距离平方和度量两图相似度，该相 似性度量函数的表达式为：

$F_{similarity}\left({\mathrm{\delta P}}_{m,n,o}\right)=\underset{\Omega }{\int \int \int }{(X-L({\mathrm{\delta P}}_{m,n,o};X\left)\right)}^{2}dX---\left(6\right)$

X表示CARTO电解剖图中的点，L(δP_m,n,o；Χ)表示CT图像中与X相对应 的点。

E_smoothness(δP_m,n,o)为平滑约束项，用于避免出现形状扭曲，保持配准网格的 规整性，E_smoothness(δP_m,n,o)的计算公式为：

$E_{smoothness}=\underset{\Omega }{\int \int \int }\left(\right|\left|\frac{\partial \delta L({\mathrm{\delta P}}_{m,n,o};X)}{\partial x}\right|{|}^2+\left|\right|\frac{\partial \delta L({\mathrm{\delta P}}_{m,n,o};X)}{\partial y}|{|}^2+|\left|\frac{\partial \delta L({\mathrm{\delta P}}_{m,n,o};X)}{\partial z}\right|{|}^2)dX---\left(7\right)$

α是平衡E_simialrity和E_smoothness两项的一个常量，α越小，收敛越快；α越大， 则形变区域越平滑。

作为本发明的一个具体实施例，所述计算控制点形变量子单元634包括：

第一计算子单元6341，用于对控制网格上每一个控制点的形变量δP_m，n，o， 根据以下公式计算目标函数的梯度

$▿E=\frac{\partial E\left({\mathrm{\delta P}}_{m,n,o}\right)}{\partial {\mathrm{\delta P}}_{m,n,o}}---\left(8\right)$

判断子单元6342，用于判断是否成立，如果否，按照公式(10) 更新该控制点形变量δP_m,n,o，并将更新后的控制点形变量δP_m,n,o返回第一计算子 单元；如果是，输出控制点P_m,n,o的形变量δP_m,n,o；

${\mathrm{\delta P}}_{m,n,o}={\mathrm{\delta P}}_{m,n,o}+\mu \frac{▿E}{\left|\right|▿E\left|\right|}---\left(10\right)$

其中，μ表示每一步沿梯度向量迭代的步幅，ε表示阈值。

作为本发明的另一具体实施例，所述确定CT图像形变后各点位置子单元 635包括：

第二计算子单元6351，用于根据所有控制点P_m,n,o的形变量δP_m,n,o，由B样 条的自由形变模型计算得出CT图像上每一个点Y＝(x,y,z)的形变量δL(Y)， 计算公式如下：

$\delta L\left(Y\right)=\underset{k=0}{\overset{3}{\Sigma }}\underset{l=0}{\overset{3}{\Sigma }}\underset{r=0}{\overset{3}{\Sigma }}{B}_k\left(u\right)B_l\left(v\right)B_r\left(w\right){\mathrm{\delta P}}_{i+k,j+l,q+r}---\left(11\right)$

这里：

表示向下取整，B_k(u)表示B样条的第k个基函数， δP_i+k,j+l,q+r((k,l,r)∈[0,3]×[0,3]×[0,3])是覆盖在CT图像点Y＝(x,y,z)上的相邻64个 控制点的形变量；

确定子单元6352，用于根据CT图像上每一个点Y＝(x,y,z)的形变量δL(Y)， 确定形变后的CT图像各点Y的最终空间位置L(Y)，其中，L(Y)＝Y+δL(Y)。

以上所述仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限 制。虽然本发明以较佳实施例揭露如上，然而并非用以限定本发明。任何熟悉 本领域的技术人员，在不脱离本发明技术方案范围情况下，都可利用上述揭示 的方法和技术内容对本发明技术方案做出许多可能的变动和修饰，或修改为等 同变化的等效实施例。因此，凡是未脱离本发明技术方案的内容，依据本发明 的技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均仍属于本 发明技术方案保护的范围内。