基于Hessian矩阵的DSA图像冠状动脉直径的测量

王力，孟庆民，郭永新，焦青

1.山东第一医科大学（山东省医学科学院）医学信息工程学院，山东泰安271016；2.泰安市中心医院介入放射科，山东泰安271000；3.山东第一医科大学（山东省医学科学院）放射学院，山东泰安271016

前言

心血管疾病具有高患病率、高致残率和高死亡率等特点［1］。据《2017～2022年中国医疗服务市场行情动态及发展前景预测报告》显示，我国心血管疾病患病人数接近3 亿，心血管疾病死亡率居首位，占居民疾病死亡构成的40%以上，而且我国心血管疾病患病率及死亡率处于上升阶段，心血管疾病的负担日渐加重，已成为重大的公共卫生问题，防治心血管疾病刻不容缓［2］。

数字减影血管造影技术（Digital Subtraction Angiography,DSA）是一种使X射线序列图片中的血管可视化的强大技术，是血管疾病无创诊断与介入治疗手术导航的重要依据，是心血管疾病诊断的金标准，广泛应用于X射线序列成像中血管的可视化系统中［3］。冠状动脉血管壁堆积较多脂质、糖类物质或动脉粥样硬化逐渐明显，易引起部分血管腔堵塞或狭窄。医生根据狭窄处血管直径的大小选择是否置入支架，及置入支架的尺寸［4］。置入过大的支架直径，冠状动脉血管可能会出现新破口或内漏，损伤冠状动脉血管管腔内壁，导致支架两端主动脉新的病变；置入过小的支架直径，会使支架的扩张应力不足，支架固定不稳固，易出现支架移位［5］。因此精确地判断冠状动脉血管直径，是确保支架置入手术的必要前提，这对于医生的临床决策具有重要的指导意义，对于冠状动脉搭桥手术的效果评价具有重要价值［6］。

该研究基于临床实践需要，在Windows7 环境下，以MATLAB为平台，通过对DSA图像进行去噪、均衡化等处理，并通过Hessian矩阵，提取冠状动脉血管的中心线和骨架，从而实现DSA 图像冠状动脉血管直径的测量。

1 DSA图像平滑

由于人为或仪器的各种原因，DSA 图像在获取和信号传输过程中会受到噪声的干扰，从而降低DSA 图像的质量，这对DSA 图像的处理与视觉效果均产生不利的影响。通常采用平滑技术对DSA图像进行处理，以消除图像的噪声［7］。

高斯滤波是一种平滑滤波器，取滤波器窗口内像素的均值作为输出，使用高斯滤波器对DSA 图像进行滤波，其效果是降低DSA 图像灰度的“尖锐”变化，能够有效抑制噪声，平滑图像［8］。对于DSA二维图像，用高斯函数做平滑滤波器，表达式如下：

式中，x、y是整个DSA 图像对应的矩阵坐标，H是模板的系数，表示每一个像素点的值，都由邻域内像素的加权平均灰度值代替。σ是DSA图像像素的标准差，σ取值较小时，高斯滤波器的频带较窄，信号边缘的定位精度较高，但平滑作用较小，对噪声的抑制能力较差（欠平滑）；反之，σ取值较大时，对信号的平滑作用较大，去噪效果较好，但会使信号边缘变得模糊，且边缘移位现象严重（过平滑）。通过调节σ值，可在过平滑与欠平滑间取得折衷［9］。

在MATLAB 平台中读入DSA 图像（图1），使用rgb2gray函数将图像转化为灰度图像，在此基础上使用fspecial函数与imfilter函数进行高斯滤波，选取较大的fspecial函数参数值（HSIZE参数设为1 000），得到DSA 图像的背景（图2），选取较小的fspecial 函数参数值（HSIZE 参数设为10），得到消除噪声后的DSA图像（图3）。

图1 原始DSA图像Fig.1 Original digital subtraction angiography image

图2 背景图像Fig.2 Background image

图3 消除噪声后的图像Fig.3 Image after noise removal

2 提取冠状动脉血管的中心线

Hessian矩阵是由图像在各个像素点的二阶偏导构成［10］，可表示为原图与高斯模板函数二次偏导数的卷积，反映各个像素点处的局部灰度几何信息［11］：

式中，x、y是平滑后灰度DSA 图像对应的矩阵坐标，H是曲面的曲率，I是平滑后的灰度图像，是高斯函数的二阶偏导数，Ixx、Ixy、Iyy分别为与I进行卷积，即灰度图像（图3）的二阶方向导数。

通过求取Hessian 矩阵的两个特征值和特征向量，可对图像中不同特征区域所具有的结构特性进行描述［12］。将Hessian 矩阵的特征值分解，得到两个正交的方向向量，分别用ν1、ν2 表示，对应的特征值用λ1、λ2 表示。特征值λ1 与其对应的特征向量ν1分别表示曲面曲率小的强度和方向，即ν1 和局部血管轴平行；特征值λ2 与特征向量ν2 分别表示曲面曲率大的强度和方向，即ν2 和局部血管轴垂直［13］。

在MATLAB 平台中，对于图3图像，使用struct函数构建多尺度立体数组，并采用sigma函数得到当前空间的值，在此基础上对Hessian矩阵使用atan2函数，获取不同方向三维曲面的曲率，进而获得血管的走向。然后采用imerode 函数、filter 函数和bwareaopen 函数对图像进行腐蚀，消除DSA 图像中较小的冠状动脉血管，最后使用bwmorph 函数得到冠状动脉血管的中心线（图4）。

图4 血管中心线图像Fig.4 Vessel centerline image

3 冠状动脉血管直径的测量

3.1 DSA图像增强

对经过高斯滤波的DSA图像进行直方图均衡化处理，该过程是采用灰度变换函数修正图像的直方图，使其成为灰度均匀分布的均衡直方图［14］，然后用此均衡直方图修正灰度图像，即对图像进行非线性拉伸，重新分配图像象元值，使一定灰度范围内象元值的数量大致相等，此时图像的信息熵最大，图像包含的信息量最大，图像对比度增强，层次分明，图像更加清晰［15］。对灰度图像进行直方图均衡化［16］的函数表达式为：

式中，nk为图像中出现该灰度的像素数（k=0,1,…,n-1），n是图像中的像素总数，是概率论的频数，si是灰度级的概率分布。

功能实现时，分别采用motion算子与unsharp算子来消除血管周围的失真晶格，及增大血管边缘的对比度。在MATLAB中分别采用histeq函数、fspecial函数及imfilter函数，对灰度图像（图5）进行均衡化处理，并增大冠状动脉血管边缘对比度（图6）。

图5 灰度DSA图像Fig.5 Grayscale digital subtraction angiography image

图6 图像均衡化Fig.6 Image equalization

3.2 提取冠状动脉血管的骨架

由于高斯滤波模板尺寸选取较大，因此滤波后图像接近于原图像的背景［17］（图2），经平滑与直方图均衡化处理后的图像，高斯滤波模板尺寸选取较小，起到去噪的作用［18］（图3），将这两步获得的后者矩阵减去前者矩阵，即可得到层次分明的血管图像（图7）。

图7 血管图像Fig.7 Vessel image

冠状动脉血管图像在物理意义上是有相同特点的区域集合，这些区域的形状或灰度或颜色大致相同、纹理属性相近等，通常采用边缘检测算子提取血管骨架。边缘检测算子包括Prewitt、Sobel、Roberts、Canny 和拉普拉斯算子等。Prewitt 算子边缘的定位准确度不够高、检测出的边缘较宽、间断点略多；Roberts算子检测出的边缘不是很平滑、边缘比较宽；Canny算子运算很复杂；拉普拉斯算子对噪声过于敏感［19］。而Sobel算子算法可对图像进行平滑去噪，边缘检测效果良好，且算法简单易实现［20］。该文使用Sobel 算子进行图像边缘检测及提取血管骨架。在MATLAB 中使用edge 函数对血管图像（图7）处理，进行边缘提取，即可得到冠状动脉血管的骨架（图8）。

图8 血管骨架图像Fig.8 Vessel skeleton image

3.3 测量冠状动脉血管的直径

由于冠状动脉血管中心线图像是以“0”、“1”数字为单位的矩阵，“0”代表图像背景，“1”代表冠状动脉血管的中心线。为了区分冠状动脉血管的骨架与中心线，遍历整个矩阵，把冠状动脉血管中心线矩阵中“1”的位置，全部用“2”替换。替换后的冠状动脉血管中心线图像矩阵与冠状动脉血管骨架图像矩阵相加，即可得到冠状动脉中心线血管骨架图像（图9）。采用Hessian 矩阵提取冠状动脉血管中心线时，得到Hessian 矩阵的两个特征向量（图10 中的v1、v2），v1 与局部血管轴的走向一致，v2 与局部血管轴垂直［21］。以v2 表示的方向作为斜率，过冠状动脉血管中心线上的点（图10中的a点），作一条直线（图10中的直线1），直线恰好与冠状动脉血管骨架有两个交点（图10中的b、c点），该交点间的距离即冠状动脉血管的直径。

图9 骨架中心线图像Fig.9 Skeleton centerline image

图10 血管直径测量图示Fig.10 Vessel diameter measurement

在MATLAB中，使用ginput函数获取光标位置，以此点为中心，使用两个for 循环，遍历矩阵，找到位于冠状动脉血管中心线上（矩阵中的值为“2”）且距离光标位置最近的点(x1,y1)。将此点与特征向量代入直线方程中，得到直线与冠状动脉血管骨架的两个交点坐标，使用dist函数求取这两个交点的距离，即冠状动脉血管的直径。

4 讨论

本研究通过对DSA 图像进行平滑、均衡化及锐化等处理，提取冠状动脉血管的骨架和中心线，从而实现DSA图像冠状动脉血管直径的测量。本研究具有以下特点：

图像处理效果好，适用范围广泛。在众多的图像平滑处理方法中，根据实际需求，本研究选取高斯滤波方法对图像进行处理，使得处理后的心血管图像适合人眼的视觉、消除图像的干扰并且保护血管边缘。为了突出血管边缘并进一步削弱图像中的杂质，采用直方图均衡化的方法，使得图像对比度增强，层次分明。为了使冠状动脉血管的边缘更为清晰，又使用滤镜消除血管周围失真的晶格。将处理后的图像减去高斯滤波后的背景图像，进一步弱化背景的干扰，使得提取的血管骨架更加明显。本研究采取的思路与方法，极大地弱化或消除图像中的杂质，突出血管边缘，适用于绝大多数医学图像的处理。

直径测量精确度高。针对Hessian 矩阵的特点，将冠状动脉血管的中心线看作线条边缘，冠状动脉血管边界看作阶跃边缘。根据线条边缘一阶导数值为零，二阶导数幅度最大，阶跃边缘一阶导数幅度最大，二阶导数值为零的微分几何特征，消除掉较小的血管，使得冠状动脉血管中心线精确度高。冠状动脉血管中心线上每个点的血管走向可由Hessian矩阵的特征向量得到，使用直线与冠状动脉血管骨架的交点测量血管直径，结果准确可靠，精度达到亚像素级，从而为操作人员提供精准的冠状动脉血管信息。

本研究的功能实现操作简便且易于扩展。用户只需读入DSA 图像，平台就会将处理后的DSA 图像以及冠状动脉血管中心线图像予以显示。医生只需用鼠标点击选择感兴趣的血管区域，平台即可自动计算并显示该位置的血管直径，极大地提高操作人员的工作效率。将算法代码在MATLAB上生成文件包，并与Java web 相融合，即可简单高效地运行在医生工作站。也可将冠状动脉血管的骨架、中心线图像以及直径等数据存放到数据库中，供科研人员进行分析处理。

本研究可供临床医师及科研人员对冠状动脉血管的DSA 图像进行回顾性分析与处理，也可为学生展示冠状动脉血管的DSA图像，易于操作，具有一定的实用价值。