1 引言

大多数软组织生物的力学性质依赖于胶原纤维网络和细胞外基质的组成成分。在动物体中,胶原蛋白主要以长纤维形式存在,坚韧的胶原纤维束是细胞外基质的主要组成成分。而胶原纤维方向提供了生物病变的依据,准确获取胶原纤维的方向分布和数量信息,可以为诊断病变组织,评估治疗反应、局部损伤以及监测组织发展状况等提供潜在的生物标记物。已有大量的用于二维图像纤维方向检测的算法,如傅里叶变换方法^[1-3]、霍夫变换^[4-5]、主成分分析^[6-8]等。这些二维算法已经验证了利用纤维方向评估胶原组织的可行性。三维分析比二维能够更准确地表现胶原纤维的结构,最近的研究表明,对于形态的描述,三维分析更精确有效^[9-11]。

当前,虽然有许多二维胶原纤维方向分析的方法,但是只有少量的三维分析技术用于纤维图像堆栈中检测三维纤维的方向。例如:光纤跟踪方法^[12-15]通过能量最小化来定义和跟踪纤维对象,但算法复杂度与图像复杂度相关,当应用于稠密的纤维物体时,随着迭代次数增加,该算法的可行性降低,且需要一定程度的用户交互;基于傅里叶变换的方法通常用于二维图像分析,现已扩展到三维纤维的分析中^[16],但需要创建一组感兴趣区域,并求出每个感兴趣区域的整体方向,然而,纤维区域的确定需要折中算法的准确性和消耗时间;还有一些其他的算法,比如弦长变换^[17]、高斯函数方向空间^[18-19]和惯性矩^[20]等方法,然而,这些方法都存在一定的局限性。当前最有效的方法是向量加权算法,该算法能够检测二维图像中每个像素点的纤维方向。该算法主要依赖于纤维不同方向上的图像强度的变化来获取纤维的方向^[21-26],是研究人员根据纤维方向上强度变化很小而纤维直径方向上强度变化剧烈的特点提出的一种加权算法。该方法比傅里叶算法具有更高的精确度和时间复杂度,但时间复杂度与图像大小的相关性较大。

为了准确高效地提取三维纤维方向,本文提出一种基于facet模型的三维纤维方向检测方法。首先,阐述了本文算法检测二维纤维方向的原理;然后,将二维情况扩展到三维;接着,分析算法的性能和准确性;最后,通过模拟已知方向的纤维图像验证算法的复杂度,并评估算法的准确性。

2 纤维方向检测方法

2.1 二维纤维方向检测方法

定义纤维方向为图像内每个像素在其邻域窗口内强度变化最缓慢的方向。一幅图像是一个连续函数的采样阵列,使用离散逼近法计算所需要的性质,并在子像素分辨率上计算纤维方向的相关信息,能够使最终结果更准确。

一幅数字图像可以用对应的离散三维图形表示,其中纵坐标表示像素灰度值。设连续图像函数为z=f(x,y),利用数字图像的灰度值重构连续函数。在处理复杂图像时,连续强度函数可能包含x和y的超阶幂方,因而重构原始函数复杂度高。将图像模型分解为简单的分段函数,复杂任务即可分解为重构简单分段解析函数。基于facet的模型就是试图找到每一像素点邻域内最接近其邻域强度值的简单函数^[27]。连续图像强度函数可在每个像素点邻域内得到一幅局部逼近的n×m图像,因而能得到n×m个逼近函数,每个函数仅对图像中的一个特定像素有效。这些函数(非像素值)可以确定边缘位置,而许多复杂解析函数可用来逼近图像强度。对于简单图像,强度函数的逼近可使用分段常量或分段双变量线性函数;而对于复杂区域图像,则要选择二次、三次甚至高次幂函数。本文采用三次多项式建立图像邻域的facet模型,即

f(x,y)=m1+m2x+m3y+m4x2+m5xy+m6y2+m7x3+m8x2y+m9xy2+m10y3,(1)

式中:(x,y)为相对于5×5窗口中心点的坐标;m_i为拟合系数,i=1,2,…,10。利用最小二乘法,通过奇异值分解计算方程中逼近函数的系数,由于该方法的时间复杂度高,故采用模板计算系数方法。在离散整数空间{-w,-w+1,…,0,1,…,w}×{-w,-w+1,…,0,1,…,w}域中,x∈{-w,-w+1,…,0,1,…,w},y∈{-w,-w+1,…,0,1,…,w},w为正整数,考虑一组二元离散正交多项式集{1,x,y,x²-b,xy,y²-b,x³-dx,(x²-b)y,(y²-b)x,y³-dy},其中b、d是与w有关的待定系数。则有

f(x,y)=k1+k2x+k3y+k4(x2-b)+k5xy+k6(y2-b)+k7(x3-dx)+k8(x2-b)y+…+k9x(y2-b)+k10(y3-dy),(2)

且有

∑x=-ww(x2-b)×1=0,(3)

则

b=∑x=-wwx22w+1,(4)

又因为

∑x=-ww(x3-dx)·x=0,(5)

所以

d=∑x=-wwx4∑x=-wwx2,(6)

式中:k_i为拟合系数,i=1,2,…,10。

给定w之后,可求出相对应的b、d。通过正交多项式的特性求取系数k_i,即多项式中任意两个不同项的乘积在空间域的累加和为0。其中Ω为空间域{-w,-w+1,…,0,1,…,w}×{-w,-w+1,…,0,1,…,w},p_i为正交多项式集{1,x,y,x²-b,xy,y²-b,x³-dx,(x²-b)y,(y²-b)x,y³-dy}的第i项。则有

∑Ωf(x,y)pi=∑Ωkipi2,(7)ki=∑Ωf(x,y)pi∑Ωpi2,(8)

(1)式中的m为k的线性组合。所以有

m2=k2-d·k7-b·k9,(9)fx(x,y)=m2,fy(x,y)=m3,(10)

函数f(x,y)对x的偏导为m₂,对y的偏导为m₃。所求的纤维方向与梯度方向正交,则计算纤维方向的公式为

θ=arctanm3m2-90°。(11)

当w=2时,可以使用窗口为5 pixel×5 pixel的模板逼近函数系数,具体方法为利用模板和待逼近图像邻域对应位置的像素点相乘并求和。

利用方向不同、长度不同、宽度不同的多条线段模拟二维纤维网络,分别采用本文算法和向量加权算法计算误差。得到的方向灰度图(通过灰度描述纤维在空间中的分布角度0°~180°)如图1所示。向量加权算法得到的方向误差最小为4.93°(窗口为9 pixel×9 pixel);本文算法得到的方向误差最小为3.59°(窗口为5 pixel×5 pixel)。本文算法的准确性略高于向量加权算法。从图1中可以看出,两种算法的主要误差集中在线段的末端和交叉处,而一般情况下,线段的长度越长,误差越小。

图 1. 仿真结果。(a)模拟二维纤维图;(b)实际方向灰度图;(c)向量加权算法;(d)所提算法

Fig. 1. Simulation results. (a) Simulation of two-dimensional fiber; (b) grayscale in th actual direction; (c) vector weighted algorithm; (d) proposed algorithm

下载图片 查看所有图片

2.2 三维纤维方向检测方法

对于二维纤维图像,上述算法循环遍历图像即可求出纤维图像中每个像素点的方向。当图像由二维扩展到三维时,本文算法要利用每个体素在一个三维纤维图像堆栈的强度变化检测方向。通常,利用横切面上的一个角θ和倾角ϕ定义三维空间中的一个方向,如图2所示,覆盖所有三维方向范围的θ和ϕ均为0°~360°。图2中,实线代表一个感兴趣的纤维,θ为纤维在xy平面的投影和x轴的夹角,ϕ为纤维和z轴之间的夹角。本文算法定义了两个额外的角度—β和γ,β为纤维在zx平面的投影和x轴之间的夹角,γ为纤维在yz平面的投影和y轴之间的夹角。虚线为计算纤维方向。

图 2. 4个角在三维空间的示意图

Fig. 2. Diagram of four angles in three-dimensional space

下载图片 查看所有图片

在三维空间中,纤维方向由θ和ϕ确定。θ、β和γ可以通过投影求取,即将三维纤维方向的检测转换到二维平面,ϕ可以利用β和γ得到,计算公式为

tan2ϕ=1tan2β+1tan2γ。(12)

本文算法利用一个立方体窗口计算窗口中心体素的方向(计算ϕ、β和γ三个角)。每一个体素的方向即为强度变化最小的方向。将它们投射到xy、yz和xz平面上,分别求取θ、β和γ、ϕ。为了节省计算时间,不再通过计算立方体窗口内所有方向向量来确定立方体中心体素的方向,而是利用三个正交平面中的一个平面对立方体进行平均,将立方体窗口缩小为正方形窗口,然后计算正方形中心像素的二维纤维方向。这样即可将三维矩阵数据处理转化为二维矩阵数据处理,大幅减少了计算时间。

3 技术性能评价

使用MATLAB软件生成已知方向的三维纤维模拟图,对比本文算法计算出的方向与模拟设定方向,求出每个体素的平均方向误差,并与其他算法的结果进行对比,评价算法的可靠性。实验硬件为Intel i7-6700 CPU,16 GB内存,软件环境为Windows 10操作系统,MATLAB R2016b软件平台。

模拟的三维纤维提供了一个理想的纤维分布模型,以此验证算法的有效性。首先评估纤维粗细以及算法选择的窗口与误差之间的关系。模拟生成7个粗细不同的三维纤维图,纤维直径分别为3,5,7,9,11,13,15 pixel。在三维空间中利用圆柱体模拟三维纤维,圆柱体中轴线上的灰度值最大,中轴线方向即纤维方向,沿与中轴线垂直的方向,灰度值逐渐减小。圆柱体在三维空间随机散乱分布,即θ和ϕ随机分布于0°~180°,记录所有模拟方向作为纤维的实际方向。圆柱体可以相距很近,甚至只有一个体素点的间隔,但不能相交。因为如果圆柱体相交,就无法确定这些体素属于哪个纤维。7个模拟纤维图中纤维方向位置、纤维数量完全相同,区别仅在于纤维半径大小不一样。图3给出了纤维直径分别为3,9,15 pixel的模拟纤维图。

图 3. 三维模拟纤维。(a)直径为3 pixel;(b)直径为9 pixel;(c)直径为15 pixel

Fig. 3. Three-dimensional simulated fiber. (a) Diameter of 3 pixel; (b) diameter of 9 pixel; (c) diameter of 15 pixel

下载图片 查看所有图片

计算所得方向误差与窗口的选择有关。为了研究窗口大小对误差的影响,利用不同大小的窗口求取同一幅三维纤维图的方向。分别模拟直径为3,5,7,9,11,13,15 pixel的7幅纤维图,每幅模拟纤维图的纤维方向、长度和数量都保持一致,仅纤维直径不同。算法窗口大小为w×w×h,其中w为5~15,h固定为11。图4给出了不同纤维直径和不同窗口对应的方向测量误差。

图 4. 不同窗口和纤维直径所对应的方向测量误差。(a) ?的误差;(b) θ的误差

Fig. 4. Direction measurement errors corresponding to different windows and fiber diameters. (a) Error of ?; (b) error of θ

下载图片 查看所有图片

从图4中可以看出,对于不同直径的纤维,选择合适的窗口能够减小误差。对于直径为3 pixel和5 pixel的纤维,7 pixel×7 pixel的窗口得到的测量误差最小;对于直径为7 pixel的纤维,9 pixel×9 pixel的窗口得到的测量误差最小;对于9 pixel和11 pixel的纤维,分别使用11 pixel×11 pixel和13 pixel×13 pixel的窗口最合适;对于13 pixel和15 pixel的纤维,15 pixel×15 pixel的窗口可以得到最准确的测量结果。从图4还可以看出,随着纤维直径的增加,测量误差呈现上升趋势。因为随着纤维直径的增加,要选择更大的窗口才能减小误差,但是大窗口又包含了更多邻近纤维的信息,引入了更多误差,所以纤维直径越大,误差越大。

虽然图4仅给出了h=11时的结果,但是通过实验仿真可以发现,改变h不会改变曲线的走势;对于不同直径的纤维,最合适的w也不会发生改变。例如,对于直径为5 pixel的纤维,h取不同值时,依然会在w=7时得到最小的测量误差。实际上,h仍然影响测量误差,h和w可以视为两个影响测量误差的独立变量。因此,对于不同直径的纤维,可以先确定一个合适的w,然后改变h,找到该直径下测量误差最小的h。图5给出了不同纤维直径下的最小测量误差。

图 5. 所提算法与向量加权算法的测量误差对比

Fig. 5. Measurement error comparison between proposed algorithm and vector weighted algorithm

下载图片 查看所有图片

从图5可以看出:求取模拟纤维方向时,本文算法精度高于向量加权算法;随着纤维直径的增加,误差呈线性增长。表1给出了不同纤维直径对应的w和h。

为验证纤维数量与测量误差的关系,生成4幅纤维直径为7 pixel,纤维数量分别为30、50、70、90的模拟纤维图,并给出在9 pixel×9 pixel×11 pixel窗口下的测量误差,如图6所示。

图 6. 纤维数量与测量误差关系

Fig. 6. Relationship between number of fibers and measurement error

下载图片 查看所有图片

从图6可以看出,θ和ϕ误差随着纤维数量的增加而缓慢增加。由于在给定空间内随着纤维数量的增加,纤维间的间隔越来越小,窗口包含了更多邻近纤维的信息,从而导致误差增加。

除了准确性,算法效率也是评估算法性能的一个重要因素。考虑到本文算法所消耗的时间与窗口大小有关,因此,利用实验分析窗口大小和时间的关系,并和向量加权算法进行比较,结果如图7所示。

图 7. 所提算法与向量加权算法的耗时对比

Fig. 7. Consuming time comparison between vector weighted algorithm and proposed algorithm

下载图片 查看所有图片

从图7可以看出:本文算法消耗的时间随窗口的增大呈现缓慢的线性增长,而向量算法呈现指数级递增趋势;在窗口为21 pixel之前,向量加权算法的效率略优于本文算法;窗口选择大于21 pixel后,向量加权算法的效率大幅下降。实际上对于直径为3~15 pixel的纤维,向量加权算法的最优窗口为11~27 pixel,而本文算法的最优窗口为7~15 pixel,因此,综合考虑在不同纤维直径下选择最优窗口。在平均消耗时间方面,本文算法优于向量加权算法。

为进一步验证算法的准确性,创建一个包含100根方向各异纤维的200 pixel×200 pixel×200 pixel的三维纤维图,纤维直径为5~11 pixel,纤维长度各异,如图8所示。本文算法使用11 pixel×11 pixel×13 pixel窗口得到的θ和ϕ最小误差为9.41°和11.90°,消耗时间仅为88 s;而向量加权算法使用27 pixel×27 pixel×27 pixel窗口得到的θ和ϕ最小误差为11.8804°和15.3402°,所用时间为198 s,该结果进一步验证了本文算法的准确性和高效性。

图 8. 100根方向、长度、直径各异的模拟纤维

Fig. 8. 100 simulated fibers with different directions, length, and diameters

下载图片 查看所有图片

4 结果与讨论

纤维的厚度和密度是影响纤维方向检测精度的重要因素。对于大多数窗口,随着纤维厚度的增加,测量误差呈缓慢线性增长的趋势。对于特定的纤维厚度,平均误差并不随窗口增加(或减小)而单调地增加(或减小),而是存在一个极小值。例如,对于纤维直径为7 pixel的纤维,窗口为9 pixel×9 pixel×11 pixel时,能够提供最准确的纤维取向评估,窗口进一步增加或减小,误差都会增加。这种情况发生的原因是,纤维之间的邻近性也是影响方向测量的一个因素,当选择的窗口越来越大时,更多的邻近纤维被包含,从而导致测量误差增加;同理,当纤维密度增大时,邻近的纤维也更容易包含进当前窗口,进而引起测量误差增大。这些结果表明,纤维直径是影响测量准确度的主要因素,因此,所选择窗口的尺寸略大于纤维直径时才能达到最优的结果。

算法消耗时间与三维图像的大小以及体素窗口的大小有关。一个较大的窗口或者一个更大的三维图像所需的时间更多,消耗的时间随窗口增大呈缓慢线性增加的趋势。传统算法提供每个感兴趣区域(ROI)的平均方向,而不是每个体素。与传统算法相比,本文算法不仅提高了结果的精确度,而且提高了计算效率;与向量加权算法相比,尽管向量加权算法在使用小窗口时效率略高于本文算法,但随着窗口的增大,向量加权算法耗时呈指数级增长趋势,而本文算法呈线性缓慢增长。实际上向量加权算法需要大窗口才能获得最优的测量结果。例如,对于直径为3~15 pixel范围的纤维,向量加权算法需要11 pixel×11 pixel×11 pixel到27 pixel×27 pixel×27 pixel之间的窗口才能获得最优的测量结果。因此,考虑到测量一系列不同直径的纤维,本文算法所需要的平均时间相比向量加权算法大大缩短。

5 结论

提出一种基于facet模型的三维纤维方向检测方法,将三维纤维方向检测转换为二维平面检测,降低了计算复杂度,并基于facet模型进行了三维纤维方向检测,实现了不同直径纤维的准确检测。实验结果表明,相比于现有的典型算法,本文算法可快速准确地获取不同直径三维纤维的方向数据,模型窗口尺寸影响最优检测结果。如何根据生物组织三维纤维直径选择合适的窗口,达到更好的检测效果,将是下一步研究的重点。