CN101650835B - 基于局部线性嵌入法构建狗左心室传导系统三维几何结构 - Google Patents

Abstract

一种基于局部线性嵌入法构建狗左心室传导系统的三维几何结构，具体涉及一种将平面图像映射成曲面图像的方法。本发明由下述步骤完成心室传导系统的构建：首先，将狗的左心室三维几何构形运用LLE算法映射到平面，记为心室传导系统的二维线状图；然后，将真实的心室传导系统的二维平面图中的传导系统提取出来并与其整合到新的图中；最后，将整合后的新图，再一次运用LLE算法映射到左心室的三维几何构型中，从而得到左心室传导系统的三维几何结构。本发明具有保持空间结构的局部线性特征的优势，同时兼顾高效快速的特点，并利用真实的数据，反应了较为真实的狗的左心室传导系统。

Description

基于局部线性嵌入法构建狗左心室传导系统三维几何结构

技术领域

本发明涉及一种将平面图像映射成曲面图像的方法，具体涉及一种基于局部线性嵌入法构建狗左心室传导系统三维几何结构。

背景技术

心室传导系统和普通的心肌细胞都有传导性。在心室中，正常的兴奋传导主要依靠心室传导系统。该传导系统包括：左右束支和浦肯野纤维，其中在左心室的部分包括左束支和浦肯野纤维。它位于心室的内膜层，肉眼不可见。由于浦肯野纤维较细的部分直径已接近5μm，因此用断层扫描术或核磁共振的办法很难得到能够反映其几何结构的数据。目前国际上的常用方法包括以大片心内膜代替，或以分形的方法生长等，这些方法都是对浦肯野纤维的近似描述，仅在一定范围内可近似反映心室的电生理状况，但由于没能如实刻画传导系统的几何空间结构，而将其应用于心室上的生理、病理研究，对反映心室活动的真实性来说是很有限的。

作为流形学习方法中的一种，局部线性嵌入(Locally Linear Embedding，即LLE)算法于2000年在Science杂志上发表后，一直倍受关注。LLE算法是一种通过局部线性关系的联合来揭示全局非线性结构的非线性降维方法。其核心思想是在降维映射前后保持源数据的局部近邻性质，即在嵌入空间每个采样点可以用它的近邻点线性表示，在低维空间中保持每个邻域中的权值不变，重构原数据点，使重构误差最小。算法中每个点的近邻权值在平移，旋转，伸缩变换下保持不变。这一特点使LLE算法在高维空间与低维空间的转化过程中，保持点的局部位置关系。就该发明而言，当这个转化过程发生在曲面与平面之间时，也同样适用。

LLE算法是基于几何直觉的，即把高维空间数据点 

按维数映射到低维嵌入空间 即 

如图1，其具体步骤为：

1.寻找各数据点 

的邻近点：

设原始数据由N维实值向量组成，记做 

(i＝0，1，…，N)，由于数据由真正光滑的多面体取样而来，故每个数据点和它的邻近点位于或近似位于该多面体的局部线性平面上。这样 就能通过线性组合系数刻画出局部平面的几何特征。在LLE中，通过度量欧氏距离的方法可找到每个数据点的K个最近的邻近数据点。每个数据点的重构错误用成本函数来衡量：

2.计算权值W_ij，通过W_ij与邻近点构造数据点 

并由约束条件最小化(1)式成本：

权值W_ij说明第j个数据点对重构第i个数据点的所做的贡献。为了得到合适的权值，须对成本函数进行最小值计算。

3.通过权值矩阵W_ij构造向量 由非零特征值最小化方程(2)所示二次型：

每个数据点的重构权值是通过他的邻近数据点计算的(和其他数据点的权值毫无关系)，但是低维坐标系统却要通过计算一个N×N矩阵而得到(即要考虑被权值矩阵所确定的图中所有存在连接的结点对)。不同维数的嵌入空间的连续计算，仅需要通过式(2)计算出最小的几个特征值对应的几个特征向量即可(即较小的特征值对应的特征向量)。

然而，现有的心室传导系统模型与真实几何体的存在着较大偏差，想要更进一步对心室兴奋传导等电生理活动研究，仍需要很大的努力。

发明内容

为解决上述技术问题，本发明旨在提供一种基于局部线性嵌入法构建狗左心室传导系统三维几何结构，以真实的数据为基础，将流形学习的方法应用于构建心室传导系统的三维曲面上，以获得由真实数据生成的心室传导系统的几何构型，使心室传导系统的三维信息得以直观地再现，不再仅由简单的分形生长等方法粗略描述，从而更进一步对心室兴奋传导和电生理活动的研究奠定良好的空间结构基础。

本发明解决上述技术问题的技术方案是：基于局部线性嵌入法构建狗左心室传导系统的三维几何结构：首先将狗的左心室传导系统(心内膜层的左束支和浦肯野纤维网)展成平面，固定拍摄成图片。虽然经过特殊染色处理，该传导系统与心内膜背景颜色仍不易区分，这时需要将拍好的图片通过半自动过程提取图片中的心室传导系统的线状图，该过程主要利用阈值分割的办法将图像进行分割提取子图，并将提取出来的子图拼接到一起，从而得到狗的左心室传导系统的二维平面线状图数据点的位置信息。

获取狗的左心室心肌三维结构数据，生成左心室内膜层三角网格曲面；在进行曲面到平 面的映射时，以左心室内膜层的离散点数据作为点数据，运用局部线性嵌入算法(LLE算法)将其映射到二维平面，得到心室内膜层的二维点数据。

通过调整比例尺、边界提取及特征点匹配，将内膜层的二维点数据与传导系统平面线状图整合到一起，即，将二者以相同比例大小，按照边界范围，尽可能多地匹配二者的特征点，最终将二者置于同一坐标系下，形成新图；由于所提取的线状展开图存在一定程度拉伸后的形变，因此在整合的过程中，必要时须局部进行放缩，使之可以更好地匹配；最终再将新图运用LLE算法，映射到狗心室的三维曲面结构中。

由于LLE算法得到的映射是一一映射，因此，将心室传导系统的平面线状图以心室内膜层的二维映射点线性表示；由于该二维映射点与心室的三维曲面是一一映射，因此心室传导系统也可由心室的三维曲面数据点线性表示；至此，心室传导系统得以完整地展现在心室三维结构中。

为实现以上内容，本发明的具体步骤如下：

一、进行离体的动物实验，将狗的左心室作特殊的显影技术处理，使其可以清晰地看到内膜层附近心室传导系统的线状结构，将切开的心室展成平面，固定并拍摄图片。

二、提取步骤1得到的图片中的心室传导系统，即左束支及浦肯野纤维，将提取的离散点集记为点集C，提取的过程是半自动过程，即：

1)先以灰度阈值来提取传导系统的轮廓，即颜色较深的部分；

2)再进行区域分割，将分割后的每一个子部分再以阈值范围来提取传导系统的分支部分。

三、获取狗的左心室心肌三维结构数据，生成左心室内膜层三角网格曲面，将左心室内膜层的离散点数据集F，经局部线性嵌入算法映射到二维嵌入空间P，映射后的点数据图记为图5；由于LLE算法的特点是保持各点与周围点的局部线性关系，因此，经该算法映射后的点数据仍将保持局部的位置特点，使各点之间的空间位置远近几乎不因映射而发生改变，这也为后面步骤的顺利进行奠定了良好的基础。

1)通过康奈尔大学发布公开的狗左心室的三维模型构建离散点集P的数据信息，采用基于分片投影和Delaunay三角化的三角剖分的方法建立了三维左心室内膜层的网格曲面。

算法根据点云的位置信息提取到相关的法矢和邻域信息，利用以上信息自动地在特征点附近对点云进行分块，然后对每个块进行投影三角剖分，最后利用点云的相关信息对无效的三角剖分进行删除，最终得到了很好的剖分结果。

2)利用高低维空间中相邻点间具有同样的相邻或者密切相关位置关系，应用局部线性嵌入(LLE)算法，将左心室的内膜层构建点集F映射到二维平面点集P；

a.计算出每个三维模型构建点F_i的k个近邻点，该近邻点数目k，由最近邻域问题求出：通常是指在笛卡尔坐标系中欧几里德距离下，包含M个点的数据集 中，找到一个点q ∈R^D的k个最近点(0＜k＜M)，把相对于所求点F_i距离最近的k个点规定为所求点的k个近邻点。k的取值需要根据后面步骤再做调整，在该步骤中可先在取值范围内任取一值。

b.计算出点集F的局部重建权值矩阵W_ij。

c.将所有的数据点F_i映射到在流形中表示内部全局坐标的低维向量P_i上。

3)获取P中的边界点集Pb，利用三维数据点集F与二维数据点集P具有一一对应关系，获得与Pb对应的狗的左心室内膜层的三维边界点集Fb，具体步骤如下：

a.建立极坐标，确定新坐标原点。

b.将点P_i(P_i∈P)极坐标化。

c.由于α_i∈[-π/2，π/2]，将区间[-π/2，π/2]进行细分，取每个小区间内极半径最大的点作为边界点。

d.设置区间个数评判数k，当区间点数k_i＜k时，将这一小区间与下一个区间合并直至满足判断条件。

e.因为 

，即F与P一一对应，因此同时获得了二维数据点集边界Pb对应的三维数据点的边界Fb。

四、优化匹配点集P与C并将C嵌入进P，从而得到H：

匹配点集P与C，即调整图4，使之可以相同比例尺整合到图5中，并将其整合到图5中，生成新图。这里的整合，是指将图4和图5中的点置于同一坐标系中，并将点集P和点集C中的点置于新的图中。

五、再一次运用LLE算法的思想，将H中的数据点映射回三维结构中。

此时，心室传导系统的线状图得以完整地展现在心室三维结构中。这时的映射，只是图C中部分点的映射，即图C中图4中的点，经过LLE算法进行映射。

本发明的优点在于，成功地将曲面与平面相互映射，使得在平面上的数据完整地再现于曲面的对应几何位置。利用狗心室传导系统真实的图像数据，以及狗心室的三维数据模型，将流形学习的方法应用于图像的曲面三维点集与平面二维点集之间的映射，并利用局部线性 嵌入算法可保持点信息的局部线性关系这一特点，构建狗心室传导系统的三维几何，对心室传导系统的真实刻画和再现，因此本发明较之以往的心室传导系统模型，有着本质的不同和明显的优势。

附图说明

图1是背景技术中所述局部线性嵌入算法(即LLE算法)的结构示意图；

图2是本发明的流程示意图；

图3是经特殊显影技术处理后的狗的左心室传导系统的线状图；

图4是提取的心室传导系统的二维线状图；

图5是心室内膜点的二维映射点；

图6是基点切平面投影点的初始三角化；

图7是基点切平面投影点三角剖分的检测(不合理的三角剖分)；

图8是基点切平面投影点三角剖分的修正(修改后的三角剖分)；

图9不同基点切平面投影点三角剖分的合成；

图10是狗左心室三角剖分曲面图；

图11是将图4与图5整合后的新图；

图12最终的心室传导系统的三维结构图。

具体实施方式

结合附图1～12及实施例，对本发明作进一步的说明：

一、进行离体的动物实验，将狗的左心室作特殊的显影技术处理，使其可以用肉眼看到内部的心室传导系统的线状结构(如图3所示)。

为后面的映射做准备，选取一条可以在映射后的二维数据图中找到对应的线或选取一些位置或形状上有特征的点，这些特征点及在心室内膜中对应的点须作好标记，以便后面步骤在曲面与平面之间相互映射过程中，定位心室传导系统线状结构。

二、提取该步骤一中的心室传导系统线状的图，即左束支及浦肯野纤维数据点，该点集记为C(如图4所示)，提取线状图的过程可分为以下两步骤：

1)先以灰度阈值来提取传导系统的轮廓部分，并以阈值范围进行区域分割。

2)将分割后的每一个子部分再以阈值范围来提取，直到将分支部分也提取出来；当分支部分达到末端附近时，分支的颜色几乎与心室肌颜色相同，再用自动的方式几乎不能有效地提取数据，这时需要更改分支部分灰度值，使之区别于背景色并将其提取出来，由于提取数据主要是提取点的位置信息而不是颜色信息，因此适当地修改线状图中的RGB颜色不会对结果有任何影响。

三、狗左心室三维数据点集P到二维空间F的嵌入，分为如下步骤：

1)对于狗左心室的三维模型构建离散点集P的数据信息，采用基于分片投影和Felaunay三角化的三角剖分的方法建立三维左心室内膜层的网格曲面，其具体步骤包括：形状索引、切平面的选取、法向一致化处理、点云预处理、基点和投影点的选择、三角化。

a.形状索引

为了估计曲面每个点的曲率，用二次曲面：

f(x，y)＝ax²+by²+cxy+dx+ey+f    (3)

来拟合点和它的邻域，通过最小二乘方法来估计二次曲面的参数，然后通过微分几何来计算曲面法线、高斯曲率、平均曲率和主曲率，通过遍历，可以得到其它点的形状索引，形状索引是指曲面点的形状量化测量，其定义如下：

$S_i\left(p\right)=\frac{1}{2}-\frac{1}{x}{\tan }^{-1}\frac{k_1\left(p\right)+k_2\left(p\right)}{k_1\left(p\right)-k_2\left(p\right)}---\left(4\right)$

其中，k1、k2分别为曲面的最大和最小主曲率。根据定义，所有的形状被映射到区间[0，1]，大的形状索引对应凸曲面，小的形状索引对应凹曲面。

b.切平面的选取

对于点云数据集P∈R^D的任意一点，可以通过其k邻域得到其切平面，问题可以描述为：对于点V的K邻域S＝{p_i|i＝1，2，...，K}，确定一个平面T_p＝(V，n)，使得邻域点到平面的距离和最小，得到近似的法线矢量n，这个平面就是点V的近似切平面，法线矢量n可以通过主成份分析方法得到，考察协方差矩阵：

$\mathrm{CVM}=\left[\begin{array}{ccc}\mathrm{\Σ}{(p_{\mathrm{ix}}-V_x)}^2& \mathrm{\Σ}(p_{\mathrm{ix}}-V_x)(p_{\mathrm{iy}}-V_y)& \mathrm{\Σ}(p_{\mathrm{ix}}-V_x)(p_{\mathrm{iz}}-V_z)\\ \mathrm{\Σ}(p_{\mathrm{iy}}-V_y)(p_{\mathrm{ix}}-V_x)& \mathrm{\Σ}{(p_{\mathrm{iy}}-V_y)}^2& \mathrm{\Σ}(p_{\mathrm{iy}}-V_y)(p_{\mathrm{iz}}-V_z)\\ \mathrm{\Σ}(p_{\mathrm{iz}}-V_z)(p_{\mathrm{iy}}-V_y)& \mathrm{\Σ}(p_{\mathrm{iz}}-V_z)(p_{\mathrm{iy}}-V_y)& \mathrm{\Σ}{(p_{\mathrm{iz}}-V_z)}^2\end{array}\right]---\left(5\right)$

由于协方差矩阵CVM是3×3实对称矩阵，特征向量必然存在，可以通过雅可比(Jacobi) 方法求得全部特征值λ₁，λ₂，λ₃和单位特征向量v₁，v₂，v₃，其中最小特征值所对应特征向量与CVM表示的微切平面法向量方向平行，因此可以通过计算协方差矩阵CVM的最小特征值来确定微切平面的法向量n，最后要统一微切平面的主向量方向。

c.法向一致化处理

在得到的每个点的法线矢量中，有的指向物体内部，有的指向物体外部，因此必须对得到的法线矢量进行一致化处理，使其都指向物体的外部。

设测点x_i和x_j是曲面上距离很近的两点，且x_i，x_j∈X，则相应法线矢量的点积n_i·n_j≈+1，否则n_i或者n_j应当被反向，法线矢量调整采用从一点开始，进行法线矢量方向传播的方法。

d.点云预处理

由于获取的点云集是散乱点，预处理是从点云的位置信息提取出相关的法线和邻域信息，点云的三角化是利用预处理得到的信息对点云进行三角剖分。

I.求得每个点的准确的k邻域，把点云数据分割为不同的长方体，采用哈希表提高搜索的速度；

II.根据公式(4)，得到每个点的形状索引值；

III.根据每个点的k邻域，利用公式(2)得到每个点的法线矢量，并对法线矢量进行一致化处理。

e.基点的选择

基点是点云分片中重要的点，点云的分片就是以基点和基点的法线作为基础进行处理的，它主要分为初始基点和其它基点；在处理中首先得到所有的初始基点，然后对初始基点进行处理，而其它基点的选择则是每得到一个基点则进行处理，处理后再判断是否有符合条件的基点，若有，则要进一步处理，直到没有这样的点存在。

f.投影点的选择算法

给定一个基点，使用以下两个条件决定哪些点将被投影到该基点的切平面上。

I.条件1

II.大部分的曲面都是局部连续的，假设选择了一个基点，它的法线矢量为n_b，对于其它点法线矢量n_i，当它们满足θ＜θ_t时，这些点将有可能被投影到切平面上；

其中，θ＝arccos(n_b·n_i)。θ为基点法线矢量和其它点法线矢量的夹角，θ_t为决定投影点 的角度阈值，其值小于π/2，当三维的点投影到二维平面上时，位置关系不发生改变；

III.条件2

假设基点到其它点的方向矢量为n_bi，w为n_bi和n_b之间的夹角，通过设置w的阈值ω_t，可以进一步确定投影点，它的值介于π/2和π之间，

w＜ω_t，其中，w＝arccos(n_b，n_bi)。

g.三角化

在把点投影到基点的切平面上后，可以把三维三角剖分转换为二维三角剖分：

首先，对所有点以不出现相交线的原则，进行粗略的三角剖分：

求出投影的到基点切平面上所有点(集合J)的外围边界点(集合K)，并将所有外围边界点依次连接得边界集L；

I.由集合K出发，求出集合K与集合J-K距离最小的点对。{(s，t)|s∈K，t∈J-K}，验证线段st是否与边界集L中已有的线段相交，若不相交则将点t扩充到集合K中K＝K∪s，并连接(s，t)扩充至L中，否则取直至满足判别条件的距离次小点s；

II.重复执行b，直至J＝K；

III.对边界集L进行检测，是否存在四变形(五边形)，若存在则将其分别化为两个(三个)三角形；

其次，判断边界集L中的三角形是否合理，合理的三角形是其外界圆不包含集合B内的其他定点(如图5所示)，如若三角形的外接圆包含一个点其他点，则重新对这四个点进行三角化(如图6所示)。

最后，对不同基点的投影图进行拼接(如图7所示)。

在二维得到的三角剖分结果都是有效的，然而在映射回三维曲面后，由于点选择的关系，并不是每个三角形都是有效的，就会出现冗余的三角剖分，因此必需对这些三角形进行删除处理，从而使最终的三角剖分结果接近最优。

该过程根据点的位置信息提取相关的法向矢量和邻域信息，并自动地在特征点附近对点进行分块，然后对各块进行投影三角剖分，之后删除无效的三角剖分，最终得到了很好的剖分结果(如图8所示)。

2)利用高低维空间中相邻点间具有同样的相邻或者密切相关位置关系，应用局部线性嵌入(LLE)算法，将左心室的内膜层构建点集P映射到二维平面点集F。

a.寻找相对于所求点的欧氏距离最近的k个样本点规定为所求点的近邻点，k不妨在取 值范围内预先设定一个值；

b.由每个样本点的近邻点计算出该样本点的局部重建权值矩阵；

c.由该样本点的局部重建权值矩阵及其近邻点计算出该样本点的输出向量；

具体操作如下：

定义一个误差函数：

$\min \mathrm{\ϵ}\left(W\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{|x_i-\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{w}_j^i{x}_{\mathrm{ij}}|}^2---\left(6\right)$

式中：x_ij(j＝1，2，L，k)为x_i的k个近邻点；w_j为x_i与x_i之间的权值，且要满足 $\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{w}_j^i=1$ 条件，在求取局部重建权值矩阵W时应使用式(6)值达到最小，即构造一个k×k的局部协方差矩阵

$Q_{\mathrm{jm}}^i={(x_i-x_{\mathrm{ij}})}^T(x_i-x_{\mathrm{im}})---\left(7\right)$

用于求取W矩阵，对于式(6)，结合 $\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{w}_j^i=1,$ 并采取拉尔朗日乘子法，既可求出局部最优化重建权值矩阵：

$w_j^i=\underset{m=1}{\overset{k}{\mathrm{\Σ}}}{\left(Q^i\right)}_{\mathrm{jm}}^{-1}/\underset{p=1}{\overset{k}{\mathrm{\Σ}}}\underset{q=1}{\overset{k}{\mathrm{\Σ}}}{\left(Q^i\right)}_{\mathrm{pq}}^{-1}---\left(8\right)$

通常，式(8)中的Qⁱ是一个奇异矩阵，必须对Qⁱ进行正则化，引入一个正则化参数r，则

Qⁱ＝Qⁱ+rI    (9)

式中，I为k×k的单位矩阵，将所有的样本点映射到地位空间的数据，并使输出数据在低维空间中保持原有的拓扑结构，为此，需构造一个损失函数，映射过程中必须使损失函数值达到最小，即

$\min \mathrm{\ϵ}\left(Y\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{|y_i-\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{w}_j^i{y}_{\mathrm{ij}}|}^2---\left(10\right)$

式中：y_i为x_i的输出向量；y_ij(j＝1，2，L，k)为y_i的k个邻近点，且满足条件： $\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{y}_i=0,$ $\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{y}_i{y}_i^T=I,$ 其中I为m×m的单位矩阵；然后求取的最优解y_i应使式(10)值达到最小， w_i ^j(i＝1，2，L，M)可以存储在M×M的稀疏矩阵W中，当x_j是x_i的近邻点时， $W_{\mathrm{ij}}=w_j^i;$ 否则，W_i ^j＝0，则式(10)可重写为

$\min \mathrm{\ϵ}\left(Y\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{G}_{\mathrm{ij}}{y}_i^T{y}_j---\left(11\right)$

式中，G为一个M×M得对称矩阵，其表达式为

G＝(I-W)^T(I-W)    (12)

由式(11)可知，要是损失函数值达到最小，则取Y为M的最小m个非零特征值的特征向量。

LLE算法中共有3个参数需要设置，即近邻点的个数k、输出维数m和正则化参数r，k的选取在算法中起到关键因素，如果k取值太大，LLE不能体现局部特性；反之，，LLE不能保持样本点在低维空间中的拓扑结构，而且，k的取值决定着二维数据点集F的空间分布，因此需要在后面的工作对k的取值作进一步的修正。

3)获取P中的边界点Pb(Pb＝{p_i|p_i∈P，p_i为点集P_{所构成图形的边界点}})，利用三维数据点集F与二维数据点集P具有一一对应关系，获得与Pb对应的狗左心室内膜层的三维边界点集Fb。

LLE算法在将数据点由三维空间映射到二维空间时，由于该算法的特点是可以保持数据点的局部线性关系，因此在三维数据点的局部位置关系也被映射到二维数据点中；这时求出二维边界点，一方面为步骤4的优化匹配点集做参照；另一方面，可以由一一映射关系得到三维边界点，从而为进行动物实验时如何解剖狗的左心室提供边界。

a.建立极坐标，确定新坐标原点；

由于映射到二维空间的点集P在空间上分布不均匀，因此为更快速准确地找到二维边界，不妨将二维点集生成图形的几何中心作为坐标原点：

$(x_o,y_o)|x_o=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{x}_i/M,y_o=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{y}_i/M,$ 其中M为点集中点的个数(13)

b.将直角坐标转换为极坐标。以上一步求得的(x_o，y_o)为新的坐标原点，将点集中的各点P_i(P_i∈P)相对于新坐标原点极坐标化，获得各点对应的辐角α_i和半径r_i如下：

α_i＝arctan((y_i-y_o)/(x_i-x_o))，    (14)

$r_i=\sqrt{{(y_i-y_o)}^2+{(x_i-x_o)}^2}---\left(15\right)$

c.由于二维点集P分布在(x_o，y_o)的四周，对于任意的P_i∈P，满足α_i∈[-π，π]，将区间[-π，π]进行n等分，对于第i个小区间 $A_i=\{P_j\∈P|-\mathrm{\π}+\frac{2(i-1)}{n}\mathrm{\π}\≤{\mathrm{\α}}_j\≤-\mathrm{\π}+\frac{2i}{n}\}$ 内的点数为 m_i，则取每个小区间内极半径最大的点作为该区间的边界点，记为：

${\mathrm{Pb}}_i=\underset{1\≤j\≤m_i}{\max }\left\{r_j\right|P_j\∈A_i\}---\left(16\right)$

d.由于二维数据点集P的分布具有不规则性，使得一部分小区间包含的点过少，从而，并不是C中获得的所有小区间的边界点都是点集真正的边界点，故设置区间个数评判数t，当区间点数m_i＜t时，将这一小区间A_i与下一个区间A_i+1合并：

A_i＝A_i∪A_i+1，      (17)

m_i＝m_i+m_i+1         (18)

且直至满足判断条件m_i≥t；

e.因为LLE方法将三维数据点集F映射成二维数据点集P，L：F→P，满足 $F\↔P,$ 且有L(F_i)＝P_i，i＝1，...，M，因此获得了二维数据点集边界Pb同时得到了三维数据点集F的边界Fb。

4.优化匹配点集P与C并将C嵌入进P得到H；

1)首先将点集C二维向量[x y]^T进行三维表示[x y 1]^T以便于数学变换。通过缩放变换K、旋转变换R、平移变换T后得到

C^*＝ToRoK(C)，      (19)

并对C^*提取二维信息，即，将任一点由[x y 1]^T变回[x y]^T，处理完后仍记为C^*；

$T=\left[\begin{array}{ccc}1& 0& x_d\\ 0& 1& y_d\\ 0& 0& 1\end{array}\right],---\left(20\right)$

$K=\left[\begin{array}{ccc}{k}_1& 0& 0\\ 0& k_2& 0\\ 0& 0& 1\end{array}\right],---\left(21\right)$

$R=\left[\begin{array}{ccc}\cos \mathrm{\θ}& -\sin \mathrm{\θ}& 0\\ \sin \mathrm{\θ}& \cos \mathrm{\θ}& 0\\ 0& 0& 1\end{array}\right],---\left(22\right)$

其中k₁，k₂分别为x，y方向的缩放系数，θ为偏转角，x_d，y_d分别为x，y方向的平移位移；

2)分别在点集P和点集C^*选取用于匹配的特征点，其特征点集Z的选取须至少满足以下条件之一：

a.点集Z为C^*与P具有相同特征的点： 

b.由于所获得数据点集F仅是狗左心室内膜层的构建信息，并不包含传导系统的构建信息，所以任何在C^*或P中特征点z_i，应是具有结构代表性的点，即该点处于C^*或P中的特定位置，或该点与邻近的点构成特殊形状。算法中采用特殊边界点Fb和凹陷的边界点作为特征点集Z。

3)将由点集P构成的分布图与点集C^*构成的分布图进行匹配，匹配的过程分为特征点匹配与分布图像匹配：

a.特征点集Z匹配：

由点集 

的定义，特征点的匹配须满足

$\min \underset{i=1}{\overset{K}{\mathrm{\Σ}}}{({\stackrel{r}{a}}_i-{\stackrel{r}{b}}_i)}^2$

b.分布图像匹配，其判别标准如下：

I.点集C^*的二维分布图形不超出点集P边界Pb；

II.点集C^*的二维分布图形占据点集P分布图形的较大比例；

III.均匀占据，不能有一侧或一部分较大的空缺。

由于分布图形的判别标准的程序实现较困难，这一步的匹配采用人眼的识别方法，尝试改变不同的局部邻近点个数k，进行特征点匹配，当以上评判标准被极大地满足后，将获得最优的邻近点个数k。

4)将C^*与嵌入点集P构建点集H，形成新的二维图，这里的嵌入，是指将C^*和P中的点置于同一坐标系中，并将原C^*和原P中的点置于新图中(如图9所示)。

5.再一次运用LLE算法的思想，通过数据点间的局部线性关系，将H中的数据点映射回到三维结构中，由于图5中的点就是三维结构中的点，因此图5中的点不需要再做映射，该处的映射只是H里面点集C中的点，即点集经过LLE算法进行的映射(如图10所示)。具体的步骤如下：

1)将图4中的点作为点数据 其邻近点集为图5中的点。即以心室内膜点的二维映射点为邻近点，找到心室传导系统二维线状图的每个点的邻近点。

2)计算权值W_ij，通过W_ij与邻近点构造数据点 

首先，在新图中，将图4中每个点在图5中的邻近点权值矩阵W_ij求出，即用图5中的 部分点分别表示出图4中的每个点；

然后，由于在3中经过LLE算法处理，心室的内膜点都有在二维平面图5与之一一对应的点，因此又可以将图5中的点用心室的内膜点表示。

这样，就用心室的内膜点表示了图4中的每个点，即用心室的内膜点表示了心室传导系统二维线状图上的点，于是，心室传导系统的二维线状图被表示在了心室的内膜层三维几何曲面中(如图11如示)，从而完成了心室传导系统三维几何构型的构建(如图2如示)，对心室传导系统真实刻画、再现(如图12如示)。

以上实施例仅供说明本发明之用，而非对本发明的限制，本技术领域的普通技术人员，在不脱离本发明的精神及范围的情况下，作出各种等同变换或变化的技术方案属于本发明的保护范畴，由各项权利要求限定。

Claims (4)

1.一种基于局部线性嵌入法构建狗左心室传导系统的三维几何结构，其特征在于由下述步骤完成心室传导系统的构建：

一、将狗的左心室作特殊的显影技术处理，使其可以清晰地看到内膜层附近心室传导系统的线状结构，固定并拍摄图片；

二、提取步骤一得到的图片中的心室传导系统线状图，即左束支及浦肯野纤维，将提取的离散点集记为点集C，提取的过程是半自动过程，步骤如下：

1)先以灰度阈值来提取传导系统的轮廓，即颜色较深的部分；

2)再进行区域分割，将分割后的每一个子部分再以阈值范围来提取传导系统的分支部分；

三、获取狗的左心室心肌三维结构数据，生成左心室内膜层三角网格曲面，将左心室内膜层的离散点数据集F，经局部线性嵌入算法映射到二维嵌入空间P，映射后的点数据仍将保持局部的位置特点；

1)通过狗左心室的三维模型构建离散点集P的数据信息，采用基于分片投影和Delaunay三角化的三角剖分的方法建立三维左心室内膜层的网格曲面；

2)利用高低维空间中相邻点间具有同样的相邻或者密切相关位置关系，应用局部线性嵌入算法，将左心室的内膜层构建点集F映射到二维平面点集P；

a.计算出每个三维模型构建点F_i的k个近邻点，该近邻点数目k，由最近邻域问题求出：通常是指在笛卡尔坐标系中欧几里德距离下，包含M个点的数据集中，找到一个点q∈R^D的k个最近点(0＜k＜M)，把相对于所求点F_i距离最近的k个点规定为所求点的k个近邻点；k的取值需要根据后面步骤再做调整，在该步骤中可先在取值范围内任取任一值；

b.计算出点集F的局部重建权值矩阵W_ij；

c.将所有的数据点F_i映射到在流形中表示内部全局坐标的低维向量P_i上；

3)获取P中的边界点集Pb，利用三维数据点集F与二维数据点集P具有一一对应关系，获得与Pb对应的狗的左心室内膜层的三维边界点集Fb，具体步骤如下：

a.建立极坐标，确定新坐标原点；

b.将点P_i(P_i∈P)极坐标化；

c.由于α_i∈[-π/2，π/2]，将区间[-π/2，π/2]进行细分，取每个小区间内极半径最大的点作为边界点；

d.设置区间个数评判数k，当区间点数k_i＜k时，将这一小区间与下一个区间合并直至满足判断条件；

e.因为

，即F与P一一对应，因此同时获得了二维数据点集边界Pb对应的三维数据点的边界Fb；

四、优化匹配点集P与C并将C嵌入进P，从而得到H：匹配点集P与C，调整心室传导系统的二维线状图，将心室传导系统的二维线状图和心室内膜点的二维映射点置于同一坐标系中，并将点集P和点集C中的点生成新图；

五、将H中的数据点映射回三维结构中。

2.根据权利要求1所述的一种基于局部线性嵌入法构建狗左心室传导系统的三维几何结构，其特征在于所述提取心室传导系统线状图可分为以下两个步聚：

1)先以灰度阈值来提取传导系统的轮廓部分，并以阈值范围进行区域分割；

2)将分割后的每一个子部分再以阈值范围来提取，直到将分支部分也提取出来。

3.根据权利要求1所述的一种基于局部线性嵌入法构建狗左心室传导系统的三维几何结构，其特征在于所述狗的左心室三维数据点集F到二维空间P的嵌入分为如下步骤：

1)对于狗左心室的三维模型构建离散点集P的数据信息，采用基于分片投影和Delaunay三角化的三角剖分的方法建立了三维左心室内膜层的网格曲面，其具体步骤包括：形状索引、切平面的选取，法向一致化处理，点云预处理，基点和投影点的选择或三角化；

2)利用高低维空间中相邻点间具有同样的相邻或者密切相关位置关系，应用局部线性嵌入(LLE)算法，将左心室的内膜层构建点集F映射到二维平面点集P；

3)获取P中的边界点Pb(Pb＝{p_i|p_i∈P，p_i为点集P所构成图形的边界点})，并利用三维数据点集F与二维数据点集P具有一一对应关系，获得与Pb对应的狗左心室内膜层的三维边界点集Fb。

4.根据权利要求1所述的一种基于局部线性嵌入法构建狗左心室传导系统的三维几何结构，其特征在于所述优化匹配点集P与C并将C嵌入进P得到H；

1)首先将点集C二维向量[x y]^T进行三维表示[x y 1]^T以便于数学变换。通过缩放变换K、旋转变换R、平移变换T后得到

C^*＝ToRoK(C)，                        (19)

并对C^*提取二维信息，即，将任一点由[x y 1]^T变回[x y]^T，处理完后仍记为C^*。

$T=\left[\begin{array}{ccc}1& 0& x_d\\ 0& 1& y_d\\ 0& 0& 1\end{array}\right],---\left(20\right)$

Images