利用自动图像分析测试机能微循环几何结构和速度分布

摘要

本发明提供了对微脉管视频序列进行量化估计的分析算法，该视频序列提供每个脉管段的脉管厚度、脉管长度和血流速度。本发明还提供了用于计算分布在视觉上具有不同厚度的脉管上的机能微脉管密度和血流速度的方法。

说明书

技术领域

本发明涉及分析微脉管视频的方法。更具体地，本发明涉及测量机能微循环几何结构和血液速度分布。

背景技术

用于研究体内人体微循环的技术在近几十年已快速发展。毛细管显微镜检查是使用活体显微镜对微脉管床成像的技术。长时间以来，其为研究易进入组织(例如，皮肤或指甲下方的毛细管床)的微循环的唯一可用技术。然而，基于显微镜的结构使毛细管显微镜检查系统不能进入深处的组织。最近的研究描述了正交偏振光谱(OrthogonalPolarization Spectral，OPS)成像[3，4，18，20，23，25，30]或测流暗视野(Side-stream Dark Field，SDF)成像[13]。这些技术基本上是具有改进的目标照射的手持式显微镜，从而其获得的对比度比毛细管显微镜检查系统获得的对比度更高[23]。这些系统的另一优点在于长的探针，其使得尖端可达活体显微镜不能到达的器官表面。OPS成像系统已被证明在研究例如人体甲皱[23]、败血症[30]期间的舌下组织[20]、动脉瘤手术[25]期间的脑组织、鼠结肠[3]、鼠皮瓣[18]和仓鼠皮瓣[10]的微循环中具有价值。一些OPS研究者使用半定量方式的分析，其中将血流(没有流动、间歇流动、缓慢流动、持续流动)记为[4，30]三个脉管类型：小(10-25μm)、中(25-50μm)和大(50-100μm)。

除了用于描述微循环的定量几何参数之外，越来越关注微循环血液流速的量化。已描述了几个方法来根据视频帧序列来确定速度[12，14，19，20]。空间相关技术[20]从参考帧中选择片段并在接下来的帧中跟踪该片段。在微循环中，这种片段经受形态改变。这是由于脉管壁附近的细胞以低于中心处的细胞的速度移动、并且例如由于脉管曲率而观察到径向运动和垂直于焦平面的运动而造成的。光流技术[12]依赖于可追踪对象的强度随时间变化保持恒定的事实。细胞或细胞群也在垂直于焦平面的方向移动。这使得对象的强度在其它细胞重叠时会及时改变，从而不满足持续对象强度的约束。使用各向异性扩散滤波可克服这一问题，但需要存在相对大的等离子体间隙以检测大的红血球栓。

用于速度估计的时空图[14]已在微循环的许多研究中[7，15，16，17，19，23，32]使用。这些图绘制了直脉管段随时间变化的纵向强度曲线。时空图的对角带表示穿过脉管移动的对象。当等离子体间隙或白血球通过脉管时观察到亮带，而暗带则表示存在红血球。该技术的一大优势在于其包括了用于速度估计的全部可用的空间和时间数据，并且研究者从时空图中出现的线接收流动类型的即时光学反馈。

Klyscz和coworkers[17]描述了展示量化分析微脉管视频的技术的计算机程序。用屏幕上测径器确定局部脉管宽度，用允许交互追踪脉管的制图工具获得脉管长度。其提供了机能毛细管密度(FCD)，该密度定义为脉管总长度(L)(由用户追踪得到)与感兴趣的图像面积(A)的比值。其还利用时空图[14，17，19]估计血流速度。该方法要求用户在脉管的中心线处绘制直线。该直线指示获取图像数据以生成时空图的位置。尽管该图在其区域内是唯一的，但其需要大量的用户交互，从而增加了观察者偏差并使分析耗时。

发明内容

本发明的目的在于提供分析微脉管视频以确定血流速度的方法，该方法需要较少的用户交互并减小了观察者偏差。

该目的通过提供分析微脉管视频以确定血流速度的方法而实现，该方法包括：

接收微脉管视频图像序列，每个视频图像均显示脉管结构；

在脉管结构中识别脉管中心线像素；

通过对中心线像素分组而定义脉管段；

为每个脉管段自动生成时空图，该时空图在与每个脉管段相关联的中心线像素处生成；

利用时空图自动估计每个脉管段中的血流速度。

根据本发明的一个方面，提供了一种计算机程序产品，其被设置为在加载到计算机上时使计算机能够执行上述方法。

根据本发明的另一方面，提供了包括这种计算机程序产品的数据载体。

附图说明

下面参照附图利用大量示例性实施方式对本发明进行更详细的描述。

图1A和1B示出了指示根据本发明实施方式的处理步骤的流程图；

图2示出了预定的捕获区域，在该区域中搜寻下一中心线像素；

图3A示出了正交像素栅格中描述实际脉管中心线的假想线；

图3B和3C示出了局部脉管方向的两个实例；

图4A示出了原始的脉管图像；

图4B示出了图4A的脉管的拉直版本(中心线在每个水平线的中途)；

图4C示出了图4A的原始图像中的中心线以及脉管壁检测的早期结果(示出了脉管壁被误解的位置)；

图4D示出了填充间隙和平滑脉管壁的结果；

图5示出了时空图的实例；

图6A示出了弯曲的段中的中心线的一部分的实例；

图6B示出了由于脉管方向改变导致的时空图中的方向(即，速度方向)的可能的错误改变；

图7A是在一个灰度级得到的、作为和ρ的函数的霍夫计数的图示；

图7B是至少为图7A的最大值的90％的单元的霍夫计数的图示；

图7C是图7B中的霍夫计数在显示全部灰度级(或选定灰度级)的总和的图示；

图8A是作为的函数的霍夫计数的图示；

图8B是显示的定义的相应时空图；

图9示出了不同方向的模拟脉管的平均长度估计相对于实际脉管长度和宽度(以像素为单位)的误差，误差条指示由脉管方向引起的误差；

图10示出了具有l＝100像素且的模拟脉管在不同的分析尺度(σ)下长度估计误差；

图11示出了不同分析尺度下的厚度估计误差(尺度为sigma＝3.0)的图示；

图11B示出了不同方向的模拟脉管的平均厚度估计相对于实际脉管厚度的误差的图示；

图12A示出了速度增加的十一个时空图；

图12B示出了利用结构张量、灰度级霍夫变换和交互地追踪图12A的时空图中的五条线进行速度估计的精度的图示；

图12C示出了交互估计的速度误差的图示，其中实线表示平均速度误差，虚线表示伴随的速度误差范围；

图13A示出了健康男性个体的舌下视频序列的平均帧；

图13B示出了叠加有分析结果的图13A的图像；

图14A示出了健康男性个体的机能微循环密度分布；

图14B示出了健康男性个体的速度分布的图示；

图15A示出了显示脱位过程之前的舌下微循环的视频序列中的平均帧；

图15B示出了与图15A相同的数据，但其叠加有分析结果；

图15C示出了在脱位过程中获得的帧的平均；

图15D显示了叠加在图15C上的分析结果；

图16示出了由图15A中的箭头标记的小静脉的时空图；

图17A示出了心脏脱位之前和之中的机能微循环密度分布；以及

图17B示出了心脏脱位之前和之中的速度分布的图示。

具体实施方式

利用当前可用的成像技术(例如，毛细管显微镜检查、OPS或SDF成像)，在存在包含血红蛋白的红血球的情况下仅观察到“脉管”，与背景介质不同，血红蛋白高度吸收入射波长。在不存在红血球的情况下，毛细管本身基本是这些成像技术不可见的。微循环的视频仅显示出由脉管壁描绘并因此被称为“脉管”的红血球栓。

利用时空图的脉管分段和血流速度估计需要多个图像处理步骤。图1A和1B示出了根据本发明实施方式的处理步骤的流程图。这些处理步骤可由配置为执行这些处理步骤中的一个或多个的计算机执行。在第一步骤2中，从成像系统接收微循环视频序列。在产生该视频序列时，轻微的尖端移动与手持式微脉管的高效光学放大相结合可导致图像间移动(inter-image translation)。该图像内移动会妨碍速度测量，并且，根据一个实施方式，利用2D互相关(见Altman[2])在图像登记步骤4对图像内移动进行补偿。在图像登记步骤4，可通过两个方式加强图像对比度。其一，可通过减去与图像最佳匹配的二次多项式表面并加上原始图像的平均图像强度，来减小每一帧的背景中的强度变化；其二，可通过利用通常所说的传递函数以Pries[26]描述的方式将输入图像的每个灰度级映射到输出图像的灰度级，来操作图像灰度级直方图(8比特范围)，从而实现对比度增强。在接下来的步骤6，对图像帧进行时间平均，以填补由于等离子体间隙或白血球的存在而导致的毛细管中断。该平均使得毛细管将被检测为整体，而不管中断的细胞流。该平均还减小了噪声的作用，这对于分段处理是有益的。在下一步骤8中，如Steger[31]所描述地，对预处理的图像进行中心线检测。该方法是基于计算海赛(Hessian)矩阵的本征向量，并得到指向最大表面曲率方向(即，垂直于脉管方向)的向量[29]以及正交方向(即，脉管方向)的向量。如果由本征值λ_n表示的、方向n上的表面曲率明显高于由λ_t表示的垂直方向t(沿脉管的切向)上的表面曲率，图像像素则被认为是候选的中心线像素。通过评估|λ_n|/|λ_n|+|λ_t|)≥ε_thr，对该条件进行测试，其中ε_thr是给定的阈值(表1)。如果交叉强度曲线(即，方向n上的强度曲线)为其极值[31]，候选像素则实际上记为中心线像素。

中心线检测得到二进制图像，其中中心线像素被设置为1。在步骤10，将各中心线像素连接为线段，用于确定脉管长度。由于可参考“脉管”而不是像素，因而减轻了用户交互。脉管方向影响每个像素的长度贡献。该方向由海赛分析得到，并在步骤12中用于计算每个像素的长度贡献和包括大量像素的脉管段的总长度。在下一步骤14，对每个中心线像素检测脉管壁，并将检测到的脉管壁用点标记，其中交叉强度曲线显示其最大陡度。这提供了对局部脉管厚度和平均脉管厚度的估计。

如果微循环图像是未聚焦的，图像则出现模糊并从而对小脉管中的脉管厚度过度估计。因此，在步骤16中，将未聚焦的脉管排除。在一个实施方式中，使用脉管的全部边缘点处的平均梯度量来确定每个脉管的聚焦分数F。通过将该聚焦分数相对于每个边缘像素局部(200像素×200像素)的背景强度标准化，使得该聚焦分数对背景变化较不敏感。该聚焦分数在[6]中定义。在一个实施方式中，用户能够仅计入具有超出给定限制的聚焦分数的脉管。

血流在分叉处分流为新的支路，这使血流速度改变。为了精确估计血流速度，需要根据由分叉描绘的像素确定时空图。在步骤18，通过确定当前脉管的端部和相邻脉管的壁之间的距离，自动执行在分叉处切割脉管的处理。如果该距离小于3/4×相邻脉管厚度，则在最接近当前脉管端部的点处将相邻脉管切割为二。3/4×因子(即，1¹/₂×相邻脉管的半径)允许切割接近相同交叉点的脉管段，但防止了由邻近的脉管(例如，平行的脉管)产生的切割。对全部可用的脉管段均执行步骤18。

在步骤20，用户能够通过删除、切割或连接脉管段而与中间结果进行交互。可通过以可选的检测尺度进行的局部分析将未检测到的脉管段“交互地”添加。大尺度分析集中于大的结构细节，而小尺度则集中于小的结构。交互分析包括与自动分析相同的分析技术，但是，其以可选择的尺度进行并且在由用户识别的局部区域中进行。分析区域基于包围用户追踪线的端坐标的矩形。该区域被稍微放大，以将希望的目标脉管检测为整体。通常在这种矩形区域中会发现许多脉管段，但仅选择与用户追踪线最佳匹配的一个脉管段。接下来，对于每个脉管段，在步骤22估计血流速度。最终结果由报告生成器显示，示出脉管密度和速度分布，见步骤24。

利用图1B更详细地解释步骤22。图1B示出了计算脉管中的速度所采用的步骤的流程图。对于感兴趣区域内的每个脉管段，重复图1B的步骤。根据一个实施方式，在每个脉管的中心线处自动获取时空图，见图1B的步骤26。在步骤28，对时空图中的脉管曲率进行校正。在步骤30，利用图像直方图均衡，例如见[27]，自动改善时空图中的线结构的可见性。

在步骤32，通过估计时空图中的线的方向，计算血流速度。现今，线的方向是通过在时空图中交互地对线进行追踪或通过利用相关技术[17]而实现的。根据本发明的一个实施方式，通过利用通常所说的结构张量自动计算线的方向。在另一实施方式中，通过利用灰度级霍夫(Hough)变换进行该估计。自动估计的优势在于，在该步骤中不需要用户交互来追踪时空图中的线。在图1B所示的下一步骤34中，用户可接受或拒绝自动或交互速度分析的结果(验证)。最后，将时空图中的结构方向转换为每个脉管段中的实际速度值。

下面，更详细地解释上述某些方面。本申请中描述的分析算法的性能依赖于表1中列出的大量参数。

脉管分段

分析尺度

上述步骤8的中心线检测需要确定图像导数。如果将这些图像导数计算为相邻像素之间的差异，其对噪声是非常敏感的。因此，在一个实施方式中，使用高斯导数，其包括高斯核的工作距离内的图像数据。通过将图像与相应的高斯导数卷积，获得以上使用的导数的高斯分量：

$I_i(x,y)=I(x,y)\otimes G_i(x,y)${i＝x，y}    (1a)

$I_{i,j}(x,y)=I(x,y)\otimes G_{i,j}(x,y)${i＝x，y；j＝x，y}    (1b)

$G(x,y)=\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{1}{2}\frac{x^2+y^2}{{\sigma }^2}},$$G_i(x,y)=\frac{\partial G}{\partial i},$$G_{i,j}(x,y)\frac{{\partial }^{2}G}{\partial i\partial j}---\left(1c\right)$

将高斯滤波的标准偏差σ用作尺度参数。以更大的尺度参数进行分析增加了分析的空间范围，并聚焦于更大的图像特征。尺度参数σ用于中心线检测。其它与距离相关的分析参数基于表1所示的尺度参数。

将中心线点分组为脉管段

步骤10中对中心线像素的分组需要将可能属于相同脉管段的像素连接起来。由于噪声而存在的像素、或者交叉或平行延伸的脉管的像素被分配给其它的脉管段。海赛分析提供每个中心线像素x处的切向向量t和法向向量n，见[31]。

该连接从具有最大曲率(最大二阶导数λ_n)的中心线像素处开始。由于脉管可在像素的两侧扩张，因此，以当前像素x作为参考，在两个相反方向(朝t或-t)继续追踪相邻的像素。该方法在预定的捕获区域内搜索下一个像素p＝(p_x，p_y)，该捕获区域最好由三角形描绘，该三角形由张角α和给定搜索范围(r_max)内距离为r处的中垂线(朝向n)限定，该三角形的平分线指向切向向量t，如图2所示。候选像素p如下计算：

p＝x+rt+in    (i＝0，±1，±2，...，±w；r＝1，2，...，r_max)    (2a)

$w=r\tan \left(\frac{\alpha }{2}\right)---\left(2b\right)$

可使用多个条件在属于相同或不同脉管段的中心线像素之间求导[29，31]。下面列出其中的两个条件：

1)由于脉管中心线的方向在从一个像素到另一像素时不会突然改变，因此在候选像素位置处的切向向量t应该与参考像素(x)的切向向量类似。如果两个切向向量的内积(c1)较大，则满足该条件。

2)候选像素应和参考像素属于相同的脉管，而不应属于与之平行延伸的脉管。通过评估参考像素的切向向量与参考向量到候选向量的距离向量的内积(c2)，对该条件进行测试。

将距离参考向量最近(最小的r)且具有最大c1+c2值的候选向量，如[31]，接受为脉管段中的下一个中心线像素。通过线性内插获得参考像素与最佳相邻像素之间丢失的像素。重复上述处理，从而对仍然保持在其它脉管段中的全部各中心线像素进行分组。

每个中心线像素均对总脉管段长度具有贡献。中心线检测处理的有利性质在于，其对线产生单一的响应[31]。中心线像素可由图像噪声而得到，尤其是当以小尺度检测脉管时。因为由于噪声的存在而存在的脉管段通常非常小，因此将具有有限长度(s_min)的全部段均去除。

脉管长度估计，见图1的步骤12

脉管方向影响每个像素的长度贡献。因此，在估计脉管长度时，包括了局部脉管方向。图3A在正交像素栅格中示出了假想线(虚线)。虚线中心线由最佳描述实际脉管中心线的像素(正方形)表示。实线段示出了每个像素的长度贡献，其和构成总的脉管长度。在图3B中，示出了局部方向的两个实例。局部方向由从海赛分析[31]得到的切向向量(t)导出。在已知时，可计算出像素-长度贡献Δl。如果主方向(principle direction)是x方向，即且Δx等于一个像素在水平方向上的长度贡献，那么，否则，其中Δy表示一个像素在垂直方向上的长度贡献。

脉管壁检测

现在参照图4A、4B、4C和4D，讨论检测单一脉管的壁的实施例。图4A示出了原始脉管图像。图4B显示了和图4A相同的图像，但通过二次抽样将该图像在法向拉直。中心线位于每个水平图像线的中途。图4C将早期脉管壁检测的结果叠加在图4A的原始图像上。图4D示出了脉管壁检测的最终结果。

利用每个中心线像素处的切向向量t和法向向量n，可通过在垂直于切向向量t的方向在子像素级(通过线性内插)对图像进行采样，来确定截面强度曲线。对每个中心线像素重复该过程，以获得拉直的脉管(见图4B所示的实例)。如果等离子体间隙或白血球妨碍了红血球的持续流动，在拉直的脉管中则可出现间隙。为此，使用各向异性扩散内核[11]，其具有主要在图4B的垂直方向(σ_along)上延伸的高斯响应。其有效地将中断闭合，并将脉管检测为整体。在水平方向，将高斯滤波内核的一阶导数用作最大梯度检测器，其在横向(σ_cross)具有小的范围以保持良好定位的边缘检测[5]。

与边缘检测内核的卷积可使脉管厚度估计具有偏差，尤其是对于小脉管。滤波器自身的脉冲响应在x＝±σ_cross处显示出最高梯度(G”(x)＝0))。因此，脉管被检测为至少具有2σ_cross个像素的宽度。

在出现大的等离子体间隙的情况下，或者当脉管临时移出焦平面从而减小了与噪声比的对比度时，在例如分叉点处，边缘点不能被正确地检测，见图4C。当具有更高强度梯度的脉管平行延伸、或当脉管部分重合时，边缘点也可能被误解。被误解的边缘点(伪像)基本上偏离其到中心线的平均距离。通过在每个重复过程中将超出平均距离两个标准偏差的最远边缘点去除，将该性质用于重复地去除伪像。重复该过程，直到全部剩余的距离采样均为与平均距离的偏离在两个标准偏离之内。得到的平均距离被分配至所有的伪像位置。最后，执行高斯滤波(σ_edge)以平滑从一个脉管端到另一脉管端之间的边缘-中心线距离，如图4D所示。对相对的脉管壁重复上述处理，并生成对局部和平均脉管厚度的估计。

速度估计

在一个实施方式中，根据时空图估计血流速度，该时空图是通过绘制中心线强度随时间t的变化而自动创建的，如图5所示。时空图中的线的斜率β是血流速度的测量，其计算为v＝Δs/Δt＝tanβ，其中Δs是时间段Δt内的位移。

与最小毛细管中的血流不同，较厚脉管中的空间速度曲线将不相同，并且类似于Poiseuille抛物线流，其中，速度在脉管中心线处具有最大值，而在脉管壁处为零。在这种脉管类型中，通常将中心线处或其附近的速度用作速度指示[7，17，19，20，23，32]。现今，经常交互地估计线结构的方向，例如，通过交互地追踪时空图中的线和通过计算平均方向。根据本发明，提供了一种自动确定时空图中的线的方向估计的方法。

曲率校正

时空图的时间轴是帧间隔的倍数，其在CCD相机中是非常精确的。另一方面，空间轴不是均匀分布的，因为每个像素的长度贡献(Δl_i)依赖于局部脉管方向。图6A和6B示出了这一问题。如果暗对象移动穿过在水平和垂直方向等间距的像素栅格上的弯曲的脉管(图6A)，那么，该对象在水平截面(Δl_B)穿过相同数量的像素比在对角截面(Δl_A)穿过相同数量的像素快倍。如果时空图的空间轴仅为中心线像素的堆叠，那么，当脉管方向改变时，速度线将向不同的方向弯曲(图6)，从而错误地暗示速度改变。通过利用线性内插将随机间隔的中心线像素映射到时空图的等距间隔，可避免该偏转误差(deflection error)。在一个实施方式中，距离采样的数量等于描绘脉管的中心线像素的数量。

利用结构张量的自动速度估计

在一个实施方式中，通过用结构张量估计线的方向，根据时空图自动计算血流速度。令I(x)表示2D线结构图像，x＝(s，t)。若u＝(Cos(β)，Sin(β))^T为指向线结构方向β的单位向量，则I(x+u)～I(x)。在垂直于u的方向，即，梯度方向强度的变化最大。因此，通过下式局部地得到方向u：

$u·▿I=0---\left(3\right)$

可通过使E最小化来估计平均方向，该E值为时空图中的全部像素根据式(3)得到的偏离的平方和，该E值由下式得到：

$E=\underset{x}{\Sigma }{\{I_t\left(x\right)\mathrm{Cos}\left(\beta \right)+I_s\left(x\right)\mathrm{Sin}\left(\beta \right)\}}^2---\left(4\right)$

最佳方向通过对下式求解得到：

Ju＝0    (5a)

$J=\left(\begin{array}{cc}\Sigma I_t^2& \Sigma I_t{I}_s\\ \Sigma I_t{I}_s& \Sigma I_s^2\end{array}\right)---\left(5b\right)$

其中，J是包含高斯图像导数的(尺度为σ_ST)的结构张量。

可见，最小化问题等价于本征值问题Ju＝λu，例如见[14]。具有最小本征值的本征向量产生线的方向。该方法的优势在于，得到提供结构细节的两个本征值。假设λ_t和λ_n是属于线方向及其法向的本征向量的本征值，当线在单一方向上出现时，则λ_n＞0且λ_t≈0。如果观察到多个方向，则λ_n＞λ_v＞0。最后，如果难以发现优势方向，则λ_n≈λ_t≈0。如果|λ_n|/(|λ_n|+|λ_t|)＞μ_thr，这些性质则用于接受自动方向估计，并因而接受速度估计。

利用灰度级霍夫变换的自动速度估计

时空图有时会显示与全局线结构相比具有不同方向的小的线伪像。这是由于利用具有有限分辨率的正交栅格在脉管中心线处进行强度采样而导致的。这种小的线伪像可对平均线结构具有大的贡献，并且不集中于感兴趣的全局线结构。霍夫变换[9]提供了一种方式以在时空图中计入线的长度，并能排除小的线伪像。霍夫变换基本上是点-曲线的变换，其检测图像中直线的参数。该技术考虑了线的极坐标表示：

其中，(x_i，y_i)是每个线像素在时空图中的坐标，是从原点开始并垂直于该线的向量的方向，ρ是该向量的长度，该长度等于到原点的线距离。每个线像素均映射到参数空间的正弦曲线参数空间的离散图像由累加器单元构成，并且对通过该单元的每个正弦曲线递增。通过将时空图中的全部线像素均转换为正弦曲线，将累加器单元加到线长度(L，像素单位)上。因此，累加器单元生成线的特征参数通过将时空图作为输入图像，希望在指定方向对距离原点不同距离(ρ)处的多个线具有高响应。

上述的传统霍夫变换需要时空图的二进制图像，该图像可通过例如阈值设定或线检测而获得。由于时空图的有限的图像质量，这些技术难以揭示线结构，并会使霍夫变换提供不精确的结果。通过为每个灰度级考虑二进制图像[21](如果该灰度级存在则将像素设为“1”，否则将像素设为“0”)，包含了更多的信息，该图像识别其自身的线集从而在霍夫空间中显示其特征峰值。通过对每个灰度级重复霍夫变换，全部像素和灰度级都包含在变换中并能对响应进行累加。对于每个灰度级霍夫变换，见图7A，维持至少为最大值90％的累加器单元，而将其它均设为零。该结果如图7B所示。这仅包括了每个灰度级的最长的线，并丢弃了小的伪像。然后，通过求和将“每个灰度级”的累加器图像相结合，并产生“长线”霍夫空间，如图7C所示。

由于仅对全局方向感兴趣，因此确定具有相同方向(但表示与原点相距不同的线距离ρ)的累加器单元的总数量，从而得到每个方向的“分数”，如图8A所示。最后，利用高斯滤波核(σ_H)对该曲线进行平滑。滤波后的曲线的最高峰值表示速度估计范围内的速度(如下文所述)，提供对全局线方向的最佳估计，见图8B，并代表不同灰度级的相同方向的长线。

速度估计范围

速度估计的物理上限依赖于脉管长度(L，单位为μm)和帧频(f)。根据慢采样场景的速度测量可能受到混叠的影响。理论上，如果对象以恒定的速度移动并仅在两个连续的帧可见，则可根据时空图计算血流速度。然而，不可能确定地判断出第一帧中的单元对象是否与第二帧中的单元对象相同。如果从其它视频帧观察到对象在连续帧之间以相当恒定的位移移动，该对象“更可能”是同一对象。因此，选择最少三个帧间隔来确定最大物理速度极限(v_max)：

$v_{\max }=\frac{L*f}{3}[\mathrm{\mu m}/s]---\left(7\right)$

通常，在时空图中用水平线标记例如由于相交的脉管产生的固定位置处的暗点(见图8A中处的峰值)，或由垂直线标记例如由于照明变化而产生的暗点。因此，拒绝与如上所述的物理极限以上或给定下限(v_min)以下的速度相对应的方向，并需要手动估计，即，通过交互地追踪时空图中的线。

分布

在脉管长度、脉管直径和血流速度可得到的情况下，可生成这样的分布，该分布显示：1)脉管类型在视图方面的贡献；以及2)在具有不同脉管直径的脉管上血流速度是如何分布的。如果脉管段被认为是中心线像素的级联，每个中心线像素则均对脉管长度提供贡献，并提供局部脉管直径。如果全部脉管段被链接为一个大的假想脉管，该脉管则将显示由长度-厚度分布提供的厚度变化。该长度-厚度分布是通过将分布箱(distribution bin)上的每个中心线像素的直径贡献进行细分并将其长度作为加权因子而获得的。

以类似的方式对面积-厚度分布(密度分布)进行估计。在该过程中，脉管段被细分为片，每个片具有一个中心线像素。这种片的投影面积是依赖于方向的长度和局部脉管直径的乘积。面积-厚度分布通过将厚度箱上的片的直径贡献细分并将投影面积用作加权因子而获得。

时空方法提供了对分配给每个中心线像素的每个脉管的血流速度的测量。通过在厚度箱上对每个片进行细分、并使用平均速度和片长度的乘积作为加权因子，确定速度-厚度分布。

实验

验征

为了验证脉管长度和厚度确定的性能，创建模拟视频(500×500像素)。每一帧包含具有高斯横截面曲线(具有标准偏差σ_l)、长度不同(50、100、150、200和250个像素)的5条线。模拟的脉管的“壁”由具有最大梯度的点标记，即，位于±σ_l处，从而得到d＝2σ_l，其中d为线厚度。线厚度在连续的帧中是改变的，其范围为[1，20]像素，并以0.5个像素递增。将背景强度和中心线强度设为200au和50au(au表示任意单位)。最后，添加高斯噪声σ_noise＝10au，其近似为典型SDF图像的高斯噪声的两倍。通过包括具有[0，90]°范围内以15°递增的不同线方向的帧，研究脉管方向的影响。

创建另一模拟视频(250×250像素)用于验证速度估计。每个视频帧均用高斯截面强度曲线(σ_cell＝3像素)显示包含为圆块的“单元”的模拟脉管。将这些单元(每个单元大约具有5个像素的脉管长度、背景强度为200au、中心强度为50au)绘制在随机的位置处，但位于延伸至每帧边缘的、10个像素宽的假想脉管的边界内。通过绘制距离脉管壁的距离为σ_cell的块来界定假想脉管壁。再次将幅度为10au的高斯噪声添加至图像。通过将单元移出脉管并将新的单元移到连续帧内，对灌注进行模拟。在方向为0°、速度范围为[2.5，2000]像素/秒(即，在25帧/秒的情况下为[0.1，80]像素/帧)的脉管中测量交互自动速度估计的精度。首先通过交互地追踪时空图中多达5条可用的线来获得速度结果。将交互获得的速度结果用作参考，用于利用灰度级霍夫变换和结构张量方法确定自动时空图分析中的误差。在自动分析过程中，与交互估计相比高达20％的速度误差水平被接受。

通过创建[0，90]°范围内以15°递增的不同方向显示脉管的类似模拟，以100像素/秒的适中速度测试旋转相关性。在此实验中，交互地确定速度。

在分段和中心线检测之前，每个速度段均覆盖被首先平均的100帧。在验证实验中，估计的精度以像素/秒为单位，因此与光学放大率无关。

临床应用

舌下视频记录通过具有标准的5倍光学放大的MicroScan SDF系统[13，22](荷兰阿姆斯特丹市的Micro Vision Medical公司)而得到，该系统得到像素间隔大约为1.5×1.4μm(h×v)的微循环图像。该硬件展现点传播，该点传播类似于在x方向具有小于1.4像素、在y方向小于1.0像素的标准偏差的高斯分布。因此，直径约为4-5μm的毛细管在标准MicroScan图像中大约为3-4像素宽。

舌下微循环通过外部颈动脉的舌动脉分支接收血液，其为中枢体循环的一部分。选择健康的男性个体的舌下视频记录用于其高对比度的适中的血流速度，这允许我们评估自动分析时空图的可行性。在用非体外循环冠状动脉旁路移植术(OPCAB)进行了心脏搭桥手术的患者的心脏脱位时，获得其它的舌下记录。在OPCAB过程中，血红蛋白浓度和血细胞比容不明显变化。在心脏脱位时，平均动脉压明显降低(＜60mmHg)。该视频序列选自一系列临床实验，因为其显示了实际同一舌下位置在脱位过程之前和之中的微循环。

分析参数

算法被配置用于利用表1给定的设置对估计和临床图像序列进行分析。通过局部分析以可选的尺度(即，σ＝1.5、3.0、6.0或12.0个像素)“交互地”添加未检测到的脉管段。

(*)对交互估计设置为2×σ

(**)对交互设置为0

表1用于自动微脉管分析的参数设置

验证实验

脉管长度

在d＝1像素的两种情况下，由于噪声的存在而使脉管长度出现错误。在这种情况下，交互地对脉管进行调整。图9示出了像素单元中观察到的长度偏差。图中的条示出了在不同方向获得的平均长度确定。误差条指示由于线方向和噪声而产生的小的误差范围(平均小于±2％)。该图示出了，长度估计的精度与线长度无关，但强烈依赖于模拟的脉管的厚度。

长度估计还依赖于分析的尺度(σ)，如图10所示。该图示出了对l＝100像素且的模拟脉管的长度估计误差。如图所示，当分析尺度较小(σ＝1.5)时，脉管端部的误差是有限的。然而，在小的尺度下，检测对噪声更敏感，且对宽度大于12个像素的脉管检测到参差不齐的脉管壁。

脉管厚度

模拟视频还用于在四个不同尺度(线检测尺度σ＝1.5、3.0、6.0、12.0个像素，即，σ_cross＝0.5、1.0、2.0、4.0个像素)下进行脉管厚度验证。在该实验中，交互地选择脉管段。将“平均”厚度用作厚度参数。通过包括具有[0，90]°范围内以15°递增的不同脉管方向的帧，再次研究脉管方向的影响。

图11A示出了小脉管的厚度被过度估计为具有增大的尺度。这是由滤波器自身的脉冲响应引起的，该脉冲响应使脉管被检测为具有2σ_cross或更宽的宽度，如“脉管壁检测”部分中所讨论的。因此，选择与可用脉管类型一致的检测尺度非常重要。图11b进一步示出了使用默认的分析尺度(σ＝3)且边缘检测核σ_cross＝1时，模拟脉管的相对的厚度-估计误差。条指示平均厚度，误差条指示作为脉管方向和误差的结果的范围。厚度为[2×σ_cross，13]个像素范围内的脉管显示高达1个像素的绝对厚度误差，即，在每个壁的检测中具有0.5个像素的误差。对于宽度大于3个像素的脉管，该小误差导致降至20％以下的相对误差(图11B)。

速度

利用图12a所示的时空图，通过交互的自动方向估计对模拟视频中的血流速度进行估计。如图所示，在低速条件下该线结构清楚可见，而在高速条件下图像则变为有噪声的。在该环境中约为2000像素/秒物理检测极限附近，时空图几乎不能显示线结构。然而，对于该模拟视频，在一些实践之后仍然能用肉眼发现线方向。图12b示出了利用以下三种方法获得的速度估计中的偏差：1)利用结构张量；2)利用灰度级霍夫变换；3)通过交互地追踪时空图中的五条线。利用结构张量的速度估计适用于速度高达100像素/秒的情况(误差小于5％)。高速由时空图中的局部伪像标记，这些伪像导致大的速度误差(＞20％)。利用霍夫变换的方法适用于速度高达750像素/秒的情况(误差小于5％)。在更高的速度(大于1250像素/秒)处，该方法则失效并会选择通常导致大的速度误差的可选方向。交互地追踪时空图中的线提供了最佳结果，并对高达该模拟实验中可接受的物理速度极限2000像素/秒的情况都是可行的(误差小于5％)。图12C显示了通过交互地追踪时空图中的任意五条线获得的平均速度误差(v_err)和速度误差范围(虚线)。对于高达1000像素/秒的速度，显示了小于4％的相当小的速度误差。由于追踪陡线段的不精确性，更高的速度显示了更大的误差。脉管旋转的影响在交互和霍夫确定的速度估计的误差上增加了多达1.1％(未示出)。利用结构张量，由于旋转而导致的的最大误差为4.8％。

临床应用

健康的个体

在对2秒时间间隔内的帧进行平均之后，对健康个体的视频记录进行分析，如图13a所示。在图13b中叠加了脉管分段的结果。在该实验中，总脉管长度的31％需要手动交互。对该图像中的脉管进行颜色编码，并标记[5(红色)，650(白色)]μm范围内的血流速度。将具有由于图像对比度不够而不能被分析的时空图的脉管分段和具有表现出大于物理极限(式7)的速度的线结构的脉管分段标记为黑色。由于用于自动检测的尺度，不对背景中的厚脉管自动分段。

图14a给出了机能微循环密度分布。如图所示，大部分图像区域由5-10μm范围内的毛细管占据。

总共对207个脉管分段进行了分析。在99个(48％)段中，时空图未显示可见的线结构。在其它情况下，脉管段太短而不能允许式7界定的可用范围内的速度分析。其余108个(52％)脉管段的时空图显示了交互追踪的线结构。图14b的速度分布显示出，在给定厚度范围d＜60μm下，血流速度处于相同的量级。然后对这108个脉管中的血流速度进行自动分析。利用灰度级霍夫变换，29个段(27％)落入了20％的误差接受级别内。在本实验中，张量方法成功地分析了15个段。

心脏脱位

图15A和15C显示了在心脏脱位前和脱位中的、舌下视频记录中的250个帧(10秒)的平均。右侧的图(图15B和15D)示出了叠加有分析结果的相同的视频数据。在这两个实验中，总脉管长度的95％(脱位前)和80％(脱位时)被自动分段，其余脉管被交互地添加。图16示出了由图15A中的箭头标记的脉管的时空图，并示出了心脏脱位前和脱位时的区间。其清楚示出了血流速度在脱位过程中逐渐减小。在每个时空图中对多达20条线进行追踪，以获得对每10秒区间内的平均速度的印象。时空图示出了脱位前44％的脉管分段的清楚的线结构，以及脱位中48％的脉管分段的清楚的线结构。这代表了两种情况下大约75％的分段的脉管长度。再次如先前实验中所述将图15B和15d中的脉管颜色编码。可清楚地看到，血流速度在心脏脱位过程中减小。图17B提供的速度分布显示了在整个厚度范围内，平均血流速度减小为原始速度的大约三分之一。

在脱位实施例中观察到的由脉管占据的图像区域由17.1％降至14.6％，减少了15％。通过研究图15所示图像，确认了这一发现。如图所示，在心脏脱位过程中，某些小脉管是不可见的，即，红血球缺乏或减少。如所期望的，除了45-50μm范围内的脉管之外，密度分布的形状大致相同(见图17A)。这可能由于脱位过程中观察到的红血球稍微减少以及分布的离散特性而造成，观察到的红血球的减小导致观察到的脉管厚度减小。在单一厚脉管的分段可发现另一误差源，该脉管从左下向上延伸且厚度处于可用范围内。该脉管被具有更高对比度的另一脉管重叠，妨碍了图15b所示的厚脉管的分段。

本发明提供了先进的图像分析技术，其将单一尺度的自动分析与四个可选尺度(例如，σ＝1.5、3.0、6.0、12.0)中任意尺度的交互分析相结合。自由改变局部脉管分段的尺度增加了检测范围。

长度验证实验示出了主要由于阶梯边缘脉管端部而引起的误差，该阶梯边缘脉管端部不指示实际脉管并且在对实际微循环图像进行分析之后不会被观察到。考虑到实际脉管的长度(平均长度约为100像素)，模拟脉管的绝对长度误差产生很小的相对误差(对于宽度高达5像素的毛细管而言，该误差小于5％)。长度估计还依赖于中心线敏感度参数λ_thr。减小阈值参数将增大线检测敏感度，并导致对脉管长度的过度估计。当脉管临时移动到焦平面外时，这可为有利的。实验显示，长度估计还受到分析尺度、线检测阈值和方向的影响。在将分析尺度设置为σ＝1.5像素时，模拟毛细管脉管段长度(宽度高达15像素)被估计为5个像素。在该小分析尺度下，边缘检测主要受图像噪声的影响。为此，最好以更大的分析尺度交互地对更厚的脉管(宽度大于12像素)进行分析。MicroScan图像中，厚脉管的片段相对较小(在长度上)，并因此要求有限的用户交互。

对于宽度为[1，20]像素的期望范围内的脉管，线厚度被检测为具有约为±1像素范围内的绝对误差。脉管方向的作用使该误差增加了1-2像素，尤其是对于宽度大于13像素的脉管。对于宽度大于例如标准MicroScan图像中的3像素的脉管而言，这产生了小于20％的相对厚度估计误差。自由选择交互分析的尺度还引入了不利条件。在大尺度下分析小脉管导致对脉管厚度的过度估计(图11a)。另一方面，利用小尺度工具分析厚脉管导致对脉管中心线及其边界的参差不齐的检测，并产生对脉管长度的过度估计。对分析结果的视觉检查和操作者的理解对分析结果的估计也是重要的。

交互速度估计非常精确，在[0，1000]m/s的情况下，平均速度偏差小于4％。脉管方向使该速度误差增大1％。该速度误差还依赖于操作者的技术和追踪速度线的精确性。在将来的研究中测试用户对速度误差的贡献是重要的。对时空图的自动分析是优选的，因为其进一步减小了观察者偏差。模拟视频显示，利用霍夫变换自动分析时空图对高达750像素/秒(±5％)的速度是可行的。对于高达100像素/秒的速度且具有相似精度(±5％)的情况，执行结构张量的效果次于霍夫检测。如果视频记录对比度大、血流为基本恒定的(对于毛细管而言通常如此)、且具有足够的等离子体间隙以在时空图中生成清晰的线，那么，对从临床获得的时空图的自动分析将是成功的。对于将来的研究，建议进一步增大现有SDF系统的图像对比度，这将改善对时空图的自动估计。脉动血流生成弯曲的时空图，其不能利用上述的灰度级霍夫方法进行分析，该方法依赖于直线。结构张量方法的局部范围产生对速度的平均测量，并可适用于这些类型的血流。

大于v_max(式7)的血流速度不能用时空图检测出。在短脉管段和低的视频帧速率下，该范围尤其受限。因此，速度检测将受益于高的帧速率。本文显示了临床采样中的速度检测适用于具有高对比度的视频记录，即使是在足够长的脉管段中使用标准视频帧速率(25Hz)。

到目前为止，未显示出由患者的心脏脱位引起的阻塞性休克对组织的适当灌注有害。利用SDF成像，对心脏脱位对舌下微循环的微脉管血流动力学的直接影响进行离线观测和分析。在由心脏脱位引起的具有严重低血压的休克情况发生期间，在给定的实验中，舌下微循环灌注减少至原始血流速度的三分之一。SDF视频图像显示了某些毛细管脱落(不灌注红血球)，而其它图像示出了血流的减小，一些图像甚至示出血流在心脏脱位期间完全停止。

可以总结，新描述的上述脉管图像分析方法是有价值的工具，其进一步减少了用户交互和观察者偏差，并使我们能够确定通过其它方式不能获得的脉管参数分布。

以上参照大量示例性实施方式对本发明进行了解释。本领域技术人员将理解，可对其进行各种变形和修改，而不偏离由权利要求限定的本发明的范围。

参考文献

1.Acton ST，Wethmar K，Ley K，“Automatic tracking of rollingleukocytes in vivo(体内滚动白细胞的自动跟踪)”，MicrovascularResearch(2002)，63：139-148.

2.Altman DG，“Practical statistics for medical research(医学研究的实际统计)”，Chapman&Hall，USA 1999，ISBN：0-412-27630-5，p611.

3.Biberthaler P，Langer S，“Comparison of the new OPS imagingtechnique with intravital microscopy：Analysis of the colonmicrocirculation(新的OPS成像技术与活体显微镜方法的比较：结肠微循环的分析)”，European Surgery Research(2002)，24：124-128.

4.Boerma EC，Mathura KR，Van der Voort PHJ，Spronk PE，Ince C，“Quantifying bedside-derived imaging of microcirculatory abnormalitiesin septic patients：a prospective validation study(败血症患者的微循环异常的床旁成像的量化)”，Critical Care(2005)，9(6)：R601-R606.

5.Canny J，“A computational approach to edge detection(边界检测的计算方法)”，IEEE Transactions on Pattern Analysis and MachineIntelligence(1986)，PAMI-8(6)：679-698.

6.Dobbe JGG，Streekstra GJ，Hardeman MR，Ince C，GrimbergenCA，“Measurement of the distribution of red blood cell deformabilityusing an automated rheoscope(利用自动验电器测量红血球变形性的分布)”，Clincal Cytometry(2002)，50：313-325.

7.Ellis CG，Ellsworth ML，Pittman RN，Burgess WL，“Applicationof image analysis for evaluation of red blood cell dynamics in capillaries(在毛细管中进行估计红血球动态的图像分析)”，MicrovascularResearch(1992)，44：214-225.

8.GeusebroekJM，“Robust autofocusing in microscopy(显微镜检查中的鲁棒自动聚焦)”，Cytometry(2000)，39：1-9.

9.Gonzalez RC，Woods RE，“Digital image processing(数字图像处理)”，Addison-Wesley publishing company，Massachusetts 1992，ISBN：0-201-50803-6.

10.Groner W，Winkelman JW，Harris A，Ince C，Bouma GJ，MessmerK，Nadeau RG，“Orthogonal polarization spectral imaging：A new methodfor study of the microcirculation(正交极化光谱成像：一种微循环研究的新方法)”，Nature America(1999)，5(10)：1209-1213.

11.Ter Haar Romeny BM，“Front-end vision and multi-scale imageanalysis：Multi-scale computer vision theory and applications，written inMathematica(前端视觉和多尺度成像分析：数学上的多尺度计算机视觉理论和应用)”，Kluwer Academic Publishers，Dordrecht TheNetherlands 2003，ISBN：1-4020-1507-0，p466.

12.Horn BKP，Schunck BG，“Determining optical flow(确定光流)”，Artificial Intelligence(1981)，17：185-203.

13.Ince C，“The microcirculation is the motor of sepsis(review)(微循环是败血症的马达)”，Critical Care(2005)，9(suppl 4)：S13-S19.

14.B，“Digital image processing(数字图像处理)”，6threvised and extended edition，ISBN 3-540-24035-7，Springer，BerlinHeidelberg New York，p607.

15.Japee SA，Ellis CG，Pittman RN，“Flow visualization tools forimage analysis of capillary networks(用于毛细管网络的图像分析的流视觉化工具)”，Microcirculation(2004)，11：39-54.

16.Japee SA，Pittman RN，Ellis CG，“A new video image analysissystem to study red blood cell dynamics and oxygenation in capillarynetworks(研究毛细管网络中的红血球动态和氧化的新的视频图像分析系统)”，Microcirculation(2005)，12：489-506.

17.Klyscz T，Jünger M，Jung F，Zeintl H，“Cap Image-einneuartiges computerunterstütztes Videobildanalysesystem für diedynamische Kapillarmikroskopie(Cap Image：一种分析动态毛细管的新型的计算机视频图像分析系统)”，Biomedizinische Technik(1997)，Band 42 Heft 6：168-175.

18.Langer S，Biberthaler P，Harris AG，Steinau HU，Messmer K，“Invivo monitoring of microvessels in skin flaps：Introduction of a noveltechnique(皮瓣中微脉管的体内监控：一种新技术的介绍)”，Microsurgery(2001)，21：317-324.

19.Lentner A，Berger F，Wienert V，“Das“Spatial Shift Alignment(SSA)”-eine neue Methode zur Bestimmung derBlutflussgeschwindigkeit in de Video-Kapillarmikroskopie(“空间移位对齐(SSA)”：测定毛细管视频中血流速度的新方法)”，BiomedizinischeTechniek(1994)，Band 39 heft 7-8：170-175.

20.Lindert J，Werner J，Redlin M，Kuppe H，Habazettl H，Pries AR，“OPS imaging op human microcirculation：a short technical report(对人体微循环的OPS成像：技术简报)”，Journal of vascular research(2002)，39：368-372.

21.Lo RC，Tsai WH，“Gray-scale Hough transform for thick llinedetection in gray-scale images(灰度级图像中用于厚线检测的灰度级霍夫变换)”，Pattern Recognition(1995)，28(5)：647-661.

22.MicroVision Medical，Meibergdreef 45，1105 BA Amsterdam，The Netherlands.www.microvisionmedical.nl.

23.Mathura KR，Vollebregt KC，Boer K，De Graaff JC，Ubbink DT，Ince C，“Comparison of OPS imaging and conventional capillarymicroscopy to study the human microcirculation(OPS成像和传统毛细管显微镜方法在研究人体微循环中的比较)”，Journal of AppliedPhysiology(2001)，91：74-78.

24.Nolte D，Zeintl H，Steinbauer M，Pickelmann S，Messmer K，“Functional capillary density：an indicator of tissue perfusion？(机能毛细管密度：组织灌注的指示？)”，International Journal of Microcirculation(1995)，15：244-249.

25.Pennings FA，Bouma GJ，Ince C，“Direct observation of thehuman cerebral microcirculation during aneurysm surgery revealsincreased arteriolar contractility(动脉瘤手术中直接观察人体大脑微循环展示了增加的小动脉收缩性)”，Stroke(2004)，35：1284-1288.

26.Pries AR，“A versatile video image analysis system formicrocirculatory research(用于微循环研究的多用途视频图像分析系统)”，Int J Microcirc Clin Exp(1988)，7：327-345.

27.Russ JC，“The image processing handbook(图像处理手册)”，第四版，ISBN 0-8493-1142-X，CRC Press LLC，Florida(2002)，p732.

28.Slaaf，DW，Jeurens TJM，Tangelder GJ，Reneman RS，Aarts T，“Methods to measure blood flow velocity of red blood cells in vivo at themicroscopic level(以显微镜级别测量体内红血球血流速度的方法)”，Annual review of Biomedical Engineering(1986)，14(2)：175-186.

29.Staal J，Abràmoff MD，Viergever MA，Van Ginneken B，“Ridge-based vessel segmentation in color images of the retina(视网膜彩色图像中基于棱线的脉管分段)”，IEEE Transactions on MedicalImaging(2004)，23(4)：501-509.

30.Spronk PE，Ince C，Gardien MJ，Mathura KR，Oudemans-Vanstraaten HM，Zandstra DF，“Nitroglycerin in septic shock afterintravascular volume resuscitation(脉管内体积复苏后败血病性休克的硝化甘油)”，The lancet(2002)，360：1395-1396.

31.Steger C，“An unbiased detector of curvilinear structures(曲线结构的无偏检测器)”，IEEE Transactions on pattern analysis and machineintelligence(1998)，20(2)：113-125.

32.De Vriese AS，Verbeuren RJ，Vallez MO，Lameire MH，DeBuyzere M，Vanhoutte PM，“Off-line analysis of red blood cell velocity inrenal arterioles(肾微动脉中红血球速度的离线分析)”，Journal ofVascular Research(2002)，37：26-31.