CN103268630A - 一种基于血管内超声影像的血管三维可视化方法 - Google Patents

Abstract

一种基于血管内超声影像的血管三维可视化方法，涉及计算机医学图像分析领域，其特征在于，首先，结合多图像平均去噪、中值滤波和小波降软阈值噪方法对图像序列进行降噪处理，该方法能减少图像噪声，很好的保留图像的重要细节信息，并且图像降噪效率高；其次，利用二次多项式拟合图像形变，实现图像配准，以补偿图像序列采集过程中产生的变形；再次，利用光线投射算法绘制出三维血管模型；最后利用切片重组方法实现对三维血管模型的任意角度平面剖切，显示血管内部结构信息，为病变分析创造了条件。

Description

一种基于血管内超声影像的血管三维可视化方法

技术领域

本发明涉及计算机医学图像分析领域，特别涉及一种基于血管内超声影像的血管三维可视化方法。

背景技术

血管内超声图像是一种断层切片图像，能够显示当前位置血管的横截面图，详细描述血管壁、内腔和斑块的组织成分，进而对血管腔径、横截面积进行计算，而且可根据斑块声学特征对其进行组织学分析，发现早期粥样硬化斑，在斑块病变诊断上具有冠状动脉造影无法比拟的优势。临床中，对采集到的某段血管的血管内超声(Intravascular Ultrasound，IVUS)图像序列，医生不仅限于单独观察每一帧图像，还需要主观重构血管的三维视图，以便帮助理解血管及病变的空间毗邻关系，或是比较手术段和非手术段的差异。从而，不可避免地造成诊断结果的主观性，同时也给医生的工作带来一定困难。

基于IVUS图像序列的血管三维可视化，是采用合适的方法由血管内超声序列重建出可以从任意视角进行观察的血管三维投影图像，并通过对重建模型进行平面剖切或是强化图像中的细节，清晰地显示血管的复杂特征和空间定位关系，能够帮助医生理解血管及病变的空间毗邻关系，做出正确的医疗诊断方案，对临床应用具有很大的价值。

目前，医学图像三维重建方法大致可以分为面绘制和体绘制两大类。面绘制能得到三维表面的细节描述，但对于血管等亮度变化小、形状不明显的重建数据，绘制的效果并不理想。而体绘制利用的是全部体数据，能保留每一个细节，对于形状特征模糊不清的组织和器官，如血管等软组织具有较好的三维显示效果。对于由体绘制得到的三维模型，虽然能显示完整的三维信息，却缺乏形象、清晰的局部信息。利用切片重组方法，则能够实现对三维模型进行任意角度平面剖切，显示医生感兴趣的局部切面信息，为诊断病情提供直观、准确的依据。

由于成像设备和成像过程的特殊性，超声图像在采集过程中，易受到加性随机噪声的污染。另外，在IVUS图像序列获取的过程中，由于血管壁部位的自发生理运动以及患者的移动都将会使内部组织的位置、形状和大小发生变化。因此，在血管三维重建之前，应对IVUS图像序列进行降噪和配准，一方面可以降低图像噪声对三维重建效果的影响，另一方面可以利用非刚性变换来补偿图像的变形。

血管三维可视化方法，以血管内超声影像为基础，根据IVUS图像序列的特点，进行图像序列降噪、配准，血管三维重建，以及对血管模型的任意角度平面剖切，能够得到较好的血管三维可视化效果。这一方法在不增加附加设备的情况下，充分利用现有血管内超声设备提供的超声图像信息和三维可视化方法，能够得到较好的血管三维显示效果。

发明内容

本发明的目的在于，通过提供一种基于血管内超声影像的血管三维可视化方法，以获得直观、形象的血管三维信息，本发明的特征如下：

步骤(1)利用血管内超声仪，以0.5mm/s的速度匀速回拉导管，获得人体冠状动脉的血管内超声视频影像；

步骤(2)将步骤(1)得到的血管内超声视频影像导入计算机，从视频中截取连续的900帧血管内超声图像作为实验图像，图像分辨率为384*384，以下简称为超声图像；

步骤(3)对上述超声图像依次按以下步骤进行超声图像序列降噪和超声图像噪声平滑：

步骤(3.1)取连续10帧超声图像组成一组序列，表示为I_k(x,y),1≤k≤10，k表示超声图像序号，对每帧超声图像进行[3,3]的中值滤波，得到中值滤波图像序列

$I_k^{\′}(x,y)=\underset{\left[\mathrm{3,3}\right]}{\mathrm{med}}\left\{I_k(x,y)\right\},$ 符号med{·}表示中值滤波，

步骤(3.2)按下式求中值滤波图像序列

的平均图像 ${\stackrel{\‾}{I}}^{\′}(x,y)=\frac{1}{10}\underset{k=1}{\overset{10}{\mathrm{\Σ}}}\left[I_k^{\′}(x,y)\right],$ 1≤k≤10，

步骤(3.3)对平均图像

按下式进行小波去噪，得到小波变换系数矩阵WT：

其中ψ^*(x,y)是小波函数ψ(x,y)的共轭， $\mathrm{\ψ}(x,y)=a_0^{-j/2}\mathrm{\ψ}(\frac{x-m{a}_0^j{a}_1}{a_0^j},\frac{y-n{a}_0^j{a}_2}{a_0^j}),$ 其中a₀，a₁，a₂是设定值，a₀＝a₁＝a₂＝2，称为扩展步长，1≤j≤3，是分辨索引，m≥1，n≥1，表示水平和垂直方向的有限平移，

步骤(3.4)对小波系数WT进行软阈值化： $W_{\mathrm{\&zeta;}}=\left\{\begin{array}{cc}\mathrm{sgn}\left(\mathrm{WT}\right)\left(\right|\mathrm{WT}|-\mathrm{\&zeta;}),& \left|\mathrm{WT}\right|\&GreaterEqual;\mathrm{\&zeta;}\\ 0,& \left|\mathrm{WT}\right|<\mathrm{\&zeta;}\end{array}\right.$ 其中sgn()为符号函数，若WT＞0，sgn(WT)＝1，若WT＝0，sgn(WT)＝0，若WT＜0，sgn(WT)＝-1，ζ为阈值，

式中L为信号长度，σ为噪声强度，设图像噪声为高斯白噪声，且σ＝1，

步骤(3.5)按下式得到小波重构后的小波软阈值降噪图像

${\stackrel{\&OverBar;}{I}}^{\&prime;\&prime;}(x,y)=\underset{\&OverBar;\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\underset{\&OverBar;\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&zeta;}}\mathrm{\&psi;}(x,y),$

步骤(3.6)按以下步骤得到最终降噪图像

$I_k^{\&prime;\&prime;\&prime;}(x,y)=I_k^{\&prime;\&prime;}(x,y)+{\stackrel{-}{I}}^{\&prime;\&prime;}(x,y),1\&le;k\&le;10,$ 其中为中值滤波图像

与平均图像

的差值：

步骤(4)按以下步骤实现图像序列的配准以补偿图像形变：

步骤(4.1)把步骤(3.6)得到的最终降噪图像作为待配准图像G(x,y)，与预先设定的基准图像F(x,y)组成一个图像组合[F(x,y),G(x,y)]，从中选取12对控制点，分别记为f(x_i,y_i)和g(x_s,y_s)，1≤i≤12，1≤s≤12，符合如下关系：f(x_i,y_i)＝H^-1[g(x_s,y_s)]，其中H^-1为形变关系， $\left\{\begin{array}{c}{H_{x_i}}^{-1}=b_1{x}_s+b_2{y}_s+b_3{x_s}^2+b_4{y_s}^2+b_5{x}_s{y}_s+b_6\\ {H_{y_i}}^{-1}=c_1{x}_s+c_2{y}_s+c_3{x_s}^2+c_4{y_s}^2+c_5{x}_s{y}_s+c_6\end{array}\right.,$ b₁…b₆，c₁…c₂为形变系数，共12个，

步骤(4.2)根据选定的12个控制点f(x_i,y_i)和g(x_s,y_s)求12个形变系数，公式如下 $\left\{\begin{array}{c}{x}_i=b_1{x}_s+b_2{y}_s+b_3{x_s}^2+b_4{y_s}^2+b_5{x}_s{y}_s+b_6\\ y_i=c_1{x}_s+c_2{y}_s+c_3{x_s}^2+c_4{y_s}^2+c_5{x}_s{y}_s+c_6\end{array}\right.,$

步骤(4.3)按下式求得配准的超声图像F′(x,y)，F′(x,y)＝H^-1[G(x,y)]；

步骤(5)把步骤(4.2)得到的配准的超声图像F′(x,y)，按以下步骤利用光线投射算法，重建血管三维模型：

步骤(5.1)利用配准的超声图像构建体数据场，为体数据场设定0,50,200,255四个灰度阈值，将体数据场分为Q₁＝[0,50]，Q₂＝[50,200]，Q₃＝[200,255]共三个阈值区间，

步骤(5.2)按下式将体数据场中各数据点的灰度值映射成为直接用于绘制的不透明度值

$A_{Q_q}\left(\mathrm{\&zeta;}\right)=\left\{\begin{array}{cc}0& 0\&le;\mathrm{\&zeta;}\&le;(c-\frac{w}{2})\\ \frac{\mathrm{\&zeta;}-c}{w\text{'}}-\frac{w}{2w\text{'}}& (c-\frac{w}{2})<\mathrm{\&zeta;}\&le;(c-\frac{w}{2}+w\text{'})\\ 1& (c-\frac{w}{2}+w\text{'})<\mathrm{\&zeta;}<(c+\frac{w}{2}+w\text{'})\\ \frac{c-\mathrm{\&zeta;}}{w\text{'}}+\frac{w}{2w\text{'}}& (c+\frac{w}{2}-w\text{'})<\mathrm{\&zeta;}\&le;(c+\frac{w}{2})\\ 0& (c+\frac{w}{2})<\mathrm{\&zeta;}\&le;255\end{array}\right.,$ 其中q为阈值区间Q₁或Q₂或Q₃的序号，ζ为数据点的灰度值，感兴趣物质的灰度值范围取决于数据中心c和宽度w两个变量，w′表示斜坡的宽度，斜坡表示的是线性的增加或降低，不透明度值用

表示，取值范围在0到1之间，

步骤(5.3)按下式给体数据场的三个阈值区间Q₁＝[0,50]，Q₂＝[50,200]，Q₃＝[200,255]赋不同颜色值

$C_{Q_q}\left(\mathrm{\&zeta;}\right)=\left\{\begin{array}{cc}{C}_1& 0\&le;\mathrm{\&zeta;}\&le;50\\ C_2& 50<\mathrm{\&zeta;}\&le;200\\ C_3& 200<\mathrm{\&zeta;}\&le;255\end{array}\right.,$ 其中q为阈值区间Q₁，Q₂，Q₃的序号，C表示颜色值域，C₁，C₂，C₃是设定的颜色值，

步骤(5.4)为数据场建立X,Y,Z三维坐标轴，将数据场中单位体积的立方体看作一个体素，

步骤(5.5)光源光线透过数据场，在屏幕上形成一个成像平面，从成像平面的每个像素点发出一条穿过数据场的光线，沿着光线选择设定的有限的K个等间距的采样点，采用下式计算各采样点的不透明度值：

$\begin{array}{c}{A}_e=A_1+x(1-y)(1-z)(A_2-A_1)+x(1-y)z(A_3-A_4)+(1-y)z(A_4-A_1)\\ +y(A_5-A_1)+\mathrm{xy}(1-z)A(A_6-A_5)+\mathrm{xyz}(A_7-A_8)+\mathrm{yz}(A_8-A_5)\end{array},$ 其中A_e代表当前采样点的不透明度值，A₁，A₂，…，A₇，A₈代表距离采样点最近的8个数据点的不透明度值，

步骤(5.6)按下式计算采样点的颜色值：

$\begin{array}{c}{C}_l^{\&prime;}=C_1^{\&prime;}+x(1-y)(1-z)(C_2^{\&prime;}-C_1^{\&prime;})+x(1-y)z(C_3^{\&prime;}-C_4^{\&prime;})+(1-y)z(C_4^{\&prime;}-C_1^{\&prime;})\\ +y(C_5^{\&prime;}-C_1^{\&prime;})+\mathrm{xy}(1-z)A(C_6^{\&prime;}-C_5^{\&prime;})+\mathrm{xyz}(C_7^{\&prime;}-C_8^{\&prime;})+\mathrm{yz}(C_8^{\&prime;}-C_5^{\&prime;})\end{array},,$ 其中

代表当前采样点的颜色值，代表距离采样点最近的8个数据点的颜色值，

步骤(5.7)对光线上的采样点进行累加，直到不透明度值增加到1，结束累加，此时的颜色值就是成像平面上像素的最终颜色，即得到三维血管模型，公式如下： $\begin{array}{c}{C}_{\mathrm{out}}{A}_{\mathrm{out}}=C_{\mathrm{in}}{A}_{\mathrm{in}}+C_{\mathrm{now}}{A}_{\mathrm{now}}(1-A_{\mathrm{in}})\\ A_{\mathrm{out}}=A_{\mathrm{in}}+A_{\mathrm{now}}(1-A_{\mathrm{in}})\end{array},$ 其中C_out、A_out分别为经过第u个采样点后的颜色值、不透明度值，C_now、A_now为第u个采样点的颜色值、不透明度值，C_in、A_in为已合成的前u-1个采样点的颜色值、不透明度值，u＝1,2,...,u,...,K，K为采样点个数；

步骤(6)按以下步骤对所述三维血管模型进行任意方向平面裁剪，获得血管内部信息：

步骤(6.1)在由X,Y轴组成的水平面X-Y上，取相邻的4个坐标点：x，x+1，y，y+1构成正方形的4个顶点：(x,y)，(x,y+1)，(x+1,y)和((x+1),(y+1))，形成一个水平的剖切平面，

步骤(6.2)在水平面X-Y的高度Z方向，建立一个三维坐标空间，在有限个数的z坐标值上建立P个在Z轴方向上相互平行的空间剖切平面p，构成一个体素，每个剖切平面p与三维血管模型的4条棱线共有4个交点：z(x_p,y_p)，z(x_p,(y+1)_p)，z((x+1)_p,y_p)和z((x+1)_p,(y+1)_p)，p=1,2,…,P，从剖切平面p与三维血管模型棱线的交点中，任意选择4个能构成斜切平面的点，就可以实现任意斜面剖切，

步骤(6.3)每一个剖切平面p的中心点即为采样点O_p，坐标为

其中： ${\stackrel{\&OverBar;}{x}}_p=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}{x}_p,{\stackrel{\&OverBar;}{y}}_p=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}{y}_p,{\stackrel{\&OverBar;}{z}}_p=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}{z}_p,$

步骤(6.4)每一个空间剖切平面p的各顶点z(x_p,y_p)，z(x_p,(y+1)_p)，z((x+1)_p,y_p)和z((x+1)_p,(y+1)_p)到对应采样点的距离d_p由下式得到： ${d_p}^2=(x_{p,f}-{\stackrel{\&OverBar;}{x}}_{p,O})+(y_{p,f}-{\stackrel{\&OverBar;}{y}}_{p,O})+(z_{p,f}-{\stackrel{\&OverBar;}{z}}_{p,O}),$ 其中f是体素的顶点，上下共8个，f＝0,1,2,…,7，x_p,f是第p个空间剖切平面的第f个顶点的横坐标，

是第p个空间剖切平面中心点O的横坐标，其它类推，

步骤(6.5)按下式计算某一个剖切平面p的中心点O_p的灰度值

其中：

d_p为所述第p个空间剖切平面的8个顶点到中心点距离之和，

d_p,f为所述第p个空间剖切平面上某一个顶点f到中心点的距离，

h_f为入射光线在所述某个顶点f处发出的光线强度。

本发明的效果为：

实验利用血管内超声成像仪，以0.5mm/s的速度匀速回拉导管，获得图像尺寸大小为384*384的超声图像序列，取两组超声图像序列作为实验图像，每组超声图像序列包含连续的900帧超声图像。超声图像序列降噪效果见图2，超声图像配准结果见图3，血管三维可视化效果见图4。可以看出，本文提出的三维可视化方法，可以综合血管内超声图像的特征，较好的显示血管三维信息，为病变分析创造了条件。

附图说明

图1是原血管内超声图像；

图2是连续10帧血管内超声图像降噪效果图，图2.1到2.10分别为第一帧到第十帧图像的降噪效果；

图3是超声图像配准结果图，图3.1是基准图像，图3.2是待配准图像，图3.3是配准结果图像；

图4是血管三维可视化效果图，图4.1为由第一组超声图像序列得到的血管三维可视化效果图，图4.2为由第二组超声图像序列得到的血管三维可视化效果图；

图5是本发明方法的主程序流程图。

具体实施方式

本发明是采用以下技术手段实现的：

一种基于血管内超声影像的血管三维可视化方法。首先，结合图像序列平均、中值滤波和小波软阈值降噪方法对血管内超声图像序列进行降噪处理，然后，利用二次多项式拟合图像形变，实现超声图像配准，补偿图像序列采集过程中产生的变形，最后，利用光线投射算法和切片重组方法实现血管三维模型的重建及任意角度平面剖切，以获得直观、形象的血管三维可视化效果。

上述基于血管内超声影像的血管三维可视化方法，包括下述步骤：

步骤1、利用血管内超声仪，以0.5mm/s的速度匀速回拉导管，获得人体冠状动脉的血管内超声视频影像；

步骤2、将步骤(1)得到的血管内超声视频影像导入计算机，从视频中截取连续的900帧血管内超声图像作为实验图像，图像分辨率为384*384，以下简称为超声图像；

步骤3、图像序列降噪，平滑图像噪声，并尽可能保留原图像中的细节信息；

取相邻10帧IVUS图像组成一个图像序列，表示为I_k(x,y),1≤k≤10，对每帧图像做[3,3]中值滤波：

$I_k^{\&prime;}(x,y)=\underset{\left[\mathrm{3,3}\right]}{\mathrm{med}}\left\{I_k(x,y)\right\}---\left(1\right)$

式中med{·}表示中值滤波运算符。求图像序列的平均图像：

${\stackrel{\&OverBar;}{I}}^{\&prime;}(x,y)=\frac{1}{10}\underset{k=1}{\overset{10}{\mathrm{\&Sigma;}}}\left[I_k^{\&prime;}(x,y)\right]---\left(2\right)$

对平均图像

进行小波分解，求得图像小波变换系数矩阵WT：

$\begin{array}{c}\mathrm{WT}=<{\stackrel{\&OverBar;}{I}}^{\&prime;}(x,y),\mathrm{\&psi;}(x,y)>\\ =\&Integral;\&Integral;{\stackrel{\&OverBar;}{I}}^{\&prime;}(x,y){\mathrm{\&psi;}}^{*}(x,y)\mathrm{dxdy}\end{array}---\left(3\right)$

式中，ψ^*(x,y)是小波函数ψ(x,y)的共轭。ψ(x,y)由下式求得：

$\mathrm{\&psi;}(x,y)=a_0^{-j/2}\mathrm{\&psi;}(\frac{x-m{a}_0^j{b}_1}{a_0^j},\frac{y-n{a}_0^j{b}_2}{a_0^j})---\left(4\right)$

其中a₀，a₁，a₂是设定值，做a₀＝a₁＝a₂＝2，称为扩展步长，1≤j≤3，是分辨索引，m≥1，n≥1，表示水平和垂直方向的有限平移。

对小波系数矩阵WT进行软阈值化：

$W_{\mathrm{\&zeta;}}=\left\{\begin{array}{cc}\mathrm{sgn}\left(\mathrm{WT}\right)\left(\right|\mathrm{WT}|-\mathrm{\&zeta;}),& \left|\mathrm{WT}\right|\&GreaterEqual;\mathrm{\&zeta;}\\ 0,& \left|\mathrm{WT}\right|<\mathrm{\&zeta;}\end{array}\right.----\left(5\right)$

其中，其中sgn()为符号函数，若WT＞0，sgn(WT)＝1，若WT＝0，sgn(WT)＝0，若WT＜0，sgn(WT)＝-1，ζ为阈值，

式中L为信号长度，σ为噪声强度，设图像噪声为高斯白噪声，且σ＝1。

由矩阵W_t，进行小波重构，得到小波软阈值降噪结果

${\stackrel{\&OverBar;}{I}}^{\&prime;\&prime;}(x,y)=\underset{\&OverBar;\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\underset{\&OverBar;\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&zeta;}}\mathrm{\&psi;}(x,y)---\left(6\right)$

计算各帧图像

与

的差值：

$I_k^{\&prime;\&prime;}(x,y)=I_k^{\&prime;}(x,y)-{\stackrel{\&OverBar;}{I}}^{\&prime;}(x,y),1\&le;k\&le;10---\left(7\right)$

将与

相加，得到最终降噪图像

$I_k^{\&prime;\&prime;\&prime;}(x,y)=I_k^{\&prime;\&prime;}(x,y)+{\stackrel{\&OverBar;}{I}}^{\&prime;\&prime;}(x,y),1\&le;k\&le;10---\left(8\right)$

综上所述，图像序列去噪的步骤总结如下：

①取相邻10帧IVUS图像组成一个图像序列，记为I_k(x,y),1≤k≤10；

②对每帧图像进行[3,3]的中值滤波，结果记为1≤k≤10；

③求图像序列

(1≤k≤10)的平均图像，记为I′(x,y)；

④计算图像

与平均图像I′(x,y)的差值，记为

⑤对平均图像I′(x,y)进行小波软阈值降噪得到图像I′′(x,y)；

⑥将与I′′(x,y)相加，得到最终降噪图像

步骤4、实现图像序列配准，补偿图像形变；

把最终降噪图像

作为待配准图像G(x,y)，与预先设定的基准图像F(x,y)组成一个图像组合??F(x,y),G(x,y)??，从中选取12对控制点记为分别f(x_i,y_i)和g(x_s,y_s)，1≤i≤12，1≤s≤12，它们之间存在如下关系：

f(x_i,y_i)＝H^-1[g(x_s,y_s)](9)其中H^-1为形变关系，表达式为：

$\left\{\begin{array}{c}{H_{x_i}}^{-1}=b_1{x}_s+b_2{y}_s+b_3{x_s}^2+b_4{y_s}^2+b_5{x}_s{y}_s+b_6\\ {H_{y_i}}^{-1}=c_1{x}_s+c_x{y}_s+c_3{x_s}^2+c_4{y_s}^2+c_5{x}_s{y}_s+c_6\end{array}\right.---\left(10\right)$

式中b₁…b₆，c₁…c₂为形变系数，共12个。

将控制点值代入(10)式得：

$\left\{\begin{array}{c}{x}_i=b_1{x}_s+b_2{y}_s+b_3{x_s}^2+b_4{y_s}^2+b_5{x}_s{y}_s+b_6\\ y_i=c_1{x}_s+c_x{y}_s+c_3{x_s}^2+c_4{y_s}^2+c_5{x}_s{y}_s+c_6\end{array}\right.---\left(11\right)$

解方程(11)式便可计算出形变系数b₁…b₆，c₁…c₂进而能得到形变关系H^-1。

由待配准图像G(x,y)和形变关系H^-1求得配准结果为：

F′(x,y)＝H^-1[G(x,y)]              (12)

综上所述的图像序列配准过程，包括以下步骤：

①确定基准图像F(x,y)和待配准图像G(x,y)，在两幅图像上选择12对控制点f(x_i,y_i)和g(x_s,y_s)；

②将控制点值代入形变方程式，求解形变系数；

③利用形变关系求解配准结果F′(x,y)＝H^-1[G(x,y)]。

步骤5、载入血管内超声图像序列，利用光线投射算法，重建血管三维模型；

根据血管内超声图像的灰度直方图特性，为体数据场确定0,50,200,255四个阈值，将数据场分为Q₁＝[0,50]，Q₂＝[50,200]，Q₃＝[200,255]共三个阈值区间。

按下式将体数据场中各数据点的灰度值映射成为直接用于绘制的不透明度值

$A_{Q_q}\left(\mathrm{\&zeta;}\right)=\left\{\begin{array}{cc}0& 0\&le;\mathrm{\&zeta;}\&le;(c-\frac{w}{2})\\ \frac{\mathrm{\&zeta;}-c}{w\text{'}}-\frac{w}{2w\text{'}}& (c-\frac{w}{2})<\mathrm{\&zeta;}\&le;(c-\frac{w}{2}+w\text{'})\\ 1& (c-\frac{w}{2}+w\text{'})<\mathrm{\&zeta;}<(c+\frac{w}{2}+w\text{'})\\ \frac{c-\mathrm{\&zeta;}}{w\text{'}}+\frac{w}{2w\text{'}}& (c+\frac{w}{2}-w\text{'})<\mathrm{\&zeta;}\&le;(c+\frac{w}{2})\\ 0& (c+\frac{w}{2})<\mathrm{\&zeta;}\&le;255\end{array}\right.---\left(13\right)$

其中q为阈值区间Q₁或Q₂或Q₃的序号，ζ为数据点的灰度值，感兴趣物质的灰度值范围取决于数据中心c和宽度w两个变量，w′表示斜坡的宽度，斜坡表示的是线性的增加或降低，不透明度值用表示，取值范围在0到1之间。

按下式给体数据场的三个阈值区间Q₁＝[0,50]，Q₂＝[50,200]，Q₃＝[200,255]赋不同颜色值

$C_{Q_q}\left(\mathrm{\&zeta;}\right)=\left\{\begin{array}{cc}{C}_1& 0\&le;\mathrm{\&zeta;}\&le;50\\ C_2& 50<\mathrm{\&zeta;}\&le;200\\ C_3& 200<\mathrm{\&zeta;}\&le;255\end{array}\right.--\left(14\right)$

其中q为阈值区间Q₁，Q₂，Q₃的序号，C表示颜色值域，C₁，C₂，C₃是设定的颜色值。

为数据场建立X,Y,Z三维坐标轴，将数据场中单位体积的立方体看作一个体素。光源光线透过数据场，在屏幕上形成一个成像平面，从成像平面的每个像素点发出一条穿过数据场的光线，沿着光线选择设定的有限的K个等间距的采样点，采用下式计算各采样点的不透明度值：

其中A_e代表当前采样点的不透明度值，A₁，A₂，…，A₇，A₈代表距离采样点最近的8个数据点的不透明度值。

按下式计算采样点的颜色值：

$\begin{array}{c}{C}_l^{\&prime;}=C_1^{\&prime;}+x(1-y)(1-z)(C_2^{\&prime;}-C_1^{\&prime;})+x(1-y)z(C_3^{\&prime;}-C_4^{\&prime;})+(1-y)z(C_4^{\&prime;}-C_1^{\&prime;})\\ +y(C_5^{\&prime;}-C_1^{\&prime;})+\mathrm{xy}(1-z)A(C_6^{\&prime;}-C_5^{\&prime;})+\mathrm{xyz}(C_7^{\&prime;}-C_8^{\&prime;})+\mathrm{yz}(C_8^{\&prime;}-C_5^{\&prime;})\end{array},---\left(16\right)$

其中代表当前采样点的颜色值，

代表距离采样点最近的8个数据点的颜色值。

对光线上的采样点进行累加，直到不透明度值增加到1，结束累加，此时的颜色值就是成像平面上像素的最终颜色，即得到三维血管模型，公式如下：

$\begin{array}{c}{C}_{\mathrm{out}}{A}_{\mathrm{out}}=C_{\mathrm{in}}{A}_{\mathrm{in}}+C_{\mathrm{now}}{A}_{\mathrm{now}}(1-A_{\mathrm{in}})\\ A_{\mathrm{out}}=A_{\mathrm{in}}+A_{\mathrm{now}}(1-A_{\mathrm{in}})\end{array}---\left(17\right)$

其中C_out、A_out分别为经过第u个采样点后的颜色值、不透明度值，C_now、A_now为第u个采样点的颜色值、不透明度值，C_in、A_in为已合成的前u-1个采样点的颜色值、不透明度值，u＝1,2,…,u,…,K，K为采样点个数。

综上所述，血管三维重建的步骤总结如下：

①设定阈值，将体数据场归类为若干阈值区间；

②设计传递函数给数据点赋予不同的不透明度值和颜色值；

③计算光线上采样点的不透明度值和颜色值；

④完成图像合成，得到三维重建效果。

步骤6、实现对血管模型的任意角度平面剖切，获得形象、清晰的局部切面信息；

在由X,Y轴组成的水平面X-Y上，取相邻的4个坐标点：x，x+1，y，y+1构成正方形的4个顶点：(x,y)，(x,y+1)，(x+1,y)和((x+1),(y+1))，形成一个水平的剖切平面。

在水平面X-Y的高度Z方向，建立一个三维坐标空间，在有限个数的z坐标值上建立P个在Z轴方向上相互平行的空间剖切平面p，构成一个体素，每个剖切平面p与三维血管模型的4条棱线共有4个交点：z(x_p,y_p)，z(x_p,(y+1)_p)，z((x+1)_p,y_p)和z((x+1)_p,(y+1)_p)，p=1,2,…,P，从剖切平面p与三维血管模型棱线的交点中，任意选择4个能构成斜切平面的点，就可以实现任意斜面剖切。

每一个剖切平面p的中心点即为采样点O_p，坐标为

${\stackrel{\&OverBar;}{x}}_p=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}{x}_p,{\stackrel{\&OverBar;}{y}}_p=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}{y}_p,{\stackrel{\&OverBar;}{z}}_p=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}{z}_p---\left(18\right)$

每一个空间剖切平面p的各顶点z(x_p,y_p)，z(x_p,(y+1)_p)，z((x+1)_p,y_p)和

z((x+1)_p,(y+1)_p)到对应采样点

的距离d_p由下式得到：

${d_p}^2=(x_{p,f}-{\stackrel{\&OverBar;}{x}}_{p,O})+(y_{p,f}-{\stackrel{\&OverBar;}{y}}_{p,O})+(z_{p,f}-{\stackrel{\&OverBar;}{z}}_{p,O})---\left(19\right)$

其中f是体素的顶点，上下共8个，f＝0,1,2,…,7…，x_p,f是第p个空间剖切平面的第f个顶点的横坐标，

是第m个空间剖切平面中心点O的横坐标，其它类推。

按下式计算某一个剖切平面p的中心点O_p的灰度值

其中：

d_p为所述第p个空间剖切平面的8个顶点到中心点距离之和，

d_p,f为所述第p个空间剖切平面上某一个顶点f到中心点的距离，

h_f为入射光线在所述某个顶点f处发出的光线强度；

综上所述，利用切片重组，实现对血管三维模型的交互操作的步骤总结如下：

①求剖切平面与体素的交点；

②求出剖切平面中心点的坐标值；

③求出剖切平面中心点的灰度值。

最后应说明的是：以上实施例仅用以说明本发明而并非限制本发明所描述的技术方案；因此，尽管本说明书参照上述的各个实施例对本发明已进行了详细的说明，但是，本领域的普通技术人员应当理解，仍然可以对本发明进行修改或等同替换；而一切不脱离发明的精神和范围的技术方案及其改进，其均应涵盖在本发明的权利要求范围当中。

Claims (1)

1.一种基于血管内超声影像的血管三维可视化方法，其特征在于，是在连接着血管内超声仪的计算机中依次按以下步骤进行的：

步骤(1)利用血管内超声仪，以0.5mm/s的速度匀速回拉导管，获得人体冠状动脉的血管内超声视频影像；

步骤(2)将步骤(1)得到的血管内超声视频影像导入计算机，从视频中截取连续的900帧血管内超声图像作为实验图像，图像分辨率为384*384，以下简称为超声图像；

步骤(3)对上述超声图像依次按以下步骤进行超声图像序列降噪和超声图像噪声平滑：

步骤(3.1)取连续10帧超声图像组成一组序列，表示为I_k(x,y),1≤k≤10，k表示超声图像序号，对每帧超声图像进行[3,3]的中值滤波，得到中值滤波图像序列

$I_k^{\&prime;}(x,y)=\underset{\left[\mathrm{3,3}\right]}{\mathrm{med}}\left\{I_k(x,y)\right\},$ 符号med{·}表示中值滤波，

步骤(3.2)按下式求中值滤波图像序列的平均图像

${\stackrel{\&OverBar;}{I}}^{\&prime;}(x,y)=\frac{1}{10}\underset{k=1}{\overset{10}{\mathrm{\&Sigma;}}}\left[I_k^{\&prime;}(x,y)\right],1\&le;k\&le;10,$

步骤(3.3)对平均图像

按下式进行小波去噪，得到小波变换系数矩阵WT：

其中ψ^*(x,y)是小波函数ψ(x,y)的共轭， $\mathrm{\&psi;}(x,y)=a_0^{-j/2}\mathrm{\&psi;}(\frac{x-m{a}_0^j{a}_1}{a_0^j},\frac{y-n{a}_0^j{a}_2}{a_0^j}),$ 其中a₀，a₁，a₂是设定值，a₀＝a₁＝a₂＝2，称为扩展步长，1≤j≤3，是分辨索引，m≥1，n≥1，表示水平和垂直方向的有限平移，

步骤(3.4)对小波系数WT进行软阈值化： $W_{\mathrm{\&zeta;}}=\left\{\begin{array}{cc}\mathrm{sgn}\left(\mathrm{WT}\right)\left(\right|\mathrm{WT}|-\mathrm{\&zeta;}),& \left|\mathrm{WT}\right|\&GreaterEqual;\mathrm{\&zeta;}\\ 0,& \left|\mathrm{WT}\right|<\mathrm{\&zeta;}\end{array}\right.,$ 其中sgn()为符号函数，若WT＞0，sgn(WT)＝1，若WT＝0，sgn(WT)＝0，若WT＜0，sgn(WT)＝-1，ξ为阈值，

式中L为信号长度，σ为噪声强度，设图像噪声为高斯白噪声，且σ＝1，

步骤(3.5)按下式得到小波重构后的小波软阈值降噪图像

${\stackrel{\&OverBar;}{I}}^{\&prime;\&prime;}(x,y)=\underset{\&OverBar;\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\underset{\&OverBar;\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&zeta;}}\mathrm{\&psi;}(x,y),$

步骤(3.6)按以下步骤得到最终降噪图像 $I_k^{\&prime;\&prime;\&prime;}(x,y)=I_k^{\&prime;\&prime;}(x,y)+{\stackrel{\&OverBar;}{I}}^{\&prime;\&prime;}(x,y),1\&le;k\&le;10,$ 其中

为中值滤波图像

与平均图像

的差值：

步骤(4)按以下步骤实现图像序列的配准以补偿图像形变：

步骤(4.1)把步骤(3.6)得到的最终降噪图像

作为待配准图像G(x,y)，与预先设定的基准图像F(x,y)组成一个图像组合[F(x,y),G(x,y)]，从中选取12对控制点，分别记为f(x_i,y_i)和g(x_s,y_s)，1≤i≤12，1≤s≤12，符合如下关系：f(x_i,y_i)＝H^-1[g(x_s,y_s)]，其中H^-1为形变关系， $\left\{\begin{array}{c}{H_{x_i}}^{-1}=b_1{x}_s+b_2{y}_s+b_3{x_s}^2+b_4{y_s}^2+b_5{x}_s{y}_s+b_6\\ {H_{y_i}}^{-1}=c_1{x}_s+c_2{y}_s+c_3{x_s}^2+c_4{y_s}^2+c_5{x}_s{y}_s+c_6\end{array}\right.,$ b₁…b₆，c₁…c₂为形变系数，共12个，

步骤(4.2)根据选定的12个控制点f(x_i,y_i)和g(x_s,y_s)求12个形变系数，公式如下

$\left\{\begin{array}{c}{x}_i=b_1{x}_s+b_2{y}_s+b_3{x_s}^2+b_4{y_s}^2+b_5{x}_s{y}_s+b_6\\ y_i=c_1{x}_s+c_2{y}_s+c_3{x_s}^2+c_4{y_s}^2+c_5{x}_s{y}_s+c_6\end{array}\right.,$

步骤(4.3)按下式求得配准的超声图像F′(x,y)，F′(x,y)＝H^-1[G(x,y)]；

步骤(5)把步骤(4.2)得到的配准的超声图像F′(x,y)，按以下步骤利用光线投射算法，重建血管三维模型：

步骤(5.1)利用配准的超声图像构建体数据场，为体数据场设定0,50,200,255四个灰度阈值，将体数据场分为Q₁＝[0,50]，Q₂＝[50,200]，Q₃＝[200,255]共三个阈值区间，

步骤(5.2)按下式将体数据场中各数据点的灰度值映射成为直接用于绘制的不透明度值

$A_{Q_q}\left(\mathrm{\&zeta;}\right)=\left\{\begin{array}{cc}0& 0\&le;\mathrm{\&zeta;}\&le;(c-\frac{w}{2})\\ \frac{\mathrm{\&zeta;}-c}{w\text{'}}-\frac{w}{2w\text{'}}& (c-\frac{w}{2})<\mathrm{\&zeta;}\&le;(c-\frac{w}{2}+w\text{'})\\ 1& (c-\frac{w}{2}+w\text{'})<\mathrm{\&zeta;}<(c+\frac{w}{2}-w\text{'})\\ \frac{c-\mathrm{\&zeta;}}{w\text{'}}+\frac{w}{2w\text{'}}& (c-\frac{w}{2}-w\text{'})<\mathrm{\&zeta;}\&le;(c+\frac{w}{2})\\ 0& (c+\frac{w}{2})<\mathrm{\&zeta;}\&le;255\end{array}\right.,$ 其中q为阈值区间Q₁或Q₂或Q₃的序号，ζ为数据点的灰度值，感兴趣物质的灰度值范围取决于数据中心c和宽度w两个变量，w′表示斜坡的宽度，斜坡表示的是线性的增加或降低，不透明度值用

表示，取值范围在0到1之间，

步骤(5.3)按下式给体数据场的三个阈值区间Q₁＝[0,50]，Q₂＝[50,200]，Q₃＝[200,255]赋不同颜色值

$C_{Q_q}\left(\mathrm{\&zeta;}\right)=\left\{\begin{array}{cc}{C}_1& 0\&le;\mathrm{\&zeta;}\&le;50\\ C_2& 50<\mathrm{\&zeta;}\&le;200\\ C_3& 200<\mathrm{\&zeta;}\&le;255\end{array}\right.,$ 其中q为阈值区间Q₁，Q₂，Q₃的序号，C表示颜色值域，C₁，C₂，C₃是设定的颜色值，

步骤(5.4)为数据场建立X,Y,Z三维坐标轴，将数据场中单位体积的立方体看作一个体素，

步骤(5.5)光源光线透过数据场，在屏幕上形成一个成像平面，从成像平面的每个像素点发出一条穿过数据场的光线，沿着光线选择设定的有限的K个等间距的采样点，采用下式计算各采样点的不透明度值：

$\begin{array}{c}{A}_e=A_1+x(1-y)(1-z)(A_2-A_1)+x(1-y)z(A_3-A_4)+(1-y)z(A_4-A_1)\\ +y(A_5-A_1)+\mathrm{xy}(1-z)A(A_6-A_5)+\mathrm{xyz}(A_7-A_8)+\mathrm{yz}(A_8-A_5)\end{array},$ 其中A_e代表当前采样点的不透明度值，A₁，A₂，…，A₇，A₈代表距离采样点最近的8个数据点的不透明度值，

步骤(5.6)按下式计算采样点的颜色值：

$\begin{array}{c}{C}_l^{\&prime;}=C_1^{\&prime;}+x(1-y)(1-z)(C_2^{\&prime;}-C_1^{\&prime;})+x(1-y)z(C_3^{\&prime;}-C_4^{\&prime;})+(1-y)z(C_4^{\&prime;}-C_1^{\&prime;})\\ +y(C_5^{\&prime;}-C_1^{\&prime;})+\mathrm{xy}(1-z)A(C_6^{\&prime;}-C_5^{\&prime;})+\mathrm{xyz}(C_7^{\&prime;}-C_8^{\&prime;})+\mathrm{yz}(C_8^{\&prime;}-C_5^{\&prime;})\end{array},,$ 其中

代表当前采样点的颜色值，

代表距离采样点最近的8个数据点的颜色值，

步骤(5.7)对光线上的采样点进行累加，直到不透明度值增加到1，结束累加，此时的颜色值就是成像平面上像素的最终颜色，即得到三维血管模型，公式如下： $\begin{array}{c}{C}_{\mathrm{out}}{A}_{\mathrm{out}}=C_{\mathrm{in}}{A}_{\mathrm{in}}+C_{\mathrm{now}}{A}_{\mathrm{now}}(1-A_{\mathrm{in}})\\ A_{\mathrm{out}}=A_{\mathrm{in}}+A_{\mathrm{mow}}(1-A_{\mathrm{in}})\end{array},$ 其中C_out、A_out分别为经过第u个采样点后的颜色值、不透明度值，C_now、A_now为第u个采样点的颜色值、不透明度值，C_in、A_in为已合成的前u-1个采样点的颜色值、不透明度值，u＝1,2,…,u,…,K，K为采样点个数；

步骤(6)按以下步骤对所述三维血管模型进行任意方向平面裁剪，获得血管内部信息：

步骤(6.1)在由X,Y轴组成的水平面X-Y上，取相邻的4个坐标点：x，x+1，y，y+1构成正方形的4个顶点：(x,y)，(x,y+1)，(x+1,y)和((x+1),(y+1))，形成一个水平的剖切平面，

步骤(6.2)在水平面X-Y的高度Z方向，建立一个三维坐标空间，在有限个数的z坐标值上建立P个在Z轴方向上相互平行的空间剖切平面p，构成一个体素，每个剖切平面p与三维血管模型的4条棱线共有4个交点：z(x_p,y_p)，z(x_p,(y+1)_p)，z((x+1)_p,y_p)和z((x+1)_p,(y+1)_p)，p=1,2,…,P，从剖切平面p与三维血管模型棱线的交点中，任意选择4个能构成斜切平面的点，就可以实现任意斜面剖切，

步骤(6.3)每一个剖切平面p的中心点即为采样点O_p，坐标为

其中： ${\stackrel{\&OverBar;}{x}}_p=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}{x}_p,{\stackrel{\_}{y}}_p=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}{y}_p,{\stackrel{\&OverBar;}{z}}_p=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}{z}_p,$

步骤(6.4)每一个空间剖切平面p的各顶点z(x_p,y_p)，z(x_p,(y+1)_p)，z((x+1)_p,y_p)和z((x+1)_p,(y+1)_p)到对应采样点

的距离d_p由下式得到： ${d_p}^2=(x_{p,f}-{\stackrel{\&OverBar;}{x}}_{p,O})+(y_{p,f}-{\stackrel{\&OverBar;}{y}}_{p,O})+(z_{p,f}-{\stackrel{\&OverBar;}{z}}_{p,O}),$ 其中f是体素的顶点，上下共8个，f＝0,1,2,,7?，x_p,f是第p个空间剖切平面的第f个顶点的横坐标，是第p个空间剖切平面中心点O的横坐标，其它类推，

步骤(6.5)按下式计算某一个剖切平面p的中心点O_p的灰度值h_p：

其中：

d_p为所述第p个空间剖切平面的8个顶点到中心点距离之和，

d_p,f为所述第p个空间剖切平面上某一个顶点f到中心点的距离，

h_f为入射光线在所述某个顶点f处发出的光线强度。

Images