0引言分形最早由Mandelbrot(1977)，Mandelbrot和Wheeler(1983)提出，是具有高度几何复杂性的数学集合，可以模拟许多自然现象，在图像处理(Novianto等，2003)、物理学(Zborovský，2018)、医学(Marchette等，1997；Fortin等，1992；Jiang等, 2008)、地质学(Ai等，2014；Cai等，2018；Cai，2014)等领域应用广泛，成为描述和定量研究复杂的、不规则几何特征的重要工具。Tolle等人(2003)总结了理论分形维数值(fractal dimension，FD)的计算方法，主要包括计盒法、小波变换法、功率谱法、赫斯特系数法、布利冈—闵可夫斯基法、变分法和容量维数法等。这些方法中，由于计盒法效率高，FD感应范围广，逐渐成为最常用的图像分形维数计算方法之一。计盒法计算为 1 ${D_b} = \mathop {{\rm{lim}}}\limits_{r \to 0} \frac{{{\rm{ln}}\left({{N_r}} \right)}}{{{\rm{ln}}(1/r)}}$ 式中，$D_{b}$为盒维数，$r$为盒子尺度大小，$N_{r}$表示尺度$r$的盒子覆盖整个图像所需最少的盒子数量。计盒方法最早由Gagnepain和Roques-Carmes(1986)提出，随着计盒法的广泛使用，研究人员对计盒法进行了大量的分析与研究，Voss(1986)通过引入概率论对计盒方法进行了改进，Keller等人(1989)通过线性插值对计盒方法进行了改进。Sarkar和Chaudhuri(1994)提出了差分计盒法(differential box-counting, DBC)计算分形维数；Jin等人(1995)对比多种分形维数计算方法，提出了相对计盒法(the relative differential box-counting, RDBC)，精度相对于DBC更高，效率与DBC相当；Sandau和Kurz(1997)利用分形的最大不变性，提出了扩展盒维数法(the extended counting method)，用局部最大盒维数来反映整个图像的复杂度。Buczkowski等人(1998)对拟合点序列进行了改进，提高了拟合点数量，提高了计盒法精度；在此之后，Foroutan-pour等人(1999)确定了普遍适用的盒子尺度范围，并利用线的骨架进行$D_{b}$的计算，进一步提高了计盒法精度；Tolle等人(2003)通过模糊聚类算法及遗传算法确定了次优的盒子覆盖数量，该方法不需要拟合直线，相比BCM方法精度有所提高；Chen等人(2003)提出了SDBC(shifting differential box-counting)与SBC(scanning box-counting)方法；Li等人(2009)、Liu等人(2014)、Kaewaramsri和Woraratpanya(2015)、Lai等人(2016)以及Panigrahy等人(2019)分别在DBC基础上提出了改进的DBC算法。需要强调的是，分形维数总是在极限条件下进行定义的，这意味着物体结构可以在任意高分辨率的条件下进行分析，在极限的情况下，FD值应与所选方法无关。然而在实践中，数字图像的结构是由有限的数据集组成。有界信息意味着测量总是被限制在有限的放大范围内，这与分形维数的定义不相符，导致现有分形维数计算方法都具有经验性，方法的优劣需要大量的实际经验来检验，并根据实际经验来优化各个参数。目前，已经对计盒法的尺度范围、拟合点选取、盒子覆盖方式以及线的厚度等进行了研究，并取得了积极进展(Buczkowski等，1998；Foroutan-pour等，1999)。在这些影响参数中，图像旋转这一问题并没有得到妥善解决，Sandau和Kurz(1997)首先研究了图像的旋转对FD的影响，根据分形维数的定义，满足分形条件的集合应具有旋转与平移不变性，但根据对Cantor尘埃的24幅图像(15°)计算分形维数，却得到了相差较大的结果；随后，Foroutan-pour等人(1999)认为：在用网格盒子进行覆盖时，应使网格外框的面积最小，即应取最小包容矩形，但文章中并没有提出位图旋转的方法。以往的方法无一例外都是在位图的基础上进行FD的计算，本文提出了一个将位图转化为矢量图，并在矢量图的基础上进行盒子覆盖，通过在矢量图的基础上旋转图形，避免了直接对位图进行旋转的插值问题，接下来采用多尺度的盒子覆盖矢量图形得到多个拟合点，最后按最小二乘法对拟合点进行拟合得到$D_{b}$。本文对3种类型的图像进行了对比分析，结果显示，与BCM方法相比，本文方法在3幅图的平均误差与变化幅度均有所降低。该方法充分利用了矢量图易旋转、无厚度的特点，减小了图像旋转对盒维数的影响，评价结果优于BCM方法。1BCM方法及其缺点BCM将二值位图看做2维空间中的一个平面$z=f(x, y)$，其中$x$，$y$分别对应图像各像素的位置，$z$对应图像像素的灰度值(0或1)，设物体结构对应灰度值为1，设图像大小为$M×N$，将$x-y$图像平面分为多个$r×r$大小的网格格子，格子坐标设为$(i, j)$，即 2 ${n_r}\left({i, j} \right) = \left\{ {\begin{array}{*{20}{c}}1&{{z_{{\rm{max}}}}\left({i, j} \right) = 1}\\0&{{z_{{\rm{max}}}}\left({i, j} \right) = 0}\end{array}} \right.{\rm{ }}$ 式中，$(i, j)$为网格格子的坐标，$n_{r}(i, j)$为单个格子的覆盖盒子数，$z_\text{max}(i, j)$为单个格子中像素灰度值的最大值。设$cl=\text{ceil}(M/r)$，$rw=\text{ceil}(N/r)$，其中ceil为向上取整，则$i≤cl$，$j≤rw$，接着计算总的覆盖盒子数量为 3 ${N_r} = \sum\limits_{i = 1}^{cl} {\sum\limits_{j = 1}^{rw} {{n_r}\left({i, j} \right)} } $ 式中，$N_{r}$为覆盖图像的盒子的总数量。改变盒子尺度$r$，选取合适的无标度区(尺度范围)，得到一组$(-\text{ln}\ r, \text{ln}\ N_{r})$序列，在无标度区内用最小二乘法计算$(-\text{ln}\ r, \text{ln}\ N_{r})$的斜率，该斜率就是图像的分形维数。接下来探讨图像旋转对FD值的影响，本文定义图像中物体结构的旋转角度$θ$为：取图像中物体结构的最小包容矩形，图像的最小包容矩形旋转角度$θ$后，能使最小包容矩形的各边与$x-y$图像的坐标轴平行，其中最小的锐角$θ$称为图像的旋转角度，$θ∈[π/4]$。用${\bf{R}}^N$来表示$N$维空间，设$\boldsymbol{X}_{1}$、$\boldsymbol{X}_{2}$分别为${\bf{R}}^N$的一个非空有界子集，因为分形维数是一种增量映射，若$\boldsymbol{X}_{1}⊆\boldsymbol{X}_{2}$，则$FD(\boldsymbol{X}_{1})≤FD(\boldsymbol{X}_{2})$(Sandau和Kurz，1997)。图 1显示了拍照成像的过程，图中左下方的蓝色矢量三角形，设其为集合$\boldsymbol{A}$，在无旋转的情况下用相机给物体拍照，图 1右图显示了拍照后所成的二值图像，若将单个像素看做一个方形面，记覆盖后生成的二值图像中物体结构的像素集合为$\boldsymbol{A}′$。用同一相机从另一角度对图像拍照(如图 2)，得到二值图像，设该图像结构的像素所组成的集合为$\boldsymbol{A}′′$。则存在$\boldsymbol{A}⊆\boldsymbol{A}′$与$\boldsymbol{A}⊆\boldsymbol{A}′′$，即 4 $\left\{ \begin{array}{l}FD\left(\mathit{\boldsymbol{A}} \right) \le FD\left({\mathit{\boldsymbol{A}}^{\prime }} \right)\\{\rm{ }}FD\left(\mathit{\boldsymbol{A}} \right) \le FD\left({\mathit{\boldsymbol{A}}^{\prime \prime} } \right)\end{array} \right.$ 图1 旋转角度$θ=0°$时，相机的成像过程 The imaging process of the camera when the rotation angle $θ=0°$Fig 1 图2 旋转角度$θ=9°$时，相机的成像过程 The imaging process of the camera when the rotation angle $θ=9°$Fig 2式中，$FD(\boldsymbol{A})，FD(\boldsymbol{A}′)，FD(\boldsymbol{A}′′)$分别为集合$\boldsymbol{A}, \boldsymbol{A}′, \boldsymbol{A}′′$的理论分形维数。通常情况下$\boldsymbol{A}′≠\boldsymbol{A}′′$，因此在使用计盒法时，不同旋转角度$θ$的图像具有不同的$D_{b}$值，这也解释了图像分辨率对计盒法的影响，当分辨率越低，覆盖图像所用的盒子就越大，像素点也越少，$D_{b}$值与真实值的误差也越大。为了缩小图像$\boldsymbol{A}$与原物体$\boldsymbol{A}′$之间的差异，应使原物体所占的像素点数量越少越好。以图 1为例，对于同一相机，当图像中物体结构的分辨率(单个像素点所代表真实范围)一定，即相机与物体在$Z$方向距离相同时，只能调节相机(或物体)的$X、Y$方向的位置和旋转角度，设单个像素点所代表的真实范围$σ$的长度为$l$(假设为一个方形区域)，如果保持$θ$一定，通过移动相机位置对原物体拍照，根据相机的位置不同，所成图像将分别在$X$、$Y$方向呈周期为$l′=scale×l$的变化，其中$scale$为生成的图像相对于真实物体的缩小或放大倍率，由于$l′$通常都极小，要在$σ$范围内找到合适的位置使得物体结构所占像素点最少是相当困难的，幸运的是，现在相机的分辨率普遍较高，相机位置对于$D_{b}$值的影响很小，通常情况下可以忽略不计。接下来探究$θ$对于$D_{b}$值的影响，图 1和图 2显示了同一物体不同旋转角度下所成的二值图像，它们具有不同的形状结构，这也意味着它们具有不同的$D_{b}$值。为了进一步探究旋转对于$D_{b}$的影响，将5次迭代的科赫海岸线(矢量图，图 3(a))看做原物体，从不同旋转角度对图像进行拍照，结果如图 4所示，然后改用不同分辨率的相机从多个角度对科赫海岸线进行拍照，即采用不同尺度$r$的盒子对矢量图进行覆盖，记录每次覆盖的盒子数量$N_{r}$，结果(如图 5所示)表明：当图像旋转后，虽然$N_{r}$变大，按照分形理论$D_{b}$值也应增大，但不同尺度$r$条件下$N_{r}$增量不相等，在大尺度条件下的$ \text{ln}\ N_{r}$增加幅度大于小尺度条件下的增加幅度，反而导致了$D_{b}$的下降；在旋转角度为0°时，$N_{r}$通常最小，这也印证了Foroutan-pour等人(1999)认为的应使图像旋转角度为0°的观点。然而，在拍照时很难保证旋转角度为0°，对于已经发生旋转的二值图像，如果直接在二值图像(设结构所占集合为$\boldsymbol{P})$上对结构进行旋转，设新生成的二值图像中结构的像素集合为$\boldsymbol{P}′$，由于二值图像的旋转算法不可避免地会使用插值，则通常$\boldsymbol{P}≠\boldsymbol{P}′$，这会导致$FD$值的误差进一步扩大，因此，解决图像旋转的问题对于计盒法准确性的提高具有重要意义。 图3 5次迭代的科赫海岸线 Koch coastline of 5 iterations((a) original vector image; (b) binary image; (c) skeleton of binary image)Fig 3 图4 相同分辨率(512×512像素)条件下，不同旋转角度的覆盖盒子数量 Number of cover boxes with different rotation angles at the same resolution (512×512 pixels)Fig 4((a) $N_{r}=340$；(b) $N_{r}=362$；(c) $N_{r}=358$；(d) $N_{r}=345$；(e) $N_{r}=354$；(f) $N_{r}=349$) 图5 不同角度、不同盒子尺度条件下5次迭代的科赫海岸线的覆盖盒子数量 Number of cover boxes of koch coastline with five iterations under different angles and box sizesFig 52旋转骨架法为了解决二值图像的旋转问题，提出一种新的分形维数计算方法，称之为旋转骨架法(skeleton rotating method，SRM)。首先提取图像的骨架，如图 6所示，图像大小为15×15像素的字符$“6”$的二值图像及其提取的骨架，图像左边显示了骨架提取的判定顺序、搜索范围等信息，右边为提取的骨架。设原物体、二值图像的结构(将单个像素点看做一个方形面)、骨架的集合分别$\boldsymbol{S}_{1}$、$\boldsymbol{S}_{2}$、$\boldsymbol{S}_{3}$。结合式(4)可知$\boldsymbol{S}_{1}⊆\boldsymbol{S}_{2}$，$\boldsymbol{S}_{3}⊆\boldsymbol{S}_{2}$，事实上，可以提取任意多种图像的骨架，本文采用的骨架是较为简单的一种，通常$\boldsymbol{S}_{1}≠\boldsymbol{S}_{3}$，则$FD(\boldsymbol{S}_{1})≠FD(\boldsymbol{S}_{3})$。提取骨架的过程也就是反推原物体结构的过程，虽然骨架与原物体会有一定差距，但在合适的尺度范围内却更能反映原物体的FD值，同时也能够避免二值图像旋转产生的插值问题。旋转骨架法的伪代码如下： 图6 二值图像骨架提取流程图 Skeleton extraction flowchart of binary imageFig 6$i=1$；#起始判定点横坐标$j=1$；#起始判定点纵坐标$M, N=\text{crop}(structure)$；#裁剪图像为结构所占大小，并将结构长宽赋值给$M$、$N$$\boldsymbol{lines}=set()$；#创建空集合，保存骨架线段While $(i≤M，j≤N)$ do begin  If $z(i, j)=1$   $count=0$；#统计满足条件的线段个数   For each point $(i′, j′)$ around $(i, j)$    If $z(i′, j′)=z(i, j)$ and $i′≤M$、$j′≤N$     $\boldsymbol{lines}.add([i-0.5, j-0.5, i′-0.5, j′-0.5])$；#满足条件的线段加入集合     $count+=1$；   End For   If $count=0$ #若范围内无满足条件的线段     $\boldsymbol{lines}.add([i-1, j-0.5, i′-1, j′-0.5])$；#孤立点提取骨架End While$\boldsymbol{lines}=Unique(\boldsymbol{lines})$；#去除多余线段该算法考虑到提取骨架连续性和复杂性等问题，同一个物体在观测尺度通常都是连续的，若只对像素点的上下左右4个方向进行搜索，不能保证图像中结构的连续性，造成图像结构的离散化，出现很多孤立点，降低原本结构的复杂度，影响计算精度。图 3显示了原物体(矢量图)、拍照所得二值图像以及二值图像的骨架(矢量图)，得到骨架之后，需要将其旋转一定的角度，使骨架的旋转角度为0°，其过程如图 7所示。首先确定骨架的最小包容矩形，使用python编程软件的inspyred库中的遗传算法进行最小包容矩形的计算，设置初始子代数量100，突变概率0.1，共繁殖10 000代，适应度为包容矩形的面积，求得最小包容矩形与图像的旋转角度$θ$(该设置条件下$θ$的变化幅度为±0.002°，满足精度条件)，然后旋转$θ$角度得到最终的骨架。最后，使用不同尺度$r$的盒子对骨架进行覆盖，记录每次覆盖的盒子数量$N_{r}$，得到一组$(-\text{ln}\ r, \text{ln}N_{r})$序列，在无标度区内用最小二乘法计算$(-\text{ln}\ r, \text{ln}N_{r})$的斜率，该斜率即为该骨架的分形维数。 图7 骨架的旋转示意图 Schematic diagram of skeleton rotation((a) extracted skeleton; (b) determining minimum containment rectangle of the skeleton; (c) rotating skeleton; (d) final skeleton)Fig 73图像分析及结果根据分形理论，分形可以分为两大类：规则分形与无规则分形。规则分形是具有严格自相似的几何图像，例如科赫曲线、康托集和谢尔宾斯基三角垫片等，除此以外的称为无规则分形。为了更好地比较SRM方法的实用性，通过BCM方法与旋转骨架法(SRM)对自相似分形图像、文字扫描图像以及植物图像进行了对比计算分析，其中自相似分形图像代表了规则分形图像，文字扫描图像与植物图像代表非规则分形图像，而文字扫描图像可更具体地代表人造物体的图像，植物图像代表自然物体的图像。3.1自相似分形图像分析自相似的分形图像是图 8中3幅带旋转的自相似分形图像D01(理论值1.585 0)、D02(理论值1.465 0)和D03(理论值1.261 9)，每一种图像有0°至44°旋转角度的图像，共计45幅图像，单幅图像的分辨率为2 048×2 048像素，图像生成方式采用python编程方式，通过特定大小的盒子覆盖原矢量图像得到。结合表 1和图 9的计算结果可以看出： 图8 不同旋转角度的自相似分形图像 Self-similar fractal images with different rotation angle((a) D01;(b) D02;(c) D03)Fig 8 表1 带旋转的自相似分形图像的分形维数计算结果 图像名称 理论分形维数 BCM 本文(SRM) 平均误差/% 变化幅度 平均误差/% 变化幅度 D01 1.585 0 2.615 1 0.069 9 1.889 9 0.012 6 D02 1.465 0 6.306 8 0.101 3 3.246 3 0.013 0 D03 1.261 9 3.407 2 0.095 7 1.108 7 0.009 8 Fractal dimension calculation results of self-similar fractal images with rotationTable 1 加粗字体表示最优结果。 图9 旋转角度0°~44°的D01、D02、D03图像分形维数计算结果 Fractal dimension calculation results of D01, D02 and D03 images with rotation angle of 0°~44°((a)D01;(b)D02;(c)D03)Fig 91) BCM计算的分形维数值具有较大偏差和变化幅度，相较之下，SRM的误差和变化幅度更小，$D_{b}$值更为准确；2) 旋转骨架法减少了旋转对图像$D_{b}$值的影响，但当$θ 2°$时，采用该方法可能会增大$D_{b}$值的误差，这是由于当$θ$角度过小，此时提取骨架的过程对$D_{b}$值的影响大于旋转角度$θ$对$D_{b}$值的影响，说明旋转骨架法适用于$θ$角度大于一定旋转角度的图像。3.2文字扫描图像分析除了自然景观以外，人造物体在日常生活中也随处可见，大都呈现较规则的几何形态，例如住房、汽车、飞机和文字等。图 10为4幅分辨率为512×512像素的文字扫描图像，其中T01是旋转角度为0.503°的图像，因其旋转角度小，将T01图像采用BCM求取的$D_{b}$值作为参考值，从表 2可以看出： 图10 不同旋转角度的文字扫描图像 Characters scanning images with different rotation angles((a) T01;(b) T02;(c) T03;(d) T04)Fig 10 表2 不同旋转角度的文字扫描图像计算结果 图像名称 BCM 本文(SRM) 分形维数 拟合误差 分形维数 拟合误差 T01 1.195 91 0.006 07 1.208 74 0.005 69 T02 1.175 89 0.006 36 1.201 14 0.005 05 T03 1.169 39 0.006 05 1.194 97 0.004 84 T04 1.186 00 0.007 07 1.198 37 0.004 85 Calculation results of text scanning images with different rotation anglesTable 2 加粗字体表示最优结果。1) 采用BCM时，同样也出现了图像旋转导致的偏差问题；2) 相对于参考值1.195 91，BCM方法的平均变化幅度为0.026 52，平均拟合误差为0.006 39，SRM的平均变化幅度为0.013 77，平均拟合误差为0.005 11。SRM方法计算出的分形维数与参考值最接近，同时分形维数的变化幅度与拟合误差也相对小，具有很高的稳定性。3.3植物图像分析自然物体是分形理论应用的重点，因此，对自然物体的图像进行分析具有重要意义，常见的自然物体包括河流、山脉、云、动物和植物等，甚至人的血管、大脑和虹膜都呈现出高度的分形特征。本文从叶子的正上方以不同距离(分辨率不同)、不同旋转角度对同一片叶子进行拍照，拍照结果如图 11所示，将图 11中叶子提取轮廓后如图 12所示。图 11中的F01与图 12中的F01*对应，其他依次对应。图 12中单幅图像的右下角显示了图像的分辨率。对图 12中的轮廓图像进行$D_{b}$值计算，由于图像分辨率不同，因此需要选择合适的盒子尺度来使不同图像所使用的覆盖盒子的感受范围相同(即单个盒子所能代表的真实范围相同)。由于拍摄的是同一片叶子，所以各图像中的物体存在成比例的放大或缩小关系，它们的最小包容矩形也存在这种比例关系，通过计算各图像中物体的最小包围矩形，根据其长边(或距边)比例设置对应比例大小的盒子，例如有两幅图像中物体的最小包容矩形的长边分别为$l_{1}=512、l_{2}=1 024$，则所采用的盒子的大小$r_{1}/r_{2}=l_{1}/l_{2}=0.5$。最后的计算结果见表 3，BCM方法的$D_{b}$值变化幅度为0.044 52，平均拟合误差为0.008 01，SRM方法为0.027 48，平均拟合误差为0.007 51。说明旋转骨架法对不同分辨率及旋转角度的图像具有较高的稳定性，拟合误差也相对较低。 图11 不同分辨率与旋转角度下的植物叶子图像 Plant leaf images with different resolution and rotation anglesFig 11 图12 不同分辨率与旋转角度下的植物叶子轮廓图像 Outline drawing of plant leaf at different resolutions and rotation anglesFig 12 表3 不同分辨率与旋转角度的叶子轮廓图像的计算结果 图像名称 BCM 本文(SRM) 分形维数 拟合误差 分形维数 拟合误差 F01* 1.023 24 0.008 02 1.030 59 0.007 50 F02* 0.998 77 0.007 12 1.040 90 0.006 31 F03* 1.022 97 0.007 51 1.037 78 0.007 39 F04* 1.043 29 0.014 95 1.055 16 0.013 44 F05* 1.034 11 0.010 15 1.056 96 0.007 20 F06* 1.019 38 0.005 89 1.029 48 0.005 42 F07* 1.026 19 0.002 42 1.034 24 0.005 36 Calculation results of leaf countour images with different resolution and rotation anglesTable 3 加粗字体表示最优结果。4结论本文提出了一种新的二值图像的分形维数计算方法——旋转骨架法(SRM)，将二值图像转化为矢量图骨架，并将骨架旋转$θ$角度使其旋转角度变为0，最后对旋转后的骨架进行分形维数求取。通过对已知自相似分形图像、文字扫描图像和植物图像的分析，对比实验证明了旋转骨架法的高稳定性与准确度，相比于BCM方法，自相似分形图像的平均误差分别降低了0.725 2%、3.060 5%、2.298 5%，变化幅度分别降低了0.057 3、0.088 3、0.085 9；文字扫描图像的变化幅度降低了0.012 75，平均拟合误差降低了0.001 28；植物图像的变化幅度降低了0.017 04，平均拟合误差降低了0.000 5。同时，该方法可以解决尺度$r$非整数的问题，尤其适用于$θ≥2°$的二值图像，为图像分形维数的研究提供了新的思路。本文算法在计算二值图像分形维数时有很好的稳定性与准确度，但是当旋转角度$θ 2°$时，SRM方法的精度有所下降，该方法适用于旋转角度$θ≥2°$的二值图像。同时，使用矢量图会增加计算时长，特别当图像结构的像素点数量较多时，计算时长会大幅增加。将来的研究方向可以围绕矢量图的计算参数进行优化，比如优化盒子尺度范围、盒子覆盖方式等参数，以提高计算速度。