CN101953696B - 一种icus图像序列中血管的三维形态参数测量方法 - Google Patents

Abstract

一种ICUS图像序列中血管的三维形态参数测量方法，用于提高冠状动脉血管形态参数的测量精度。其技术方案是：首先利用连续回撤超声导管获取的、覆盖多个心动周期的冠状动脉内超声图像序列和在超声导管回撤路径起点采集的、两个近似垂直方位上的X射线冠状动脉造影图像之间的交叉信息完成血管的三维重建，然后采用该三维血管模型，用几何方法对临床重要的血管形态参数进行测量和计算。实验证明，本发明对血管形态参数的测量结果比传统方法更为精确，从而为冠心病的可视化诊治和对介入治疗效果的评价提供了可靠依据。

Description

一种ICUS图像序列中血管的三维形态参数测量方法

技术领域

本发明涉及一种根据连续回撤超声导管获取的、覆盖多个心动周期的冠状动脉内超声(Intracoronary Ultrasound，ICUS)图像序列定量测量冠状动脉血管的三维形态结构参数的方法，属保健技术领域。

背景技术

X射线冠状动脉造影和冠状动脉内超声(Intracoronary Ultrasound，ICUS)是目前临床普遍采用的诊断冠状动脉病变的介入影像手段。客观准确地测量冠状动脉血管及可能存在的病变(包括粥样硬化斑块和狭窄等)的形态结构参数是诊断冠心病的关键，也是制定介入治疗方案的重要依据。传统的诊断过程中，医生需要利用解剖、病理等专门知识和临床经验，根据X射线冠状动脉造影图像想象血管的三维形态，估计其形态和结构参数，因此结果不够客观和准确。基于两个近似正交角度的X射线冠状动脉造影图像的血管三维重建(S.J.Chen，J.D.Carroll.3Dreconstruction of coronary arterial tree to optimize angiographic visualization.IEEETransactions on Medical Imaging.2000，vol.19，no.4，pp.318-336)可为医生提供形象、直观的三维血管图像，同时还可对血管的形态参数(如长度、分支夹角、曲率等)进行定量测量。但是，该重建结果仅描述了血管内腔的空间位置和形态，无法得到血管腔截面的形态和斑块的弥漫程度。且重建时一般假设血管腔横截面为椭圆，事实上当发生狭窄时血管腔的形状复杂多样，狭窄多呈偏心型和不规则型，因此这种假设是不准确的。与造影相比，ICUS可以显示血管腔的横断面、管壁的形态结构、斑块的组织形态学特征等。传统的ICUS图像序列中血管的三维重建方法(P.Schoenhagen.3D Intravascular ultrasound(IVUS)and IVUS-Palpography：insightsinto the mechanica]behavior of the coronary vessel wall.The International Journal ofCardiovascular Imaging.2006，vol.22，pp.153-155)是把一系列的血管内超声图像按采集顺序叠加起来形成一个三维直血管段，完全不考虑血管本身的弯曲和扭曲、导管在回撤过程中的扭曲所造成的图像旋转、以及由于心脏运动和搏动血流所导致的运动伪影，因而此类方法的结果是不准确的，进而也无法保证在此基础上进行的定量测量的精度。

发明内容

本发明的目的是克服现有技术的不足、提供一种ICUS图像序列中血管的三维形态参数测量方法，以提高测量的精度、为冠心病的诊治提供可靠依据。

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

一种ICUS图像序列中血管的三维形态参数测量方法，它首先利用连续回撤超声导管获取的、覆盖多个心动周期的冠状动脉内超声图像序列和在超声导管回撤路径起点采集的、两个近似垂直方位上的X射线冠状动脉造影图像之间的交叉信息完成血管的三维重建，然后采用该三维血管模型，用几何方法对临床需要的血管形态参数进行测量和计算，这些血管形态参数是斑块体积、血管段长度、血管曲率、血管挠率和血管段容积，具体测量和计算方法如下：

a、对原始ICUS图像序列和X射线冠状动脉造影图像进行包括去噪、分割和边缘提取的预处理；

b、从在超声导管回撤路径起点采集的、一对近似垂直方位上的X射线冠状动脉造影图像中三维重建出超声导管的回撤路径；

c、对ICUS图像序列中存在的由周期性心脏运动所致的运动伪影进行补偿；

d、从完成运动伪影抑制的ICUS图像序列中选择出在相同心脏相位处采集的各帧图像；

e、将分割后的、在相同心脏相位处采集的各帧ICUS图像按照采集顺序沿三维导管路径顺序排列；

f、拟合血管腔的内外表面；

g、测量斑块体积：

用经过血管段上下截面重心连线的平面将血管段分割成多个六面体网格体元，按照逆时针顺序分别用V₀、V₁、V₂、V₃表示六面体网格体元一个表面的四个顶点，由此可得两个三角形V₀V₁V₂和V₀V₂V₃，则斑块边界与外膜边界之间的第k个体元的体积V_k为：

$V_k=\left|\underset{j=1}{\overset{6}{\mathrm{\Σ}}}\underset{i=1}{\overset{2}{\mathrm{\Σ}}}\frac{1}{6}\right[V_{j,0}\·(V_{j,i}\×V_{j,i+1})\left]\right|$

其中，V_j，i表示坐标原点O与第j个表面上的第i个三角形的顶点之间的向量，

斑块体积为： $V=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{V}_k$

其中，N是斑块边界与外膜边界之间的六面体网格体元的数目。

上述ICUS图像序列中血管的三维形态参数测量方法，按照采集顺序，对在相同心脏相位处采集的各帧ICUS图像中管腔横截面轮廓重心进行B样条曲线拟合，即可得到三维管腔轴线，利用B样条曲线的积分计算三维管腔轴线的长度，即得到血管段的长度。

上述ICUS图像序列中血管的三维形态参数测量方法，设用B样条曲线表示的三维管腔轴线的曲线方程为c(s)＝(x(s)，y(s)，z(s))(s∈[0，1])，则血管曲率κ和血管挠率τ的计算公式如下：

$\mathrm{\κ}=\sqrt{\frac{{(y^{\′}{z}^{\′\′}-y^{\′\′}{z}^{\′})}^2+{(z^{\′}{x}^{\′\′}-z^{\′\′}{x}^{\′})}^2+{(x^{\′}{y}^{\′\′}-x^{\′\′}{y}^{\′})}^2}{{(x^{\′2}+y^{\′2}+z^{\′2})}^3}}$

$\mathrm{\τ}=\frac{\left|\begin{array}{ccc}{x}^{\′}& y^{\′}& z^{\′}\\ x^{\′\′}& y^{\′\′}& z^{\′\′}\\ x^{\′\′\′}& y^{\′\′\′}& z^{\′\′\′}\end{array}\right|}{{(y^{\′}{z}^{\′\′}-y^{\′\′}{z}^{\′})}^2+{(z^{\′}{x}^{\′\′}-z^{\′\′}{x}^{\′})}^2+{(x^{\′}{y}^{\′\′}-x^{\′\′}{y}^{\′})}^2}$

上述ICUS图像序列中血管的三维形态参数测量方法，对从各帧ICUS图像中提取出的血管壁内膜离散轮廓点进行B样条曲线拟合，得到封闭的、用连续参数曲线表示的内膜轮廓曲线。以该曲线的重心作为坐标原点建立局部坐标系，通过在该坐标系中进行曲线积分求出该曲线所包围的面积，就是该帧ICUS图像中血管壁内膜的横截面积；设血管段上下截面(第i帧和第j帧ICUS图像)的法向矢量分别为n_i和n_j，n₀＝n_i+n_j，上下截面血管壁内膜的横截面积分别为G_i和G_j，上下截面重心之间的矢量为d，n₀和n_i之间的夹角为μ_i，n₀和n_j之间的夹角为μ_j，n₀和d之间的夹角为φ，则上下截面血管壁内膜的横截面在矢量n₀方向上的投影面积分别为：G_i′＝G_i cos μ_i和G_j′＝G_j cos μ_j；血管段上下截面之间的距离为：h＝|d|·cosφ，该段血管的容积V_i为： $V_i=h\·({G_i}^{\′}+\sqrt{{G_i}^{\′}{G_j}^{\′}}+{G_j}^{\′})/3.$

本发明在根据冠状动脉内超声图像序列重建三维血管模型的基础上，采用几何方法，对该三维血管模型的具有重要临床价值的形态结构参数进行定量测量(包括斑块体积、血管段长度、容积、曲率和挠率)，保证了测量结果的客观性和准确性。与现有的方法相比，该方法不仅简化了图像采集过程(即无需ECG门控图像装置，而是连续回撤超声导管采集ICUS图像)，保证了ICUS图像序列数据集合的完整性(无需抛弃有用帧)，而且抑制了ICUS图像序列中存在的、由周期性心脏运动所致的运动伪影，因而重建结果可全面反映冠状动脉血管在整个心动周期中的形态。

实验证明本发明的测量结果较仅采用X射线冠状动脉造影图像重建出的血管模型或传统的从ICUS图像序列中三维重建血管的测量结果更为精确，从而为冠心病的可视化诊治和对介入治疗效果的评价提供了可靠依据。

附图说明

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

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

图2是根据本发明方法的相邻帧ICUS图像的运动参数示意图；

图3是回撤路径上各帧ICUS图像采集点处的单位切矢；

图4是ICUS图像平面垂直于其采集点处的单位切矢；

图5是根据本发明方法的血管段上下截面示意图；

图6是网格化的三维血管段示意图；

图7是六面体网格体元的构造方法示意图；

图8是六面体网格体元体积计算方法示意图。

图中各符号为：γ₁、γ₂、时刻t₁和t₂的血管壁外膜边界曲线；C₁、C₂、γ₁和γ₂的重心；Δx、Δy、C₁和C₂之间分别在x和y方向上的位移；Δα、γ₁和γ₂之间的旋转角；(x₀，y₀)、导管中心的坐标；c(s)、表示超声导管三维回撤路径的B样条曲线；c(s₀)、曲线c(s)上当s＝s₀时的点；c(s₁)、曲线c(s)上当s＝s₁时的点；c(s₂)、曲线c(s)上当s＝s₂时的点；c(s_i)、曲线c(s)上当s＝s_i时的点；c(s_i+1)、曲线c(s)上当s＝s_i+1时的点；t₀、曲线c(s)在点c(s₀)处的单位切矢量；t₁、曲线c(s)在点c(s₁)处的单位切矢量；t₂、曲线c(s)在点c(s₂)处的单位切矢量；t_i、曲线c(s)在点c(s_i)处的单位切矢量；t_i+1、曲线c(s)在点c(s_i+1)处的单位切矢量；t、曲线c(s)上某帧ICUS图像采集点处的单位切矢量；n_i、n_j、血管段上下截面(即两帧ICUS图像)的法向矢量；d、血管段上下截面(即两帧ICUS图像)重心之间的矢量；C_L1、C_L2、血管段上下截面的重心；A、B、C、D、E、F、G、H、六面体网格体元的顶点；V₀、V₁、V₂、V₃、六面体网格体元一个表面的四个顶点。

文中所用符号：γ₁、γ₂、时刻t₁和t₂的血管壁外膜边界曲线；C₁(xc₁，yc₁)、C₂(xc₂，yc₂)、γ₁和γ₂的重心；Δx、Δy、C₁和C₂之间分别在x和y方向上的位移；Δα、γ₁和γ₂之间的旋转角；(x₀，y₀)、导管中心的坐标；Δx_d、Δx的运动分量；Δx_g、Δx的几何分量；Δy_d、Δy的运动分量；Δy_g、Δy的几何分量；Δα_d、Δα的运动分量；Δα_g、Δα的几何分量；M、ICUS图像序列的总帧数；(Δx_k，Δy_k)、第k帧和第k+1帧图像中管腔边界曲线重心之间的位移；Δα_k、第k帧和第k+1帧图像中管腔边界曲线重心之间的旋转角；ΔX_k、Δx_k的幅度谱；ΔY_k、Δy_k的幅度谱；ΔA_k、Δα_k的幅度谱；ΔX_k，d、Δx_k的运动分量的幅度谱；ΔY_k，d、Δy_k的运动分量的幅度谱；ΔA_k，d、Δα_k的运动分量的幅度谱；I_k(x，y)、ICUS图像序列中的第k帧图像；I_k(x′，y′)、补偿运动伪影后的第k帧ICUS图像；c(s)、表示超声导管三维回撤路径的B样条曲线；t、曲线c(s)上某帧图像采集点处的单位切矢量；c′(s)、c(s)的一阶导数；x(s)、y(s)、z(s)、c(s)分别在x、y和z方向上的分量；x′、x″、x″′、x(s)的一阶、二阶和三阶导数；y′、y″、y″′、y(s)的一阶、二阶和三阶导数；z′、z″、z″′、z(s)的一阶、二阶和三阶导数；κ、曲率；τ、挠率；G_i、G_j、第i帧和第j帧ICUS图像中血管壁内膜的横截面积；n_i、n_j、第i帧和第j帧ICUS图像的法向矢量；n₀、血管段的方向矢量；G_i′、G_j′、第i帧和第j帧ICUS图像中的血管壁内膜的横截面在矢量n₀方向上的投影面积；μ_i、n₀和n_i之间的夹角；μ_j、n₀和n_j之间的夹角；d、血管段上下截面重心之间的矢量；h、血管段上下截面之间的距离；φ、n₀和d之间的夹角；V_i、血管段的体积；V_k、六面体网格体元的体积；V、目标区域(斑块)的体积；N、斑块边界与外膜边界之间的六面体网格体元的数目；C_L1、C_L2、血管段上下截面的重心；A、B、C、D、E、F、G、H、六面体网格体元的顶点；V₀、V₁、V₂、V₃、六面体网格体元一个表面的四个顶点；O、坐标原点；V_j，i、坐标原点O与第j个表面上的第i个三角形的顶点之间的向量。

具体实施方式

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

(1)对原始X射线冠状动脉造影图像和ICUS图像进行预处理：

对于在超声导管回撤路径起点采集的、一对近似垂直方位上的X射线冠状动脉造影图像，首先应用Top-Hat变换消除图像中尺寸较大的组织结构噪声，然后再采用旋转一维高斯模板法(XuZhi，Xie Hongbo，Chen Xiaodong，et al.Detection of vessel boundary in coronary angiogram.Proceedings of SPIE International Conference on Electronic Imaging and MultimediaTechnology III.2002，vol.4925，pp.558-562)对图像进行模板匹配增强，消除尺寸较小的噪声。

对于ICUS图像序列，采用snake模型(M.Kass，A.Witkin，D.Terzopoulos.Snakes：activecontour models.International Journal of Computer Vision.1987，vol.1，no.4，pp.321-331)法从各帧图像中提取出血管壁的管腔-斑块边界和外膜边界。

(2)从在超声导管回撤路径起点采集的、一对近似垂直方位上的X射线冠状动脉造影图像中三维重建出超声导管的回撤路径：

本发明采用国家发明专利，一种冠状动脉血管轴线的四维重建方法(专利号：ZL200810055038.3.)的方法，完成一对近似垂直方位上的X射线冠状动脉造影图像中超声导管回撤路径的三维重建。方法的基本原理是：运用snake模型技术(M.Kass，A.Witkin，D.Terzopoulos.Snakes：activecontour models.International Journal of Computer Vision.1987，vol.1，no.4，pp.321-331)，表示导管路径的snake曲线在约束力的作用下直接在三维空间中变形，外部约束力保证曲线在左右成像平面上的投影位于相应的导管投影，内部约束力保证曲线的连续和光滑。通过求解能量函数的最小化问题，使模型变形，直接获得导管路径的三维坐标，结果是用三次B样条曲线表示的三维曲线。

(3)补偿ICUS图像序列中存在的由周期性心脏运动所致的运动伪影：

如附图2所示，相邻帧ICUS图像之间血管腔横截面的平移和旋转可分别用血管壁外膜边界曲线的重心的位移和边界曲线之间的旋转角来表示。设时刻t₁和t₂的外膜边界曲线分别为γ₁和γ₂，其重心分别为C₁(xc₁，yc₁)和C₂(xc₂，yc₂)，γ₁和γ₂之间的位移为：

$\left\{\begin{array}{c}\mathrm{\Δx}={\mathrm{xc}}_2-{\mathrm{xc}}_1\\ \mathrm{\Δy}={\mathrm{yc}}_2-{\mathrm{yc}}_1\end{array}\right.$ ①

γ₁和γ₂之间的旋转角为Δα，旋转中心为超声导管中心(x₀，y₀)：

Δα＝arctg(yc₂/xc₂)-arctg(yc₁/xc₁)    ②

位移(Δx，Δy)和旋转角Δα分别由两部分组成：

$\left\{\begin{array}{c}\mathrm{\Δx}={\mathrm{\Δx}}_d+{\mathrm{\Δx}}_g\\ \mathrm{\Δy}=\mathrm{\Δ}{y}_d+\mathrm{\Δ}{y}_g\\ \mathrm{\Δ\α}={\mathrm{\Δ\α}}_d+{\mathrm{\Δ\α}}_g\end{array}\right.$ ③

其中脚标d和g分别表示运动分量和几何分量。本发明方法采用谱分析的方法完成运动和几何分量的分离，即对计算出的一系列(Δx_k，Δy_k，Δα_k)(k＝1，2，...，M，M为ICUS图像序列的总帧数)进行傅立叶变换，得到其幅度谱ΔX_k、ΔY_k和ΔA_k，则运动和几何分量分别对应于其中的高频和低频分量。然后采用高通滤波器分别对ΔX_k、ΔY_k和ΔA_k进行滤波，滤波器的输出即是运动分量的幅度谱ΔX_k，d、ΔY_k，d和ΔA_k，d，再对其进行逆傅立叶变换即可得到Δx_k，d、Δy_k，d和Δα_k，d。由于运动分量主要由周期性心脏运动引起，因此本发明方法将高通滤波器的通带截止频率设定为病人的心率值(单位：次/秒，即Hz)，该值可以从ICUS图像序列中估计出来，具体方法参见(国家发明专利，一种提高冠状动脉内超声图像视觉效果的方法.申请号：200910075133.4)。

估计出ICUS图像序列中的第k帧图像I_k(x，y)(k＝2，3，...，M)的(Δx_k，d，Δy_k，d，Δα_k，d)之后，将I_k(x，y)中、基于以导管中心为坐标原点的坐标系的各像素点的坐标先反向平移，再反向旋转，即得到补偿运动伪影后的图像I_k(x′，y′)：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=\left[\begin{array}{cc}\cos (-\underset{i=2}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\Δ\α}}_{i,d})& \sin (-\underset{i=2}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\Δ\α}}_{i,d})\\ -\sin (-\underset{i=2}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\Δ\α}}_{i,d})& \cos (-\underset{i=2}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\Δ\α}}_{i,d})\end{array}\right]\left[\begin{array}{c}x-\underset{i=2}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\Δx}}_{i,d}\\ y-\underset{i=2}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\Δy}}_{i,d}\end{array}\right]$ ④

(4)从完成运动伪影抑制的ICUS图像序列中选择出在相同心脏相位处采集的各帧图像：

本发明采用(国家发明专利，一种提高冠状动脉内超声图像视觉效果的方法.申请号：200910075133.4)中的方法，对于完成运动伪影抑制的ICUS图像序列，找到相邻心动周期中各帧ICUS图像之间的对应关系，并选择出在相同心脏相位处采集的各帧图像。方法的基本原理是：首先分析各帧ICUS图像之间灰度特征的差异度，构造ICUS图像序列的差异矩阵；然后采用动态规划算法(Geiger D，Gupta A，Vlontzos JA，Vlontzos J.Dynamic programming for detecting，tracking and matching deformable contours.IEEE Transactions on Pattern Analysis andMachine Intelligence.1995，vol.17，no.3，pp.294-302)，通过在差异矩阵中寻找一条累计差异度最小的最优路径，为各帧ICUS图像找到其在相邻心动周期中的对应帧。

(5)将分割后的、在相同心脏相位处采集的各帧ICUS图像按照采集顺序沿三维导管路径顺序排列：

如附图3所示，设c(s)(0≤s≤1)是表示超声导管三维回撤路径的B样条曲线，为了沿c(s)准确排列在相同心脏相位处采集的各帧ICUS图像，计算c(s)上各帧图像采集点处的单位切矢量t＝c′(s)/|c′(s)|，然后如附图4所示，在沿导管路径按照采集顺序依次排列各帧ICUS图像时，令各帧图像平面垂直于其采集点处的单位切矢量即可。

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

采用NURBS(非均匀有理B样条)曲面，对沿导管路径排列的管腔横截面中的血管壁边界点进行拟合，得到光滑、连续的三维血管腔内外表面。

(7)根据重建出的三维血管模型定量测量血管的形态结构参数，包括血管段长度、容积、曲率、挠率和斑块体积的测量：

(7.1)血管段长度：

如附图4所示，按照采集顺序，对在相同心脏相位处采集的各帧ICUS图像中管腔横截面轮廓重心进行B样条曲线拟合，即可得到三维管腔轴线。利用B样条曲线的积分计算三维管腔轴线的长度，即得到血管段的长度。

(7.2)血管段曲率和挠率：

在得到用B样条曲线表示的三维管腔轴线之后，利用微分几何中的相关公式即可计算出管腔轴线上各点处的曲率和挠率。设曲线方程为c(s)＝(x(s)，y(s)，z(s))(s∈[0，1])，曲率κ和挠率τ的计算公式如下：

$\mathrm{\κ}=\sqrt{\frac{{(y^{\′}{z}^{\′\′}-y^{\′\′}{z}^{\′})}^2+{(z^{\′}{x}^{\′\′}-z^{\′\′}{x}^{\′})}^2+{(x^{\′}{y}^{\′\′}-x^{\′\′}{y}^{\′})}^2}{{(x^{\′2}+y^{\′2}+z^{\′2})}^3}}$ ⑤

$\mathrm{\τ}=\frac{\left|\begin{array}{ccc}{x}^{\′}& y^{\′}& z^{\′}\\ x^{\′\′}& y^{\′\′}& z^{\′\′}\\ x^{\′\′\′}& y^{\′\′\′}& z^{\′\′\′}\end{array}\right|}{{(y^{\′}{z}^{\′\′}-y^{\′\′}{z}^{\′})}^2+{(z^{\′}{x}^{\′\′}-z^{\′\′}{x}^{\′})}^2+{(x^{\′}{y}^{\′\′}-x^{\′\′}{y}^{\′})}^2}$ ⑥

(7.3)血管段容积：

对从各帧ICUS图像中提取出的血管壁内膜离散轮廓点进行B样条曲线拟合，得到封闭的、用连续参数曲线表示的内膜轮廓曲线。以该曲线的重心作为坐标原点建立局部坐标系，通过在该坐标系中进行曲线积分求出该曲线所包围的面积，就是该帧ICUS图像中血管壁内膜的横截面积。设G_i和G_j分别是第i帧和第j帧ICUS图像中血管壁内膜的横截面积。如附图5所示，设血管段上下截面(即第i帧和第j帧ICUS图像)的法向矢量分别为n_i和n_j，则该段血管的方向矢量n₀近似为：

n₀＝n_i+n_j           ⑦

第i帧ICUS图像中的血管壁内膜的横截面在矢量n₀方向上的投影面积为：

G_i′＝G_i cos μ_i    ⑧

其中μ_i是n₀和n_i之间的夹角。第j帧ICUS图像中的血管壁内膜的横截面在矢量n₀方向上的投影面积为：

G_j′＝G_j cos μ_j    ⑨

其中μ_j是n₀和n_j之间的夹角。该段血管上下截面之间的距离为：

h＝|d|·cos φ      ⑩

其中矢量d是血管段上下截面重心之间的矢量；φ是n₀和d之间的夹角。则该段血管的体积V_i可用圆台体积公式近似计算得到：

$V_i=h\·({G_i}^{\′}+\sqrt{{G_i}^{\′}{G_j}^{\′}}+{G_j}^{\′})/3---\left(11\right)$

(7.4)斑块体积：

如附图6所示，先计算内腔-斑块边界与外膜边界之间的各个六面体网格体元的体积V_k，那么目标区域(斑块)的体积V就是这些网格体元的体积之和：

$V=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{V}_k---\left(12\right)$

其中N是斑块边界与外膜边界之间的六面体网格体元的数目。六面体网格体元的构成方法如附图7所示，其中C_L1和C_L2是血管段上下截面的重心，则ABCDEFGH就构成了一个六面体网格体元。如附图8所示，体元的每个表面有四个顶点，按照逆时针顺序分别用V₀、V₁、V₂、V₃表示。由此可得到两个三角形V₀V₁V₂和V₀V₂V₃，每个三角形和坐标原点O都构成一个四面体。对于体元的每个表面，计算出各四面体的体积，然后相加即可求出体元体积V_k：

$V_k=\left|\underset{j=1}{\overset{6}{\mathrm{\Σ}}}\underset{i=1}{\overset{2}{\mathrm{\Σ}}}\frac{1}{6}\right[V_{j,0}\·(V_{j,i}\×V_{j,i+1})\left]\right|---\left(13\right)$

其中V_j，i表示坐标原点O与第j个表面上的第i个三角形的顶点之间的向量。计算结果与坐标原点的选取无关。

Claims (4)

1.一种ICUS图像序列中血管的三维形态参数测量方法，其特征是，首先利用连续回撤超声导管获取的、覆盖多个心动周期的冠状动脉内超声图像序列和在超声导管回撤路径起点采集的、两个近似垂直方位上的X射线冠状动脉造影图像之间的交叉信息重建血管的三维血管模型，然后采用该三维血管模型，用几何方法对临床需要的血管形态参数进行测量和计算，这些形态参数是斑块体积、血管段长度、血管曲率、血管挠率和血管段容积，具体测量和计算方法如下：

a、对原始ICUS图像序列和X射线冠状动脉造影图像进行去噪、分割和边缘提取的预处理；

b、从在超声导管回撤路径起点采集的、一对近似垂直方位上的X射线冠状动脉造影图像中三维重建出超声导管的回撤路径；

c、对ICUS图像序列中存在的由周期性心脏运动所致的运动伪影进行补偿；

d、从完成运动伪影抑制的ICUS图像序列中选择出在相同心脏相位处采集的各帧图像；

e、将分割后的、在相同心脏相位处采集的各帧ICUS图像按照采集顺序沿三维导管路径顺序排列；

f、拟合血管腔的内外表面；

g、测量斑块体积：

用经过血管段上下截面重心连线的平面将血管段分割成多个六面体网格体元，按照逆时针顺序分别用V₀、V₁、V₂、V₃表示六面体体元一个表面的四个顶点，由此可得两个三角形V₀V₁V₂和V₀V₂V₃，则斑块边界与外膜边界之间的第k个体元的体积V_k为：

$V_k=\left|\underset{j=1}{\overset{6}{\mathrm{\Σ}}}\underset{i=1}{\overset{2}{\mathrm{\Σ}}}\frac{1}{6}\right[V_{j,0}\·(v_{j,i}\×V_{j,i+1})\left]\right|$

其中，V_j，i表示坐标原点O与第j个表面上的第i个三角形的顶点之间的向量，

斑块体积为： $V=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{V}_k$

其中，N是斑块边界与外膜边界之间的六面体网格体元的数目。

2.根据权利要求1所述ICUS图像序列中血管的三维形态参数测量方法，其特征是，按照采集顺序，对在相同心脏相位处采集的各帧ICUS图像中管腔横截面轮廓重心进行B样条曲线拟合，即可得到三维管腔轴线；利用B样条曲线的积分计算三维管腔轴线的长度，即得到血管段的长度。

3.根据权利要求1或2所述ICUS图像序列中血管的三维形态参数测量方法，其特征是，设用B样条曲线表示的三维管腔轴线的曲线方程为c(s)＝(x(s)，y(s)，z(s))(s∈[0，1])，则血管曲率κ和血管挠率τ的计算公式如下：

$\mathrm{\κ}=\sqrt{\frac{{(y^{\′}{z}^{\′\′}-y^{\′\′}{z}^{\′})}^2+{(z^{\′}{x}^{\′\′}-z^{\′\′}{x}^{\′})}^2+{(x^{\′}{y}^{\′\′}-x^{\′\′}{y}^{\′})}^2}{{(x^{\′2}+y^{\′2}+z^{\′2})}^3}}$

$\mathrm{\τ}=\frac{\left|\begin{array}{ccc}{x}^{\′}& y^{\′}& z^{\′}\\ x^{\′\′}& y^{\′\′}& z^{\′\′}\\ x^{\′\′\′}& y^{\′\′\′}& z^{\′\′\′}\end{array}\right|}{{(y^{\′}{z}^{\′\′}-y^{\′\′}{z}^{\′})}^2+{(z^{\′}{x}^{\′\′}-z^{\′\′}{x}^{\′})}^2+{(x^{\′}{y}^{\′\′}-x^{\′\′}{y}^{\′})}^2};$

式中，x′、x″、x′″分别是x(s)的一阶、二阶和三阶导数；y′、y″、y′″分别是y(s)的一阶、二阶和三阶导数；z′、z″、z′″分别是z(s)的一阶、二阶和三阶导数。

4.根据权利要求1或2所述ICUS图像序列中血管的三维形态参数测量方法，其特征是，对从各帧ICUS图像中提取出的血管壁内膜离散轮廓点进行B样条曲线拟合，得到封闭的、用连续参数曲线表示的内膜轮廓曲线；以该曲线的重心作为坐标原点建立局部坐标系，通过在该坐标系中进行曲线积分求出该曲线所包围的面积，就是该帧ICUS图像中血管壁内膜的横截面积；设血管段上截面第i帧、下截面第j帧ICUS图像的法向矢量分别为n_i和n_j，n₀＝n_i+n_j，上下截面血管壁内膜的横截面积分别为G_i和G_j，上下截面重心之间的矢量为d，n₀和n_i之间的夹角为μ_i，n₀和n_j之间的夹角为μ_j，n₀和d之间的夹角为φ，则上下截面血管壁内膜的横截面在矢量n₀方向上的投影面积分别为：G_i′＝G_i cos μ_i和G_j′＝G_j cos μ_j；血管段上下截面之间的距离为：h＝|d|·cosφ，所述血管段的容积V_i为：

Images