CN103077550B - 一种非门控icus图像序列中血管的四维重建方法 - Google Patents

Abstract

一种非门控ICUS图像序列中血管的四维重建方法，所述方法利用与ICUS图像同步采集的一对近似正交的CAG图像序列精确定位超声导管回撤路径，结合从ICUS图像中提取的血管腔横截面信息，对靶血管段（包含可能存在的斑块）的形态结构进行四维重建，再现冠脉血管在心动周期中各时相的真实形态。本发明实现了冠状动脉靶血管段（包含可能存在的斑块）的形态结构的四维重建，全面和真实地反映靶血管段在心动周期不同时相的真实形态，为冠状动脉粥样硬化性病变的临床诊治和研究冠脉血管组织的生物力学特性等提供更为可靠的依据。所述方法不需要心电门控采集装置，也无需利用ECG信号，大大缩短了介入检查的时间，降低了对原始图像数据的要求。

Description

一种非门控ICUS图像序列中血管的四维重建方法

技术领域

本发明涉及一种从非门控的冠状动脉内超声(intracoronaryultrasound，ICUS)图像序列中四维(三维+时间)重建血管的方法。属于医学成像技术领域。

背景技术

冠状动脉内超声(intracoronaryultrasound，ICUS)是目前临床常用的诊治冠状动脉粥样硬化性病变的介入影像手段。由于二维图像的种种局限性，对ICUS图像数据予以三维重建，从而展示血管的整体观，对于辅助冠心病的诊断、制定介入手术规划、实现手术的精确定位等具有很好的临床意义和应用前景。

冠状动脉附着在心外膜表面上，随着心脏有节律地收缩和舒张，因而在心动周期的不同时相，血管的形态结构和空间位置都会发生很大变化。目前国内外对ICUS图像序列中冠脉血管三维重建的研究仅限于重建某一个时相(通常是舒张末期)的血管形态，其基本原理是采用心电(ECG)门控的图像采集方式，仅采集相同心脏时相(通常是R波，即舒张末期)的等间距ICUS图像序列，并将其按照采集顺序，沿从同步采集的X射线冠状动脉造影(coronaryangiography,CAG)图像中重建出的超声导管回撤路径排列，完成该时相血管的重建。其代表性的方法包括“A.Wahle,G.P.M.Prause,S.C.DeJong,etal.Geometricallycorrect3-Dreconstructionofintravascularultrasoundimagesbyfusionwithbiplaneangiography—methodsandvalidation.IEEETransactionsonMedicalImaging.1999,18(8):686-699”、“C.J.Slager,J.J.Wentzel,J.C.H.Schuurbiers,etal.True3-dimensionalreconstructionofcoronaryarteriesinpatientsbyfusionofangiographyandICUS(ANGUS)anditsquantitativevalidation.Circulation.2000,102(5):511-516”和“一种血管三维模型的重建方法，中国发明专利，专利号：ZL200810055037.9”。目前的重建方法仅能获得一个特定心脏时相(通常是舒张末期)的重建结果，而研究冠脉血管组织的生物力学特性、分析血管壁受到的由搏动血流引入的剪应力分布情况、建立血管弹性图等，都需要分析冠脉血管在心动周期中不同时刻的变形，因而限制了ICUS在此方面的应用。同时，心电门控的图像采集方式需要专门的ECG门控采集装置，而目前临床采用的多数血管内超声成像系统并不包含此功能。并且由于每个心动周期只采集一帧，与连续回撤超声导管的非门控采集方式相比，图像采集时间延长至少三倍，这不仅增加了病人的痛苦，也延长了医患暴露在X射线下的时间。

总之，目前还没有一种理想的、对非门控ICUS图像序列中的冠脉血管进行四维重建的方法，能够准确重建出冠脉血管在心动周期中不同时刻的解剖结构，获得对血管及其病变的全面了解。

发明内容

本发明的目的在于针对现有技术之弊端、提供一种对连续回撤超声导管采集的非门控ICUS图像序列中的靶血管段进行四维(即三维+时间)重建的方法，再现冠状动脉血管(包含可能存在的斑块)在心动周期中各时相的真实形态，从而为定量测量血流动力学参数、探讨粥样硬化斑块的分布和发展规律、研究冠脉血管组织的生物力学特性、制定介入治疗计划以及评价介入治疗效果等提供依据。

本发明所述问题是以下述技术方案实现的：

一种非门控ICUS图像序列中血管的四维重建方法，所述方法利用与ICUS图像同步采集的一对近似正交的CAG图像序列精确定位超声导管回撤路径，结合从ICUS图像中提取的血管腔横截面信息，对靶血管段(包含可能存在的斑块)的形态结构进行四维重建，再现冠脉血管在心动周期中各时相的真实形态，具体步骤如下：

a、采集靶血管段的ICUS图像序列和一对近似正交的CAG图像序列，其中ICUS图像采用连续回撤超声导管的非心电门控采集方式；

b、从CAG图像序列中三维重建出各时相的超声导管回撤路径；

c、从各帧ICUS图像中提取出血管壁的内外膜轮廓；

d、补偿ICUS图像序列中存在的平面外运动伪像，对图像序列进行分组，得到在各心脏时相处采集的子序列；

通过分析ICUS图像序列灰度特征的周期性变化，为各帧图像找到其在相邻心动周期的相同时相采集的对应帧，将整个图像序列划分成在心动周期不同时相采集的子序列，以抑制ICUS纵向视图中的锯齿效应、补偿平面外运动伪像；

具体步骤如下：

①计算两帧图像之间灰度特征的差异值：

对于由n帧组成的ICUS图像序列{I₁,I₂,…,I_n}，进行逐帧比较，计算两帧图像之间灰度特征的差异值d_i,j：

$d_{i,j}=1-\frac{\underset{x=0}{\overset{N-1}{\Σ}}\underset{y=0}{\overset{M-1}{\Σ}}|I_i(x,y)-{\mathrm{\μ}}_i|\·|I_j(x,y)-{\mathrm{\μ}}_j|}{\sqrt{\[\underset{x=0}{\overset{N-1}{\Σ}}\underset{y=0}{\overset{M-1}{\Σ}}{\[I_i(x,y)-{\mathrm{\μ}}_i\]}^2\]{\mathrm{\μ}}_i\[\underset{x=0}{\overset{N-1}{\Σ}}\underset{y=0}{\overset{M-1}{\Σ}}{\[I_j(x,y)-{\mathrm{\μ}}_j\]}^2\]}},$

式中，i,j＝1,2,…,n；I_i和I_j分别表示图像序列中的第i帧和第j帧，其图像大小均为N×M像素，平均灰度值分别为μ_i和μ_j；I_i(x,y)和I_j(x,y)分别表示第i帧和第j帧在像素(x,y)处的灰度值；

②确定首帧I₁在第二个心动周期的对应帧：

分析图像序列的第一帧I₁和其它帧的差异值{d_1,1,d_1,2,…,d_1,n}，将该函数除d_1,1＝0之外的第一个局部极小值所对应的i值作为第二个心动周期中与I₁在相同时相采集的帧的序号，记为F；

③估计平均心率的近似值和心动周期长度：

计算第m帧和第m+k帧(m＝1,2,…,n-k)之间的平均差异值 $\stackrel{\‾}{D}\left(k\right)=\frac{1}{n-k}\underset{m=1}{\overset{n-k}{\Σ}}{d}_{m,m+k},$ 将函数 $\begin{array}{cc}\stackrel{\‾}{D}\left(k\right)& (k=\mathrm{0,1,2},...,n-1)\end{array}$ 的幅度谱曲线峰值所对应的频率作为平均心率的近似值H_R(单位：次/分钟)，心动周期长度的近似值为C_C＝(60×F_R)/H_R(单位：帧)，其中F_R为ICUS图像的采集速率(单位：帧/秒)；

④为后续帧找到其在相邻心动周期中的对应帧：

在矩阵D＝[d_i,j](i,j＝1,2,…,n)中搜索一条具有最小累计差异值的路径，该路径连接互为对应帧的两帧图像的差异值，该路径的起点为d_1,F(即首帧与第F帧之间的差异值)，搜索范围Δ(即相邻心动周期对应帧序号之差的搜索范围)设定为Δ∈[C_C/2,2C_C]；

e、对于各时相的子序列，分别补偿其平面内运动伪像，确定各帧图像的空间方向：

对于由步骤d得到的各时相的ICUS子序列，在以导管中心为坐标原点的坐标系中，分别计算各帧图像中血管壁内膜轮廓的几何中心，作为对其重心的近似，并计算相邻帧之间内膜轮廓重心的位移(Δx_p,Δy_p)和旋转角Δα_p(p＝2,3,…,n)，采用中心频率为平均心率的近似值H_R、通带区间为[H_R-H_R/2,H_R+H_R/2]的带通滤波器分别对序列{Δx_p}、{Δy_p}和{Δα_p}进行滤波，从中分离出运动分量{Δx_p,d}、{Δy_p,d}、{Δα_p,d}，

然后，对各帧图像中的平面内运动进行补偿，具体方法是：对于图像序列中的第p帧(p＝2,3,…,n)，将其血管壁内外膜轮廓上各像素点的坐标(基于以导管中心为坐标原点的坐标系)先反向平移再反向旋转即得到补偿平面内运动伪像之后的血管壁轮廓；

f、对于补偿了平面内运动伪像之后的各时相的子序列，按照采集顺序，沿对应时相的三维导管路径，确定各帧ICUS图像的轴向位置：

对于由步骤d得到的各时相的ICUS子序列，根据采集图像时记录的相邻帧之间的切面间距，确定在对应时相的三维导管路径上各帧ICUS图像的采集点，然后，计算导管路径上各采集点处的单位切矢量，使各帧ICUS图像平面垂直于其采集点处的单位切矢，并且图像中心与采集点重合，将其等间隔地排列于对应时相的导管路径上；

g、对于沿导管路径准确排列的各时相的子序列，采用曲面拟合技术，得到光滑连续的血管腔内外表面。

上述非门控ICUS图像序列中血管的四维重建方法，从CAG图像序列中三维重建出各时相的超声导管回撤路径的具体步骤是：

a、对原始CAG图像进行初步的滤波、去噪、畸变校正等处理，增强视觉效果；

b、建立X射线血管造影系统在两个角度的投影模型，推导两幅不同角度的造影图像之间的几何变换关系：

设点s₁和s₂为两个角度的X射线源焦点的位置，分别以s₁和s₂为原点，建立空间坐标系s₁x₁y₁z₁和s₂x₂y₂z₂；坐标系U₁V₁O₁和U₂V₂O₂为成像平面坐标系；D₁和D₂分别是s₁和s₂到各自成像平面的垂直距离，空间血管上的点P在成像平面A和B上的投影点分别为p₁(u₁,v₁)和p₂(u₂,v₂)。根据透视投影成像理论，空间点的三维坐标与其二维投影坐标之间的关系为：

[x₁y₁z₁]^T＝z₁·[ξ₁η₁1]^T

[x₂y₂z₂]^T＝z₂·[ξ₂η₂1]^T

其中， ${\mathrm{\ξ}}_1=\frac{u_1}{D_1}=\frac{x_1}{z_1},{\mathrm{\η}}_1=\frac{v_1}{D_1}=\frac{y_1}{z_1},{\mathrm{\ξ}}_2=\frac{u_2}{D_2}=\frac{x_2}{z_2},{\mathrm{\η}}_2=\frac{v_2}{D_2}=\frac{y_2}{z_2}$

s₁x₁y₁z₁和s₂x₂y₂z₂之间的几何变换关系为：

[x₂y₂z₂]^T＝R([x₁y₁z₁]^T-t)

其中，R为3×3的旋转矩阵：R＝R_Y(α₂)·R_X(β₂)·R_X(-β₁)·R_Y(α₁)

t为平移向量：t＝R_Y(-α₁)·R_X(β₁)·[00L₁-L₂]^T

L₁和L₂为X射线源s₁和s₂到旋转中心的距离；(α₁,β₁)和(α₂,β₂)分别是A和B的造影角度，(x₁,y₁,z₁)和(x₂,y₂,z₂)分别表示空间点P在坐标系s₁x₁y₁z₁和s₂x₂y₂z₂中的坐标；R_Y(α₁)、R_Y(-α₁)和R_Y(α₂)分别表示绕Y轴旋转α₁、-α₁、α₂角度的旋转矩阵；R_X(β₁)、R_X(-β₁)和R_X(β₂)分别表示绕Y轴旋转β₁、-β₁、β₂角度的旋转矩阵。

c、在序列第一时刻的图像中，通过人工取点获得超声导管在左右投影中的近似中心线，用折线表示，并利用外极约束得到两个角度间对应点的匹配；

d、对手动选取的采样点进行三维重建；

e、连接各三维点，所得折线作为snake模型的初始位置，通过使预先设定的能量函数最小，snake模型在空间中变形，最终得到具有最小能量的最优位置，就是第一时刻的三维超声导管路径；

f、后续时刻超声导管路径的三维重建：

对于序列中的后续时刻，以前一时刻的三维超声导管路径作为当前时刻snake模型的初始位置，通过snake模型的变形，完成整个图像序列中各时刻超声导管路径的三维重建。

上述非门控ICUS图像序列中血管的四维重建方法，从各帧ICUS图像中提取出血管壁的内外膜轮廓的具体步骤是：

a、对原始ICUS图像进行预处理，包括滤波去噪和去除环晕伪像：

首先采用中值滤波和高斯平滑两种通用预处理方法，减少ICUS图像中的椒盐噪声和随机噪声，然后对各帧ICUS图像进行极坐标变换，得到其极坐标视图，再按照下式去除极坐标视图中的环晕伪像：

$I^{\′}(r,\mathrm{\θ})=\left\{\begin{array}{cc}{I}_{cathetet},& r\≤\mathrm{\ϵ}\×\mathrm{Im}ageWidth/2\\ I(r,\mathrm{\θ}),& else\end{array}\right.$

其中，I(r,θ)和I'(r,θ)分别为原极坐标视图和去除环晕伪像后的极坐标视图中像素点(r,θ)处的灰度值；I_catheter是原ICUS图像中导管区域的像素灰度值；r为像素点的极径；ImageWidth为以像素为单位的图像宽度；ε为权重参数；最后经过极坐标逆变换，即可得到直角坐标系下去除环晕伪像后的ICUS图像；

b、分别取ICUS图像序列的沿血管长轴方向的四个纵向截面视图，即垂直切面A、水平切面B、左对角线切面C和右对角线切面D；

c、从ICUS纵向视图中提取出血管内腔和中-外膜边界：

从纵向视图的中轴线开始，分别向左和向右逐行遍历四个纵向视图中的各像素，用I(i,j)表示坐标为(i,j)的当前像素的灰度值，若I(i,j+1)-I(i,j)≥η，η为阈值，则为目标边界点，否则不是；将每行左右两部分的像素中，第一个符合上述条件的像素记为内腔边界点，第二个记为中-外膜边界点；

d、将四个纵向视图中的边界曲线映射到各帧ICUS图像中，纵向视图的轮廓线对应到横向视图中为轮廓点，依次连接各轮廓点，得到各帧ICUS图像中的血管内腔和中-外膜初始轮廓；

e、初始轮廓通过snake模型变形获得各帧ICUS图像中的血管内腔边界和中-外膜边界：

将各帧ICUS图像中血管壁的初始轮廓作为snake模型的初始形状，将初始轮廓离散成由N个点组成的有序点集，则能量函数的离散表达式为：

$E=\underset{i=0}{\overset{N-1}{\Σ}}\[E_{\mathrm{int}}\left(i\right)+E_{ext}\left(i\right)\],$

$E_{\mathrm{int}}\left(i\right)=\mathrm{\α}\frac{\left|\stackrel{\‾}{d}-\left|c_i-c_{i-1}\right|\right|}{\max d}+\mathrm{\β}\frac{\left|c_{i-1}-2c_i+c_{i+1}\right|}{2\max d},$

其中，E_int是内部能量；E_ext是外部能量；c_i(x_i,y_i)(i＝1,2,…,N-2)是第i个snake点；和maxd分别是相邻snake点之间的平均距离和最大距离；α和β是权重参数，取值区间是[0,1]；I(x_i,y_i)和分别是像素(x_i,y_i)的灰度和灰度梯度值；γ,λ∈[0,1]是权重参数。

通过使能量函数E最小，snake模型从初始形状开始不断变形，最终停留在能量函数取得全局最小值的最优位置，即为目标轮廓。

上述非门控ICUS图像序列中血管的四维重建方法，所述冠状动脉造影图像序列的两个采集角度之间夹角的取值范围为60°至90°。

上述非门控ICUS图像序列中血管的四维重建方法，在第一时刻图像的感兴趣血管段中选取的三维重建采样点包括起点、终点和3～6个中间点。

上述非门控ICUS图像序列中血管的四维重建方法，所述权重参数ε的取值范围为[0.1,0.35]，阈值η的取值区间为[10,20]。

上述非门控ICUS图像序列中血管的四维重建方法，能量函数的全局最优化采用Williams贪婪算法完成。

本发明利用连续回撤超声导管采集的ICUS图像序列及与ICUS图像同步采集的一对近似正交的CAG图像序列对靶血管段(包含可能存在的斑块)的形态结构进行四维重建，全面地反映靶血管段在心动周期不同时相的真实形态，为冠状动脉粥样硬化性病变的临床诊治和研究冠脉血管组织的生物力学特性等提供更为可靠的依据。所述方法不需要心电门控采集装置，也无需利用ECG信号，大大缩短了介入检查的时间，降低了对原始图像数据的要求。

附图说明

下面结合附图对本发明作进一步说明。

图1是本发明方法的流程图；

图2是ICUS图像序列中存在的平面外运动示意图；

图3是ICUS图像序列中存在的平面内运动示意图；

图4是从X射线血管造影图像序列中三维重建各时相的管腔轴线流程图。

图5是X射线血管造影系统在两个角度的成像示意图；

图6是X射线血管造影成像系统的造影角度示意图；

图7是对ICUS图像序列进行三维分割的流程图；

图8是一帧ICUS图像的环晕伪像去除过程和结果。

图9是获取ICUS图像序列纵向视图示意图。

图中或说明书中所用符号清单为：n、ICUS图像序列的总帧数；{I₁,I₂,…,I_n}、由n帧组成的ICUS图像序列；I_i和I_j、ICUS图像序列中的第i帧和第j帧；d_i,j、I_i和I_j之间灰度特征的差异值；μ_i、I_i的平均灰度值；μ_j、I_j的平均灰度值；{d_1,1,d_1,2,…,d_1,n}、首帧I₁和其它帧的差异值；F、第二个心动周期中与I₁在相同时相采集的帧的序号；第m帧和第m+k帧之间的平均差异值；H_R、平均心率的近似值(单位：次/分钟)；C_C、心动周期长度的近似值(单位：帧)；F_R、ICUS图像的采集速率(单位：帧/秒)；D、由{d_i,j}(i,j＝1,2,…,n)组成的矩阵；d_1,F、首帧与第F帧之间的差异值；Δ、相邻心动周期对应帧序号之差的搜素范围；(Δx_p,Δy_p)和Δα_p、第p-1帧和第p帧之间血管壁内膜轮廓重心的位移和旋转角；{Δx_p,d}、{Δy_p,d}、{Δα_p,d}、序列{Δx_p}、{Δy_p}和{Δα_p}(p＝2,3,…,n)的心脏运动分量；{Δx_p,g}、{Δy_p,g}、{Δα_p,g}、序列{Δx_p}、{Δy_p}和{Δα_p}(p＝2,3,…,n)的血管几何分量；图中各符号为：I(r,θ)、I'(r,θ)、原极坐标视图和去除环晕伪像后的极坐标视图中像素点(r,θ)处的灰度值；I_catheter、原ICUS图像中导管区域的像素灰度值；r、像素点的极径；ImageWidth、以像素为单位的图像宽度；ε、权重参数；(i,j)、当前像素的坐标；I(i,j)、坐标为(i,j)的当前像素的灰度值；η、阈值；maxd、相邻snake点之间的平均距离和最大距离；α、β、权重参数；E、能量函数；E_ext、外部能量；I(x_i,y_i)、像素(x_i,y_i)的灰度和灰度梯度值；γ、λ、权重参数；A、B、成像平面；s₁、s₂、两次造影过程中X射线源焦点的位置；s₁x₁y₁z₁、以s₁为原点的空间坐标系；s₂x₂y₂z₂、以s₂为原点的空间坐标系；OXYZ、以旋转中心为原点的成像系统坐标系；U₁V₁O₁、成像平面A上的直角坐标系；U₂V₂O₂、成像平面B上的直角坐标系；D₁、s₁到成像平面A的垂直距离；D₂、s₂到成像平面B的垂直距离；P、空间血管上的点；p₁、P点在成像平面A上的投影；p₂、P点在成像平面B上的投影；u₁、p₁在坐标系U₁V₁O₁内的横坐标；v₁、p₁在坐标系U₁V₁O₁内的纵坐标；u₂、p₂在坐标系U₂V₂O₂内的横坐标；v₂、p₂在坐标系U₂V₂O₂内的纵坐标；K₁、K₂、外极线；(α₁,β₁)、成像平面A的造影角度；(α₂,β₂)、成像平面B的造影角度；L₁、X射线源s₁到旋转中心的距离；L₂、X射线源s₂到旋转中心的距离。

具体实施方式

下面结合附图详细说明本发明的处理步骤：

(1)图像采集：

ICUS和CAG成像是同时进行的，ICUS图像序列用血管内超声成像仪，采用连续回撤超声导管的非心电门控方式采集，可缩短介入检查的时间。采用C型臂X射线血管造影机采集一对近似正交的CAG图像序列，要求两个造影角度之间的夹角在60°至120°之间，记录成像系统参数(造影角度，X射线源到成像平面的距离)。

(2)超声导管回撤路径的四维重建：

超声导管回撤路径的四维重建步骤包括：

(1)图像预处理：

对原始CAG图像进行必要的预处理，主要包括畸变校正、均衡对比度、去除噪声和图像增强等，提高图像的视觉效果，为后续处理奠定基础。

(2)建立X射线血管造影系统在两个近似正交角度的成像模型：

如附图5所示，点s₁和s₂表示两个角度X射线源焦点的位置。分别以s₁和s₂为原点，建立空间坐标系s₁x₁y₁z₁和s₂x₂y₂z₂；坐标系U₁V₁O₁和U₂V₂O₂为成像平面坐标系；D₁和D₂分别是s₁和s₂到各自成像平面的垂直距离，随成像面的移动而改变；空间血管上的点P在成像平面A和B上的投影点分别为p₁(u₁,v₁)和p₂(u₂,v₂)。

从坐标系s₁x₁y₁z₁和s₂x₂y₂z₂的变换运动是三维空间中的刚体运动，根据刚体运动理论，推导出坐标系s₁x₁y₁z₁和s₂x₂y₂z₂之间的几何变换关系：

[x₂y₂z₂]^T＝R([x₁y₁z₁]^T-t)

(1)

其中R为3×3的旋转矩阵：

R＝R_Y(α₂)·R_X(β₂)·R_X(-β₁)·R_Y(α₁)

(2)

t为平移向量：

t＝R_Y(-α₁)·R_X(β₁)·[00L₁-L₂]^T

(3)

式中，L₁和L₂为X射线源s₁和s₂到旋转中心的距离；(α₁,β₁)和(α₂,β₂)分别是图像A和B的造影角度(如附图6所示)。距离和角度值均可从造影机上得到。R_Y(α₁)、R_Y(-α₁)和R_Y(α₂)分别表示绕Y轴旋转α₁、-α₁、α₂角度的旋转矩阵；R_X(β₁)、R_X(-β₁)和R_X(β₂)分别表示绕Y轴旋转β₁、-β₁、β₂角度的旋转矩阵。

根据透视投影成像的几何关系可知：

[x₁y₁z₁]^T＝z₁·[ξ₁η₁1]^T

[x₂y₂z₂]^T＝z₂·[ξ₂η₂1]^T

(4)

其中

${\mathrm{\ξ}}_1=\frac{u_1}{D_1}=\frac{x_1}{z_1},{\mathrm{\η}}_1=\frac{v_1}{D_1}=\frac{y_1}{z_1},{\mathrm{\ξ}}_2=\frac{u_2}{D_2}=\frac{x_2}{z_2},{\mathrm{\η}}_2=\frac{v_2}{D_2}=\frac{y_2}{z_2}---\left(5\right)$

(x₁,y₁,z₁)和(x₂,y₂,z₂)分别表示空间点P在坐标系s₁x₁y₁z₁和s₂x₂y₂z₂中的坐标。因此根据式(1)、(4)和(5)，由点p₁和p₂的坐标可反算出点P的三维坐标。

(3)第一时刻超声导管路径的三维重建：

本发明采用手动取点的方法，首先从一个角度的图像中手动选取超声导管投影上的若干采样点(一般选取超声导管投影的起点、终点和3～4个中间点即可)。然后采用外极约束得到各点在另一个角度图像上的对应点。如附图5所示，空间点P与X射线源焦点s₁和s₂构成外极平面，该平面与图像平面A和B相交形成左右两条外极线：K₁和K₂。根据外极约束原理，点P在图像A上的投影p₁在图像B中的对应点p₂一定位于K₂上；P在图像B上的投影p₂在图像A中的对应点p₁一定位于K₁上。由于几何变换矩阵和其它计算的误差，匹配点可能不会准确地出现在对应的外极线上，而在外极线的邻域内。本发明方法在该邻域内搜索和外极线距离最近的点作为匹配点。由这几组对应点分别求出它们的三维坐标。用直线段连接这些三维点，所得折线作为snake模型的初始位置。

snake模型(M.Kass,A.Witkin,D.Terzopoulos.Snakes:activecontourmodels.InternationalJournalofComputerVision.1987,vol.1,no.4,pp.321-331)是一条可变形的参数曲线：c(s)＝(x(s),y(s),z(s)),s∈[0,1]。将该曲线离散化为N个点c_i(x_i,y_i,z_i)(i＝1,2,…,N)，模型能量函数的离散表达式为

$E=\underset{i=1}{\overset{N}{\Σ}}\[E_{\mathrm{int}}\left(i\right)+E_{ext}\left(i\right)\]---\left(6\right),$

其中，E_int是内部能量：

$E_{\mathrm{int}}\left(i\right)=\mathrm{\α}(\stackrel{\‾}{d}-|c_i-c_{i-1}\left|\right)+\mathrm{\β}\left|c_{i-1}-2c_i+c_{i+1}\right|---\left(7\right)$

保证曲线变形过程中的连续性和光滑性。式中第一项约束使点平均分布，既能满足连续性的要求，又不会产生曲线收缩的现象。在每次迭代结束时，点间的平均距离的值被更新；第二项是曲率，使曲线保持平滑。α和β是权值。

外部能量E_ext是保证模型收敛的外部力，吸引表示超声导管的空间曲线移动到这样的位置：3D曲线在左右两个图像平面上的投影位于图像上灰度最小、且灰度梯度最小的区域，即超声导管的投影。根据式(1)和(5)，空间点P在图像A和B上的投影坐标(u₁,v₁)和(u₂,v₂)都可用该点的三维坐标c＝(x₁,y₁,z₁)来表示：

[u₁v₁]^T＝T_L(c)，[u₂v₂]^T＝T_R(c)(8)

其中T_L和T_R都是以c为自变量的函数。E_ext包含两部分，分别对应于两个角度的二维投影：

$E_{ext}=\mathrm{\γ}\[I_L\left(T_L\right(c_i\left)\right)+I_R\left(T_R\right(c_i\left)\right)\]+\mathrm{\λ}\[|\▿I_L\left(T_L\right(c_i\left)\right)|+|\▿I_R\left(T_R\right(c_i\left)\right)|\]---\left(9\right)$

式中I_L(T_L(c_i))＝I_L(u_1i,v_1i)和I_R(T_R(c_i))＝I_R(u_2i,v_2i)分别是左右图像点的灰度值； $\▿I_L\left(T_L\right(c_i\left)\right)=\▿I_L(u_{1i},v_{1i})$ 和分别是左右图像点的灰度梯度。由于造影图像中，血管的灰度值比背景小，所以这里权重参数γ取正值。而且感兴趣的是位于血管腔内部的超声导管，不是管腔边缘，所以λ也取正值。

(4)后续时刻超声导管路径的三维重建：

得到第一时刻的三维导管路径之后，对于其后的序列图像，根据血管形状不突变的特点，把前一时刻的三维导管路径作为当前时刻snake初始位置，实现对整个序列中各时刻导管路径的三维重建。

在这一步骤中，本发明方法在snake能量函数中嵌入相似度匹配，通过对相邻时刻之间的目标进行配准，实现对超声导管的跟踪。假设在足够短的时间间隔内，运动前后血管图像的灰度不发生变化，那么得到第一时刻的超声导管之后，在开始对后续图像进行运动跟踪时，外部能量函数采用如下的形式：

$\begin{array}{c}{E}_{ext}=\mathrm{\γ}\[I_L^t\left(T_L\right(c_i^t\left)\right)+I_R^t\left(T_R\right(c_i^t\left)\right)\]+\mathrm{\λ}\left(\right|\▿I_L^t\left(T_L\right(c_i^t\left)\right)|+|\▿I_R^t\left(T_R\right(c_i^t\left)\right)\left|\right)\\ +\mathrm{\η}\[{\Σ}_{\mathrm{\Δ}w}|I_L^t\left(T_L\right(c_i^t)+\mathrm{\Δ}w)-I_L^{t-1}\left(T_L\right(c_i^{t-1})+\mathrm{\Δ}w)|\\ +{\Σ}_{\mathrm{\Δ}w}|I_R^t\left(T_L\right(c_i^t)+\mathrm{\Δ}w)-I_R^{t-1}\left(T_L\right(c_i^{t-1})+\mathrm{\Δ}w)|\]\end{array}$

$\left(10\right)$

式中和分别表示时刻t左右图像点的灰度值；和分别表示时刻t-1左右图像点的灰度值；Δw是邻域窗口的大小。(10)式中的前两项与(9)式相同，表示吸引曲线向着图像中灰度最小、灰度梯度最小的区域移动的图像力。第三项是使得运动前后(时刻t-1和t)血管上对应点的Δw×Δw邻域内的灰度变化最小的外部约束力。

对于求解snake能量函数最小化的问题，在计算机视觉发展的最初阶段，一般都是采用变分法，迭代求解E-L偏微分方程。所述方法不仅计算复杂，运算时间长，而且很容易收敛于局部最小值。采用动态规划(“Usingdynamicprogrammingforsolvingvariationalproblemsinvision,”IEEETransactionsonpatternanalysisandmachineintelligence.vol.12,no.9,pp.855-867,1990)可避免上述问题，由于迭代向着总能量减少的方向进行，当总能量不再变化时，则迭代终止，因此保证了结果的收敛性。但是所述方法的缺点是计算速度慢，需要较大的存贮空间。针对上述问题，本发明采用Williams(“Afastalgorithmforactivecontoursandcurvatureestimation,”ComputerVision,GraphicsandImageProcessing,vol.55,no.1,pp.14-26,1992)提出的贪婪算法完成对最优snake的搜索，该算法具有动态规划方法的所有优点，并且计算速度快，所需存贮空间小。

对于重建出的各时刻的三维超声导管路径，根据式(8)即可算出其在左右两个成像平面上的二维投影。因此本发明方法将二维提取、三维重建和运动跟踪集成到一个框架中完成，提高了运算精度和速度。

(3)提取各帧ICUS图像中的血管壁内外膜轮廓：

如附图7所示，本发明方法的步骤包括：

(1)对原始ICUS图像进行预处理，包括滤波去噪和去除环晕伪像：

首先，采用中值滤波和高斯平滑两种通用预处理方法，减少ICUS图像中的椒盐噪声和随机噪声。

然后采用极坐标变换法去除环晕伪像，具体方法如下：

首先，对各帧ICUS图像进行极坐标变换，得到其极坐标视图，如附图8所示。可见，环晕伪像固定的位于极坐标视图的上部。然后，按照下式去除极坐标视图中的环晕伪像：

$I^{\′}(r,\mathrm{\θ})=\{\begin{array}{cc}{I}_{catheter},& r\≤\mathrm{\ϵ}\×\mathrm{Im}ageWidth/2\\ I(r,\mathrm{\θ}),& else\end{array}---\left(11\right)$

其中，I(r,θ)和I'(r,θ)分别为原极坐标视图和去除环晕伪像后的极坐标视图中像素点(r,θ)处的灰度值；I_catheter是原ICUS图像中导管区域的像素灰度值；r为像素点的极径；ImageWidth为以像素为单位的图像宽度；ε为权重参数，经实验确定其取值范围为[0.1,0.35]。在极坐标视图中消除环晕伪像后，再经过极坐标逆变换，即可得到直角坐标系下去除环晕伪像后的ICUS图像。

(2)获取ICUS图像序列的四个纵向视图：

如附图9所示，分别取ICUS图像序列的沿血管长轴方向的四个纵向截面视图，即垂直切面A、水平切面B、左对角线切面C和右对角线切面D。

(3)从纵向视图中提取出血管内腔和中-外膜边界：

分别遍历四个纵向视图中的各像素，判断识别血管壁边界点。由于目标边界分布于纵向视图的中部，所以从纵向视图的中轴线开始，分别向左和向右进行逐行遍历。对于坐标为(i,j)的当前像素，其灰度值为I(i,j)，若I(i,j+1)-I(i,j)≥η，η为阈值，则为目标边界点，否则不是。将每行左右两部分的像素中，第一个符合上述条件的像素记为内腔边界点，第二个记为中-外膜边界点。经实验确定η的取值区间为[10,20]。

(4)将四个纵向视图中的边界曲线映射到各帧ICUS图像中，得到各横向视图中血管壁的初始轮廓：

如附图9所示，纵向视图的轮廓线对应到横向视图中为轮廓点，四条纵向视图轮廓线映射到横向视图中为四个轮廓点。依次连接各轮廓点，即可获得各帧ICUS图像中的血管内腔和中-外膜初始轮廓。

(5)通过使预先设定的能量函数最小，初始轮廓在内外力的共同作用下不断变形，最终得到各帧ICUS图像中的血管内腔边界和中-外膜边界：

将各帧ICUS图像中血管壁的初始轮廓作为snake模型的初始形状，通过使预先设定的能量函数最小，snake模型不断变形，当能量函数取得全局最小值时，snake模型即停留在目标轮廓处，从而完成对各帧图像中血管内腔和中-外膜边界的并行提取。

能量函数的构成如下：

将初始轮廓离散成由N个点组成的有序点集，则能量函数的离散

表达式为： $E=\underset{i=0}{\overset{N-1}{\Σ}}\[E_{\mathrm{int}}\left(i\right)+E_{ext}\left(i\right)\]$

$\left(12\right)$

其中内部能量E_int的归一化表达式为：

$E_{\mathrm{int}}\left(i\right)=\mathrm{\α}\frac{\left|\stackrel{\‾}{d}-\left|c_i-c_{i-1}\right|\right|}{\max d}+\mathrm{\β}\frac{\left|c_{i-1}-2c_i+c_{i+1}\right|}{2\max d}---\left(13\right)$

其中c_i(x_i,y_i)(i＝1,2,…,N-2)是第i个snake点；和maxd分别是相邻snake点之间的平均距离和最大距离，在每次迭代结束时，它们的值都被更新；α和β是权重参数。(13)式中的第一项保证变形过程中snake的连续性，使snake点均匀分布，不致产生收缩的现象；第二项是曲线二阶导数的离散形式，保证变形过程中snake的光滑性。由于进行了归一化，故(13)式中两项的取值范围都在[0,1]区间内，权重α和β的取值区间也是[0,1]。

外部能量E_ext定义为：

$E_{ext}\left(i\right)=\mathrm{\γ}\frac{I(x_i,y_i)}{255}-\mathrm{\λ}\frac{|\▿I(x_i,y_i)|}{255\sqrt{2}}---\left(14\right)$

其中I(x_i,y_i)和分别是像素(x_i,y_i)的灰度和灰度梯度值。由于8位灰度图像的灰度和灰度梯度的取值范围分别为[0,255]和故本发明对I(x_i,y_i)和分别进行了归一化，使其取值范围为[0,1]区间。γ,λ∈[0,1]是权重参数。

通过使能量函数E最小，snake模型从初始形状开始不断变形，最终停留在能量函数取得全局最小值的最优位置，即为目标轮廓。能量函数的全局最优化采用Williams贪婪算法完成，不仅计算速度快，所需存贮空间小，而且非常适合于编程实现。

对临床采集的ICUS图像序列的实验证明，在获得与现有逐帧处理的二维串行分割方法相类似的分割精度的情况下，采用本发明的三维分割方法可大大提高处理效率，缩短处理时间。

(4)确定各帧ICUS图像的空间方向：

在采集ICUS图像序列的过程中，伴随着周期性心脏运动和冠脉管腔内血流的搏动，超声导管在血管腔内的运动导致ICUS图像序列中存在两种运动伪像：平面内运动和平面外运动。前者即导管和血管腔之间的相对平移和旋转，导致切片图像序列相邻帧之间的错位；后者即导管沿管腔轴线的前向或后向的纵向位移，导致纵向视图中的血管壁边缘呈现锯齿形(如附图2所示)。本发明分别对这两种运动伪像进行补偿，从而确定各帧ICUS图像的空间方向。

(4.1)补偿平面外运动伪像：

本发明通过分析ICUS图像序列灰度特征的周期性变化，为各帧图像找到其在相邻心动周期的相同时相采集的对应帧，将整个图像序列划分成在心动周期不同时相采集的子序列，从而抑制ICUS纵向视图中的锯齿效应，补偿平面外运动伪像。具体步骤如下：

(4.1.1)计算两帧图像之间灰度特征的差异值：

对于由n帧组成的ICUS图像序列{I₁,I₂,…,I_n}，进行逐帧比较，计算两帧图像之间灰度特征的差异值d_i,j：

$d_{i,j}=1-\frac{\underset{x=0}{\overset{N-1}{\Σ}}\underset{y=0}{\overset{M-1}{\Σ}}|I_i(x,y)-{\mathrm{\μ}}_i|\·|I_j(x,y)-{\mathrm{\μ}}_j|}{\sqrt{\[\underset{x=0}{\overset{N-1}{\Σ}}\underset{y=0}{\overset{M-1}{\Σ}}{\[I_i(x,y)-{\mathrm{\μ}}_i\]}^2\]{\mathrm{\μ}}_i\[\underset{x=0}{\overset{N-1}{\Σ}}\underset{y=0}{\overset{M-1}{\Σ}}{\[I_j(x,y)-{\mathrm{\μ}}_j\]}^2\]}}---\left(15\right)$

式中，i,j＝1,2,…,n；I_i和I_j分别表示图像序列中的第i帧和第j帧，其图像大小均为N×M像素，平均灰度值分别为μ_i和μ_j；I_i(x,y)和I_j(x,y)分别表示第i帧和第j帧在像素(x,y)处的灰度值。

(4.1.2)确定首帧I₁在第二个心动周期的对应帧：

分析图像序列的第一帧I₁和其它帧的差异值{d_1,1,d_1,2,…,d_1,n}，由于心脏的周期性运动，使得函数d_1,i(i)(i＝1,2,…,n)具有近似周期性，且随着i的增大，该函数的各局部极小值呈递增趋势。本发明将该函数除d_1,1＝0之外的第一个局部极小值所对应的i值作为第二个心动周期中与I₁在相同时相采集的帧的序号，记为F。

(4.1.3)估计平均心率的近似值和心动周期长度：

计算第m帧和第m+k帧(m＝1,2,…,n-k)之间的平均差异值 $\stackrel{\‾}{D}\left(k\right)=\frac{1}{n-k}\underset{m=1}{\overset{n-k}{\Σ}}{d}_{m,m+k}.$ 本发明将函数 $\begin{array}{cc}\stackrel{\‾}{D}\left(k\right)& (k=\mathrm{0,1,2},...,n-1)\end{array}$ 的幅度谱曲线峰值所对应的频率作为平均心率的近似值H_R(单位：次/分钟)，心动周期长度的近似值为C_C＝(60×F_R)/H_R(单位：帧)，其中F_R为ICUS图像的采集速率(单位：帧/秒)。

(4.1.4)为后续帧找到其在相邻心动周期中的对应帧：

本发明在矩阵D＝[d_i,j](i,j＝1,2,…,n)中搜索一条具有最小累计差异值的路径，该路径连接互为对应帧的两帧图像的差异值，从而为后续帧{I₂,I₃,…,I_n}找到其在相邻心动周期的相同时相采集的对应帧。该路径的起点为d_1,F(即首帧与第F帧之间的差异值)，搜索范围Δ(即相邻心动周期对应帧序号之差的搜素范围)设定为Δ∈[C_C/2,2C_C]。

(4.2)补偿平面内运动伪像：

对于由步骤(4.1)得到的各时相的ICUS子序列，在以导管中心为坐标原点的坐标系中，分别计算各帧ICUS图像中血管壁内膜轮廓的几何中心，作为对其重心的近似，并计算相邻帧之间内膜轮廓重心的位移(Δx_p,Δy_p)和旋转角Δα_p(p＝2,3,…,n)，如附图3所示。由于相邻帧之间管腔横截面重心的位移和空间方向的改变主要由周期性心脏运动和血管腔的不规则几何形状两方面因素造成，因此序列{Δx_p}、{Δy_p}和{Δα_p}(p＝2,3,…,n)包含心脏运动分量{Δx_p,d}、{Δy_p,d}、{Δα_p,d}和血管几何分量{Δx_p,g}、{Δy_p,g}、{Δα_p,g}两部分。本发明采用中心频率为平均心率的近似值H_R、通带区间为[H_R-H_R/2,H_R+H_R/2]的带通滤波器分别对序列{Δx_p}、{Δy_p}和{Δα_p}进行滤波，从中分离出运动分量{Δx_p,d}、{Δy_p,d}、{Δα_p,d}。

然后，对各帧图像中的平面内运动进行补偿。具体方法是：对于图像序列中的第p帧(p＝2,3,…,n)，将其血管壁内外膜轮廓上各像素点的坐标(基于以导管中心为坐标原点的坐标系)先反向平移再反向旋转即得到补偿平面内运动伪像之后的血管壁轮廓。

(5)沿对应时相的三维导管路径排列各帧ICUS图像：

对于由步骤(4.2)得到的各时相的ICUS子序列，根据采集图像时记录的相邻帧之间的切面间距，确定在对应时相的三维导管路径上各帧ICUS图像的采集点。然后，计算导管路径上各采集点处的单位切矢量，使各帧ICUS图像平面垂直于其采集点处的单位切矢，并且图像中心与采集点重合，将其等间隔地排列于对应时相的导管路径上。

(6)拟合血管内外表面：

对于沿三维导管路径准确排列的各时相的ICUS子序列，分别采用曲面拟合技术，对各帧图像中的血管壁内外膜轮廓进行曲面拟合，得到光滑连续的血管腔内外表面，完成对各时相冠脉血管的三维重建。

Claims (7)

1.一种非门控ICUS图像序列中血管的四维重建方法，其特征是，所述方法利用与ICUS图像同步采集的一对近似正交的CAG图像序列精确定位超声导管回撤路径，结合从ICUS图像中提取的血管腔横截面信息，对包含可能存在斑块的靶血管段的形态结构进行四维重建，再现冠脉血管在心动周期中各时相的真实形态，具体步骤如下：

a、采集靶血管段的ICUS图像序列和一对近似正交的CAG图像序列，其中ICUS图像采用连续回撤超声导管的非心电门控采集方式；

b、从CAG图像序列中三维重建出各时相的超声导管回撤路径；

c、从各帧ICUS图像中提取出血管壁的内外膜轮廓；

d、补偿ICUS图像序列中存在的平面外运动伪像，对图像序列进行分组，得到在各心脏时相处采集的子序列；

分析ICUS图像序列灰度特征的周期性变化，为各帧图像找到其在相邻心动周期的相同时相采集的对应帧，将整个图像序列划分成在心动周期不同时相采集的子序列，抑制ICUS纵向视图中的锯齿效应，补偿平面外运动伪像；

具体步骤如下：

①计算两帧图像之间灰度特征的差异值：

对于由n帧组成的ICUS图像序列{I₁,I₂,…,I_n}，进行逐帧比较，计算两帧图像之间灰度特征的差异值d_i,j：

$d_{i,j}=1-\frac{\underset{x=0}{\overset{N-1}{\Σ}}\underset{y=0}{\overset{M-1}{\Σ}}|I_i(x,y)-{\mathrm{\μ}}_i|\·|I_j(x,y)-{\mathrm{\μ}}_j|}{\sqrt{\[\underset{x=0}{\overset{N-1}{\Σ}}\underset{y=0}{\overset{M-1}{\Σ}}{\[I_i(x,y)-{\mathrm{\μ}}_i\]}^2\]\[\underset{x=0}{\overset{N-1}{\Σ}}\underset{y=0}{\overset{M-1}{\Σ}}{\[I_j(x,y)-{\mathrm{\μ}}_j\]}^2\]}}$

式中，i,j＝1,2,…,n；I_i和I_j分别表示图像序列中的第i帧和第j帧，其图像大小均为N×M像素，平均灰度值分别为μ_i和μ_j；I_i(x,y)和I_j(x,y)分别表示第i帧和第j帧在像素(x,y)处的灰度值；

②确定首帧I₁在第二个心动周期的对应帧：

分析图像序列的第一帧I₁和其它帧的差异值{d_1,1,d_1,2,…,d_1,n}，将所述图像序列中除d_1,1＝0之外的第一个局部极小值所对应的i值作为第二个心动周期中与I₁在相同时相采集的帧的序号，记为F；

③估计平均心率的近似值和心动周期长度：

计算第m帧和第m+k帧之间的平均差异值其中，m＝1,2,…,n-k，将函数的幅度谱曲线峰值所对应的频率作为平均心率的近似值H_R，H_R的单位为：次/分钟，k＝0,1,2,…,n-1；心动周期长度的近似值为C_C＝(60×F_R)/H_R，C_C的单位为：帧，其中F_R为ICUS图像的采集速率，F_R的单位为：帧/秒；

④为后续帧找到其在相邻心动周期中的对应帧：

在矩阵D＝[d_i,j]中搜索一条具有最小累计差异值的路径，其中i和j均＝1,2,…,n，该路径连接互为对应帧的两帧图像的差异值，该路径的起点为d_1,F，即首帧与第F帧之间的差异值，搜索范围Δ，即相邻心动周期对应帧序号之差的搜索范围，设定为Δ∈[C_C/2,2C_C]；

e、对于各时相的子序列，分别补偿其平面内运动伪像，确定各帧图像的空间方向：

对于由步骤d得到的各时相的ICUS子序列，在以导管中心为坐标原点的坐标系中，分别计算各帧图像中血管壁内膜轮廓的几何中心，作为对其重心的近似，并计算相邻帧之间内膜轮廓重心的位移(Δx_p,Δy_p)和旋转角Δα_p，此处p＝2,3,…,n，采用中心频率为平均心率的近似值H_R、通带区间为[H_R-H_R/2,H_R+H_R/2]的带通滤波器分别对序列{Δx_p}、{Δy_p}和{Δα_p}进行滤波，从中分离出运动分量{Δx_p,d}、{Δy_p,d}、{Δα_p,d}；

然后，对各帧图像中的平面内运动进行补偿，具体方法是：对于图像序列中的第p帧，此处p＝2,3,…,n，将其血管壁内外膜轮廓上各像素点的坐标，所述坐标为基于以导管中心为坐标原点的坐标系，先反向平移再反向旋转即得到补偿平面内运动伪像之后的血管壁轮廓；

f、对于补偿了平面内运动伪像之后的各时相的子序列，按照采集顺序，沿对应时相的三维导管路径，确定各帧ICUS图像的轴向位置：

对于由步骤d得到的各时相的ICUS子序列，根据采集图像时记录的相邻帧之间的切面间距，确定在对应时相的三维导管路径上各帧ICUS图像的采集点，然后，计算导管路径上各采集点处的单位切矢量，使各帧ICUS图像平面垂直于其采集点处的单位切矢，并且图像中心与采集点重合，将其等间隔地排列于对应时相的导管路径上；

g、对于沿导管路径准确排列的各时相的子序列，采用曲面拟合技术，得到光滑连续的血管腔内外表面。

2.根据权利要求1所述的非门控ICUS图像序列中血管的四维重建方法，其特征是，从CAG图像序列中三维重建出各时相的超声导管回撤路径的具体步骤是：

a、对原始CAG图像进行初步的滤波、去噪、畸变校正的处理，增强视觉效果；

b、建立X射线血管造影系统在两个角度的投影模型，推导两幅不同角度的造影图像之间的几何变换关系：

设点s₁和s₂为两个角度的X射线源焦点的位置，分别以s₁和s₂为原点，建立空间坐标系s₁x₁y₁z₁和s₂x₂y₂z₂；坐标系U₁V₁O₁和U₂V₂O₂为成像平面坐标系；D₁和D₂分别是s₁和s₂到各自成像平面的垂直距离，空间血管上的点P在成像平面A和B上的投影点分别为p₁(u₁,v₁)和p₂(u₂,v₂)；根据透视投影成像理论，空间点的三维坐标与其二维投影坐标之间的关系为：

[x₁y₁z₁]^T＝z₁·[ξ₁η₁1]^T，

[x₂y₂z₂]^T＝z₂·[ξ₂η₂1]^T，

其中， ${\mathrm{\ξ}}_1=\frac{u_1}{D_1}=\frac{x_1}{z_1},$ ${\mathrm{\η}}_1=\frac{v_1}{D_1}=\frac{y_1}{z_1},$ ${\mathrm{\ξ}}_2=\frac{u_2}{D_2}=\frac{x_2}{z_2},$ ${\mathrm{\η}}_2=\frac{v_2}{D_2}=\frac{y_2}{z_2},$

s₁x₁y₁z₁和s₂x₂y₂z₂之间的几何变换关系为：

[x₂y₂z₂]^T＝R([x₁y₁z₁]^T-t)

其中，R为3×3的旋转矩阵：R＝R_Y(α₂)·R_X(β₂)·R_X(-β₁)·R_Y(α₁)，

t为平移向量：t＝R_Y(-α₁)·R_X(β₁)·[00L₁-L₂]^T，

L₁和L₂为X射线源s₁和s₂到旋转中心的距离；(α₁,β₁)和(α₂,β₂)分别是A和B的造影角度，(x₁,y₁,z₁)和(x₂,y₂,z₂)分别表示空间点P在坐标系s₁x₁y₁z₁和s₂x₂y₂z₂中的坐标；R_Y(α₁)、R_Y(-α₁)和R_Y(α₂)分别表示绕Y轴旋转α₁、-α₁、α₂角度的旋转矩阵；R_X(β₁)、R_X(-β₁)和R_X(β₂)分别表示绕Y轴旋转β₁、-β₁、β₂角度的旋转矩阵；

c、在序列中第一时刻的图像中，通过人工取点获得超声导管在左右投影中的近似中心线，用折线表示，并利用外极约束得到两个角度间对应点的匹配；

d、对手动选取的采样点进行三维重建；

e、连接各三维点，所得折线作为snake模型的初始位置，通过使预先设定的能量函数最小，snake模型在空间中变形，最终得到具有最小能量的最优位置，就是第一时刻的三维超声导管路径；

f、后续时刻超声导管路径的三维重建：

对于序列中的后续时刻，以前一时刻的三维超声导管路径作为当前时刻snake模型的初始位置，通过snake模型的变形，完成整个图像序列中各时刻超声导管路径的三维重建。

3.根据权利要求1或2所述的一种非门控ICUS图像序列中血管的四维重建方法，其特征是，从各帧ICUS图像中提取出血管壁的内外膜轮廓的具体步骤是：

a、对原始ICUS图像进行预处理，包括滤波去噪和去除环晕伪像：

首先采用中值滤波和高斯平滑两种通用预处理方法，减少ICUS图像中的椒盐噪声和随机噪声，然后对各帧ICUS图像进行极坐标变换，得到其极坐标视图，再按照下式去除极坐标视图中的环晕伪像：

$I^{\′}(r,\mathrm{\θ})=\left\{\begin{array}{cc}{I}_{catheter},& r\≤\mathrm{\ϵ}\×\mathrm{Im}ageWidth/2\\ I\left(r,\mathrm{\θ}\right),& else\end{array}\right.,$

其中，I(r,θ)和I'(r,θ)分别为原极坐标视图和去除环晕伪像后的极坐标视图中像素点(r,θ)处的灰度值；I_catheter是原ICUS图像中导管区域的像素灰度值；r为像素点的极径；ImageWidth为以像素为单位的图像宽度；ε为权重参数；最后经过极坐标逆变换，即可得到直角坐标系下去除环晕伪像后的ICUS图像；

b、分别取ICUS图像序列的沿血管长轴方向的四个纵向截面视图，即垂直切面A、水平切面B、左对角线切面C和右对角线切面D；

c、从ICUS纵向视图中提取出血管内腔和中-外膜边界：

从纵向视图的中轴线开始，分别向左和向右逐行遍历四个ICUS纵向视图中的各像素，用I(i,j)表示坐标为(i,j)的当前像素的灰度值，若I(i,j+1)-I(i,j)≥η，η为阈值，则为目标边界点，否则不是；将每行左右两部分的像素中，第一个符合上述条件的像素记为内腔边界点，第二个记为中-外膜边界点；

d、将四个纵向视图中的边界曲线映射到各帧ICUS图像中，纵向视图的轮廓线对应到横向视图中为轮廓点，依次连接各轮廓点，得到各帧ICUS图像中的血管内腔和中-外膜初始轮廓；

e、初始轮廓通过snake模型变形获得各帧ICUS图像中的血管内腔边界和中-外膜边界：

将各帧ICUS图像中血管壁的初始轮廓作为snake模型的初始形状，将初始轮廓离散成由N个点组成的有序点集，则能量函数的离散表达式为：

$E=\underset{i=0}{\overset{N-1}{\Σ}}\[E_{\mathrm{int}}\left(i\right)+E_{ext}\left(i\right)\],$

$E_{\mathrm{int}}\left(i\right)=\mathrm{\α}\frac{|\stackrel{\‾}{d}-|c_i-c_{i-1}\left|\right|}{\max d}+\mathrm{\β}\frac{|c_{i-1}-2c_i+c_{i+1}|}{2\max d},$

$E_{ext}\left(i\right)=\mathrm{\γ}\frac{I\left(x_i,y_i\right)}{255}-\mathrm{\λ}\frac{|\▿I\left(x_i,y_i\right)|}{255\sqrt{2}},$

其中，E_int是内部能量；E_ext是外部能量；c_i(x_i,y_i)(i＝1,2,…,N-2)是第i个snake点；和maxd分别是相邻snake点之间的平均距离和最大距离；α和β是权重参数，取值区间也是[0,1]；I(x_i,y_i)和▽I(x_i,y_i)分别是像素(x_i,y_i)的灰度和灰度梯度值；γ,λ∈[0,1]是权重参数，

通过使能量函数E最小，snake模型从初始形状开始不断变形，最终停留在能量函数取得全局最小值的最优位置，即为目标轮廓。

4.根据权利要求3所述的非门控ICUS图像序列中血管的四维重建方法，其特征是，冠状动脉造影图像序列的两个采集角度之间夹角的取值范围为60°至90°。

5.根据权利要求4所述的非门控ICUS图像序列中血管的四维重建方法，其特征是，在第一时刻图像的超声导管回撤路径中选取的三维重建采样点包括起点、终点和3～6个中间点。

6.根据权利要求5所述的非门控ICUS图像序列中血管的四维重建方法，其特征是，所述权重参数ε的取值范围为[0.1,0.35]，阈值η的取值区间为[10,20]。

7.根据权利要求6所述的非门控ICUS图像序列中血管的四维重建方法，其特征是，snake模型能量函数的全局最优化采用Williams贪婪算法完成。