CN113838153A - 一种基于血管造影图像的血管管径分析系统 - Google Patents

Abstract

本发明公开了一种基于血管造影图像的血管管径分析系统，属于计算机视觉及机器学习领域。系统由四个模块组成。模块1：利用血管造影图像自动分割血管区域；模块2：依据血管分割结果进行距离变换分析，得到血管内像素距离血管边界的距离图像；模块3：点选局部血管段起始点和终止点，设计路径检测优化目标估计出血管内穿行的内切球的运动轨迹；模块4：获取运动内切球的直径曲线作为血管管径曲线。本发明的优势在于解释性强、计算简单，避免了传统2D图像重建策略的多步骤误差累计问题，也规避了传统分析方法不合理的管腔形态学假设，更适用于血管管径的自动分析。

Description

一种基于血管造影图像的血管管径分析系统

技术领域

本发明涉及计算机视觉及机器学习领域，更具体的，涉及一种基于血管造影图像的血管管径自动分析系统。

背景技术

动脉造影成像是一种血管成像方法，该成像方法能够高亮地凸显动脉血管。基于血管造影图像分析血管的管径曲线能够直观可视化地展示血管的管径变化规律，是一种图像信息可视化技术。

前人开发的主流血管管径分析系统采用的是基于2D图像重建的多步分析策略。这一策略通常包括以下步骤：血管路径检测、2D图像重建、基于2D图像的血管管径定量。传统2D图像重建策略的基本思想是模拟人工在2D图像上测量管径的过程。在实施过程中，这类系统首先进行血管路径的检测，再基于路径轨迹进行2D图像重建，生成展示血管纵截面的2D曲面重建图或展示血管横截面的横截面重建图，随后基于纵截面重建图或横截面重建图分割血管区域或检测血管边界(也有部分工作在血管路径检测步骤前进行血管分割)，最终进行血管管径的自动定量。在管径定量上，部分工作遵循了最短内径定量原则；也有部分工作定量的是血管横截面平均管径；还有部分工作在假设血管截面为正圆形条件下，根据圆形物体的面积与半径的转换公式粗略估计管径数值。目前最受领域内认同的是最短内径定量原则，即选择最短内径作为特定位置的管径数值。

目前主流的基于2D图像重建策略的管径分析系统主要有以下几点不足：1、血管2D截面的自动重建极其依赖于血管路径的精确检索，路径走形的丝毫偏差会导致截面重建的错误，从而影响管径定量准确性，导致较差的鲁棒性。正是由于重建图像对血管轨迹定量准确性的高度依赖，前人的方法常涉及非常复杂繁琐的血管中心线矫正的步骤。2、在基于2D截面图像的血管管径测量步骤中，通常会有一些先验的假设(如假设血管横截面形态为椭圆形)，但是这些假设在病变的血管位置往往并不成立，这类先验假设也会引入测量的误差。例如在发生斑块的位置，斑块的形状通常较为多变，此时血管腔常常不满足圆形等先验假设。3、管径定量依赖于有限截面的选择，选择不同的截面图像重建的位置或角度，会导致测量出不同的管径数值。4、不少前人研究采用面积、平均管径长度、或者根据截面积估计的管径数值生成管径曲线，仅有部分工作采用业内常用的最短内径标准生成管径曲线。由于整体步骤繁多，并且上述部分子步骤对前处理步骤的精度要求很高，误差的累计效应导致这类模拟人工测量的2D图像重建策略的精确度非常受限。

前人采用的另一类常见的狭窄分析系统采用的是：血管追踪、3D柱形拟合、管径定量的分步策略。这类柱形拟合策略通过血管追踪得到血管的路径，然后在血管形态服从3D柱形先验假设情况下进行管径定量。目前上述两类策略的研究工作基本上都需要一定的人工交互干预，例如血管起始点和终止点的手工点选。

柱形拟合策略的主要局限性在于：1、几何学模型建模过程繁琐，适应不同血管结构变异的能力较差，例如对血管分叉位置和斑块位置需要额外的建模过程；几何学模型的参数拟合计算过程也较为复杂。2、基于先验几何学假设(如柱形假设)常难以有效描述结构复杂的血管形态，这也会影响定量精度。这里值得指出的是，虽然前人研究习惯用血管中心线一词来描述血管路径，但是这些方法都未能从建模的角度严格地论证算法所寻找的路径真正位于局部血管的中心。

因此，目前亟需一种不依赖于2D图像重建以及血管形态假设进行管径定量方法。

发明内容

为了克服现有技术的不足，避免血管2D截面重建所带来的误差累计问题以及不合理的血管形态先验假设，本发明提出一种运动球囊法的血管造影图像分析方法，用于自动定量血管管径。本发明另辟蹊径，使用血管内最大内切球直径估计血管内径。可以理解，最大内切球的直径恰好对应局部血管的最短内径，符合最短内径定量标准。此外，所述方法在估计管径时不依赖血管形态的先验假设，不必对血管的边界形态进行约束。本发明进一步通过模拟一个运动的可变直径的最大内切球穿行血管的过程，将内切球的直径变化曲线视作是管径曲线，这一策略也能够避免2D图像重建步骤，规避血管2D截面重建所引入的多步骤累计定量误差。这些优势使得本发明提出的新策略更适合自动化地定量血管管径。

本发明提出的管径分析系统主要由以下模块组成：血管分割模块M1、基于血管分割结果的距离变换变换模块M2、基于路径检测模块M3求解运动球囊的球心轨迹、以及进一步的管径曲线定量模块M4。四个模块的整体流程图示如图1所示，下面进一步描述以上四个模块的具体构成：

血管分割模块M1首先估计出的血管亮度阈值范围[T_L，T_U]，然后再使用阈值范围[T_L，T_U]对血管造影图像I进行分析得到二值化的血管分割结果V。

阈值范围[T_L，T_U]以下T11、T12、T13、T14步骤分析得到：

T11：手工点选血管内种子点记为u(可选用下述模块M3中的血管段起始点s或终止点e代替)；并记录u点的信号值大小，记为Q。

T12：以u为中心，以R为半径，计算局部区域内信号亮度的统计直方图H_R；R的取值范围为1，2，3，…，L。

T13：计算加权统计直方图

记为H；并对H进行高斯滤波，得到滤波后的统计直方图

T14：在统计直方图

中检测最靠近信号值Q的信号峰，并记录峰值位置M_c和计算半峰宽M_w；选择[M_c-τM_w，M_c+τM_w]为自适应的血管亮度阈值范围[T_L，T_U]，τ为人工设定的常数。

在根据阈值范围[T_L，T_U]分割血管区域时，获取包含种子点u的连通区域∏，使得连通区域∏内所有像素的亮度均位于[T_L，T_U]范围内，而∏局部领域以内、∏以外的所有其他像素的亮度均位于[T_L，T_U]范围之外；将血管的连通区域∏标记为一，背景区域标记为零，得到二值化的血管分割结果V。

距离变换模块M2根据二值化的血管分割结果V和血管造影图像I进行二值图像的二值距离变换分析或者灰度图像的模糊距离变换分析，得到血管内像素距离血管边界的距离图像，记为图像D；图像D中血管区域内每个像素点的数值大小表示该像素点距离血管边界的最短距离，图像D中血管区域外每个像素点的数值为零。如采用二值距离变换分析，直接对二值化的血管分割结果V进行二值距离变换可得到距离图像D。如采用模糊距离变换分析，首先根据血管造影图像I和血管分割结果V定义归属度函数μ，非血管区域的归属度函数值为零，血管区域内每个像素点a的归属度数值μ(a)表示人工定义的属于血管区域的先验概率；然后定义血管内任意像素点b到血管边界的路径集合为ξ_b，最后计算μ关于ξ_b中所有路径的路径积分的最小值作为模糊距离D(b)，得到距离图像D。

路径检测模块M3由以下步骤组成：

S31：人工点选局部待分析血管段在图像中的起始点s和终止点e。

S32：连接s点和e点的轨迹可以有无穷多种可能性，即穿过血管内的从s点运动到e点的可变直径球体的运动轨迹有无穷多种可能性。为了保证球体的直径曲线恰好对应管径的变化曲线，需要要求球心轨迹p中每个点都经过血管的正中心，即让球体的半径刚好对应局部截面的最大内切球半径，则要保证上述运动球体球心在血管微元长度内的取值为极大值，即球心轨迹p中的每一个点在D图像中血管局部横截面内取值最大。此外，连接起始点s和终止点e的血管轨迹必定只经过每个微元截面一次，即轨迹经过微元截面的次数加和取值最小。上述的微元法的轨迹检测问题事实上有两点限制：1、球心轨迹p经过所有微元的次数和最小。2、对于图像D，球心轨迹p中每个像素点在相应血管微元截面中取值最大(最大值约束)。如果构造一个单调递减的一元函数F，则第2条的最大值约束可以变成最小值约束，此时轨迹检测问题可以转化为如下的优化目标：

利用上述优化目标求取球心轨迹p及对应的路径长度K；公式中i，K为整数；p(i)表示路径上第i个像素点的空间坐标；函数D(p(i))表示在坐标p(i)处D的取值；s.t.表示约束条件；F为在定义域范围内单调递减的一元函数。将球心轨迹作为血管路径，将球体直径作为管径数值，管径曲线定量问题转化为根据上述的优化目标，求取连接s点和e点的最佳轨迹曲线问题。

S33：为了对S32中的优化目标进行求解，本步骤进行等价问题转换将S32中的优化问题转化为图上的最短路径检索问题，具体实施过程为：依据图像D构造图空间G，根据如下公式赋予图G中像素点的连接边权重：

w_ab＝w(G_a，G_b)＝γ(F(D_a)+F(D_b))

公式中w_ab表示连接空间任意像素位置点a和邻域内像素位置点b的连接边的权重，γ为预设常数。对于输入的血管段起始点s和终止点e，在图空间G中计算连接s和e的最短路径作为球心轨迹p。

管径曲线定量模块M4为：获取运动球体球心轨迹p之后，根据公式

{C＝(C₁，C₂，C₃，...，C_K)|C_i＝D(p(i))，i＝1，2，...，K}

得到血管管径曲线C。

本发明的有益效果为：一种基于血管造影的血管管径分析系统，通过运动球囊法的管径定量策略，在不依赖于管腔形态学假设的情况下，快速分析血管管径曲线；新策略同样避免了主流的2D图像重建策略的多步骤误差累计局限性，具有解释性强、计算简单的优点。新策略仅需用户指定待分析血管段的起始点和终止点，便能实现全自动的血管管径曲线定量。

附图说明

图1：本发明整体流程图示。

图2：基于血管局部区域统计直方图分析确定血管亮度阈值范围；(A)手工点选种子点局部邻域内的统计直方图；(B)加权统计直方图以及高斯滤波后的加权统计直方图；(C)基于加权统计直方图确定血管亮度阈值范围。

图3：运动球囊法管径定量结果图，(A)(B)(C)(D)(E)图示中箭头指示特定位置处运动球囊的形态和大小。

图4：根据球心轨迹重建的2D投影图和管径曲线；(A)矢状位投影图；(B)冠状位投影图；(C与D)血管段3D重建图；(E1)矢状位投影图局部，图中虚线表示球心轨迹；(E2)冠状位投影图局部，图中虚线表示球心轨迹；(E3)血管段管径曲线定量结果。

具体实施方式

以下通过实施例对本发明中部分步骤的具体实施过程做进一步说明，以便更好地理解本发明的技术方案。

本实施例以三维CTA(Computed Tomography Angiography)影像中的颈动脉分叉段血管为分析对象，获取数例头颈CTA图像数据。成像分辨率为：轴位图像中像素距离的范围为0.461至0.625mm，层厚为0.5mm或0.8mm；轴位图像的图像分辨率为512×512，层数范围为300至700层。不同的血管造影成像类型、扫描机器、成像参数和收集图像数量、选取的待分析血管段不会影响本发明实施的有效性，只会在一定程度上影响本发明专利实施的效果。

在血管分割步骤之前，实施例首先通过三次插值将原始三维图像插值为XYZ轴等分辨率的图像，称之为XYZ轴等分辨插值。进行图像插值放大的一个原因是帮助后续距离变换步骤获得更加准确的距离变换数值。另一个原因是分析方法中的球心轨迹和半径曲线均以像素为基本单位，所以在图像插值放大的过程中，需要保证放大后的图像的XYZ轴是等分辨率的。假设原始图像中XYZ轴单位像素物理长度分别为r_x、r_y和r_z(单位：mm/像素)，假设对图像做N倍放大，则插值后图像的不同轴的单位像素物理长度应同为min{r_x，r_y，r_z}/N。本实施例中N的设置为3。下面详细描述发明内容中各模块的具体实施过程。

血管分割模块M1：本发明使用区域生长算法分割血管区域，通过手工点选局部血管段中的起始点作为种子点u，根据种子点u进行区域生长，区域生长方法通过不断地将邻域内像素亮度位于血管亮度阈值范围[T_L，T_U]内的像素添加入分割区域得到最终的血管分割结果。在选取区域生长自适应阈值范围[T_L，T_U]的T11、T12、T13、T14步骤中，L的取值为20，α_R取值为1/20；高斯滤波的卷积核

中的参数σ取值为8；τ取值为2，即选择距离峰值中心2倍半峰宽的信号范围作为估计出的血管亮度阈值范围。图2展示了典型的基于局部统计直方图分析估计血管亮度阈值范围的过程与结果。图2中(A)子图展示了选择不同半径R时，邻域内像素信号的统计直方图，图2中(B)展示了加权统计直方图，图2中(C)子图展示了自适应选择的血管亮度阈值范围。

在加权统计直方图中，血管的信号分布与背景的分布区分得极其明显，例如子图(B)中的加权统计直方图清晰展示了分布的两个峰值，低CT值的峰展示了背景像素亮度分布，高CT值的峰凸显了血管亮度分布。而在未加权的统计直方图中，血管与背景组织的CT值分布的区分并不明显。例如在子图(A)中，当R＝1时，以种子点为中心的1邻域内都是血管区域，CT值的分布都是血管信号的范围；当R＝5或R＝10时，局部邻域范围内同时包括血管和非血管像素，但是二者信号峰区分不明显；当R＝20时，分布中展示的0至200HU范围内的信号峰在非血管背景组织的信号范围内，由于血管所占体积较小，血管CT值范围内未有明显信号峰的出现。但是在子图(B)中展示的加权统计直方图中，对应的血管和背景组织的信号峰区分非常明显。

距离变换模块M2：使用数字图像处理中的二值距离变换方法对二值化的血管分割结果V进行分析得到距离图像D，二值距离变换中的距离选择为欧式距离。在距离图像D中，血管边界像素的取值为1，越靠近血管中心的像素数值越高，数值的含义表示像素距离血管边界的最短距离。

路径检测模块M3中各步骤的实施过程：

S31：首先人工点选颈动脉分叉局部血管段在图像中的起始点s₀和终止点e₀，人工点选的s₀与e₀均位于血管内。由于人工点选的起始和终止点常会偏离血管的正中心位置，本文根据模块M1中估计出的血管亮度阈值范围[T_L，T_U]对s₀与e₀所在的轴状位图像进行二值化和距离变换，并选择最靠近起始点s₀的局部最大值点替代手工点选的起始点，矫正之后的起始点记为s；选择最靠近终止点e₀的局部最大值点替代手工点选的终止点，矫正之后的终止点记为e。

S32：在计算连接s点和e点的运动内切球的球心轨迹时，单调递减的一元函数F设置为F(x)＝e^-x，此时球心轨迹检测问题可以转化为如下的优化目标：

s.t.p(1)＝s；p(K)＝e；

根据优化目标求取连接s点和e点的球心轨迹p及对应的路径长度K。

S33：对S32中的优化目标进行求解时，本步骤根据图像D构造图空间G，将像素点数视为图节点。在图G的连接边构造过程中，将所有单个像素点与三维图像局部26-邻域内所有像素点相连，任意位置点a和邻域内位置点b的连接边权重定义为：

w_ab＝w(G_a，G_b)＝γ(F(D_a)+F(D_b))

常数γ设置为0.5。本实施例利用Dijkstra算法计算出连接起始点s和终止点e的最短路径作为球心轨迹p。

根据上述实施例描述的分析步骤对一例典型颈动脉分叉段血管进行分析，得到的管径曲线如图3所示。图示左侧一栏展示了颈动脉分叉段血管的管径曲线定量结果，图示右侧(A)(B)(C)(D)(E)子图展示了管径曲线中特定位置对应的运动球囊在3D血管中的形态。

进一步，根据上述实施例参数设定和方法实施过程对另一例颈动脉分叉段血管进行分析，根据管径曲线进行曲面重建得到图4中的轴位和冠状位投影图像。可以观察到在轴位和冠状位投影图像中，血管边界清晰且连续，并且血管关于球心轨迹左右对称良好，这说明本文方法定量出的球心轨迹位于血管的中心，能够有效地代表血管中心线。

虽然本发明已以较佳实施例披露如上，然而并非用以限定本发明。任何熟悉本领域的技术人员，在不脱离本发明技术方案范围情况下，都可利用上述揭示的方法和技术内容对本发明技术方案作出许多可能的变动和修饰，或修改为等同变化的等效实施例。因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均仍属于本发明技术方案保护的范围内。

Claims (4)

1.一种基于血管造影图像的血管管径分析系统，其特征在于，该系统由以下模块组成：

血管分割模块M1，对采集的血管造影图像I进行分析得到二值化的血管分割结果V；

距离变换模块M2，通过对V和I进行分析得到血管内任意像素距离血管边界的最短距离，并生成距离图像D，对于在血管区域内任意像素，距离图像D记录相应的像素距离边界的最短距离数值，血管区域以外像素在距离图像D中取值为零；

路径检测模块M3，根据输入的血管段起始点s和终止点e，通过以下管径定量优化目标

求取血管路径p及对应的路径长度K；公式中i，K为整数，p(i)表示路径上第i个像素点所对应的空间坐标，D(p(i))表示在坐标p(i)处D的取值，s.t.表示约束条件，F为在定义域范围内单调递减的一元函数；

管径曲线定量模块M4，根据公式{C＝(C₁，C₂，C₃，...，C_K)|C_i＝D(p(i))，i＝1，2，...，K}得到血管管径曲线C。

2.如权利要求1所述的路径检测模块M3，其特征在于，路径检测模块中的管径定量优化目标通过如下步骤转换为最短路径检索问题：

S331：根据图像D构造图空间G，将像素点作为图节点，并根据如下公式赋予图空间G中像素节点的连接边权重：

w_ab＝w(G_a，G_b)＝γ(F(D_a)+F(D_b))

公式中w_ab表示连接空间任意像素节点a和邻域内像素节点b的连接边的权重，γ为预设常数；

S322：对于输入的血管段起始点s和终止点e，在图空间G中检索连接起始点s和终止点e的最短路径作为血管路径p，并将最短路径的长度作为K。

3.如权利要求1所述的距离变换模块M2，其特征在于，通过对图像V进行二值距离变换或者根据V和I进行模糊距离变换，得到距离图像D。

4.如权利要求1所述的血管分割模块M1，其特征在于，通过对人工输入的血管内种子点u的局部邻域的灰度统计分布直方图进行分析得到血管亮度阈值范围[T_L，T_V]，根据血管亮度阈值范围[T_L，T_U]提取血管造影图像I中的血管区域得到二值化的血管分割结果V。

Images