一种三维图像欧拉数快速计算方法

姚 斌， 何立风，2， 康世英， 赵 晓

(1.陕西科技大学 电子信息与人工智能学院， 陕西 西安 710021； 2.日本爱知县立大学 信息科学学院， 日本爱知县长久手市 4801198； 3.咸阳师范学院 计算机学院， 陕西 咸阳 712000)

0 引言

图像的欧拉特征是图像重要的拓扑性质之一，用来描述图像的结构，它与图像的几何形状无关，并且不会因为图像的扩大、缩小、旋转及变形而改变，有很强的鲁棒性.二维图像的欧拉特征称为欧拉数(Euler Number)，它被定义为该图像中连通域数量与孔洞数量之差[1]，被广泛应用于医学诊断、字符识别以及地质勘探中土壤内部孔洞连通性分析等方面[2-4].三维图像的欧拉特征也被称为三维图像的欧拉-庞加莱特征(Euler-Poincare Characteristic)或简称三维欧拉数.三维欧拉数分为三维表面欧拉数和三维体欧拉数，对于确定三维图像的内部结构，三维体欧拉数会更有效(在后续篇幅中，三维图像欧拉数特指体欧拉数).

三维图像欧拉数经常被用来描述物体的基本特征，例如用来描述结构的连接性等[5]，被广泛用于解决实际问题.比如，H.J.Vogel等[6]利用三维图像欧拉数确定土壤孔隙结构；A.Pierret等[7]利用三维图像欧拉数来重建和量化macropores.P.Lehmann等[8]分析了几何形态对水流和流体分布的影响，他们利用三维图像欧拉数来预测多孔介质中的流体相分布.在工业生产方面，A.Velichko等[9]利用三维图像欧拉数区分铸铁中的石墨形态，使人们更容易理解铸铁增长机制和铸铁性能.在医学方面，L.Kubínová等[10]用三维图像欧拉数来测量毛细血管的长度.由于欧拉数多用于物体的快速识别与分类，欧拉数的计算速度越快，就越有利于提高相关实时图像处理系统的整体性能.

1 三维图像像素的邻接性及三维图像欧拉数

1.1 三维图像像素的邻接性

在本文中，图像特指二值图像，按照惯例，用“1”表示二值图像中的前景像素(目标像素)，而用“0”表示背景像素(非目标像素).

在对三维图像的研究中，像素之间的邻接性可以分为6-邻域连接和26-邻域连接[11].对于两个不同的像素X(x1,x2,x3)和Y(y1,y2,y3)，当|x1-y1|+|x2-y2|+|x3-y3|=1时，称像素X与像素Y为6-邻域连接.如图1所示，像素p的6-邻域连接像素有p1、p3、p5、p7、p17和p26.而当max(|x1-y1|,|x2-y2|,|x3-y3|)=1时，称像素X与像素Y为26-邻域连接.在图1中，像素p的26-邻域连接像素为p1，p2，…，p26.

1.2 三维图像欧拉数

三维图像的欧拉数被定义为：E=O-H+C.其中,O、H和C分别为图像的连通域、孔洞和空腔的数量.三维图像中“孔洞”定义具有多义性与不确定性，在文献[12]中，C.N.Lee将“孔洞”解释为不产生分离的最大切割数，也被称为“tunnel”，而“空腔”则被称为“cavity”.由于直接利用定义来计算三维图像的欧拉数需要获得图像中连通域、孔洞和空腔的数量，也就需要对图像有全面了解，无法由图像的局部性质计算得到，因此不能将欧拉数的计算分解成每个网格片的子任务，无法进行局部计算.

图1 三维图像中像素的邻接性示意图

除了利用三维图像欧拉数定义中的公式，三维图像欧拉数也可以通过局部计算的方法得到.由局部计算三维图像欧拉数的算法最早是在1971年由C.M.Park等[11]的论文提出的，但其算法只对三维图像像素间的6-邻域连接关系有效，无法计算图像像素间是26-邻域连接关系时的图像欧拉数.1980年，D.G.Morgenthaler[13]提出一种对三维图像像素间26-邻域连接关系有效的欧拉数局部算法，该算法需要检查图像中每个2×2×2立方体(称为“三维立方体模式”，有256种可能的情况)的像素组成情况，而且被后续的研究所沿用[14].这种方法计算欧拉数需要的像素读取次数为图像像素总数的8倍.在后续篇幅中为了表示方便，把该算法称为DGM算法.

1992年，国内学者杨敬安等[15]提出了利用微分几何及代数拓扑原理计算三维图像欧拉数的方法并给出了相应算法.1996年P.K.Saha等[16]提出了计算三维图像欧拉数的新方法，该方法是通过计算图像中3×3×3邻域内的连通域、孔洞和空腔在数量上的变化来实现的.由于需要同时考虑到图像中连通域、孔洞和空腔的数量，效率比单纯计算欧拉数低一些.2010年，国内学者林小竹等[17]提出了计算三维图像欧拉数的新公式，利用图像中图段数及相邻图段数实现了三维图像欧拉数的计算.这种方法将二维图像中的方法推广应用到三维图像中，对于密度较低的图像效率较高.2013年，H.Snchez Cruz等[18]提出了一种新的计算三维图像欧拉数算法，利用图像几何方面的表面连接性，通过计算2×2×2三维基本立方体和2×2×1、2×1×2、1×2×2立方体的数量来实现计算.与文献[13]中提出的算法类似，采用这种方法计算给定三维图像欧拉数需要的像素存取次数也是该图像中像素总数的8倍.

2 DGM算法

在DGM算法中，通过统计给定三维图像中包含的特定2×2×2的三维立方体模式的数量来计算图像欧拉数，如图2所示.

图2 DGM算法用到的立方体模式

如果用qn表示给定三维图像中含有n个前景像素的2×2×2三维立方体的个数，例如，q3是含有3个前景像素的三维基本立方体的个数，即图2中(6)、(7)和(8)所示三种模式数量总和，则给定三维图像在26-邻域连接情况的下的欧拉数E26的计算公式为：E26=q1-q2+q3-q4+q5-q6+q7-q8.由于图2中列出的三维立方体模式是根据其中前景像素的个数及对称性分类的，在处理图像时需要从各个方向考虑，实际计算中使用起来并不方便.

为了便于计算，D.G.Morgenthaler通过对256种2×2×2的三维基本立方体模式的立体特征进行研究、分析，对其中八个像素a、b、c、d、e、f、g和h(如图3所示)的取值情况从00000000到11111111进行了分析，计算出256种模式中每种模式对于图像欧拉数的增量，并给出了相应的查找表.对照查找表进行分析，可以看出，在256种模式中，有207种模式的欧拉数增量为0，30种为-1，17种为1，其余两种为-2[8].对于图像欧拉数的增量不为0的49种三维立方体模式的编号及欧拉数增量△E如表1所示.

表1 DGM算法中对图像欧拉数有贡献的模式

在实际进行三维图像欧拉数的计算时，只需要对图像进行光栅扫描，对于扫描到的像素，只需要检查、判断该像素所在的三维立方体模式属于256种模式中的哪一种，然后查表即可得到该立方体模式对图像欧拉数的增量值.所有像素扫描完成后，三维图像中所包含的所有三维立方体模式都被检查过，所有的三维立方体模式的欧拉数增量的总和即是给定三维图像的欧拉数.

如图3所示，用DGM算法计算三维图像欧拉数时，需要对整个图像的像素进行扫描.对扫描到的每个像素，需要检查其所对应的三维立方体模式中的其他七个像素才能得到当前三维基本立方体对应的欧拉数增量.例如，对于当前扫描到的像素a，需要检查像素a所在的三维立方体模式{a,b,c,d,e,f,g,h}，这样，处理每个像素一共需要检查八个像素.

图3 DGM算法中检查三维立方体的像素

3 算法的改进

在DGM算法中，为了计算给定三维图像的欧拉数，需要检查图像中所包含的所有立方体模式.尽管需要统计的立方体模式只有49种，但对当前检查的立方体模式，还需要检查8个像素才能确定该立方体模式是不是需要统计.对于一个大小为X×Y×Z像素的图像，为了计算图像欧拉数，像素存取次数为8×X×Y×Z.下面，本文从两个方面对DGM算法进行改进.

3.1 改变像素检查顺序

DGM算法中，尽管只有一部分立方体模式对图像欧拉数有贡献，但是为了找出对图像欧拉数有贡献的立方体模式，还是需要对图像中包含的所有立方体模式进行检查.如图3所示，检查每一个立方体模式所包含的a、b、c、d、e、f、g和h八个像素的值.另一方面，既然只有49种立方体模式对图像欧拉数有贡献，可以把这些立方体模式的每一个像素值列出来，进行分析，寻找其中规律，尽可能减少不必要的像素存取以提高算法效率.

将表1中所有三维立方体模式编号用二进制数表示后(从高位到低位依次为a、b、c、d、e、f、g和h)，可以得到表2.

表2 DGM算法中对图像欧拉数有贡献的模式

在49种对图像欧拉数有影响的立方体模式中，本文对所有立方体模式编号的二进制编码进行分析，发现编码中数字“1”出现的位置是有差异的.

经过仔细分析后，发现数字“1”在所有立方体模式编号的二进制编码中“b”的位置都没有出现.根据立方体模式二进制编码的这个特点，在检查当前立方体模式的八个像素值时，可以先检查像素“b”的值，如果“b”的值为1，由表2得知，当前立方体模式对图像欧拉数没有影响，因此可以跳过该立方体模式的检查，直接检查下一个立方体模式.这样，对于前景像素较多的图像，在特定情况下(像素b为前景像素)只需要检查一个像素就可以确定该立方体模式是不是对图像欧拉数有影响，这样可以节省很大一部分时间，算法效率也随之提高.

因此，在算法实现时，本文将检查像素的顺序由原来的“abcdefgh”改变为“bacdefgh”.尽管只改变了两个像素的检查顺序，在检查图像中包含的立方体模式时，部分立方体模式只需要检查一个像素就可以确定该模式是否对图像欧拉数有影响.

3.2 合并部分立方体模式

对表1进行进一步分析，发现在对图像欧拉数有影响的49种立方体模式中，部分立方体模式的编号是连续的，如编号24-27，36-39，44-47，56-59，148-151，156-159，180-183和188-191的立方体模式，这些编号连续的立方体模式还有个共同点，就是对图像的欧拉数增量相同，如表3所示.

表3 可以合并的立方体模式

由表3可以得到，对于编号为24-27的四种立方体模式，这些立方体模式中a、b、c、d、e、f、g和h八个像素形成的二进制编码的值为00011000-00011011.不难发现，这四个编码的八位数值中只有最后两位不同(00-11)，前六位数值是相同的(000110).与此同时，可以发现这些立方体模式对图像欧拉数的增量相同，都是“-1”.因此，在算法实现时，可以把这四种立方体模式合并，检查立方体模式的像素时，只需要检查前六个像素的值就可以确定该模式对图像欧拉数的影响.

对于当前正在处理的立方体模式，当检查其像素串“abcdef”的值为“000110”时，就可以确定当前立方体模式是编号为24-27的模式中的某一个，因此可以直接将图像欧拉数增加“-1”，不需要再检查最后两个像素的值.同理，当检查立方体模式时，如果二进制编码“abcdef”的值为“001001”、“001011”或“001110”时，图像欧拉数增加“-1”.

而对于正在处理的立方体模式，当检查像素串“abcdef”的值为“100101”、“100111”、“101101”或“101111”时，同样可以不用检查最后两个像素的值，直接可以将图像欧拉数增加1.这样，同样可以检查六个像素就可以判断一个立方体模式的类型，与DGM算法相比，在处理部分立方体模式时可以减少像素存取次数，同样可以提高算法效率.

算法在实际实现时，将改变像素检查顺序和合并立方体模式两种改进方法相结合，对于当前正在处理的立方体模式，首先将像素检查顺序更改为“bacdefgh”，如果像素“b”的值为“1”，直接跳过该立方体模式，检查下一个；如果正在处理的立方体模式中像素“b”的值为“0”；再检查其他像素的值.如果前六个像素的值是表3所列的八组之一，则由前六个像素的值就能确定正在处理的立方体模式对图像欧拉数的影响，尽可能减少像素存取次数，提高算法效率.

4 实验结果及分析

为了测试提出的算法的正确性和有效性，本文采用结构较为复杂的噪声类型的图像对算法进行测试，并把本文提出算法(Ours)的运行结果与DGM算法进行对比.

本文采用41幅噪声图像作两种算法的图像密度与程序执行时间的关系对比实验.噪声图像的密度指图像中前景像素的比例，密度越大，前景像素所占比例越大.当图像密度为1时，所有像素均为前景像素；反之，当图像密度为0时，所有像素均为背景像素.不同密度的三维图像截面如图4所示.

41幅噪声图像采用门限法随机生成，阈值从0到1 000，步长值为25，大小为128×128×128像素.测试平台为台式机(Intel Core i7-6700 CPU@3.40 GHz，8 GB内存，操作系统Windows 7)，编译器为GNU C compiler(Version 4.6.1).为保证实验结果的准确性，所有计算图像的欧拉数的程序均被运行1 000次，取平均值作为实验结果.

图4 不同密度的三维图像截面

为了验证本文改进算法的正确性，分别采用两种算法计算给定图像的欧拉数.实验结果表明用两种算法计算给定图像的欧拉数结果是相同的.表4给出了用两种算法计算时部分图像的欧拉数值.

表4 两种算法计算图像欧拉数的结果

算法效率实验结果如图5所示.从图5可以看出，几乎对于所有图像，本文提出的算法的效率都高于DGM算法，特别是在图像密度值在0.4～0.8之间的图像，本文提出算法的优势更为明显.在本文提出的算法中，对于当前正在检查的立方体模式，首先检查像素“b”的值，如果该像素的值为“0”，然后再检查其他像素的值；如果像素“b”的值为“1”，则可以跳过当前立方体模式，直接检查下一个.在密度比较小的图像中，图像中包含的“1”的个数较少，因此，本文提出算法的效率和DGM算法相比优势不是很明显.从图5可以看出，对于密度接近0的图像，两种算法计算欧拉数的时间几乎相等.随着图像中前景像素的增加，本文提出算法的优势逐渐增大.另一方面，计算图像欧拉数所需要的时间除了对像素进行读取、判断之外，还需要进行一些计算，前景像素越多，需要的计算量也越多，因此当图像中前景像素的数量增加到一定程度后两种算法效率之间的差异反而逐渐缩小了.

图5 两种算法的运行时间对比图

5 结论

为了解决三维二值图像欧拉数的快速计算问题，本文在现有算法的基础上，通过改变立方体模式的像素检查顺序、合并部分立方体模式两种方法，尽可能减少处理一个立方体模式所需要检查的像素个数，提高算法效率.实验结果表明，本文提出的算法的计算结果与现有算法得到的结果相同，对于绝大多数被测图像来说，本文提出算法的性能优于现有算法.

在算法实现过程中，发现在对图像进行逐行扫描、检查每个三维立方体模式来计算图像欧拉数时，仍然存在像素重复检查的现象.所以在接下来的工作中，笔者将对算法进行优化和改进，以期得到性能更好的二值图像欧拉数算法.