CN104282036B - 由稀疏样本对心脏进行的基于模型的重构 - Google Patents

Abstract

本发明提供了一种表示心脏的一部分的参数模型，所述参数模型使用从所述部分的其他实例的数据集的所述形状的统计先验来构建。通过使用标测电极，在受检者的心脏的所述部分的多个位置中获取电数据。将所述参数模型拟合到所述电数据和所述统计先验以产生心脏的所述部分的等值面以及其形状的重构。

Description

由稀疏样本对心脏进行的基于模型的重构

相关申请的交叉引用

本申请要求提交于2013年7月9目的美国临时申请No.61/844,024的权益，该临时申请以引用方式并入本文。

背景技术

1.技术领域

本发明涉及医学成像。更具体地讲，本发明涉及三维特征例如心脏心房的成像的改进。

2.背景技术

本文所用的某些首字母缩略词和缩写词的含义在表1中给出。

表1-首字母缩略词和缩写词

FAM	快速解剖标测
		PV	肺静脉
CT	计算机断层摄影术
		3D	三维

医疗导管插入是现今常规进行的。例如，就诸如心房纤颤之类的心律失常而言，在当心脏组织区向相邻组织异常地传导电信号时，便会发生心律失常，从而扰乱正常的心动周期并造成心律不同步。用于治疗心律失常的手术包括通过手术扰乱造成心律失常的信号源，以及扰乱用于这种信号的传导通路。通过经由导管施加能量(如射频能量)来选择性地消融心脏组织，有时可以终止或更改不利电信号从心脏一部分传播到另一部分。消融方法通过形成不传导的消融灶来破坏不需要的电通路。

左心房是复杂的三维特征，虽然所有左心房具有相同的基本形状，但左心房壁的尺寸因人而异。左心房可分成多个子结构，例如肺静脉、僧帽瓣或二尖瓣和隔膜，它们在概念上易于识别。所述子结构通常也因人而异，但对于整个左心房，每个子结构具有相同的基本形状。此外，给定的子结构具有与心脏其他子结构的相同关系，而不论子结构形状的个体差异。

发明内容

根据本发明的实施例提供了一种方法，所述方法通过如下步骤进行：定义表示心脏的一部分的形状的参数模型，以及由所述部分的其他实例的数据集来构建形状的统计先验。所述方法还通过如下步骤进行：将探针插入活体受检者体中，将所述探针的标测电极推进到与所述受检者心脏的所述部分的多个位置中的组织处于接触关系，获取来自相应位置的电数据，将所述参数模型拟合到所述电数据和统计先验以产生所述受检者的心脏的所述部分的等值面，以及重构所述受检者心脏的所述部分的所述形状，其中至少一个以上步骤在计算机硬件或包含在非暂态计算机可读存储介质中的计算机软件中执行。

根据所述方法的另一个方面，所述参数模型具有内坐标，并且定义参数模型包括将所述形状表示为在边界域内的点处定义的场函数，以及将所述点变换为所述内坐标。

在所述方法的另外一个方面，计算所述参数模型包括计算对于场函数的值和径向导数的边界条件。

在所述方法的又一个方面，计算所述参数模型包括通过加上新幂和新系数扩展拉普拉斯方程的解。

根据所述方法的一个方面，所述边界域包括单位球面。

根据所述方法的另外一个方面，变换所述点包括应用偏斜变换。

根据所述方法的一个方面，变换所述点包括应用球面投影变换。

根据所述方法的另外一个方面，变换所述点包括应用拉伸变换。

根据所述方法的又一个方面，变换的点对应于参数模型中的管和椭圆体，并且场函数包括管场公式和椭圆体场公式，其中拟合所述参数模型包括将所述管场公式和所述椭圆体场公式分别应用于所述管和椭圆体。

根据所述方法的又一个方面，拟合所述参数模型还包括将混合算子应用于所述管和椭圆体。

根据所述方法的另外一个方面，构建统计先验包括由心脏计算机断层摄影扫描来制备分段数据网格。

根据所述方法的另一个方面，拟合所述参数模型包括由所述数据网格来计算解剖特征。

根据所述方法的一个方面，所述解剖特征包括管中心线、管取向、管面积、管椭圆形程度和脊点中的至少一者。

根据所述方法的另外一个方面，拟合所述参数模型包括计算解剖特征中的不同解剖特征之间的相关系数。

根据所述方法的又一个方面，使所述电数据与所述拟合参数模型相关联包括使目标函数最小化，所述目标函数描述所述参数模型相对于所述电数据的估计误差。

根据所述方法的又一个方面，使目标函数最小化包括将来自所述统计先验的约束施加于所述目标函数。

根据所述方法的另外一个方面，所述目标函数包括代价函数。

根据所述方法的另一个方面，使目标函数最小化通过如下步骤进行：将相应权重赋给所述参数模型的参数，以及通过根据优化调度在所述目标函数相应迭代中改变相应权重来对所述目标函数进行迭代。

根据所述方法的又一个方面，使目标函数最小化包括计算所述目标函数的相对于参数模型的参数的导数。

根据所述方法的一个方面，拟合所述参数模型通过基于模型分量的加权来进行。

根据所述方法的又一个方面，拟合所述参数模是通过曲率加权来进行。

根据所述方法的另外一个方面，拟合所述参数模型通过基于骨架的拟合来进行。

附图说明

为更好地理解本发明，以举例的方式提供本发明的具体实施方式。要结合以下附图来阅读具体实施方式，附图中类似的元件用类似的附图标号来表示，并且其中：

图1是用于根据本发明的公开实施例施行医疗手术的系统的图示；

图2是可根据本发明实施例使用的快速解剖标测程序的图示；

图3是根据本发明实施例示出心脏中所测量的位置与形状模型的拟合的框图；

图4是根据本发明实施例示出偏斜变换对管状结构的效应的示意图；

图5是根据本发明实施例显示拉伸参数对管状结构的效应的示意图；

图6是根据本发明实施例示出膨胀参数的效应的三个管系列；

图7是根据本发明实施例的偏斜椭圆体的例子；

图8是根据本发明实施例示出混合算子的效应的两个心脏等值面系列；

图9是根据本发明实施例示出左心房网格构建中的某个阶段的表面；

图10根据本发明实施例示出了点的稀疏集合以及基准心房表面；

图11是根据本发明实施例示出拟合法结果的心房等值面；

图12是根据本发明实施例示出场函数方位的示意图；

图13是根据本发明实施例示出对于场值的边界条件的示意图；

图14是根据本发明实施例示出场径向梯度上的边界条件的示意图；

图15是根据本发明实施例示出斜坐标变换的示意图；

图16是根据本发明实施例示出心房网孔分析的示意图；

图17是根据本发明实施例示出兰勃特投影计算的定义及角度和向量的示意图；以及

图18是根据本发明实施例示出膨胀骨架的数据表示程序的示意图。

具体实施方式

为了能够全面理解本发明的各种原理，在以下说明中阐述了许多具体细节。然而对于本领域的技术人员将显而易见的是，并非所有这些细节始终都是实施本发明所必需的。在这种情况下，为了不使主要概念不必要地模糊，未详细示出熟知的电路、控制逻辑以及用于常规算法和进程的计算机程序指令细节。

系统概述

现在转到附图，首先参见图1，其为用于在活体受检者的心脏12上施行示例性导管插入手术的系统10的图示，该系统根据本发明所公开的实施例来构造和操作。该系统包括导管14，由操作者16将该导管14经由皮肤穿过患者的血管系统插入到心脏12的心室或血管结构中。操作者16(通常为医师)将导管的远侧末端18在消融靶点处与心壁接触。然后可根据美国专利No.6,226,542和6,301,496及共同转让的美国专利No.6,892,091中公开的方法，使用位于控制台24中的处理器22制备电激活图、解剖位置信息(即，导管远侧部分的信息)和其他功能图像，所述专利的公开内容以引用方式并入本文。体现为系统10的元件的一种商业产品可以商品名3系统购自Biosense Webster，Inc.(3333DiamondCanyon Road，Diamond Bar，CA91765)，其能够根据需要产生心脏的电解剖标测图。该系统可经本领域技术人员修改而体现本文所述的本发明原理，以用于使用如本文更详细描述的建模技术对结构例如左心房进行重构。

可以通过施加热能对例如通过电激活图评价确定为异常的区域进行消融，例如通过将来自射频(RF)发生器40的射频电流通过导管中的金属线传导至远侧末端18处的一个或多个电极，这些电极将射频能量施加至心肌。能量被吸收在组织中，从而将组织加热到一定点(通常约50℃)，在该温度下组织会永久性失去其电兴奋性。当手术成功后，在心脏组织中产生非传导性的消融灶，这些消融灶可中断导致心律失常的异常电通路。

导管14通常包括柄部20，在柄部上具有合适的控制器以使操作者16能够按消融手术所需对导管的远端进行操纵、定位和定向。为了辅助操作者16，导管14的远侧部分包含位置传感器(未示出)，其为处于控制台24中的定位处理器22提供信号。

可使消融能量和电信号经由缆线34穿过位于远侧末端18处或附近的导管末端和/或消融电极32，在心脏12和控制台24之间来回传送。也可以通过缆线34和消融电极32将起搏信号和其他控制器信号从控制台24传送到心脏12。同样连接至控制台24的感测电极33设置在消融电极32与缆线34之间。

接线35将控制台24与体表电极30和定位子系统的其他部件连接在一起。电极32和体表电极30可用于在消融位点处测量组织阻抗，如授予Govari等人的美国专利No.7,536,218中所教导的那样，该专利以引用方式并入本文。温度传感器(未示出)，通常为热电偶或热敏电阻器，可安装在电极32的每一个上或附近。

控制台24通常包括一个或多个消融功率发生器25。导管14可适于使用射频能量将消融能量传导至心脏。共同转让的美国专利No.6,814,733、6,997,924和7,156,816中公开了此类方法，这些专利以引用方式并入本文。

定位处理器22为系统10中定位子系统的元件，其测量导管14的位置和取向坐标。在一个实施例中，定位子系统包括磁定位跟踪构造，该磁定位跟踪构造利用生成磁场的线圈28，通过以预定的工作容积生成磁场并感测导管处的这些磁场来确定导管14的位置和取向。定位子系统可采用阻抗测量，如以引用方式并入本文的美国专利No.7,756,576以及上述美国专利No.7,536,218中所教导的那样。

如上所述，导管14联接到控制台24，该控制台使得操作者16能够观察并调控导管14的功能。处理器22通常为具有合适的信号处理电路的计算机。将处理器22联接以驱动监视器29。信号处理电路通常接收、放大、过滤并数字化来自导管14的信号，这些信号包括上述传感器和在导管14内朝远侧定位的多个位置感测电极(未示出)所产生的信号。数字化信号经由缆线38接收并被控制台24和定位系统使用，以计算导管14的位置和取向，并分析来自电极的电信号以及生成期望的电解剖标测图。

系统10可包括心电图(ECG)监视器，其被联接以接收来自一个或多个体表电极的信号。

左心房重构

以下描述涉及左心房。这仅作为例子而不是限制。本发明的原理适用于其他心脏腔室、血管结构，并且实际上适用于整个身体的中空器官。控制台24中的处理器可被适当地编程为执行下述功能。所述程序可接受由控制台24中的相同或其他处理器(未示出)处理的上述ECG信号和其他电解剖数据的输入，例如以重构门控图像。

本发明的实施例通过将受检者的测量数据与身体中空隔室(如心脏的左心房)的形状和/或子结构的相应存储模型相关联来重构基础基本形状。就左心房而言，通过使用模型，可以评价难以或不能用心脏导管显示的子结构。所述模型描述左心房及其子结构的不同三维形状。所述模型可由本领域已知的任何成像系统例如心内回波心动描记法(ICE)成像、电解剖成像(如CARTO成像)以及通过计算机断层摄影术(CT)成像、磁共振成像或手动超声扫描获取的图像来制备。形状模型42可以若干方式定义，如，基于网格、基于点、基于曲线图、基于隐场，并且它们可基于心室为管(即血管)的假设。在一种方法中，考虑血管中的流体动力系统的拉普拉斯方程或其修改形式可应用于描述所需形状。拉普拉斯方程的使用是示例性的。如以下部分所示，可使用描述和拟合所述形状的其他方法。可基于图像数据库和/或其模型的统计分析(例如PCA)对形状模型进一步约束或重新参数化。

现在参见图2，该图是可根据本发明实施例使用的快速解剖标测程序的图示。使用模型拟合程序，可通过拟合通常在左心房导管插入过程中获取的相对稀疏数据以达到最终的最佳形状，从而重构受检者左心房的详细图像。可通过心脏腔室44的超声或荧光镜成像，获取稀疏数据。或者，可通过如图2所示的快速解剖标测(FAM)，获取稀疏数据46。如本领域已知，数据46描述由导管上的位置传感器报告的位置。例如，可使用CARTO3系统的FAM功能与标测导管(例如导管)协同获取稀疏数据，所述两种器械均购自Biosense Webster，Inc.(3333Diamond Canyon Road，Diamond Bar，CA91765)。通常，使用如图2所示导管50上的位置传感器48所提供的位置信息，基于解剖标志或基准标记，将稀疏数据与三维空间中的坐标相关联。

在模型拟合程序过程中，将患者左心房壁的一定噪声、可能较少噪声和可能较高噪声的所测量的位置在拟合程序中与所存储的形状模型进行比较，并与子结构之间的已知关系进行比较。现在参见图3，该图是根据本发明实施例示出心脏中所测量的位置与形状模型的拟合的框图。如图3所示，在导管插入过程中，在方框52中获得患者左心房壁的所测量的位置。将所测量的位置与存储于数据库54中的形状模型42进行比较，并与子结构之间的已知关系进行比较，这适用于方框56中的拟合程序。通常，拟合程序应用对所存储的形状的更改，同时保持它们的关系，以生成适合所测量的位置的最佳形状58。所应用的更改的类型和程度可预先确定，并且通常基于所测量的可用形状以及基于物理特性，例如子结构的弹性。

在本发明的实施例中，在受到一组约束制约的情况下，可通过使目标函数最小化，将参数模型拟合到所述数据。一般而言，该方法如下所述：

定义表示心房形状的参数模型；

通过基准左心房形状的数据集的统计分析，对所述形状的结构、它们的相互关系和/或形状模型参数构建统计先验；以及

建立模型拟合程序，所述模型拟合程序在受到约束制约的情况下相对于模型参数优化对数据的形状拟合。

形状模型、统计先验和模型拟合程序的实现在以下部分中描述。

第一实施例。

在该实施例中，将心房形状表示为在边界域内的所有点处定义的场函数的等值面。通过应用例如以下所述那些的一系列坐标变换，将每个点变换成模型的内坐标系。然后在变换的坐标上计算每个解剖部分的贡献。通过混合解剖部分的贡献，计算最终场函数。该实施例的一种实现的坐标变换和场公式将在下文更详细描述。

坐标变换。

包围球变换。

通过应用变换T边界，将患者坐标系中给定的点t变换为由单位球面界定的域

x_偏斜＝T_边界(t)

在一个实施例中，使用仿射变换。将变换参数选择为使得所有所变换的关注点均在单位球面内，||x_偏斜||＜1.

偏斜变换。

对于每个解剖部分j，定义偏斜中心然后应用坐标变换，使得原点变换为在一个实施例中，将变换定义为：

其中T_j＝||x^j||。可颠倒该变换以在给定坐标向量x_偏斜情况下计算坐标向量x^j。现在参见图4，该图是根据本发明实施例显示偏斜变换的各种参数对管状结构的效应的示意图。

球面投影变换。

对于一些解剖部分，通过应用球面投影例如赤平投影，将“非偏斜”坐标x^j投影到平面坐标系中，使得：

其中T_投影为球面投影变换。

拉伸变换。

对于一些解剖部分，通过应用拉伸变换，在z方向(垂直于投影平面)上拉伸所投影的坐标。在一个实施例中，该变换通过由α_j参数化的幂变换定义，使得：

现在参见图5，该图是根据本发明实施例显示拉伸参数对管状结构的效应的与图4类似的示意图。与最上表面62相比，改变拉伸参数α_j对管弯曲率的效应在下表面60中较明显。如同管目标位置(偏斜中心)一样，管在其开口处的位置和取向也是恒定的。

幂变换可详细阐述为具有连续值和导数的类似分段幂变换，且每个管j和分段k具有单独α_jk参数。

解剖的部分场。

通过将场公式应用于所变换的点坐标，计算每个解剖部分在任何给定点处的场贡献。在一个实施例中，使用了两种类型的解剖部分：管和椭圆体。

管场公式。

管j由定义管中心的单位向量μ^j、定义主轴方向的正交单位向量定义轴长的标量膨胀函数η_j(x)和场函数f_管(·)参数化。将管在给定点处的场贡献定义为：

其中横截面椭圆体矩阵由下式给出：

在一个实施例中，可将膨胀函数η_j(x)定义为由β_i参数化的平面幂变换，使得：

其中r＝||x||。

现在参见图6，该图是根据本发明实施例示出膨胀参数β_j的效应的三管系列64、66、68。管64的右下部逐渐变细至一点，但在管66、68中逐渐变成球状。

膨胀函数也可详细阐述为连续平滑分段函数，且每个管j和分段k具有单独β_jk参数。

场函数f_管(·)可描述标准衰减函数，例如高斯函数和洛伦兹函数。的等值面将描述在μ^j处与单位球面相交的管，所述管具有中心线方向和大致椭圆的横截面。所述管以由α_j确定的速率朝向其端点弯曲，以由膨胀函数η_j确定的速率逐渐膨胀或收缩。

椭圆体场公式。

椭圆体j由定义其中心的单位向量μ^j、定义主轴方向的正交单位向量定义轴长的标量和场函数f_椭圆体(·)参数化。将椭圆体在给定点处的场贡献定义为：

其中矩阵 的等值面将描述在μ^j处居中且具有由给定的偏斜的偏斜椭圆体。现在参见图7，该图是根据本发明实施例的偏斜椭圆体70的例子。

混合函数。

通过应用混合算子，将各个解剖点的贡献组合。在一个实施例中，这通过贡献的逐点线性组合来完成，其中权重参数w_j：

现在参见图8，该图是根据本发明实施例示出混合算子的效应的两个心脏等值面系列。该图上部的心脏72的单独部分有区别。混合算子已应用于产生心脏74，其中失去了区别性，最明显的是失去了大血管的区别性。

统计先验模型。

可通过分析患者左心房形状的数据集的解剖特征，构建统计先验。可使用适当软件，通过处理CT扫描，构建表示左心房形状的网格的数据集。在一个实施例中，将形状模型拟合到数据集中的每个网格(使用诸如上述的模型拟合程序)，并将统计分析应用于所得的形状模型和/或它们的参数。或者，可直接由数据网格计算特征。

可通过自动化程序由网格或形状来计算解剖特征，例如肺静脉位置、取向和面积。为了比较整个数据集的特征，基于解剖标志将形状配准至共用坐标系。可通过拟合到适当的多变量概率分布，估计整个数据集的归一化特征(和/或模型参数)的联合分布。所得概率分布定义左心房的解剖特征的先验。

数据集构建。

在一个实施例中，使用分割软件例如CARTO系统，由CT扫描构建左心房网格。可通过手动地将其切割使得短残段保持连接，来将心房与肺静脉树、僧帽瓣和心耳分开。现在参见图9，该图是根据本发明实施例示出左心房网格构建中的某个阶段的表面。由分离所得的孔76、78、80、82可易于基于其在CT坐标中的位置用其解剖部分识别。脊点84、86、88由图标指示。然后可使用可免费使用的网格处理工具例如得自Sourceforge.net的MeshLab，对所得网格进行平滑处理、简化和修正拓扑。

特征抽取。

在一个实施例中，使用下文所述程序，利用取自网格表面的稠密点数据将形状模型拟合到每个心房网格。然后可使用下文所述公式由所得模型计算解剖特征。

管中心线。

在ξ^j坐标系中，管中心线只是直线。因此对于所需的高度h，管中心线点可通过如下计算：

其中是定义管j的球面投影T_投影的极点的单位向量。

然后可使用以上在名称为“坐标变换”的部分中定义的坐标变换，将中心线坐标变换为任何所需坐标系。可通过选择对应于用于样本的管切割位置的高度h，定义管中心特征。

管取向。

管取向在ξⁱ坐标系中以单位向量给出。这些向量可通过乘以所需坐标变换的雅可比矩阵来变换。取向可通过如下方式表示：所述单位向量绕其平均方向正交投影，产生管方向的两参数表示。

管面积

管横截面椭圆面积A^j在任何坐标系中以A^j＝πl₂l₃给出，其中l₂，l₃由矩阵的逆特征值给出，其中J为变换的雅可比式，并且#为由管权重和场阈值计算的归一化因数。

管椭圆形程度

管横截面椭圆可作为另外一种选择通过计算投影来描述，其中代表位于椭圆平面中的单位向量的预定义集合。例如，对于管j，可将指向指定邻近管j′的向量定义为：然后可将描述椭圆剩余自由度的3个单位向量的集合定义为： 其中R(轴线，角度)代表以给定角度围绕轴的旋转矩阵。

脊点。

对于两个邻近管j和j^t，可通过将连接管中心线的向量投影到管的椭圆矩阵上，而将近似中点线x_中点(h)定义为如下：

其中中心线点和椭圆矩阵使用标准点和双线性算子变换方法变换成x_偏斜共用坐标系。

现在参见图9，该图是根据本发明实施例示出左心房网格构建中的某个阶段的表面。然后可通过对一些高度值取样并检测场函数达到阈值f＝f_阈值的点，使近似中点线与心房表面相交。该交点将出现于脊点84、86、88，如图9所示。

心房容积。

可通过对域取样并对与f＞f_阈值的所有点相关联的面积求和，计算心房容积。

次要特征。

上述解剖测量值可用于计算次要特征，例如：

管中心之间的弦长，如左下腔PV和左上腔PV之间的左弦；

左弦和右弦之间的扭转角(方位角和余纬角)；

管位置向量，其将心房中心连接至管中心，归一化到单位长度，由绕平均向量方向的正交投影表示；

管位置向量与管取向向量之间的角度；以及

管横截面面积之和。

配准。

为了比较整个心房数据集的特征，定义共用坐标系。在一个实施例中，将坐标系的该原点定义为左脊点与右脊点之间的中点。可将旋转且均匀的尺度因子定义为使得脊点变换为坐标(±1，0，0)。可将附加旋转定义为使得左弦与xy平面平行，从而完成相似变换的定义。此处所述的变换仅作为例子，可使用其他变型和/或变换族提供配准步骤的替代实现。

概率分布估计。

可将配准变换应用于解剖特征以便将其归一化为共用坐标系。一些特征(例如脊点之间的距离)可在原始物理坐标系中计算，以提供心房绝对尺寸的统计的描述。所得特征基可拟合到概率分布，例如多变量正态分布。标准特征选择和/或降维方法(例如PCA)可用于增强统计先验的稳健性。正态分布的使用是示例性的，其他更高级的统计分布模型也可用于构建先验。

可通过使用联合分布模型例如多变量正态分布，来利用特征之间的相关性。我们对心房数据集的研究发现了各种解剖特征之间的显著相关性，指示了在仅部分信息可用的情况下它们的预测能力。这些相关性的例子在表2中给出(在多重比较修正之后，所有相关性都是统计意义上显著的)：

表2

特征1	特征2	相关系数(r)
			对左PV面积求和	对右PV面积求和	0.65
对所有PV面积求和	瓣膜面积	0.54
			从心耳到瓣膜-心房中心线段的	瓣膜面积	0.56

特征1	特征2	相关系数(r)
			距离
管位置单位向量(正交投影)	管取向单位向量(正交投影)	0.26-0.64
			心耳中心z坐标	从心耳到瓣膜-心房中心线段的距离	0.39
心房中心z坐标	从心耳到瓣膜-心房中心线段的距离	0.43
			左弦与右弦之间的扭转角	瓣膜方向(正交投影y坐标)	-0.36

本文所述代价函数中的统计先验项基于特征的联合概率，而不只是它们的边际分布。相关性表明，在给出特征B(如，左PV附近的指示其面积的点)相关信息的情况下，在优化过程中使用特征的联合概率应提高我们猜测特征A(如，右Pv面积)的能力。

模型拟合程序实施例。

通过使描述模型相对于从患者获取的稀疏数据的估计误差的目标函数最小化并结合适当约束，来估计心房形状模型参数。目标函数和约束函数由多个项组成，例如下文所述的那些。可在受到约束制约的情况下，通过标准非线性规划方法例如序列二次规划，使目标函数最小化。

在一个实施例中，稀疏数据由根据约定的协议从心房表面获取的点组成。现在参见图10，该图是根据本发明实施例的点的稀疏集合以及基准心房表面。点中的一些可描述心房特定区域的特征，例如线条或环(如，线条90所示)。线条90指示点集中的所需点，所述点集也可包括从心房表面的任何部分获取的一般未标记点。

现在参见图11，该图是根据本发明实施例示出拟合法结果的心房等值面。最终模型表面区域92叠加在基准表面94上。基准表面94可通过CT扫描建立。

近似距离项。

来自模型心房表面的数据点p的符号距离可通过计算下式近似：

其中f为所述点处的模型场值，为其空间梯度，f_阈值为阈值参数，并且为在大值处饱和的S型函数。近似距离可通过损失函数

E_距离＝L_距离(d_近似)

例如方距离或者鲁棒函数例如Huber损失函数，促成目标函数。近似距离也可用作约束，如通过要求数据点位于模型外的适当公差范围内，d_近似＜d_公差。可通过对给定特征中的点求平均值或通过对给定特征中的点计算柔性最大值，组合各个数据点的贡献。

环匹配项。

对于通常通过在管状结构例如肺静脉中操纵而获取的环状点集合，可计算匹配分数以将环与模型所述的对应管横截面比较。匹配分数可基于整个环或环上定义的特定点。在一个实施例中，使用例如以下的相似性度量，在拟合到数据点的椭圆和对应模型管横截面椭圆之间计算匹配分数：

下标0，1代表数据和模型椭圆。向量代表数据点椭圆平面中的椭圆中心，并且矩阵描述平面中的椭圆。L_环(·)为损失函数，例如上述那些。

隶属度项。

对于已知属于特定解剖部分的数据点，可使用模型解剖部分在所述点处对场的相对贡献，计算隶属度项。在一个实施例中，已知属于解剖部分j的点p的隶属度分数通过下式计算：

其中L_隶属度(·)为损失函数，例如上述那些。可例如通过求平均值或柔性最大值，组合多点贡献。隶属度项可并入目标函数或通过要求其超过适当阈值而用作约束。

本征模型约束。

为了确保优化过程的稳定性，可将多种约束应用于模型参数及其特征。可应用边界和线性约束以确保模型参数保持合理值。例如，可通过将线性约束应用于所有t，t^ε∈{1，2，3}，t≠t^ε，实施κ的最大椭圆体纵横比。附加约束可以非线性方式应用。例如，可通过要求每个管的偏斜中心从其他解剖部分的场接收足够的贡献，如实施模型形状的邻接。这些约束可在优化过程中严格实施，或可通过将约束馈入损失函数并将结果加到完整目标函数来应用为软约束。

统计先验项。

为了进一步确保稳定性，尤其是噪声或缺失数据的情况，可引入心房形状模型的统计先验。先验可应用于模型特征的集合，例如如以上在名称为“统计先验模型”部分中所述的管中心、横截面面积等等。先验也可在归一化到共用坐标系后直接应用于模型参数的子集。对心房大数据库计算所选特征。特征可由基准心房形状(如，从CT扫描获取)计算和/或由在根据多个数据点使用上述程序拟合到样本后的模型心房计算。模型参数可相似地取自数据集中每个样本的模型拟合结果。先验分布可采取简单参数分布例如多变量正态分布的形式，或更复杂的统计形式例如贝叶斯网络。先验分布的参数可使用标准程序例如最大似然法由样本的所计算特征估计。

然后所得先验分布可用于约束优化过程并改进由稀疏和/或噪声数据的模型估计。在每次优化器迭代时，由当前估计的心房模型计算所需特征。这些特征的值和/或当前模型参数自身可用于计算当前模型的先验概率。先验概率的适当函数(如，对数)可从目标函数减去，或用作优化约束，以将搜索空间限制于具有高先验统计似然的模型，从而即使在具有较少、噪声或无数据点的心房区域中也可获得合意的结果。

第二实施例。

在本发明另一个实施例中，使用标准非线性优化方法将非线性参数模型拟合到数据。一般而言，该方法如下所述：

(a)定义具有小参数空间的形状的参数模型。

(b)定义基于统计研究的约束和/或参数降维公式。

(c)使用代价函数和优化调度对数据(优化)的拟合合并有对数据的形状拟合和统计似然项。

形状42可以若干方式定义，例如基于网格、基于点、基于曲线图。

在模型生成程序过程中，将患者左心房壁的一定噪声、可能较少噪声和可能较高噪声的所测量的位置在拟合程序中与所存储的形状模进行比较，并与子结构之间的已知关系进行比较。所应用的更改的类型和程度可预先确定，并且通常基于所测量的可用形状以及基于物理特性，例如子结构的弹性。

上述形状、以及可用的较少所测量的位置的随后拟合基于假设心室为血管或管。在一种实现中，拉普拉斯方程的解可被概括以描述匹配所需血管的形状和取向的隐场函数。使用FAM技术的重构即使在低信噪比的条件下，并且在数据的量非常有限时也是可能的。如图3所示，在导管插入过程中，在方框52中获得患者左心房壁的所测量的位置。将所测量的位置与使用数据库54的形状模型42进行比较，并且与子结构之间的已知关系进行比较，这适用于方框56中的拟合程序。通常，拟合程序应用对所存储的形状的更改，同时保持该关系，以生成最佳形状58。所应用的更改的类型和程度可预先确定，并且通常基于所测量的可用形状以及基于物理特性，例如子结构的弹性。

在本发明的实施例中，来自稀疏数据的左心房重构的程序采用了以三维容积定义的场函数的隐式表面。边界条件为管位置、外形和方向。用于模型的数据为心房的患者专用网格，其产生自成像模式，例如CT扫描。然而，网格还为建立统计模型和降低参数空间维度的目的而结合。

左心房模型中的参数具有直观、自然的含义：

肺静脉(PV)、瓣膜和心耳管位置；

PV的轴线和其他管。

PV场方向性影响权重。

脊深度。

整体容积(阈值)。

心房主体“膨胀”的量。

心房中心的偏斜。

包围椭圆体。

这些参数在以下更正式地描述。

边界条件。

在单位球内的所有位置x处定义场函数

因此，其边界条件在单位球面的表面上的所有位置处定义

现在参见图12，其为根据本发明实施例示出场函数方位的示意图。此处由球面96表示的进入心房(PV、瓣膜和心耳)的每个管促成对场函数和其第一径向导数的边界条件。每个管被模型化为椭圆柱98，并且通过以下参数充分描述：

γ1是单位向量，其定义管中心线与单位球面的交叉处。

δ₂，δ₃是单位向量，其定义管的椭圆轴线的方向，δ₂⊥δ₃。

λ₂，λ₃为管的椭圆轴线的长度。

在任一点ξ处的管的场函数被模型化为单位高度高斯函数，其中协方差矩阵基于管的椭圆轴线而定义

其中∑^-1≡ΔΛ^-2Δ^T并且Δ≡(δ₂，δ₃)并且

要在单位球面上的点处模型化管场的影响，将该点标测到点γ1周围的相切平面上，这使用兰勃特等面积投影：

其中I_3x3为3×3单位矩阵。

定义在单位球面上的点处的管的影响可写为：

在边界场的径向梯度上的管的影响。

肺静脉可以不同角度与心房主体混合。这由管的中心线方向δ₁≡δ₂×δ₃表示。要完全表达心房模型中的这些方向，需定义场的径向梯度上的额外边界条件。管对边界场的n阶径向梯度的贡献被模型化为管场相对于(w.r.t.)所投影的坐标ξ，在法线γ₁的方向上的n阶梯度：

要计算这些导数，函数通过完全平方重新排列：

其中Z_uv＝u^T∑^-1v。变量的改变被应用，产生：

指数项的n阶导数由熟知的公式给出：

其中He_n(b)为盖然论者的埃尔米特多项式。

导数现在可使用以下公式来计算：

代入a＝0，完整表达式现在由以下公式给出：

应注意，如果管垂直于球面，则其对所有梯度的贡献为零，因为：

还应注意：

根据需要而定。

完全边界条件。

每个管通过其场函数和影响权重促成对于场值的边界条件。

对于场值的完全边界条件通过对所有管的贡献(由j索引)与添加的基线值f₀求和来定义：

其中N_j为管的数目。相似地，场的n阶径向梯度上的边界条件被定义为所有管的贡献的加权和：

现在参见图13和图14。图13是根据本发明实施例示出对于场值的边界条件100、102、104、106的示意图。图14是根据本发明实施例示出场径向梯度上的边界条件的示意图。箭头108、110、112、114、116示出了进入球面的管(未示出)的中心线方向。

基本场函数计算。

基本场函数的形式如下定义：

其中r为距原点的距离，为指示点的角位置的单位向量，并且为实值球谐函数(SPHARM)。指数l∈{0，...，N_j}覆盖基线值和管的贡献。指数l的范围从0到最大SPHARM等级(目前为20)，并且指数m的范围从-l到+l。场的径向从属关系根据l与幂模型化为幂定律。目前，l-从属关系模型化为：

c_jkl＝a_jkl

因数α_jk为控制心房脊的深度的模型的参数。该模型通过加上新幂α_jkl和新系数A_jktm扩展拉普拉斯方程的解，其中

通过向每个管单独地施加边界条件来计算系数A_jktm。基线值f₀通过额外的系数A₀₀₀₀＝4πY_00f0来施加。一般而言，在管l∈{0，...，N_j}的场贡献的n阶导数(n＝0、1......)上施加条件导出以下每个l、m的线性方程：

其中：

为管j对场函数的n阶导数的贡献上的边界条件的SPHARM扩展。

因数由以下递归关系给出：

单位球面上的所有积分可通过使用合适的均匀的球面网格(例如正二十面体)来离散化。

在目前的实施中，施加两个边界条件：一个对于场值施加，而一个对于其径向导数施加，n∈{0，1}。因此，仅使用指数k∈{1，2}。因此，系数的表达式由每个l、m的两个线性方程的方程组的解给出：

当c₀(l)＝c₁(l)时，对于导数的边界条件被忽略。对于c_k(l)＝α_k(l)的目前选择而言，如同其平均数(0阶分量)为零那样处理导数边界条件。

心房主体的膨胀。

在计算系数之后，将额外的膨胀操作应用到场函数：

使用参数β＜0导致原点周围的心房主体的各向同性膨胀，增加了其与球面的相似之处。当计算解系数A_jklm时，发现该操作产生更合意的结果，与使用不同的1从属关系(例如c_jk(l)＝a_jkl+β_jk)相反。

阈值分割和等值面生成。

从场函数中减去最终阈值，得出：

阈值可由其值定义，或者通过指定心房应占据单位球容积的百分比，并且计算单位球内的场函数值的合适的百分率来定义。改变阈值以使心房在所有位置更厚，与膨胀相反，其优选地使靠近原点的区域变厚。

最终心房表面被定义为模型场函数的零等值面：

场函数空间梯度和海赛函数。

模型的场函数w.r.t.x的空间梯度由以下公式给出：

其中

和

为前两个向量球谐函数。

斜坐标变换。

现在参见图15，其为根据本发明实施例示出斜坐标变换的示意图。变换的目标是将心房的中心偏斜到点x₀，同时保持包围球的表面上的所有点固定在其原始位置。容积中的所有其他点的位置应平滑地变换。针对单位球中的所有点x定义的标测表示如下：

所需标测被构建(例如)以满足两个条件：

(1)单位球面的表面上的所有点保持固定：

(2)存在标测v(x)的同质性(发现简单的等级1同质性产生最合意的结果)：

要求上述两个条件产生标测：

在球面坐标系中：

逆变换由以下公式给出：

配准坐标变换。

直至该点，心房模型定义在单位球面的边界内。应用线性坐标变换在所需坐标系中将单位球面变换为椭圆体。目前，仅允许没有反射的可逆的仿射变换，以产生：

T_配准(x)＝Mx+t^配准，

其中M为具有正行列式的3×3矩阵，并且treg为平移向量。如果需要，可应用更一般的变换，例如非线性透视变换。

对于以下所讨论的精细拟合优化而言，变换由固定参数M₀和表示初始变换的以及可优化参数M₁和来表示，如下所述：

其中固定参数被设定为：

M₀＝M

并且可优化参数被初始化(在运行精细拟合优化算法之前)为：

M₁＝I

导致参数自然居中并且定标在1阶。M₁的行列式相似地被约束为正的。

模型概要。

心房模型由场函数和坐标变换组成，并且概括如下：

场函数：

场：

系数取决于模型参数：

场函数参数：

全局参数：

β＝中心周围的膨胀。

f₀＝基线(约与球混合)。

f_阈值＝等值面生成的阈值。

每管参数：

α_jk＝脊的深度。

约束：

单位向量；

正交性：

完整性：

坐标变换：

完全变换：

T＝T_配准οT_偏斜

偏斜：

配准：

坐标变换参数：

参数：

M₁＝相对线性变换矩阵

约束：

偏斜边界：

无反射：

优化框架。

通过使用下述优化框架，可自动估计已知网格的模型参数，所述模型参数表示患者心房的CT扫描。以这种方式，可通过指定上述模型参数来充分描述患者心房。在一种方法中，可使用优化框架和其所估计的模型参数来分析配准患者心房的较大数据集。参数值的数据库可用于构建参数值的联合分布的统计模型。降低的维度能够应用关于心房形状的先验知识，以在不过度拟合所得数据的情况下实现挑战性、噪声和局部FAM数据的良好解释。

粗略拟合。

粗略拟合的目标是在给出表示患者心房的网格的情况下自动计算心房参数的初始估计。初始估计应在没有任何手动参数调整的情况下产生对给定心房的定性形状的可接受的拟合，并且充当随后的精细拟合法的初始条件。

粗略拟合算法的输入为基于对拓扑误差进行校正的计算机断层摄影(CT)扫描的心房网格，并且使用可用于MeshLab中的标准方法来进行平滑。

在这样的网格中，肺静脉(PV)、心耳和瓣膜被切割为短的管状残段。网格恰好包含针对每个PV、心耳和瓣膜的各一个孔。给出了孔的中心、边界顶点和面，以及每个孔的特性(左PV、右PV、心耳和瓣膜)。在原始CT扫描的物理坐标系中给出网格。

定义配准坐标变换的包围椭圆体：

包围椭圆体的中心treg被简单地定义为心房网格的质心。该质心目前计算为其面的质心的加权平均数。当对网格进行体素化并且计算其质量中心时获得相似的结果。

根据心房标志来定义标准椭圆体轴线方向第一轴线(“左至右”方向)被定义为从平均左PV中心到平均右PV中心的方向：

其中c_R[ε]PVs为右PV管残段孔中心或左PV管残段孔中心的平均数。第三轴线(“上”方向)被定义为从心房中心到所有PV的平均中心的方向，其与正交：

第二轴线方向通过计算向量积由右手定则给出：

包围椭圆体轴线长度优化。

考虑到固定方向，包围椭圆体的轴线应尽可能短，同时仍包封心房网格中的所有点。网格中的点t在轴线方向上的投影被定义为其中每个i∈{1，2，3}。椭圆体轴线长度表示为{ε_i}_{i∈{1，2，3)}。在所有点保持在椭圆体内的约束下，通过对平方轴线长度之和最小化来计算轴线长度。通过实施以下优化来实现最小化：

在该公式中，将优化应用到因此，说明所有点都落在椭圆体内的约束变为线性的。然后使用标准优化算法容易地解决问题。

配准变换矩阵的列现在由包围椭圆体的轴线给出：

M＝(ε₁，ε₂，ε₂)，

其中

完全包围椭圆体估计：

使得：

detM＞0

网孔分析。

要定义心房模型的边界条件，必须指定单位球面上的管位置，以及其主轴方向和长度。所有这些参数通过分析给定心房网格中的孔来估计。仅有邻接孔边界的面的薄“套管”用于该分析。在将网格的点反向变换到单位球面坐标系中之后，对网格进行所有分析：

x≡M^-1(t-t^配准)

现在参见图16，其为根据本发明实施例示出心房网孔分析的示意图。示出了五个管118、120、122、124、126。在管120、124、126中清楚地看出，实心三角形128代表孔的边界面。粗箭头130示出管的方向向量细箭头132表示管的外形椭圆的轴线(在孔内绘制)。点134代表孔的边界顶点以及它们在的方向上到单位球面的投影(由平行线136连接)。

管方向估计。

继续参见图16，该过程的目标是根据管的孔的边界面(三角形128)来估计管的方向。计算每个边界面的法线。管的方向向量定义为“尽可能正交于”所有边界面的法线。因此，边界面的法线上的投影应最小化：

其中矩阵N的列包含孔j的边界面的法线，矩阵w的对角包含它们的面积，并且ρ为拉格朗日乘子，其被引入以在上实施单位范数。

通过将选择为具有最小特征值NWN^T的特征向量来获得该问题的解。被归一化到单位范数并且使得指向向外的方向。

管外形椭圆轴线估计。

通过将原始孔中心在管方向向量的方向上投影到单位球面的表面上来计算管中心

选择因数c_j≥0以使得位于单位球面上：

通过选择来实现该条件：

粗略拟合概要。

总而言之，以下参数在粗略拟合法中自动估计：

配准变换(包围椭圆体)

T_配准(x)＝Mx+t^配准；

管中心

管椭圆的主轴和

管椭圆的轴线长度

仅仅根据网孔与其外边界，在不使用任何其他关于心房形状的信息的情况下计算参数。结果已示出，在没有任何手动参数调整，没有使用所有其他参数的固定值的情况下，这些参数足以预测各种患者的心房的定性形状。粗略拟合参数被用作随后的精细拟合分析的初始条件。

精细拟合框架。

在获得心房参数的粗略初始估计之后，使用一般优化框架来获得所有模型参数的精确估计。优化过程的数学基础由目标函数E的定义，以及相对于所有可优化的模型参数p的目标函数的导数的分析计算组成，该分析计算由以下公式给出：

精细拟合框架以模块化方式实施，以允许模型参数的任何子集的优化，同时使其他参数保持固定。这允许对优化过程的改进的控制，以及对各种参数的效应的理解。

要确保合适的会聚，精细拟合法应受到适当的约束。固有模型参数约束在上文中给出。表面距离误差度量(“主要候选方法”)不随场函数上的相乘常数的变化而变化，所以至少一个参数保持固定(例如瓣膜的场强影响权重如果需要，可采用额外的约束，例如限制坐标变换T_配准，以使得包围椭圆体保持在合理的尺寸。

已分析地研究了两种误差度量：

表面对表面误差度量(主要方法)：对误差度量最小化，该误差度量反映模型心房表面与已知心房表面之间的距离，以及其相对取向。为了效率，已知表面的距离变换可预先计算。该方法在下文中详细描述。

完全场拟合(替代方法)：使模型的完全场函数与基于数据计算的相似场函数之间的误差最小化。该方法的潜在优点是效率-大多数计算可在SPHARM系数A_klm上完成，而不是在物理空间中，这显著减小了计算复杂性。然而，实施该方法需要适当地修改模型场函数，或距离变换的合适的归一化，以确保两个“场”函数表现出相当的行为。该方法在下文中有所描述。

数据表示。

数据给定为表示患者左心房的CT扫描的网格，在于RSIP处建立的网格处理管线中进行处理。将这些网格的PV、心耳和瓣膜切割成短残段。目前实施两种变型：

“封顶”-管残段在网格处理管线中被自动封顶。

“延伸”：管在按名称为“网孔分析”的部分中所述而计算的方向上延伸直至它们超过衬垫式包围立方体(参见下文)，并且随后被重新封顶。

通常以物理单位计的数据的坐标系由坐标向量t≡T(x)表示，其中x为模型的单位球面坐标系中的相应的点，并且T为变换。

数据由其距离变换D_数据表示，在体素的离散栅格上以t为单位以数字方式计算。典 型体素尺寸为Δt＝0.5mm。距离变换在包含初始模型的包围椭圆体加上一些额外垫料的初 始立方体形状区域中明确地预先计算。在内部，该容积以体素为单位在与立方体对准的t′ 坐标中表示，使得：其中的列为初始 包围椭圆体的坐标轴单位向量，t^配准为其中心，t₀为旋转坐标中的衬垫式包围立方体的原 点，并且1＝(1，1，1)^T。应注意，和t^配准取自在粗略拟合过程中计算的初始包围椭圆体，并且 在优化过程中不改变。距离变换以t为单位通过乘以Δt来给出。

空间梯度和第二导数的海赛矩阵通过使用有限差分和恒等式在该区域内相似地预先计算：

在优化过程中，通过将线性插值与合适的坐标变换一起使用来在所需位置处取样这些函数。

在优化的过程中，包围椭圆体可改变。在一些情况下，该需要可出现以外推D_数据、和超过初始区域。要促进外推，夹持坐标被定义为：

t′_夹持≡最小[最大(t′，1)，t′_sz]，

其中t′_sz代表体素中的立方体的尺寸，并且最小、最大针对向量的每个分量分别计算。对所有点而言，使用以下公式来计算距离值：

D_数据(t′)＝D_数据(t′_夹持)+||t′-t′_夹持||₁

其中||·||₁代表曼哈顿度量。当与欧几里得度量相比时，曼哈顿度量产生更平滑的外推场，特别是靠近例如具有1个夹持坐标的区域与具有2个或多个夹持坐标的区域之间的过渡。

现在空间梯度由以下公式给出：

应注意，由于t′_夹持不改变的是这些方向，所述的夹持分量为0。sign函数应用到向量t′-t′_夹持的每个分量，并且对于所有非夹持分量为0。

海赛函数仅由以下公式给出：

仅有非夹持分量促成海赛函数。

数据权重在原始心房表面上被定义为1，而在管的延伸部分和顶端上被定义为0。加权函数通过使用在心房表面上的最近邻而被定义在整个3D容积中。使用具有等于最大取样距离最大(ε_i)Δx的半宽度的高斯滤波器来使该函数平滑，其中ε_i为初始包围椭圆体轴线长度，并且Δx表征模型的x坐标系中的栅格点之间的距离。通过使用针对距离函数的上述方法来计算加权函数的空间梯度。

表面对表面误差度量。

误差度量基于模型与实际已知的心房表面之间的平均距离，以及它们的相对取向。超出目前模型的包围椭圆体的已知心房表面的部分通过向其分配距离值而受到惩罚。

误差度量被定义为向前平均距离量度与向后平均距离量度之和，计算为表面上的加权平均数。所有距离和表面积在数据坐标系中测量，t≡T(x)。从模型到数据的距离基于数据中给出的已知心房表面的预先计算的距离变换。从数据到模型的距离通过使模型场函数的梯度相对于t而归一化来近似。靠近表面，给出了近似的距离量度，如将在下文中示出。

完全误差度量由以下公式给出：

E≡W_{模型→数据}+E_{数据→模型}

模型对数据误差度量由以下公式给出：

其中为从模型表面到数据表面的在模型表面上求积分的累积加权绝对距离，为累积加权取向量度，为以距离为单位的恒定元参数，并且S_模型为在模型表面上求积分的总权重。

数据对模型误差度量由以下公式给出：

其中为从数据到模型的累积近似绝对距离，其在落入模型的包围椭圆体内的数据部分上求积分，并且为椭圆体内的取向量度。对于数据表面的外部而言，代表所受到的累积惩罚，以距离为单位。数据的总加权表面积为预先计算的常数，其由表示。

累积外部惩罚由以下公式给出：

其中代表已知心房表面，φ_外部(t)定义由已知心房表面上的单位面积引起的惩罚，该心房表面超出边界ta的加权函数。应注意，φ_外部(t)具有t坐标系中的距离单位。模型的包围椭圆体由变换的单位球给出：

并且外部区域由其补集给出。

目前，使用恒定的惩罚使得

其中为数据的总表面积，该数据落入模型的包围椭圆体内。

表面泛函和

量和通过使用通用框架来计算，该通用框架由三个泛函和定义：

表面由一些平滑场函数f的零值等值面表示。心房模型的场函数f_模型由以上名称为“场函数”的部分中所述来定义。对于已知心房数据而言，符号距离变换D_数据目前被用作数据场函数，f_数据＝D_数据。

估计的从一些源表面到一些目的表面的在边界域内的源表面上求积分的绝对距离由泛函给出，该泛函定义为：

其中源表面和目的表面通过其场函数f_src和f_dst来定义，并且权重由源加权函数W_src给出。t坐标系中的目的表面的绝对距离绝对距离_t[f_dst]将在稍后定义。该函数在t坐标系中的源表面上求积分，

并且dS_t为该坐标系中的面积元素。在当前设置中，边界域为模型的包围椭圆体，其由变换的单位球给出：

域内的源表面的加权面积由泛函给出，该泛函定义为：

目前，所有测量在模型的x坐标系中的固定球面栅格上取样(这不影响测量单位，仅影响数字取样方案)。源等值面上的面积元素通过按以下所述进行变量的改变来变换：

其中为变换T的雅可比矩阵，并且代表其行列式的绝对值。算子和分别代表空间梯度w.r.t.，t和x坐标系。

因此，表面泛函由以下公式给出：

其中为x坐标系中的边界域。

原则上，和S_模型可使用该方法通过在每次迭代中抽取模型等值面来计算。然而，等值面生成是计算密集的。此外，由于积分区域取决于目前模型参数，因此所述计算这些项w.r.t.模型参数的导数不是简单的。相似地，以这种方式计算将需要在每次迭代中重新计算在数据心房表面上的所有点处的模型场函数。然而，模型场值在球面坐标中(即在容积中)的固定栅格上最方便地计算。

因此，表面泛函通过应用余面积公式而重新用公式表示为容积积分。该公式的推论针对任何表现良好的规定，

其中δ(·)代表狄拉克δ函数。

为了能够在当前设置中应用该公式，x归一化场函数被定义为：

该函数具有x的单位，即模型的坐标系中的距离。代入并且假设靠近表面(即该表面为表现良好的二维流形)，该假设始终适用于数据的符号距离变换。对于心房模型表面而言，该假设适用于除了原点以外的任何地方，其中然而，会聚是缓慢的，并且φ在原点附近相当平坦，所以这应该不会引起结果的过分偏离。由此得出，以及源表面上的无因次量产生：

分别代入绝对距离(f_dst)和表 面泛函现在由以下公式给出：

在实践中，容积积分离散在球面栅格上，并且狄拉克δ函数通过使用以下公式在x坐标系中近似：

其中∈_x＝1.5Δx为半宽度，并且Δx表征栅格点之间的典型距离。使用该δ函数近似意味着，表面的“厚度”在模型的x坐标系中是恒定的，根据需要在大约3个点上求平均值。

要估计绝对距离_t[f_dst]，目的场函数应转换为以t为单位的近似距离函数。这通过使用t归一化场函数来实现，其被定义为：

对于数据而言，场函数目前为符号距离变换，所以几乎无处不在，产生：

其中D_数据是以t坐标系为单位的距数据心房表面的距离，其从预先计算的距离变换中取样。

对于模型场函数而言，引入饱和度以避免振动远离心房模型表面，如下所述：

恒定的确定以t为单位的饱和度值，并且通常被设定成一些较大距离，例如计算的数据距离变换的包围立方体的对角。靠近零等值面，空间梯度的量值产生具有所需单位的近似距离函数。

要靠近表面实施绝对距离_t[f_dst]的平滑度，需使用符号函数的抹平近似：

其中具有在x坐标系中的均匀厚度的符号函数近似被定义为sign(x)＝2H(x)-1，并且海维赛德阶跃函数H(x)由以下公式表示：

取向泛函由偶极-偶极相互作用的能量激发。该项的目的是确保两个表面的相应贴片具有相似的取向。该泛函由以下公式给出：

应注意，被积函数是无因次的。

路径泛函

路径泛函沿着由t_路径(s_t)参数化式地给出的路径(例如心房骨架)评估任选地归一化的标量场函数的路径积分，其中参数s_t代表t坐标系中的弧长。路径积分由以下公式给出：

其中为t坐标系中的总路径长度。

定义x_路径＝T^-1(t_路径)，路径积分可由x坐标系中的弧长表示，这通过代入：

产生：

其中为x坐标系中的路径长度。

路径积分可转换为容积表示，如下所述：

其中δ_x(·)为x坐标系中的三维δ函数，并且表示路径方向向量函数通过最近邻延伸至全容积。

δ函数δ_x(·)可在任何正交坐标系中由3个一维δ函数的乘积表示。此处，我们假装我们在使用可分离的δ函数表示，例如S.Zahedi、A.K.Tornberg所建议的高斯近似法，“Delta function approximations in level set methods by distance functionextension”，Journal of Computational Physics229(2010)2199-2219。具体地讲，可使用以下表示：

δ_x(x-x_路径)＝δ(0)δ(Δx_||)δ(Δx_⊥)

其中

为x-x_路径的分量，其平行于该路径，并且

为在位置s_x处垂直于该路径的分量。按照定义，由于局部坐标轴的选择，第三分量为零。

路径泛函现在由以下公式给出：

对于任何给定的x而言，内部路径积分沿着实现Δx_||＝0的路径筛选出点。对于该点而言，仅为x与路径之间的最近距离。因此，路径泛函由以下公式给出：

在实践中，最近邻方向向量和距离是在t坐标系中预先计算的，其中路径为固定的，并且在点t＝T(x)处取样。假设坐标变换T对路径附近的最近邻的特性具有可忽略的效应。距离可使用以下公式来近似：

或者通过以下公式来近似：

误差度量导数。

对于高效和精确的优化而言，误差度量w.r.t.的导数、模型参数应分析地计算。应注意，模型参数影响坐标变换T和场函数f_模型。分析导数计算可使用矩阵计算手动地进行，或者为了该目的，通过实施标准自动微分算法来进行。

数据对模型误差度量(替代公式)。

所有数据对模型误差度量都可作为另外一种选择表示为已知心房表面上的表面积分。其导数的表达式可借助于散度定理进行简化。

精细拟合框架概要。

误差度量。

表面泛函。

归一化场函数。

模型场函数。

模型场函数空间梯度。

模型系数

数据场函数。

数据场函数空间梯度。

兰勃特等面积投影。

混合到单位球面中的管的场所需的边界条件应成形为类似定义在球面表面上的高斯函数。这种函数的规范定义被称为Kent分布。通过使用兰勃特投影，管场函数f_管近似该分布的形状(Kent，J.T.1982：“The Fisher-Bingham distribution on the sphere”，J.Royal.Stat.Soc.，44：71-80。使用投影而不是原始Kent分布能够在无需明确计算在文件中定义的特征向量γ₂，γ₃的情况下计算场。这产生了对于在优化阶段中计算的目标函数的导数的更简单的表达式。

现在参见图17，其为根据本发明实施例示出兰勃特投影计算的定义及角度和向量的示意图。兰勃特方位角等面积投影是从球面到盘的标测，其准确地表示球面的所有区域中的面积。该盘可被认为定位在以球面表面上所选的极点为中心的相切平面上。在这种情况下，点的投影ζ由以下公式给出(Kent，J.T.，1982年)：

如图17所示，其中θ为极角，为方位角，为定义坐标系的正交单位向量。所需极点位于PV中心线与球面表面的交叉处，角θ、的所需函数可表达如下：

其中：

代入这些表达式并且使用轴线的正交性，获得以下表达式：

其中I_3x3为单位矩阵。应注意，该表达式与轴线z□2、z□3的选择无关。根据需要，靠近γ₁，变换的点近似等于原点：

如通过采用如下范数平方的极限可发现，范数||ξ||在极点处达到其最小值0，并且在对极处达到其最大值2：

替代的精细拟合方法-完全场拟合。

该方法的优点是，优化可在系数A_ktm的空间中而不在物理空间中进行。然而，将需要特殊的修改以确保完全模型场函数和数据场函数具有相似的行为。例如，目前模型场函数在原点处发散至无穷大。

目标函数被定义为模型场函数与表示已知心房数据的场函数之间的均方误差(MSE)。目前，这根据患者的心房网格的距离变换来计算，并且如在使用名称为“模型概要”的部分中所述的模型中定义的坐标变换T来取样。

通过使用球谐函数的正交性并且评估积分，目标函数可用公式表示为SPHARM指数l，m上的总和，更确切地为完全3D容积上的总和，使得能够进行高效的计算。然后目标函数由以下公式给出：

在RHS中的所有指数上进行求和，并且该项定义如下：

A_ktm为模型中定义的系数。

归一化因数：Z_kkflm≡[d_k(l)+d_kf(l)+3]^-1。

指数：d_k(l)≡c_k(l)+β。

数据积分：

数据场SPHARM扩展：

数据幂：

通过省去指数l，m并且将指数k，k’合并到矩阵向量乘积符号中，目标函数可写为：

E＝A^TZA-2A^TI+Φ_数据

要进行高效准确地优化，需要使用相对于模型参数的目标函数导数的分析计算。由于完全以分析项指定心房模型，所以还可分析地计算相对于所有模型参数的导数。

实施细节。

优化改进。

以下部分使用改进的初始线性变换估计、数据加权积分和基于骨架-的拟合来概述用于在优化中的逐步改进的程序。

改进的初始线性变换。

以下程序运行：

优化所有变换参数(变换矩阵、中心)；

找出最小-容积椭圆体，其包括所有数据点：

使得

det M＞0

心房模型加权的两个方案可被采用：

基于模型分量的加权

通过模型的管分量分解模型函数：

恒定分量(球)通过因数添加到每个管分量。

计算每个分量的函数：

探试法：表示距该分量的距离。脊为其中的(变换的)乘积为高的地方。

变换(Student-t函数：

所选择的分量乘积：

w_模型＝w(φ_右上，μ，σ_右)*w(φ_右下，μ，σ_右)+α(w(φ_左上，μ，σ_左)*w(φ_左下，μ，σ_左)+w(φ_左上，μ，σ_{left_s app})*w(φ_心耳，μ，σ_{left_s app})+w(φ_左下，μ，σ_{left_i app})*w(φ_心耳，μ，σ_{left_i app})+w(φ_左上，μ，σ_{left_s_i app})w(φ_左下，μ，σ_{left_s_i app})*w(φ_心耳，μ，σ_{left_s_i app}))

该加权方案的目标是向脊区域分配增加的权重。然而，其不像曲率加权(如下所述)一样准确。

曲率加权

隐式表面在3D中的曲率给定为：

高斯曲率

平均曲率

该方法取决于海赛矩阵，并且可以：数字地计算(在球面坐标中)或分析地计算。

基于骨架-的拟合。

原理：-起始自模型的薄的型式，并且将其拟合到数据骨架。-逐渐使数据表示和模型膨胀。始终保持模型表面在数据表面内，以防止局部最小。

数据骨架计算。

将心房中心定义为数据表面的距离变换D的最大值处。

使用迪杰斯特拉算法计算从管孔中心到心房中心的最佳路径：

能量函数：1/D。

曲线图定义：连接所有相邻体素，边缘权重＝{体素之间的距离}*(1/Dsrc+1/Ddst)/2

基于骨架-的优化框架。

使用薄的模型，其具有固定阈值，没有“球”。从数据骨架到模型表面的符号距离可如下计算：

模型表面应平行于骨架，并且可能平行于数据表面。这可使用取向泛函来进行量化：

模型表面不应为外侧数据表面：修改现存的距离量度以惩罚外侧部分。通过使符号距离泛函最大化，骨架在模型内应尽可能远。在优化过程中，可通过在代价函数或非线性约束函数中合并取向泛函来约束骨架取向。

数据表示-膨胀骨架。

现在参见图18，其为根据本发明实施例示出膨胀骨架的数据表示程序的示意图。

计算骨架臂138、140、142、144(迪杰斯特拉，1/D代价、由路径段长加权)

起始自孔的椭圆，由尺度因子进行定标。

通过变换以前的椭圆轴线，沿着骨架臂的每个点被分配了椭圆。

变换由与骨架臂相切的变化来定义。参见向量146的设定：

轴线≡t_i-1×t_i

角度＝cos^-1(t_i-1·t_i)

通过骨架上的最近邻的椭圆来标记整个容积。

中心斑点(Central blob)：将阈值定义为横跨所有中心点椭圆(每个管一个)的距外心房表面的最小距离。当尺度因子较大时，中心斑点为完整心房。

优化阶段。

初始较弱猜测-基于骨架顶端的PV位置和方向。使模型内的骨架深度最大化(仅改变管位置和方向)：

将模型拟合到具有薄瓣膜的薄的膨胀骨架。使用“完整尺寸的”瓣膜，(首先仅为全局参数和瓣膜参数，然后为所有参数)。

逐渐膨胀。

优化约束。

模型必须保持在下一阶段的数据表示内。相反地，下一阶段数据必须保持在模型外。如上所述，骨架必须保持在模型内。

使用“柔性最小值/最大值”约束泛函。

要帮助所有模型面积保持在边界内，需将“切断”指数项添加到代价函数：

在替代实施例中，可实施以下选项：

使用通式的高阶导数边界条件：

Z_uv≡U^T∑^-1V

这允许略微尖锐的管角。

对于具有短的共同口的情况，可加入额外的弱“管”。

管约束。

对于每个PV和心耳：

加权靠近管开口的面积(数据的“长切割”)。

距离的柔性最大值必须在公差范围内。

取向匹配的柔性最大值必须在公差范围内。

每管2×2约束：(数据、模型)×(距离、取向)，用于所有拟合阶段。

未约束的瓣膜(当作心房主体)。

本领域的技术人员会认识到，本发明并不限于在上文中具体示出和描述的内容。更确切地说，本发明的范围包括上文所述各种特征的组合与子组合，以及这些特征的不在现有技术范围内的变型和修改形式，这些变型和修改形式是本领域技术人员在阅读上述说明后可想到的。

Claims (22)

1.一种用于重构心脏的形状的方法，包括以下步骤：

定义表示所述心脏的一部分的所述形状的参数模型；

由所述部分的其他实例的数据集来构建所述形状的统计先验；

将所述参数模型拟合到所述统计先验；

将探针插入活体受检者体中，所述探针具有标测电极；

将所述标测电极推进到与所述受检者的所述心脏的所述部分中的多个位置中的组织处于接触关系；

获取来自相应位置的电数据；

使所述电数据与所拟合的参数模型相关联以产生所述受检者的所述心脏的所述部分的等值面；以及

响应于使所述电数据关联来重构所述受检者的所述心脏的所述部分的所述形状，其中上述步骤中的至少一个步骤在计算机硬件或体现在非暂态计算机可读存储介质中的计算机软件中执行。

2. 根据权利要求1所述的方法，其中所述参数模型具有内坐标，并且其中定义参数模型包括：

将所述形状表示为在边界域内的点处定义的场函数；以及

将所述点变换为所述内坐标以定义变换的点。

3.根据权利要求2所述的方法，其中所述场函数具有值和径向导数，还包括通过计算对于所述值和所述径向导数的边界条件来计算所述参数模型。

4.根据权利要求2所述的方法，其中所述参数模型包括具有幂和系数的球谐函数展开，还包括通过添加新幂和新系数来扩展拉普拉斯方程的解。

5.根据权利要求2所述的方法，其中所述边界域包括单位球面。

6.根据权利要求2所述的方法，其中变换所述点包括应用偏斜变换。

7.根据权利要求2所述的方法，其中变换所述点包括应用球面投影变换。

8.根据权利要求2所述的方法，其中变换所述点包括应用拉伸变换。

9.根据权利要求2所述的方法，其中所变换的点对应于所述参数模型中的管和椭圆体，并且所述场函数包括管场公式和椭圆体场公式，其中拟合所述参数模型包括将所述管场公式和所述椭圆体场公式分别应用于所述管和椭圆体。

10.根据权利要求9所述的方法，其中拟合所述参数模型还包括将混合算子应用于所述管和椭圆体。

11.根据权利要求1所述的方法，其中构建统计先验包括由心脏计算机断层摄影扫描来制备分段数据网格。

12.根据权利要求11所述的方法，其中拟合所述参数模型包括由所述数据网格来计算解剖特征。

13.根据权利要求12所述的方法，其中所述解剖特征包括管中心线、管取向、管面积、管椭圆形程度和脊点中的至少一者。

14.根据权利要求12所述的方法，其中拟合所述参数模型包括计算所述解剖特征中的不同解剖特征之间的相关系数。

15.根据权利要求1所述的方法，其中使所述电数据与所拟合的参数模型相关联包括使目标函数最小化，所述目标函数描述所述参数模型相对于所述电数据的估计误差。

16.根据权利要求15所述的方法，其中使目标函数最小化包括将来自所述统计先验的约束施加于所述目标函数。

17.根据权利要求15所述的方法，其中所述目标函数包括代价函数。

18.根据权利要求15所述的方法，其中使目标函数最小化通过以下步骤进行：将相应权重赋给所述参数模型的参数；以及通过根据优化调度在所述目标函数的相应迭代中改变所述相应权重来对所述目标函数进行迭代。

19.根据权利要求15所述的方法，其中使目标函数最小化包括计算所述目标函数的相对于所述参数模型的参数的导数。

20.根据权利要求1所述的方法，其中拟合所述参数模型包括基于模型分量的加权。

21.根据权利要求1所述的方法，其中拟合所述参数模型包括曲率加权。

22.根据权利要求1所述的方法，其中拟合所述参数模型包括基于骨架的拟合。