基于流形正则化的血管3D/2D弹性配准方法及装置与流程

本发明涉及医学图像处理的技术领域，尤其涉及一种基于流形正则化的血管3d/2d弹性配准方法，以及基于流形正则化的血管3d/2d弹性配准装置。

背景技术：

血管3d/2d配准目的是，在术前的3d图像和术中的2d图像中找到对应血管解剖结构的对应关系，并估计3d血管的刚性和弹性变换，使得3d和2d血管在对应结构上达到空间位置的一致。3d和2d血管分别从cta(计算机断层扫描血管造影)和xra(x射线血管造影)图像中提取，经3d/2d配准后，3d和2d血管在对应解剖结构上达到空间位置的一致。根据3d/2d配准任务的不同，对其进行任务分解，包括：初始化、刚性变换和弹性变换。

弹性变换中的形变模型可以分为两类：一类是基于大量训练数据统计的形变模型，如统计形状模型(ssm)、主动表观模型(aam)和可形变模型(morphablemodel)；另一类是基于平滑先验假设的形变模型，如薄板样条模型(tps)、b样条模型和相干点漂移模型(cpd)。

针对血管配准问题，许多研究者在以上这些模型的基础上做出了适用于血管模型的改进。groher等人于2009年提出了一种基于单视角的血管形变模型，在构建的能量函数中对血管长度变化和非自然弯曲进行惩罚，通过求解目标函数对3d血管点位移的偏导构建目标函数梯度场，利用bfgs优化方法求得最优的弹性变换参数。2010年，groher等人改进了方法，将配准问题表示为包含以图像信息响应为主的外能和以3d血管长度变化约束为主的内能的能量最小化问题，该方法直接利用2d造影图像的管状结构增强，无需对2d图像进行血管分割。liao等人将所有血管点视为自由点，在腹主动脉血管的图模型的基础上，构建以距离测度、长度保持项和平滑正则项为主的相似性测度，使用多尺度优化策略求取血管点的位移值。rivest-henault等人使用了与类似的内能和外能构建的弹性配准能量函数，不同的是作者利用了血管分叉点处的形变约束来构建外能，并在双视图的冠脉造影图像上求解冠脉3d中心线的位移值。liu等人利用2dtps函数实现投影后的3d血管点到2d血管中心线的非刚性配准。serradell等人针对心脏运动的先验知识，提出了一种合成的冠脉血管树形变模型。利用该先验随机的生成大量的血管中心线训练集，并使用ssm来构建血管形变模型。baka等人则在真实的临床冠脉4dcta数据上，提取出3d+t冠脉中心线，并构建真实冠脉的ssm。但是，这些改进的血管模型均只适用于特定的具有拓扑结构的点云结构。

技术实现要素：

为克服现有技术的缺陷，本发明要解决的技术问题是提供了一种基于流形正则化的血管3d/2d弹性配准方法，其能够适用于任意具有拓扑结构的点云结构，可以得到最优的血管3d/2d弹性配准结果。

本发明的技术方案是：这种基于流形正则化的血管3d/2d弹性配准方法，其包括以下步骤：

(1)在已知血管3d和2d点的匹配关系基础上，表示3d点，表示2d点，i和j分别为点集的索引值；用表示3d点集矩阵，其中为3d点坐标，用表示2d点集矩阵，其中为2d点坐标，两个点集中对应点具有相同的索引值，k为点数量；

(2)构建空间变换使得3d点的变换表示为τ(yk)为yk的位移；

(3)在血管拓扑流形上的约束描述为在拓扑上相邻的点的位移具有一致性，用流形正则化表示；

(4)获得血管3d/2d弹性配准的目标函数；

(5)基于梯度优化的算法对目标函数进行求解。

本发明提出一种基于流形正则化的血管弹性形变模型，该形变模型利用了血管点的位移在其血管拓扑上具有平滑约束的先验，因此该模型适用于任意具有拓扑结构的点云结构；采用基于梯度优化的算法对目标函数进行求解，可以得到最优的血管3d/2d弹性配准结果。

还提供了基于流形正则化的血管3d/2d弹性配准装置，其包括：

初始化模块，其在已知血管3d和2d点的匹配关系基础上，用表示3d点集矩阵，其中为3d点坐标，用表示2d点集矩阵，其中为2d点坐标，两个点集中对应点具有相同的索引值，k为点数量；

构建模块，其构建空间变换使得3d点的变换表示为τ(yk)为yk的位移；

流形正则化模块，其在血管拓扑流形上的约束描述为在拓扑上相邻的点的位移具有一致性，用流形正则化表示；

目标函数获得模块，其获得血管3d/2d弹性配准的目标函数；

求解模块，其基于梯度优化的算法对目标函数进行求解。

附图说明

图1示出了根据本发明的基于流形正则化的血管3d/2d弹性配准方法的流程图。

具体实施方式

如图1所示，这种基于流形正则化的血管3d/2d弹性配准方法，其包括以下步骤：

(1)在已知血管3d和2d点的匹配关系基础上，表示3d点，表示2d点，i和j分别为点集的索引值；用表示3d点集矩阵，其中为3d点坐标，用表示2d点集矩阵，其中为2d点坐标，两个点集中对应点具有相同的索引值，k为点数量；

(2)构建空间变换使得3d点的变换表示为τ(yk)为yk的位移；(为作用符，等同于)

(3)在血管拓扑流形上的约束描述为在拓扑上相邻的点的位移具有一致性，用流形正则化表示；

(4)获得血管3d/2d弹性配准的目标函数；

(5)基于梯度优化的算法对目标函数进行求解。

本发明提出一种基于流形正则化的血管弹性形变模型，该形变模型利用了血管点的位移在其血管拓扑上具有平滑约束的先验，因此该模型适用于任意具有拓扑结构的点云结构；采用基于梯度优化的算法对目标函数进行求解，可以得到最优的血管3d/2d弹性配准结果。

优选地，所述步骤(2)中，使用在可再生核希尔伯特空间rkhs中运行的正则化技术，令为一个rkhs，通过tikhonov正则化来最小化目标函数，利用公式(1)求得定义在的空间变换

其中第一项为距离误差，第二项则限制了空间变换的复杂度，λ为权衡两项贡献的正则化系数，为rkhs的范数；

对于定义在rkhs的空间变换利用公式(2)表达为任意可再生核函数的线性组合，

其中k为正定的核函数为任意非空集合，k满足可再生性，对于任意与核函数k的內积等于其自身<f,k(ci,·)>k＝f(ci,·)，为核函数的输入参数；在此定义下rkhs的范数用其內积表示

优选地，所述步骤(2)中，采用高斯函数φ作为核函数，空间变换表示为公式(3)

其中为核函数φ输入参数，则为输入的控制点矩阵，c为控制点数量，核函数的形式为β为高斯半径，为空间变换系数；

则为系数矩阵，φ(·,k)为矩阵的第k列，φ的各元素φik＝φ(ci,yk)。

优选地，所述步骤(2)中，空间变换的范数用公式(4)计算，

其中矩阵的各元素为uij＝φ(ci,cj)，tr(·)为矩阵的迹，空间变换的范数项约束了由控制点c和系数矩阵w描述的空间形变场的平滑程度。

优选地，所述步骤(3)中，在血管拓扑流形上的约束用公式(5)表示

优选地，所述步骤(4)中，血管3d/2d弹性配准的目标函数为公式(6)

表示3d点的位移矩阵，w为弹性形变模型的参数，λ1和λ2是权重系数。

优选地，所述步骤(5)中，对目标函数第一项进行梯度计算，令tk表示则由于得到公式(7)

因此逐项计算得到其中为对于tk的雅克比矩阵。

优选地，所述步骤(5)中，目标函数的第二项q2＝λ1tr(w^tuw)的梯度通过公式(8)计算

优选地，所述步骤(5)中，对于目标函数第三项q3＝λ2tr(τ^tlτ)的梯度通过公式(9)计算，τ＝φ^tw-y，

本领域普通技术人员可以理解，实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，所述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，包括上述实施例方法的各步骤，而所述的存储介质可以是：rom/ram、磁碟、光盘、存储卡等。因此，与本发明的方法相对应的，本发明还同时包括一种基于流形正则化的血管3d/2d弹性配准装置，该装置通常以与方法各步骤相对应的功能模块的形式表示。该装置包括：

初始化模块，其在已知血管3d和2d点的匹配关系基础上，用表示3d点集矩阵，其中为3d点坐标，用表示2d点集矩阵，其中为2d点坐标，两个点集中对应点具有相同的索引值，k为点数量；

构建模块，其构建空间变换使得3d点的变换表示为τ(yk)为yk的位移；

流形正则化模块，其在血管拓扑流形上的约束描述为在拓扑上相邻的点的位移具有一致性，用流形正则化表示；

目标函数获得模块，其获得血管3d/2d弹性配准的目标函数；

求解模块，其基于梯度优化的算法对目标函数进行求解。

以下更详细地说明本发明的内容。

在已知血管3d和2d点的匹配关系基础上，用表示3d点集矩阵，其中为3d点坐标，用表示2d点集矩阵，其中为2d点坐标，两个点集中对应点具有相同的索引值，k为点数量。弹性配准的目的是估计的空间变换使得3d点的变换可表示为τ(yk)为yk的位移。该问题是不适定的(ill-posed)，因其具有无限数量的解。解决该问题的常用方法是使用在可再生核希尔伯特空间(reproducingkernelhilbertspace，rkhs)中运行的正则化技术。令为一个rkhs，通过tikhonov正则化来最小化目标函数可求得定义在的空间变换

其中第一项为距离误差，第二项则限制了空间变换的复杂度，λ为权衡两项贡献的正则化系数，为rkhs的范数。对于定义在rkhs的空间变换可以表达为任意可再生核函数的线性组合，

其中k为正定的核函数为任意非空集合。k满足可再生性，即对于任意与核函数k的內积等于其自身<f,k(ci,·)>k＝f(ci,·)，为核函数的输入参数。在此定义下rkhs的范数可用其內积表示在本方法中，采用高斯函数φ作为核函数。因此空间变换可以表示为下式

其中为核函数φ输入参数，则为输入的控制点矩阵，c为控制点数量，核函数的形式为β为高斯半径。为空间变换系数。

可写成公式(3)

则为系数矩阵，φ(·,k)为矩阵的第k列，φ的各元素φik＝φ(ci,yk)。

此外，点集y的空间变换可写成其矩阵形式为

根据rkhs的范数定义，空间变换的范数则可用下式计算，

其中矩阵的各元素为uij＝φ(ci,cj)，tr(·)为矩阵的迹。

空间变换的范数项约束了由控制点c和系数矩阵w描述的空间形变场的平滑程度。而对于3d/2d配准，本方法对血管点的位移在其血管拓扑上进行约束，令矩阵a∈{0,1}^k×k表示3d血管点拓扑的邻接矩阵，当yi和yj之间存在直接连接时aij＝1，否则aij＝0。

因此，提出一种在血管拓扑流形上的约束可以描述为在拓扑上相邻的点的位移具有一致性，用下式表示，

拓扑流形的约束可进一步限制空间变换在血管上的平滑性。通过观察对比，可发现拓扑流形约束同样可限制血管的长度变化。根据图拉普拉斯的定义，可用流形正则化表示，令d表示邻接矩阵a的度矩阵则血管拓扑的拉普拉斯矩阵为l＝d-a。令表示3d点的位移矩阵，因此血管3d/2d弹性配准的目标函数可如下所示，

w为弹性形变模型的参数，λ1和λ2是权重系数。

本方法采用基于梯度优化的算法对公式(6)进行求解，因此需要计算目标函数对参数的梯度

首先对目标函数第一项进行梯度计算，令tk表示则由于可得到如下关系

因此逐项计算得到其中为对于tk的雅克比矩阵。

目标函数的第二项q2＝λ1tr(w^tuw)的梯度可直接通过矩阵偏导公式计算

对于目标函数第三项q3＝λ2tr(τ^tlτ)，τ＝φ^tw-y，

综合上述计算结果可得目标函数对形变参数的梯度本方法采用l-bfgs优化算法对目标函数进行优化，可得到最优的血管3d/2d弹性配准结果。

以上所述，仅是本发明的较佳实施例，并非对本发明作任何形式上的限制，凡是依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属本发明技术方案的保护范围。