CN103679801A - 一种基于多视角x光片的心血管三维重建方法 - Google Patents

Abstract

本发明提供一种基于多视角X光片的心血管三维重建方法，包括了四个步骤：图像增强，用来增强图像，加大X光片对比度；血管提取，从X光图像中提取出心血管；中心线追踪，细化提取出的血管；三维重建，结合几个视角的中心线采用最优化方法重建血管的三维结构。本发明可完全基于GPU进行重建和渲染，具有重建速度快、重建输入少、重建效果逼真的特点。

Description

一种基于多视角X光片的心血管三维重建方法

技术领域

本发明涉及一种基于多视角X光片的心血管三维重建方法。

背景技术

心血管疾病对人类健康构成了巨大威胁，其中以心血管狭窄最为常见。随着医学影像技术的发展，该技术在血管狭窄诊断上也得到了应用，其中心血管造影图像已经成为心血管疾病诊断的金标准。

传统的诊断方法中，医生如果需要查看病人的心血管的三维结构数据，需要病人做核磁共振或者旋转断层扫描，并从扫描或者共振的数据中手工分割心脏区域，才能得到对应的三维结构，整个过程耗时长、人工参与多、对病人和医生的辐射都较大，不能很好的满足医患的需求。

利用多角度二维心血管图像重建的心血管三维结构不仅能够更好地对心血管狭窄进行客观诊断，而且还能够实现病变的量化描述，对心血管狭窄防治具有重要意义。造影心血管三维重建包括两个阶段：1、造影图像心血管二维信息提取；2、心血管三维结构重建。本文工作主要是围绕心血管骨架信息提取、直径信息提取及三维重建的初步研究展开的。

为了解决上述问题，本发明提出了基于多视角X光片的心血管三维重建方法。本方法能极大的提高三维重建的效率，并能够保证三维重建的精度。与此同时，本方法只需要三个角度、每个角度一个心动周期的X光片为输入源，具有输入少、获取输入方便、重建速度高的特点。依据本发明重建三维结构，能够极大的减轻造影过程对医患的危害。

发明内容

本发明解决的技术问题是：在尽量减少将医患曝光在X射线的前提下，利用尽量少的视角的光片，以较快的速度重建心血管的三维结构。

本发明采用的技术方案为：一种基于多视角X光片的心血管三维重建方法，包括以下四个步骤：

步骤（1）、图像增强：针对X光成像片对比度低，动态范围窄，成像模糊，成像层级少的缺点，采用基于Retinex理论的图像增强方法，并对处理后的图像引入gain\offset方法进行修正，使图像能够在人眼可见的范围内显示；

步骤（2）、图像中血管结构的提取：根据（1）中增强后的图像，选取不同大小的海森矩阵对图像做卷积，并采用特殊的图像过滤器计算当前像素点成为一个血管结构的可能性，最终根据该可能性构成一张代表血管结构的二值图；

步骤（3）、血管中心线提取：根据步骤（2）中计算得到的血管，采用基于局部特征的方法，将血管外层结构层层剥离，最终只留下内部代表血管结构的单像素中心线；

步骤（4）、血管的三维重建：使用步骤（3）得到的血管中心线，结合三视角的中心线特征，根据X光的成像原理，将X光发生器和成像板之间的空间划分为均匀的多个体素，并根据每个体素在光片上的位置设定每个体素的值，最终根据最优化理论，在考虑连续性的情况下求得所有体素的最优值。

本发明的原理在于：

（1）通过基于Retinex理论的图像增强方法，不同于传统的图像增强算法，如线性、非线性变换、图像锐化等只能增强图像的某一类特征，如压缩图像的动态范围，或增强图像的边缘等，Retinex可以在动态范围压缩、边缘增强和颜色恒常三方面达到平衡，因此可以对各种不同类型的图像进行自适应性地增强。

（2）为了能够对血管的结构进行快速三维重建，需要剥离出图像中血管状结构，去除杂乱的非血管结构，通过使用不同尺度的海森矩阵，能够保留源图像上在径向增强的部分，去除非径向部分，在保留血管结构的同时能极大的去除噪声，为后续处理提供良好的基础。

（3）为了对三维血管进行重建，考虑到单视角必然会有遮挡的现象，需要多个视角的数据共同作用。本方法依据X光成像的特点，将成像空间划分为空间体素，并对所有体素根据其行为的不同划分一定的标价，并最终使用信任传递算法求得所有标价的最小和，从而取得空间血管的三维走形。

本发明与现有技术相比的优点在于：

1、所需视角少，对数据源的要求低，数据的获得方式容易，所有有能力做X光片的医院均可采用本方法进行三维重建。

2、高度个性化，由于重建的数据源来源于每个病人的真实X光片，能够真实的再现病人拍摄光片时心血管的状态和结构。

3、快速，和其他方法动则数小时的重建时间相比，使用本方法的重建只需要数分钟即可完成，其中图像增强处理阶段的效果也可作为医生更好的诊断的依据。

附图说明

图1为基于多视角X光片的心血管三维重建方法的处理流程图；

图2为图像增强之后对比图；

图3为分割后血管与原始血管对比示意图；

图4为中心线追踪结果示意图；

图5为本文方法模拟的X光成像示意图；

图6为本文方法模拟的X光机空间划分示意图；

图7为最终重建结果。

具体实施方式

图1给出了基于多视角X光片的心血管三维重建方法的总体处理流程，下面结合其他附图及具体实施方式进一步说明本发明。

本发明提供一种基于多视角X光片的心血管三维重建方法，主要步骤介绍如下：

1、多尺度Retinex图像增强方法

该方法基于人眼对外界色彩的辨识特征，对低对比度、高动态范围以及模糊的图像进行增强。Retinex这个词是由视网膜（Retina）和大脑皮层（Cortex）两个词组合构成的。Retinex理论主要包含了两个方面的内容：物体的颜色是由物体对长波、中波和短波光线的反射能力决定的，而不是由反射光强度的绝对值决定的；物体的色彩不受光照非均匀性的影响，具有一致性。

根据Retinex理论，人眼感知物体的亮度取决于环境的照明和物体表面对照射光的反射，其数学表达式为：

I(x，y)＝L(x，y)*R(x，y)    （8）

式中：I（x，y）代表被观察或照相机接收到的图像信号；L（x，y）代表环境光的照射分量；R（x，y）表示携带图像细节信息的目标物体的反射分量。

将上式两边取对数，则可抛开入射光的性质得到物体的本来面貌，即有关系式：

Log[I(x,y)]＝Log[L(x,y)]+Log[R(x,y)]    (9)

将Retinex技术运用到图像处理上，就是针对我们现在已经获得的一副图像数据I（x，y），计算出对应的R（x，y），则R（x，y）认为是增强后的图像，现在的关键是如何得到L（X，Y）。Retinex理论的提出者指出这个L（x，y）可以通过对图像数据I（x，y）进行高斯模糊而得到。

但是单尺度的Retinex增强有很多问题，如偏色等，我们采用了多尺度的带色彩恢复的Retinex增强方法。其具体步骤为：

1）分别计算出Log[R（x，y）]中R/G/B各通道数据的均值Mean和均方差Var。

2）利用类似下述公式计算各通道的Min和Max值。其中Mean为各通道数据的均值，Var为均方差:

Min＝Mean-DynamiC*Var；

Max＝Mean+Dynamic*Var；          (10)

3）对Log[R（x，y）]的每一个值Value，进行线性映射：

R(x，y)＝(Value-Min)/(Max-Min)*(255-0)    (11)

同时要注意增加一个溢出判断，即：

if(R(x，y)＞255)R(x，y)＝255；elseif(R(x，y)＜0)R(x，y)＝0    (12)经过本方法处理后的图像对比图见图2。

2、X光片血管结构的提取

对X光片而言，由于血管本身的低对比度和柔软的特性，从光片中将血管提取出来是很难的。最大的挑战是如何能够在增强或者提取出血管结构的同时，尽量减少非血管结构的影响。

我们的方法中，在取得高对比度的图像之后，我们使用了一种基于海森矩阵的方法来增强光片中的血管结构。这种方法基于血管的局部特征。

定义I的海森矩阵为：

$H=\left(\begin{array}{cc}{I}_{\mathrm{xx}}& I_{\mathrm{xy}}\\ I_{\mathrm{yx}}& I_{\mathrm{yy}}\end{array}\right)---\left(13\right)$

其中I_xx为对x的二阶偏导数，I_yy为对y的二阶偏导数，I_xy＝I_yx为x，y的混合偏导数。

对图像中所有点都计算器海森矩阵的值，该矩阵的特征值和特征向量可以用来衡量某个点成为血管的可能性。利用海森矩阵提取血管结构的步骤如下：

（1）、预计算偏导数

根据当前sigma的大小，计算高斯核的大小，并构建对应的高斯核函数。构造完成后，计算所构造的高斯函数的二阶偏导和二阶混合导数，并将计算后的值传输到GPU上以便利用CUDA进行加速。

（2）、海森矩阵计算

对图像上所有点，根据预先计算的偏导数，与当前像素进行卷积计算，并将卷积的值代替当前像素值，存入一副新的图像。依次计算三个偏导数所对应的函数与图像卷积的值，并保存。根据保存下来的三个偏导数的值，即可构建当前二维图中所有像素的海森矩阵的值。对每个像素的海森矩阵，计算其特征值及特征向量。

（3）、参数估计

我们定义以下变量来对像素点进行估计。

$R_b=\frac{{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1},\left|{\mathrm{\&lambda;}}_2\right|<<\left|{\mathrm{\&lambda;}}_1\right|---\left(14\right)$

并且，

S²＝，λ₁ ²+λ₂ ²       (15)

其中，λ₁和λ₂是当前点的两个特征值。

最终，我们定义，

$v_{\mathrm{\&sigma;}}\left(s\right)=\left\{\begin{array}{cc}0,& \mathrm{if}{\mathrm{\&lambda;}}_1<0\\ \exp (-\frac{R_b^2}{2{\mathrm{\&beta;}}^2})(1-\exp (-\frac{s^2}{2c^2})),& \mathrm{else}\end{array}\right.---\left(16\right)$

来衡量某点成为血管的可能性。其中C为常数。该式计算出的值可以作为衡量每个像素点成为血管结构的可能性。我们另外初始化了一张新的同样大小的图，并设定阈值，将阈值小于一定范围的点社区，最终生成一张二值图，该图为1的区域表示此处为血管结构，为0的区域表示该点为非血管结构。

在取得二值图以后，我们还使用一个十字模板计算整张图像的联通性，其中联通点数小于一定阈值的线我们认为其也是杂质，进而进行剔除。

经过以上步骤之后，我们极大的减少了分割出来的血管的噪声，分割的血管和原始血管的对比参见图3。

3、血管中心线提取方法

得到血管结构的二值图像之后，我们使用了基于快速行程算法的中心线提取方法。这种方法本身是用来解决雅克比方程。

我们使用了一种基于多模板的快速行程算法以取得更高的效率。通过在二维空间内对每个点8临域的模板进行计算，然后在其中选择最符合要求的一个，能够取得更高的精度。

我们将二值图上的点划分为两类，分别是冻结点和窄带点。每个循环过程中，在窄带点的点集中选择距离最小的点，将它加入冻结点的点集，并重新计算距离。在整个过程中，窄带点往外扩散，同时冻结点也跟随者他们一起往外。最终，当所有点都成为冻结点后，方法收敛。

在追踪中心线的过程中，一条直线上的所有点都被记录下来，并且所有的分支点也被记录下来，最终我们得到一张包含所有分支点和所有血管线段的分支图。最终的结果可以参考图4。

4、三维重建方法

X光成像过程是将三维的物体投影到二维空间的一个过程，在投影过程中，很多的细节信息就丢失了，并且由于遮挡，也容易造成一些让人无法区分的信息。根据X光片进行三维重建，在数据源上可以分为两类，一类是双臂造影机，这种机器有双C型臂，能同时拍摄两张同部位的光片；另一类是单臂造影机，同一时刻只能拍摄同一部位一张造影片。造影方式的不同，决定了重建方式的不同。相比较而言，单臂造影剂造价更低，维护费用更少，使用的更广泛。因此。本文申请的方法主要关注在后者的重建方法上。

X光片的成像方式和普通的相机等的成像方式不同，它的光源是X射线源，成像面试X光机的增强器。整个成像的过程是透视投影，我们使用OpenGL模拟了这一过程，参见图5。

在三维世界中，被造影的血管应当处在光心附近，造影的每个视角的血管的交集应当就是要求得血管的三维位置。这些所谓的交集其实对应的是各个视角的血管的对应关系。在三维中，一个血管树可以认为是由多个血管的部分组成，每个部分都可以看成是由数量不一的采样点组成。所有这些点应当在X光射线源和增强器之间，同时，这些点也应当能较好的投影在X光片上。

我们的方法中，以上所述的每个三维采样点都能够被赋给一个成为三维血管的概率值。这个值得大小由三方面决定：

1）这个三维点能够投影在几个视角的光片上。投影在越多的光片上，证明这个点是血管点的概率越高，因为当三维中一个点无论怎么投影，都被投影在造影片的血管上的时候，这个点必定是三维血管的一部分。

2）该三维点的投影点离对应视角的光片上血管点的距离的长短。距离越短，证明这个点越是接近血管点，越有可能是被投影的血管。

3）考虑该三维点的临近点的连续性和周围临近点的距离，这主要是考虑到血管应当是连续的，而非单一的离散的个体。

我们的重建方法的具体流程是：

1）空间分层采样。我们使用三个视角的X光片进行重建，并且将遮挡最小的视角作为索引视角。我们将X光光心和增强器之间的空间划分为均匀的层，如图6所示。

2）概率计算。每个在索引视角上的点都可以对应到空间中的一条直线，这些直线与分割的层相交，被均匀的分割为空间的点。所以索引图像上的每一个二维点和空间中的一个层相结合，就可以确定空间中的一个三维点。所有三维点都被赋予一个成为三维血管的概率。

3）最优化计算。当所有点都被赋值完毕，我们的最终目的是在所有采样点中寻找能够使能量：

$E\left(f\right)=\underset{p\&Element;P}{\mathrm{\&Sigma;}}{D}_p\left(f_p\right)+\mathrm{\&lambda;}\underset{p,q\&Element;N}{\mathrm{\&Sigma;}}{V}_{p,q}(f_p,f_q)---\left(17\right)$

最小的组合，这个问题可以被归结为一个能量最小化问题，我们使用信任传递算法来解这个问题。在我们的方法中，我们定义V_p，q(f_p，f_q)为两个点p，q间的欧拉距离，而D_p(f_p)被认为是“色彩一致性”系数，我们将其定义为：

$D_p\left(f_p\right)=\frac{1}{(n-1)}\underset{i=2}{\overset{n}{\mathrm{\&Sigma;}}}{P}_i(x,y)---\left(18\right)$

其中P_ith(x，y)为点p在第i个视角的投影片上值，该值可以按照以下进行计算：

$P_{\mathrm{ith}}(x,y)=\left\{\begin{array}{cc}{W}_h,& p(x,y)\&Element;I_{\mathrm{ith}}\\ W_l,& \&cup;(p(x,y),1)\&NotElement;I_{\mathrm{ith}}\\ W_a& \mathrm{else}\end{array}\right.$

$W_a=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{Val}}_{\mathrm{ith}}(x,y)---\left(19\right)$

其中W_h为常量，代表较高的能量值；W_l为常量，代表较低的能量值；W_a为计算得出，Val_ith为投射到第i幅图像上的值，N代表被投射的图像的个数。

我们的算法中包含两个主要的步骤，信息传递和最小化能量求解，在信息传递过程中，我们定义当前的信息能量为：

V＝V_p-1+aminD+(1-a)V_pminD    (20)

其中e_i(p_t)为能量值；e_i一1(q)为上一邻居的能量值，γ为常量，D(p_t，q)礼为两点间距离，V(q)为当前点权值。

在能量最小化的过程中，我们定义能量为：

e_i(p_i)＝min[γD(p_i，q)+(1-γ)V(q)+e_i-1(q)]    (21)

其中e_i(p_t)为能量值；e_i-1(q)为上一邻居的能量值，γ为常量，D(p_t，q)为两点间距离，V(q)为当前点权值。

最终我们计算所有成组的血管的能量的最小值，并取得最优解，重建结果参见图7。

本发明未详细阐述的技术内容属于本领域技术人员的公知技术。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (8)

1.一种基于多视角X光片的心血管三维重建方法，其特征在于包括以下四个步骤：

步骤（1）、图像增强：针对X光成像片对比度低，动态范围窄，成像模糊，成像层级少的缺点，采用基于Retinex理论的图像增强方法，并对处理后的图像引入gain\offset方法进行修正，使图像能够在人眼可见的范围内显示；

步骤（2）、图像中血管结构的提取：根据步骤（1）中增强后的图像，选取不同大小的海森矩阵对图像做卷积，并利用当前像素点的特征值计算其特殊的图像过滤器成为一个血管结构的可能性，最终根据该可能性构成一张代表血管结构的二值图；

步骤（3）、血管中心线提取：根据步骤（2）中计算得到的血管，采用基于局部特征的方法，将血管外层结构层层剥离，最终只留下内部代表血管结构的单像素中心线；

步骤（4）、血管的三维重建：使用步骤（3）得到的血管中心线，结合三视角的中心线特征，根据X光的成像原理，将X光发生器和成像板之间的空间划分为均匀的多个体素，并根据每个体素在光片上的位置设定每个体素的值，最终根据最优化理论，在考虑连续性的情况下求得所有体素的最优值。

2.根据权利要求1所述的基于多视角X光片的心血管三维重建方法，其特征在于：步骤1使用基于视网膜理论的方法，增强图像，然后使用gain/offset方法对增强后的值进行计算，使其处于预定的范围内。

3.根据权利要求2所述的基于多视角X光片的心血管三维重建方法，其特征在于：血管结构的提取根据增强后的图像，使用海森矩阵对图像卷积，并使用计算出来每点的特征值和特征向量对该点进行估计；通过设定一定的阈值，过滤掉可能性较低的点，剩下的即是血管点。

4.根据权利要求3所述的基于多视角X光片的心血管三维重建方法，其特征在于：血管中心线提取，依据Fast marching算法，将血管点分类为Frozen点和邻居点，通过计算邻居点和Frozen点之间的距离将Frozen点中离候选点最近的加入邻居点集，直到所有Frozen点都被处理完毕，邻居点即是所计算的中心线的点。

5.根据权利要求4所述的基于多视角X光片的心血管三维重建方法，其特征在于：三维重建对每个空间采样点根据其投射到成像面的多少和投射结果的好坏赋予一个能量值，并通过最优化的方法求得所有点在空间能量的最优值，即是三维血管结构。

6.根据权利要求2所述的基于多视角X光片的心血管三维重建方法：其特征在于：血管结构提取进一步包括：

（1）在GPU上以每个线程对应预处理图像的每个点，计算各点的二阶偏导数和混合偏导数；

（2）用计算的偏导组成海森矩阵与当前图像做卷积，并计算特征值和特征向量；

（3）利用求得的特征值λ₁和λ₂，和特征向量根据公式（1）对像素点进行估计；

$R_b=\frac{{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1},\left|{\mathrm{\&lambda;}}_2\right|<<\left|{\mathrm{\&lambda;}}_1\right|---\left(1\right)$

并且，

S²＝λ₁ ²+λ₂ ²

最终，定义，

$v_{\mathrm{\&sigma;}}\left(s\right)=\left\{\begin{array}{cc}0,& \mathrm{if}{\mathrm{\&lambda;}}_1<0\\ \exp (-\frac{R_b^2}{2{\mathrm{\&beta;}}^2})(1-\exp (-\frac{s^2}{2c^2})),& \mathrm{else}\end{array}\right.---\left(2\right)$

来衡量某点成为血管的可能性；其中C为常数；该式（2）计算出的值作为衡量每个像素点成为血管结构的可能性；之后，使用阈值分割将低于一定阈值的点去除。

7.根据权利要求3所述的基于多视角X光片的心血管三维重建方法：其特征在于：血管中心线追踪方法进一步包括：

（1）使用多尺度的基于多模板的快速行程算法，利用当前点8临域的模板进行计算；

（2）将二值图上的点划分为两类，分别是冻结点和窄带点，窄带点的点集中选择距离最小的点，将它加入冻结点的点集，并重新计算距离；

（3）最终，当所有点都成为冻结点后，方法收敛。

8.根据权利要求4所述的基于多视角X光片的心血管三维重建方法：其特征在于：三维重建进一步包括：

（1）将X光机光心和接收器之间均匀划分为多层，并预备在每层上进行操作；

（2）根据空间中采样点的投射到光片上的位置和投射到光片的个数，对所有采样点都赋予一个能量值；

（3）在被赋予能量值之后使用信任传递算法，在考虑“颜色一致性”的基础上进行信任值传播，在考虑“连续性”基础上进行空间结构计算，最终求得能量最小值对应的血管；

（4）定义V_p，q(f_p，f_q)为两个点p，q间的欧拉距离，而D_p(f_p)被认为是“色彩一致性”系数，将其定义为：

$D_p\left(f_p\right)=\frac{1}{(n-1)}\underset{i=2}{\overset{n}{\mathrm{\&Sigma;}}}{P}_i(x,y)---\left(3\right)$

其中n为索引视角的编号，P_ith(x，y)为点p在第i个视角的投影片上值，该值可以按照以下进行计算：

$P_{\mathrm{ith}}(x,y)=\left\{\begin{array}{cc}{W}_h,& p(x,y)\&Element;I_{\mathrm{ith}}\\ W_l,& \&cup;(p(x,y),1)\&NotElement;I_{\mathrm{ith}}\\ W_a& \mathrm{else}\end{array}\right.---\left(4\right)$

$W_a=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{Val}}_{\mathrm{ith}}(x,y)---\left(5\right)$

其中W_h为常量，代表较高的能量值；W_l为常量，代表较低的能量值；W_G为计算得出，Val_ith为投射到第i幅图像上的值，N代表被投射的图像的个数；

在信息传递过程中，我们定义当前的信息能量为：

V_P＝V_p-1+aminD+(1-a)V_PminD    (6)

其中α为常量，，minD为两点间最小距离，V_pminD代表取得两点间最小距离时点所代表的值。在能量最小化的过程中，我们定义能量为：

e_i(p_i)＝min[γD(p_i，q)+(1-γ)V(q)+e_i-1(q)]    (7)

其中e_i(p_t)为能量值；e_i-1(q)为上一邻居的能量值，γ为常量，D(p_i，q)为两点间距离，V(q)为当前点权值。

Images