CN101999885A - 一种自动分离动静脉血管的内源光学成像方法 - Google Patents

Abstract

本发明公开了一种自动分离动静脉血管的内源光学成像方法，目的是提供一种有效的自动分离生物组织动静脉血管的内源光学成像方法，技术方案是采用窄带准单色光照明，对采集到的图像序列进行频谱分析，利用图像序列的低频振荡信号的幅值分布特征，通过频谱变换，局部自适应对比度增强和自动阈值分割等操作实现了动脉血管和静脉血管的识别与自动分割。采用本发明可实现动静脉血管的自动分离，能凸显并分离很多小动脉和小静脉的血管结构，避免了使用多种中心波长的光照明时成像装置设计的复杂性。

Description

一种自动分离动静脉血管的内源光学成像方法

技术领域：

本发明属于生物医学成像方法，尤其指采用窄带准单色光照明，利用内源生理信号谱值特征的动静脉血管自动分离方法。

背景技术：

脑皮层动静脉血管的自动分离对于生物医学基础研究和血管相关疾病的临床诊断治疗具有重要的现实意义。动脉肿瘤，动脉硬化等血管疾病日益威胁着人们的生命安全，通过血管分析方法能够及时可靠的诊断和预防此类疾病。目前，对动静脉血管进行自动分离的方法的报道较为有限。根据动脉血和静脉血含氧饱和度的不同而具有的光谱吸收特性差异，可通过不同中心波长的准单色光检测目标血管的血氧饱和度的值，实现动脉血管和静脉血管的自动分离（Narasimha-Iyer, H等, 2007. Automatic identification of retinal arteries and veins from dual-wavelength images using structural and functional features （双光谱图像中的结构特征和功能特征用于视网膜动脉和静脉的自动识别）, IEEE Transactions on Biomedical Engineering, 54(8), 1427-1435），附图1给出了上述方法的流程图。但是上述方法对细小血管的分离效果很差，而且受限于预先提取血管网络结构，在现有的血管网络提取方法下，无法分离出小动脉和小静脉。以MRA成像（磁共振造影成像）为技术手段的动静脉分离方法，主要依赖于动脉和静脉血管中血液的流向相反特性 （Svensson, J等, 2002. Separation of arteries and veins using flow-induced phase effects in contrast-enhanced MRA of the lower extremities（下肢MRA成像中基于血流相位特征的动静脉分离）, Magnetic Resonance Imaging, 20, 49–57），但该特性仅存在于人体的四肢血管和部分心血管***中，在复杂的大脑血管网络中是不存在的，因此无法应用于脑皮层中动静脉的分离。生物组织中内源信号包含的频率值很丰富，在生理信号的频谱分布中，低于1Hz的信号称为低频振荡信号，这些振荡信号以0.1Hz频率为中心，包括超低频振荡信号（VLFO），低频振荡信号（LFO），其中0.1Hz低频振荡信号的存在具有普遍性，是生物体内源振荡信号的主要成分之一；通过光学成像***采集的光学信号是根据不同的血液成分对不同波长的光的吸收率和反射率不同来反映动脉血和静脉血的参数变化情况，如血流，血容，血氧浓度等，从而反映血管的动态变化过程。但目前还没有利用生物组织中低频振荡信号的谱值特征进行动静脉血管的自动分离的公开报导。

发明内容：

本发明要解决的技术问题是提供一种有效的自动分离生物组织动静脉血管的内源光学成像方法，采用窄带准单色光照明，对采集到的图像序列进行频谱分析，利用图像序列的低频振荡信号的幅值分布特征，实现动脉血管和静脉血管的识别与自动分割。

为解决上述技术问题，本发明动静脉自动分离的内源光学成像方法包括以下步骤：

第一步，将窄带准单色光照射到被测对象上，用CCD或CMOS相机通过光学成像***，以相同的曝光时间和帧间隔时间连续采集N帧被测对象反射的光学图像，每帧图像的曝光时间不高于100ms，帧间隔时间∆T不高于500ms，采集的帧数N≥20；

第二步，设单幅图像的尺寸为I_R×C，其中R为图像矩阵的行数，C为图像矩阵的列数，则该图像共有R×C个像素，相应的图像时间序列中包含R×C个像素时间序列；平均像素时间序列

中第n个数据点的平均灰度值

按公式一有：

公式一，

其中

为第n幅图像的第i个像素的灰度值，

为N帧图像所有像素的平均灰度值；

第三步，利用快速傅里叶变换算法（A.V. 奥本海姆, R.W. 谢弗（美）．离散时间信号处理．***，刘树樘 译．北京：科学出版社，1998）将第二步所得的平均像素时间序列

变换到频谱域，确定频率区间(0.05Hz, 1Hz)中具有最大幅值的频率值f _m；f _m通常位于0.1Hz附近，称为0.1-Hz低频振荡。

第四步，利用快速傅里叶变换算法将每个像素时间序列变换到频谱域，获得f _m频率处的幅值A_m(r , c)，称为特征频率幅值，其中r为像素所在位置的行，c为像素所在位置的列；以每个像素时间序列对应的特征频率幅值A_m(r , c)为灰度，组成二维的特征频率幅值图A_m；

第五步，利用快速傅里叶变换算法将每个像素时间序列变换到频谱域，计算频率区间[1Hz, 1/(2∙∆T)]中所有离散频率点的谱功率总和W _g(r , c)，其中各频率点的谱功率为相应频率幅值的平方；计算特征频率比值R_a(r , c)，R_a(r , c)为各像素序列对应的特征频率f _m的谱值(A_m(r , c))²与W _g(r , c)的比值，R_a(r , c)=(A_m(r , c))²/W _g(r , c)；以各像素序列对应的特征频率比值R_a(r , c)为灰度，组成二维的特征频率比值图R_a；

第六步，为提高动脉血管和静脉血管的灰度分布差异，利用局部对比度增强方法（C. Sinthanayothin, F.J. Boyce, H.L. Cook, 1999. Automated localisation of the optic disk, fovea, and retinal blood vessels from digital colour fundus images （眼底数字彩色图像中视盘、中央凹和视网膜血管的自动定位）. British Journal of Ophthalmology 83, 902-910）对由第四步和第五步得到的特征频率幅值图A_m和特征频率比值图R_a做局部对比度增强，具体操作如下：

6.1 以图像中第r行，第c列处的像素为中心，选取像素大小为49×49的空间滑动窗W，按公式二得到中心像素的新灰度值p(r, c)，

公式二

其中，p _min为待操作图像所有像素灰度值的最小值，p _max为待操作图像所有像素灰度值的最大值，

按公式三计算，

公式三

其中

和分别为W滑动窗所覆盖的所有像素的灰度平均值和方差；

6.2 遍历图像A_m和R_a中所有的像素，将每个像素的灰度值按步骤6.1的方法进行更新，得到局部对比度增强后的特征频率幅值图I_Am和特征频率比值图I_Ra。

第七步，使用阈值分割方法（章毓晋 著．图像分割．北京：科学出版社，2001）将第六步中所得局部对比度增强后的特征频率幅值图I_Am作二值化处理，得到二维动脉血管结构图像I_artery，其中动脉血管区域像素值为1，其他区域像素值为0；

第八步，使用阈值分割方法将第六步中所得图像I_Ra作二值化处理，得到二维静脉血管结构图像I_vein，其中静脉血管区域像素值为1，其他区域像素值为0；

第九步，为体现动静脉分离后的效果，可以将图像I_artery和图像I_vein按如下方式构建新的RGB图像C_RGB，动脉区域和静脉区域分别用R（红色）、G（绿色）、B（蓝色）三基色中的一种来表示。设置动脉区域的颜色值为红色：（1，0，0），静脉区域的颜色值为蓝色（0，0，1），其他区域的颜色值为白色（1，1，1）；从而绘制出一幅以不同颜色标记动静脉血管的RGB图像。

本发明以窄带准单色光照明，用CCD或CMOS相机通过内源光学成像***连续采集若干帧光学图像，对图像序列做频谱分析，获取0.1-Hz振荡的特征频率幅值图和比值图，以及局部对比度增强等操作，实现生物组织中动静脉血管的自动分离。与其他现有的方法相比，采用本发明可以达到以下技术效果：

1. 采用一种准单色光照明，避免了使用多种中心波长的光照明时成像装置设计的复杂性。

2. 利用低频振荡信号的幅值分布特性，实现了动静脉血管的自动分离。

3. 由于低频振荡信号能够反应动脉和静脉的功能性差异，因此采用本发明能凸显并分离很多小动脉和小静脉的血管结构，可避免造影剂注射的MRA成像方法区分动静脉血管时产生的副作用，而这些小动脉和小静脉的血管结构由于相对于背景软组织的对比度很低，仅依靠光学图像的灰度特征是很难提取和分离的。

4. 本发明有助于对生物组织中二维血管的形态和动静脉血管中血流分布，及血管中血液动力学变化进行实时、动态和高时空分辨率的监测。

附图说明

图1为已报导的视网膜动静脉分离方法流程图（Narasimha-Iyer, H等, 2007. Automatic identification of retinal arteries and veins from dual-wavelength images using structural and functional features （双光谱图像中的结构特征和功能特征用于视网膜动脉和静脉的自动识别）, IEEE Transactions on Biomedical Engineering, 54(8), 1427-1435）。

图2为本发明总体流程图。

图3为采用已报导的内源光学成像***（Optical Imaging Inc., Germantown, NY, USA）采集去除颅骨的大鼠顶叶脑皮层图像并用本发明所提出的自动分离动静脉血管的方法进行处理的结果。（a）脑皮层灰度图像；（b）平均时间序列的频谱分布；（c）各像素时间序列的f _m幅值分布图A_m；（d）各像素时间序列的f _m功率比值图R_a。（e）用不同颜色表示的动静脉血管分离结果图。

具体实施方式

图2为本发明总体流程图，本发明包括以下步骤：

1、 将窄带准单色光照射到被测对象上，用CCD或CMOS相机通过光学成像***，以相同的曝光时间和帧间隔时间连续采集N帧被测对象反射的光学图像；

2、 按公式一计算平均时间序列

中第n个点的平均灰度值

：

公式一

3、 利用快速傅里叶变换算法将第二步所得的平均像素时间序列

变换到频谱域，确定频率区间(0.05Hz, 1Hz)中具有最大幅值的频率值f _m；

4、 利用快速傅里叶变换算法将每个像素时间序列变换到频谱域，计算每个像素时间序列在特征频率f _m处的幅值，组成特征频率幅值图A_m；

5、 利用快速傅里叶变换算法将每个像素时间序列变换到频谱域，计算频率区间[1Hz, 1/(2∙∆T)]中所有离散频率点的谱功率总和W _g(r , c)，其中各频率点的谱功率为相应频率幅值的平方；计算特征频率比值R_a(r , c)，R_a(r , c)为各像素序列对应的特征频率f _m的谱值(A_m(r , c))²与W _g(r , c)的比值，R_a(r , c)=(A_m(r , c))²/W _g(r , c)；以各像素序列对应的特征频率比值R_a(r , c)为灰度，组成二维的特征频率比值图R_a；

6、 利用局部对比度增强方法对由第4步和第5步得到的特征频率幅值图A_m和特征频率比值图R_a做局部对比度增强（公式二和公式三），分别得到增强图像I_Am和I_Ra；

7、 使用阈值分割方法将第6步中所得图像I_Am作二值化处理，得到二维动脉血管结构图像I_artery；

8、 使用阈值分割方法将第6步中所得图像I_Ra作二值化处理，得到二维静脉血管结构图像I_vein；

9、 将得到的动脉和静脉以不同的颜色区别标记，以显示动静脉分离的效果。

动物实验：

如图3所示，图3-(a)是(546±10)nm准单色光照射下的单帧光学灰度图像；图中所示标尺长度为0.5mm。我们采用购买自Optical Imaging公司（Germantown, NY, USA）的光学成像与采集装置，光学采集装置置于物理隔震台上。实验对象采用Sprague-Dawley大鼠，固定于立体定位架上，以卤素灯为光源，经滤光片（546

10nm）滤波后斜入射到去除头骨后的大鼠顶叶皮层上。经大鼠顶叶皮层反射后的光线由光学成像***成像，并通过图像采集卡与计算机通信，由计算机存储采集到的图像序列。曝光时间30ms，帧时间间隔66.67ms，连续采集1500帧光学灰度图像。冷光源经滤光片滤光得到的绿色光束（波长为546

10nm）通过光纤照射到待采集区域。经过血管区域反射的光线通过CCD镜头成像并由采集***完成血管图像数据的采集和存储。图像序列采集频率为15Hz，采集过程是在静息（无刺激）状态下完成的。将采集到的图像时间序列变换到频域，构建特征频率幅值图和比值图，局部自适应对比度增强等操作，构建一幅由不同颜色标记动静脉血管的伪彩色图像。图3-(b)是按本发明第二步和第三步对1500帧光学灰度图像序列中的每一帧图像做灰度平均，得到有1500个采样点的平均灰度时间序列，再用快速傅里叶变换算法将平均灰度时间序列变换到频域，搜索0.1Hz频率附近的峰值幅度，由此确定具有最大幅值的特征频率的位置f _m；图3-(c)是按本发明第四步对每个像素位置处对应的灰度时间序列用FFT方法做频谱变换，计算f _m处的幅值作为该像素位置处的灰度值，组成的特征频率幅值图并做局部对比度增强（第六步）后的结果；图3-(d)是按本发明第五步计算每个像素位置处对应的特征频率f _m的幅值与≥1Hz频率区间内谱功率之和的比值图，并做局部对比度增强（第六步）后的结果。按本发明第八步和第九步对图3-(c) 和图3-(d)分别做二值化处理，前者得到二维动脉血管结构图，后者得到二维静脉血管结构图。由图3-(c)和图3-(d)可见，0.1-Hz振荡信号在动脉区域的幅值要高于其在静脉和皮层其他区域的幅值，而静脉区域的0.1-Hz振荡信号的幅值是最低的，因此，对图3-(b)和图3-(c)通过自动阈值分割方法作二值化处理，可分别得到动脉和静脉血管结构。按发明内容所述第九步，将分割出来的动脉区域用红色表示，静脉区域用蓝色表示，其他皮层区域用白色表示，构建一幅RGB伪彩色图像，如附图3-(e)所示。由附图3-(e)可知，所分割出的动脉血管和静脉血管的走向与其解剖学特征相吻合，说明动脉血管和静脉血管得到了很好的区分，证明了本发明所提出的自动分离动静脉血管的内源光学成像方法的有效性。

Claims (2)

1.一种自动分离动静脉血管的内源光学成像方法，其特征在于包括以下步骤：

第一步，将窄带准单色光照射到被测对象上，用CCD或CMOS相机通过光学成像***，以相同的曝光时间和帧间隔时间连续采集N帧被测对象反射的光学图像，每帧图像的曝光时间不高于100ms，帧间隔时间∆T不高于500ms，采集的帧数N≥20；

第二步，按公式一计算平均时间序列 中第n个点的平均灰度值：

公式一，

其中

为第n幅图像的第i个像素的灰度值，为N帧图像所有像素的平均灰度值， R为图像矩阵的行数，C为图像矩阵的列数，单幅图像的尺寸为I_R×C，该图像共有R×C个像素，相应的图像时间序列中包含R×C个像素时间序列；

第三步，利用快速傅里叶变换算法将第二步所得的平均像素时间序列

变换到频谱域，确定频率区间(0.05Hz, 1Hz)中具有最大幅值的频率值f _m，f _m位于0.1Hz附近，称为0.1-Hz低频振荡；

第四步，利用快速傅里叶变换算法将每个像素时间序列变换到频谱域，获得f _m频率处的特征频率幅值A_m(r , c)，其中r为像素所在位置的行，c为像素所在位置的列；以每个像素时间序列对应的特征频率幅值A_m(r , c)为灰度，组成二维的特征频率幅值图A_m；

第五步，利用快速傅里叶变换算法将每个像素时间序列变换到频谱域，计算频率区间[1Hz, 1/(2∙∆T)]中所有离散频率点的谱功率总和W _g(r , c)，其中各频率点的谱功率为相应频率幅值的平方；计算特征频率比值R_a(r , c)，R_a(r , c)为各像素序列对应的特征频率f _m的谱值(A_m(r , c))²与W _g(r , c)的比值，R_a(r , c)=(A_m(r , c))²/W _g(r , c)；以各像素序列对应的特征频率比值R_a(r , c)为灰度，组成二维的特征频率比值图R_a；

第六步，利用局部对比度增强方法对由第四步和第五步得到的特征频率幅值图A_m和特征频率比值图R_a做局部对比度增强，具体操作如下：

以图像中第r行，第c列处的像素为中心，选取像素大小为49×49的空间滑动窗W，按公式二得到中心像素的新灰度值p(r, c)，

公式二

其中，p _min为待操作图像所有像素灰度值的最小值，p _max为待操作图像所有像素灰度值的最大值，

按公式三计算，

公式三

其中

和

分别为W滑动窗所覆盖的所有像素的灰度平均值和方差；

遍历图像A_m和R_a中所有的像素，将每个像素的灰度值按步骤6.1的方法进行更新，得到局部对比度增强后的特征频率幅值图I_Am和特征频率比值图I_Ra；

第七步，使用阈值分割方法将第六步中所得局部对比度增强后的特征频率幅值图I_Am作二值化处理，得到二维动脉血管结构图像I_artery，其中动脉血管区域像素值为1，其他区域像素值为0；

第八步，使用阈值分割方法将第六步中所得图像I_Ra作二值化处理，得到二维静脉血管结构图像I_vein，其中静脉血管区域像素值为1，其他区域像素值为0。

2.如权利要求1所述的一种自动分离动静脉血管的内源光学成像方法，其特征在于将图像I_artery和图像I_vein按如下方式构建新的RGB图像C_RGB，动脉区域和静脉区域分别用红色R、绿色G、蓝色B三基色中的一种来表示，设置动脉区域的颜色值为红色：（1，0，0），静脉区域的颜色值为蓝色（0，0，1），其他区域的颜色值为白色（1，1，1）；从而绘制出一幅以不同颜色标记动静脉血管的RGB图像。

Images