在体光声微血管图像提取方法研究

张智超　赵景秀　孟静

摘要摘要：光學分辨率光声显微镜能够进行活体、3D、高分辨率成像，在生物微循环研究领域具有重要的应用价值。组织微血管的变化能够反映很多疾病的发生和发展过程，对微血管图像进行分析可以提取与疾病相关的各种参数，实现疾病早期诊断和治疗。然而，微血管图像由于血管细小，光吸收较弱，直接基于光声显微成像系统得到的图像在微血管密集或较深区域比较模糊，不能很好地提取血管脉络信息。因此，提出一种融合了Hessian特征增强、血管分割和骨架提取的微血管图像处理方法，并首次用于光声肿瘤微血管图像的处理和分析。活体小鼠耳部正常血管及黑色素瘤血管图像的实验表明，所提出的微血管图像处理方法能够有效分割出活体组织的微血管结构，最终实现对光声肿瘤微血管及其骨架的提取。

关键词关键词：光声成像；肿瘤血管；图像增强；Hessian矩阵；血管分割；骨架提取

DOIDOI：10.11907/rjdk.171495

中图分类号：TP317.4

文献标识码：A文章编号文章编号：16727800（2017）005017904

0引言

血管新生是在原来存在的血管网络上长出新的血管的一系列复杂生物学过程，而肿瘤是血管新生的依赖性疾病，肿瘤的生长、侵袭与转移均有赖于血管新生。所以，肿瘤微血管的早期诊断和评估对肿瘤疾病早期治疗具有非常重要的意义。目前，在X射线、计算机断层扫描、核磁共振和超声等成像技术领域，已经有研究人员对血管图像作了量化分析和风险评估[13]。但由于这些成像方式不能分辨直径在100μm以下的新生肿瘤血管，所以只能提供普通血管的结构特征和功能成像，很难为肿瘤的早期诊断、生长和扩散提供有效信息。

光学分辨率光声显微成像（optical-resolution photoacoustic microscopy， ORPAM）是基于光声效应的一种成像手段，兼具光学成像高分辨率和超声成像深度的优势，利用血红蛋白的内源性对比，ORPAM成像分辨率可以达到微米/亚微米级，从而实现对单根毛细血管进行成像[4]。因此，ORPAM正成为研究生理和病理相关组织微循环的有力工具[5]。但ORPAM对于较深层或非常细小的血管，由于光的散射和组织光吸收少等原因，导致成像信噪比不高，某些微弱、密集信号区域的信号模糊，血管成像不连续。因此，需要通过图像后续处理技术对这些血管进行增强，提高血管脉络的清晰度和对比度。

在传统成像技术领域，如X射线成像、光学相干断层扫描技术（OCT）[67]、核磁共振成像[8]、CT成像等，许多科研人员开展了对血管图像的处理和分析工作。最近，针对ORPAM的光声血管图像，已有相关人员开展研究[9]。然而，这些研究工作都是在普通血管上进行的，并没有对血管更为精细的肿瘤微血管进行分析。因此，针对ORPAM系统微血管图像特征，本文提出了一种基于高频强调滤波、多尺度Hessian矩阵增强和基于梯度场均衡化增强的图像处理方案，实现了活体光声微血管图像的增强和分割。为验证本文提出的方法，对小鼠耳部黑色素瘤进行活体ORPAM成像，然后对黑色素瘤及周边组织的血管脉络图像进行处理。实验结果证明，本文方法可以有效增强光声微血管图像，高精度提取普通和肿瘤微血管。

1方法

为实现ORPAM系统获取的光声血管图像的高质量分析，本文设计光声微血管图像增强和提取算法流程，如图1所示。图像处理流程如下：①用高频强调滤波增强微血管的对比度；②基于Hessian矩阵的多尺度血管增强检测血管边缘。由于矩阵的特征值和特征向量能够表示血管的强度和方向，因此可以对图像中的管状结构进行判定；③基于梯度场均衡化的图像对比度增强[10]，对Hessian图像在梯度场进行均衡化操作，使图像中较亮的区域在梯度场进行增强，然后对增强后的梯度场进行重建；④用八元数矢量积改进三维区域生长算法得到最后的血管分割图像，用AFMM方法提取血管中心线。

1.1血管增强

1.1.1高频强调滤波

本文用高频强调滤波（High-frequency emphasis filter，HFE）增强微血管的对比度，该滤波器的定义如下：

Hhfe=a+b×Hhp（u，v）（1）其中，a是相对于原点的偏移量，b是乘数，Hhp（u，v）是高通滤波器的传递函数。高频强调滤波把高通滤波器与一个大于1的常数b相乘，再加上一个偏移量a。突出高频部分的同时增加了低频部分的幅度，但是在a和b值都比较小的情况下，低频增强的影响就会低于高频增强的影响。

1.1.2基于多尺度Hessian矩阵的血管增强

高频强调滤波增强了血管的对比度，但是血管的边缘仍然不清楚，血管方向模糊。在Hessian方法中，垂直于图像特征的方向是具有最大模特征向量的方向，平行于图像特征的方向是与它垂直的方向。因此，本文将Hessian矩阵应用到血管检测中，通过设计线状增强滤波函数，将图像中的噪声去除。此外，用Hessian矩阵的特征值判断图像中密度变化剧烈的点，以此检测血管边缘，而微血管的强度和方向可由Hessian矩阵的特征值和特征向量表示。假设二维图像用I表示，为实现图像I的Hessian矩阵增强，首先构造图像的Hessian矩阵：

H=IxxIxyIyxIyy（2）

其中，Ixy是图像I的二阶偏导数，可用如下公式计算：

Ixy（x，y，σ）=σ2gxy（x，y，σ）I（x，y）（3）σ为尺度，gxy为多尺度高斯核函数的二阶偏导数。Ixx，Iyx，Iyy分别可用类似公式（3）的方式计算得到。然后，计算公式（2）中Hessian矩阵的两个实特征值λ1和λ2（|λ1|>|λ2|），在尺度σ下，λ1对应的特征向量代表最大曲率方向，λ2对应的特征向量代表最小曲率方向。在此基础上，计算血管的相似性判别函数[11]，公式如下：F（x，y，σ）=0 if λ2>0，exp-M22α21-exp-N2H2β2 （4）

其中，M=|λ2/λ1|，NH=λ21+λ22，α和β是用于调节线性滤波器对测度M和NH的灵敏度阈值，α取固定值0.5，β的取值依赖于图像的灰度值范围，一般取Hessian矩阵范数最大值的1/2，F（x，y，σ）∈[0，1]；其中，F（x，y，σ）的数值越大，则表示此处与血管有越高的相似度。根据公式（4），可以得到图像的特征图：

F（x，y，σ）=maxσmin<σ<σmaxF（x，y，σ）（5）

1.1.3基于梯度场均衡化的图像对比度增强

Hessian方法实现了信号的连接，展现出了血管的脉络，但是图像中血管与周围组织的对比度比较低，影响血管的提取。因此，在此基础上，进行基于梯度场均衡化处理[10]，提高血管的亮度和清晰度。具体步骤如下：

对Hessian图像的梯度模‖F（x，y，σ）‖进行直方图均衡化L‖F（x，y，σ）‖，得到增强后的梯度场：

G=F（x，y，σ）‖F（x，y，σ）‖·L‖F（x，y，σ）‖（6）通过最小二乘法得到目标函数：

E（f）=（‖f-G‖2）dxdy（7）其中，f表示增强结果。

式（7）的Euler-Lagrange方程为：

2f=div（G）（8）

其中，2是Laplacian算子，2f=2f/x2+2f/y2。為求解公式（8），本文使用标准的有限差分近似产生它们的线性方程系统。而对于有限差分法产生的线性方程组的求解过程，详见文献[12]。

1.2血管分割

1.2.1获取血管二值图像

本文利用八元数矢量积方法[13]获取血管二值图像。首先使用传统的区域生长算法，分割出连续的血管，通过梯度计算，找到血管的边缘，将其作为八元数区域生长的初始种子点，然后用种子点6个邻域的灰度构造八元数来表示种子点的特征，运用八元数的矢量积运算，进一步对噪声影响下的断裂血管进行有效连接。流程如图2所示。

1.2.2提取中心线

在提取二值图像的基础上，用一种增广的快速行进法（Augmented fast marching method，AFMM）提取图像中轴，测定血管的中心线。AFMM算法流程如下：

（1）初始化DT，DT为血管内的点p到血管边界最近点q的距离，初始化U作为边界。

（2）用快速行进法（Fast marching method， FMM）计算距离值 DT，DT（p）=minq∈U（dist（p，q）），dist（p，q）=‖p-q‖2。

（3）修正FMM计算八邻域点的计算偏差， DT0（p）={DT（p），{q}}， q=argminq∈U（dist（p，q））。

（4）任意两点p和q之间的距离设为Dp（q）。

（5）对所有边界点q和每一个血管内的点p，如果Dp（q）

2实验与讨论

实验数据基于ORPAM系统获取的小鼠耳部黑色素瘤及其周边正常血管图像。ORPAM系统采用聚焦激光束和聚焦的单个超声换能器来获取组织深度方向的一维光声图像，通过移动激光束与超声换能器在生物组织上进行单行扫描，就可以得到生物组织的二维断面图像，三维的ORPAM图像可以通过多行扫描得到。

2.1小鼠耳部正常血管图像提取

首先，基于小鼠耳部正常血管图像测试本文提出的血管增强和分割方法的有效性。用到的图像是一个500×500×50大小的三维数据，将其在深度方向构建500×500的三维最大振幅投影图像（projection of maximum amplitude，MAP）作为耳部正常微血管图像处理的原图，如图3（a）所示。从图3（a）中可以看到：尺寸较大的血管由于直径宽，光吸收好，血管的分布情况比较清晰。但是，一些位置较深、更为细小的微血管信号则比较微弱，血管成像模糊且不连续，不容易直接分辨。应用本文提出的图像处理方法，首先对图3（a）进行高频强调滤波处理，结果如图3（b）所示，实验中，偏移量a取6，乘数b取4。可见，高频强调滤波将变换后的低频和高频分量都得到了增强，图像边缘更加清晰，血管位置更加突出。图3（c）是基于多尺度Hessian矩阵的血管增强结果，尺度σ∈（1，11），同时设置迭代步长为1。经过这一步处理后，图3（c）已经展现出大部分的血管脉络。图3（d）是基于梯度场均衡化增强的血管梯度图像，图像直方图均衡化的动态范围是0到图像梯度模的最大值。通过这一步操作，血管与周围组织的对比度得到进一步增强。图3（e）是用八元数矢量积改进三维区域生长算法得到的血管分割二值图像。图3（f）是用AFMM方法提取的血管中心线图像。通过该实验，可以看出本文提出的图像处理方案能够有效增强模糊血管区域，高质量提取出组织血管脉络。

2.2黑色素瘤血管图像提取

为验证本文方法对真实肿瘤微血管图像处理的有效性，对小鼠黑色素瘤微血管图像做了相应的处理和分析。黑色素瘤图像是在小鼠种植肿瘤5天后，通过ORPAM系统成像的结果，实验使用的波长为633nm。图4（a）表示小鼠体内、黑色素瘤的ORPAM活体图像。肿瘤图像数据是550×550×25大小的块构建出的550×550的MAP图。从图4（a）中可以看出，肿瘤图像中的血管非常细小、密集，直接观察只能看到成片的点状信号，微血管之间的界限和脉络难以分清。因此，将本文的图像处理方案应用到该肿瘤图像的血管处理和提取上。首先，对图4（a）进行高频强调滤波处理，结果如图4（b）所示，偏移量a取6，乘数b取4。同样，高频强调滤波后的图像低频和高频分量都得到了增强，血管对比度得到提高。图4（c）是基于多尺度Hessian矩阵的血管增强图，计算中尺度范围丢失，同第一个实验数据范围相同，同时设置迭代步长为1。Hessian矩阵计算后，图4（c）已经实现了点状信号的连接，展现出大部分的血管脉络，但是微血管的对比度较低，血管不是很清晰。因此，本文在其梯度域进行了均衡化增强操作，得到的血管梯度图像如图4（d）所示。可见，梯度域处理后，微血管图像的空间分辨率得到提高，血管之间的界限更为明显。图4（e）是用八元数矢量积方法得到的血管分割后的二值图像。图4（f）是在二值血管图像上提取的肿瘤血管骨架图。为了更清楚看到信号增强和骨架提取结果，图4（g）～图4（i）给出了图4（a）、（e）、（f）中矩形框所示子图的放大效果图，从图中可以更清晰地看到肿瘤血管脉络被有效增强和分割，血管骨架被精确提取出来。

3结语

本文提出了一种用于光声图像微血管提取的图像处理方法，融合了高频强调滤波、基于多尺度的Hessian矩阵计算、梯度场均衡化增强和相应的区域生长算法。该方法成功用于ORPAM系统得到的活体小鼠耳部正常血管图像和黑色素瘤血管图像的处理和提取中。研究结果表明，该方法可以对ORPAM系统成像中的正常血管和更加精细的肿瘤微血管进行有效处理，实现微弱信号和不连续信号的提取和分割。该工作将为下一步肿瘤微血管的量化分析提供有效的图像处理方法，促进光声成像技术在肿瘤血管生成和组织微循环研究中的深入应用。

参考文献参考文献：

[1]DELIBASIS K K， KECHRINIOTIS A I， TSONOS C，et al. Automatic modelbased tracing algorithm for vessel segmentation and diameter estimation[J]. Computer Methods & Programs in Biomedicine， 2010， 100（2）： 108122.

[2]BULLITT E， GERIG G， PIZER S M， et al. Measuring tortuosity of the intracerebral vasculature from MRA images[J]. IEEE Transactions on Medical Imaging， 2003， 22（9）： 11631171.

[3]ZHOU M， YUAN F， ZHAO J. Quantitative analysis of vessel morphology using microcomputed tomography[C]. Proceeding of 4th International Conference on Biomedical Engineering and Informatics，BMEI 2011，Shanghai，China，2011：457461.

[4]MASLOV K， ZHANG H F， HU S， et al. Opticalresolution photoacoustic microscopy for in vivo imaging of single capillaries[J]. Optics letters， 2008， 33（9）： 929931.

[5]HU S， WANG L V. Photoacoustic imaging and characterization of the microvasculature[J]. Journal of Biomedical Optics， 2010， 15（1）： 101110.

[6]ENFIELD J， JONATHAN E， LEAHY M. In vivo imaging of the microcirculation of the volar forearm using correlation mapping optical coherence tomography （cmOCT）[J]. Biomedical Optics Express， 2011， 2（5）： 11841193.

[7]CHU Z， LIN J， GAO C， et al. Quantitative assessment of the retinal microvasculature using optical coherence tomography angiography[J]. Journal of Biomedical Optics， 2016， 21（6）： 66008.

[8]T BOSKAMP， DF RINCK， B KUMMERLEN， et al. New vessel analysis tool for morphometric quantification and visualization of vessels in CT and MR imaging data sets[J]. Radiographics， 2004， 24（1）： 287297.

[9]YANG Z， CHEN J， YAO J， et al. Multiparametric quantitative microvascular imaging with opticalresolution photoacoustic microscopy in vivo[J]. Optics Express， 2013， 22（2）： 15001511.

[10]朱立新， 王平安， 夏德深. 基于梯度場均衡化的图像对比度增强[J]. 计算机辅助设计与图形学学报， 2007， 19（12）：15461552.

[11]ALEJANDRO F FRANGI，WIRO J NIESSEN，KOEN L VINCKEN.Multiscale vessel enhancement filtering[C]. International Conference on Medical Image Computing and ComputerAssisted Intervention. SpringerVerlag， 1998： 130137.

[12]王忆锋，唐利斌.利用有限差分和MATLAB矩阵运算直接求解二维泊松方程[J].红外技术， 2010（4）：213216，230.

[13]吴明珠，王晓蝉，李兴民.利用八元数矢量积改进三维区域生长算法[J].中国医学影像学杂志，2016（7）：5659.

责任编辑（责任编辑：陈福时）