[发明专利]一种基于血管内超声影像的血管三维可视化方法

摘要

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

主权项

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 ′ ( x , y ) = med [ 3,3 ] { I k ( x , y ) } , ]]>符号med{·}表示中值滤波，步骤(3.2)按下式求中值滤波图像序列的平均图像 I ‾ ′ ( x , y ) = 1 10 Σ k = 1 10 [ I k ′ ( x , y ) ] , 1 ≤ k ≤ 10 , ]]>步骤(3.3)对平均图像按下式进行小波去噪，得到小波变换系数矩阵WT：其中ψ^*(x,y)是小波函数ψ(x,y)的共轭， ψ ( x , y ) = a 0 - j / 2 ψ ( x - m a 0 j a 1 a 0 j , 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 ζ = sgn ( WT ) ( | WT | - ζ ) , | WT | ≥ ζ 0 , | WT | < ζ , ]]>其中sgn()为符号函数，若WT＞0，sgn(WT)＝1，若WT＝0，sgn(WT)＝0，若WT＜0，sgn(WT)＝-1，ξ为阈值，式中L为信号长度，σ为噪声强度，设图像噪声为高斯白噪声，且σ＝1，步骤(3.5)按下式得到小波重构后的小波软阈值降噪图像 I ‾ ′ ′ ( x , y ) = Σ ‾ ∞ ∞ Σ ‾ ∞ ∞ W ζ ψ ( x , y ) , ]]>步骤(3.6)按以下步骤得到最终降噪图像 I k ′ ′ ′ ( x , y ) = I k ′ ′ ( x , y ) + I ‾ ′ ′ ( x , y ) , 1 ≤ k ≤ 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为形变关系， 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 , ]]>b₁…b₆，c₁…c₂为形变系数，共12个，步骤(4.2)根据选定的12个控制点f(x_i,y_i)和g(x_s,y_s)求12个形变系数，公式如下 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 , ]]>步骤(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 ( ζ ) = 0 0 ≤ ζ ≤ ( c - w 2 ) ζ - c w ' - w 2 w ' ( c - w 2 ) < ζ ≤ ( c - w 2 + w ' ) 1 ( c - w 2 + w ' ) < ζ < ( c + w 2 - w ' ) c - ζ w ' + w 2 w ' ( c - w 2 - w ' ) < ζ ≤ ( c + w 2 ) 0 ( c + w 2 ) < ζ ≤ 255 , ]]>其中q为阈值区间Q₁或Q₂或Q₃的序号，ζ为数据点的灰度值，感兴趣物质的灰度值范围取决于数据中心c和宽度w两个变量，w′表示斜坡的宽度，斜坡表示的是线性的增加或降低，不透明度值用表示，取值范围在0到1之间，步骤(5.3)按下式给体数据场的三个阈值区间Q₁＝[0,50]，Q₂＝[50,200]，Q₃＝[200,255]赋不同颜色值 C Q q ( ζ ) = C 1 0 ≤ ζ ≤ 50 C 2 50 < ζ ≤ 200 C 3 200 < ζ ≤ 255 , ]]>其中q为阈值区间Q₁，Q₂，Q₃的序号，C表示颜色值域，C₁，C₂，C₃是设定的颜色值，步骤(5.4)为数据场建立X,Y,Z三维坐标轴，将数据场中单位体积的立方体看作一个体素，步骤(5.5)光源光线透过数据场，在屏幕上形成一个成像平面，从成像平面的每个像素点发出一条穿过数据场的光线，沿着光线选择设定的有限的K个等间距的采样点，采用下式计算各采样点的不透明度值： 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 ) + xy ( 1 - z ) A ( A 6 - A 5 ) + xyz ( A 7 - A 8 ) + yz ( A 8 - A 5 ) , ]]>其中A_e代表当前采样点的不透明度值，A₁，A₂，…，A₇，A₈代表距离采样点最近的8个数据点的不透明度值，步骤(5.6)按下式计算采样点的颜色值： C l ′ = C 1 ′ + x ( 1 - y ) ( 1 - z ) ( C 2 ′ - C 1 ′ ) + x ( 1 - y ) z ( C 3 ′ - C 4 ′ ) + ( 1 - y ) z ( C 4 ′ - C 1 ′ ) + y ( C 5 ′ - C 1 ′ ) + xy ( 1 - z ) A ( C 6 ′ - C 5 ′ ) + xyz ( C 7 ′ - C 8 ′ ) + yz ( C 8 ′ - C 5 ′ ) , ]]>其中代表当前采样点的颜色值，代表距离采样点最近的8个数据点的颜色值，步骤(5.7)对光线上的采样点进行累加，直到不透明度值增加到1，结束累加，此时的颜色值就是成像平面上像素的最终颜色，即得到三维血管模型，公式如下： C out A out = C in A in + C now A now ( 1 - A in ) A out = A in + A mow ( 1 - A in ) , ]]>其中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，坐标为其中： x ‾ p = 1 P Σ p = 1 P x p , y _ p = 1 P Σ p = 1 P y p , z ‾ p = 1 P Σ p = 1 P 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 - x ‾ p , O ) + ( y p , f - y ‾ p , O ) + ( z p , f - 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处发出的光线强度。