基于改进图割的超声瞳孔图像分割与直径测量

周亚丽 高玉芝 周 浩 周 密 郑音飞*

1(浙江大学生物医学工程与仪器科学学院，杭州 310027)2(浙江大学医学院附属第二医院急诊科，杭州 310027)

基于改进图割的超声瞳孔图像分割与直径测量

周亚丽1高玉芝2周 浩1周 密1郑音飞1*

1(浙江大学生物医学工程与仪器科学学院，杭州 310027)2(浙江大学医学院附属第二医院急诊科，杭州 310027)

传统的瞳孔直径测量是通过医生手工标定，对于眼外伤和丧失意识的患者测量不方便。针对瞳孔直径测量的人工交互量大且测量鲁棒性不强的问题，采用图割算法分割瞳孔超声图像并测量瞳孔直径。对传统图割算法进行两个方面的改进，采用自适应阈值的区域生长代替人为种子点选取，在保证分割效果的基础上减少了图割的交互量；在能量函数的数据项部分增加图像的梯度信息，减少了原始算法分割结果中出现的小区域，增强了对弱边缘的分割。最后，对采集到的超声瞳孔图像进行自动分割、自动测量瞳孔直径，可以得到患者瞳孔的直径动态变化，给临床诊断提供依据。为了验证算法的有效性，对10位患者的动态瞳孔超声图像进行基于改进图割的瞳孔直径测量，并与医生的手动测量结果对比。结果表明，本方法的结果与医生手动测量结果的绝对误差小于0.2 mm，相关系数不小于0.83。通过改进图割算法，改善了分割效果，实现了超声瞳孔动态图像的自动直径测量，并可有效代替瞳孔直径的人工测量，减少人工交互量。

图割；超声瞳孔图像；图像分割；瞳孔直径

引言

眼睛是人体中重要且调节机制复杂的器官，其中瞳孔是眼睛内虹膜中心的小圆孔，是光线进入眼睛的通道。瞳孔的缩小与扩大受交感神经、副交感神经和脑干等多个神经中枢支配，它的大小变化有许多非常重要的生理意义。例如，瞳孔直径的波动可以反映生理觉醒[1]；瞳孔直径的变化可以间接反映心率变化，实现无接触无损伤测量心率变异性[2]，反映和评价自主神经活动[3-4]。定量分析和评价瞳孔直径具有重要的临床意义。

瞳孔直径的诊查方式主要有传统的直尺、Rosebaum瞳孔尺，现代的红外瞳孔测量技术和临床上常用的波前像差仪、角膜地形图[5]等。在医院急诊科瞳孔的诊查非常重要，很多发生眼外伤、昏迷及术后的重症患者采用以上方式进行瞳孔大小诊查[6-11]。这些测量方式由于需要拨开眼睑或特定方式测量，所以对于患者可能造成二次损伤而致失明[12-13]。然而，瞳孔的超声探查只需要将超声探头放在眼睑之外，就可以清晰地观察到瞳孔及其他眼部结构[13]。同时，超声检测具有安全、实时、成本低廉、操作简单等特点[14]，可以提供瞳孔的实时动态图像，检测瞳孔的动态变化。

现今用超声定量评价瞳孔的直径都需要医生手动标定瞳孔直径的位置从而得到直径大小[13]，这样做在定量瞳孔变化时需要太多的人工交互并且耗时。图割作为交互式图像分割技术[15]，已经广泛用于各种二维图像、三维图像以及彩色图像的分割[16]。传统的图割以及其他的图割改进基本上都是将加强能量函数的约束，约束的方法有增加形状先验、像素相似性以及SIFT特征等[17-21]。这些方法仍然需要人工交互选取种子点，本研究在原有图割的基础上改进种子点的选取，对图割种子点选取结合阈值自适应及像素数目约束区域生长进行改进，在保证图像分割效果的基础上减少图割的人工交互，基本实现自动分割和测量瞳孔直径；并且改进图割能量函数的数据项，在数据项中增加图像的梯度信息，以适应超声瞳孔图像边缘较弱的特点，从而实现更好的分割。

1 方法

瞳孔以及超声瞳孔图像如图1所示，超声探头放置在眼睑上方，所成的超声瞳孔图像是一块暗区域，通过改进的图割算法将这一暗区分割出来，在得到的分割图像自动测量瞳孔的左右径，最终得到瞳孔的直径。

图1 瞳孔示意图和超声瞳孔图像。(a)瞳孔示意图；(b)超声瞳孔图像Fig.1 Pupil schematic and ultrasound pupil image.(a)Pupil schematic;(b)Ultrasound pupil image

1.1 图割改进

图割是基于图论的图像分割方法[22-24]，对图像建立图模型，在此基础上实现图像分割。对于图像I，建立相应的图G=(如图2所示)，其中V是图像中的所有结点，E是连接结点的每一条边。

图2 图像映射为图Fig.2 Image mapped into graph

在图2中，结点p,q∈V对应图像中的像素，其中结点S与T是额外加入的终端结点；图中边n-links∈E对应像素与像素之间的邻接关系，其中像素与终端结点之间的边称为t-links，像素之间的连接边称为n-links。

1.1.1 能量函数及改进

图割能量函数是一个关于标号的能量函数[23]，其中标号代表图像中的不同区域。对与二元图像的分割，标号值只有两个，一个代表目标，另外一个代表背景，具体形式为

(1)

式中，l代表将图中某一点像素q标为标号l，Es(l)是能量函数中的平滑项，Ed(l)是数据项，λ是平衡两项的常数系数。

能量函数的具体求解如下：

(2)

(3)

(4)

对于所有的像素p∈P(P是图像中所有像素的集合)，对它进行标定，把它标定为标号lp所需要消耗的代价值为Dp(lp).这个代价值表征了标号lp与像素p的相似程度，两者越相似，所需消耗的代价就越小。也可以说，将像素标为目标或背景时的概率大小，如果像素属于目标的概率大，那么它被标为目标所需要消耗的代价就小；反之，将它标为背景所需要消耗的代价大。所以，可以用图像中像素属于目标和背景的概率来得到能量函数的数据项[16]。记像素p属于目标的概率为RO(p)，属于背景的概率为RB(p)，那么数据项可以具体表示为

(5)

(6)

具体的概率可以通过选取种子点，通过判断像素与已知种子点的相似程度来估计RO(p)、RB(p)。一般地，该相似程度通过对选取的种子点进行灰度直方图统计来得到。

本研究除应用了灰度直方图统计信息外，能量函数的数据项增加了梯度的信息。因为瞳孔的超声图像动态范围小，尤其是在要分割的瞳孔附近的动态范围非常小，单用灰度信息无法实现较好的分割，所以将梯度信息引入图割的数据项中。通过梯度直方图，可以估计整幅图像每个像素点的梯度概率；通过像素-梯度-梯度直方图的相对应关系，最后得到梯度概率估计与灰度概率估计，将两者结合起来估计总的像素概率，有

RO(p)=gray_histO(p)+grad_histO(gradp)

(7)

RB(p)=gray_histB(p)+grad_histB(gradp)

(8)

式中，gray_histO、grad_histO、gray_histB和grad_histB分别是目标种子灰度直方图、梯度直方图、背景种子灰度直方图和梯度直方图；gradp是指像素p所对应的梯度值。

平滑项表示对于邻接域N，所有邻接的像素p,q∈N两者分别标为标号lp、lq时所需要消耗的代价为Vp,q(lp,lq)，这表明当两个标号值不相同时，所需要消耗的代价可以用相邻像素之间的差异来表征。在边界处，修正标号对图像分割起到平滑的作用，一般形式可以表示为

(9)

式中，Ip、Iq分别是像素p、q的灰度值，σ是邻域内的方差，dist(p,q)是相邻像素p、q的欧式距离。本算法的邻域选择为4邻域，因此相邻像素之间的欧氏距离均为1。方差σ取整幅图像的像素灰度方差，相邻像素的灰度差异用高斯滤波来得到，分别进行横向和纵向的高斯滤波，得到横向和纵向的像素灰度差异代入公式中。

综上所述，可以说明图像分割就是通过最小化能量函数来实现对图像前景与背景的分割。如何实现该能量函数的最小化，就是通过用能量函数构建图的边权值，再通过的图论的方法实现能量函数的最小化。各个边的权值分配如表1所示。

表1 图模型边权值(其中lT、lS分别代表背景和目标标号)

Tab.1 The edge weight of graph,lTandlSrepresent background and object labels)

边类型权值(代价)条件p,q{}Vp,q(lp,lq)p,q{}∈Np,T{}λ·Dp(lT)p∈Pp,S{}λ·Dp(lS)p∈P

将能量函数以权值的形式分配给图的边，运用图论中的原理最大流-最小割算法[23-24]实现该能量函数的最小化，图割的流程如图3所示。其中，B表示背景种子点，O表示目标种子点。图割的流程主要为：以一幅3×3的图像为例，分别选取种子点B和O构建能量函数，然后建图来进行能量函数最小化，最后对标号标定实现分割。

图3 图割流程Fig.3 The process of graph-cut

1.1.2 种子点选取改进

能量函数中数据项构建的关键是要得到每个像素属于图像前景和图像背景时的概率，这个概率可以通过一定的先验知识来得到。一般通过人为地标定目标和背景种子点，用目标和背景种子点来估计概率。对于灰度图像来说，统计出的结果是灰度直方图，这样每个灰度级G就对应一个频率PG，用频率PG来估计每个灰度级(像素)的概率RO(p)、RB(p)。

由以上的描述和图3所示可以说明，种子点的选取至关重要，在一定范围内选取的种子点越多、越准确，得到的概率分布也越接近真实值，分割效果也越好；种子点选取的差异，会导致分割结果的不同、算法鲁棒性的好坏与否。再者，图割算法是需要先验信息支持来得到能量函数，属于交互式图像分割。但交互过程中描点形式的种子点选取需要大量的人工干预，尤其是在超声序列图像分割中，所以本研究采用自适应区域生长来选取种子点。

对于图4(a)所示的超声瞳孔图像来说，瞳孔的处回波信号弱、灰度小且均匀，用区域生长算法选取目标种子点的人工交互少而且选取的区域准确，所以可以用区域生长算法来选取瞳孔目标种子点。

目标种子点选取的区域生长算法见表2。由于目标灰度均匀且都偏暗，所以阈值选取的范围不大，本研究选定的阈值选取范围为1～30(灰度值单位)，并且可以控制区域生长结果的像素数以得到最优的生长结果。

表2 区域生长选取种子点算法

图4 不同种子点选取得到的分割结果对比。(a)超声瞳孔原图(中间黑色近似椭圆的区域为瞳孔)；(b)～(d)种子点选取；(e)～(g)相应种子点选取后的分割结果；其中红色方框内的像素点为目标种子点，绿色方框内的像素点为背景种子点Fig.4 Comparison of segmentation results with different seeds selection. (a)The original ultrasound pupil image,(b)～(d)Seed selection;(e)～(g)Comparison about segmentation results with different seeds selecting; The pixels in red rectangles are object seeds and green are background seeds

图5 图割的目标和背景种子点选取，其中白色区域就是目标和背景种子点。(a)原图；(b)为目标种子点；(c)为背景种子点Fig.5 The white region is object and background seeds. (a)Original image; (b)Object seeds; (c)Background seeds

背景种子点选取是在目标种子点选取完成之后，选取目标种子点的上部区域为背景种子点，由于超声瞳孔成像的方式和所成图像的特点，瞳孔上部的像素灰度较大，与目标有较好的区分，所以选取此处为背景种子点，能得到较好的分割效果。目标和背景种子点的选取如图5所示。

2 结果

2.1 能量函数改进的分割结果

通过本算法中能量函数数据项的改进，对3幅超声瞳孔图像进行分割，得到原始算法、仅使用梯度信息以及灰度和梯度信息叠加使用时的分割结果对比，如图6所示。

图6 本算法改进与原始算法对比。(a)对不同瞳孔图像的原始算法分割结果；(b)采用图像梯度信息作为图割数据项的分割结果；(c)本方法的分割结果Fig.6 The results of proposed and original algorithm, (a)Results of different pupil images with original graph-cut;(b)Results of different pupil images with gradient data portion of energy function; (c)Results of different pupil images with gradient plus gray data portion of energy function;

由以上的结果对比可以看到，本研究得到的分割结果更好，没有出现原始算法中的分割小区域，同时边界清晰。

2.2 临床瞳孔直径测量结果

标准采集瞳孔图像应该是瞳孔居图像正中，可以看见晶状体等结构，但人的眼球是会转动的，眼球上翻或侧斜都会影响整个瞳孔轮廓的探查，有时获取的瞳孔图像轮廓不清或呈椭圆，所以测量时要测瞳孔的左右径。采集超声瞳孔视频的操作如图7所示，要对视频中的逐帧图像手动测量瞳孔直径。

笔者测量瞳孔直径的方法是对瞳孔avi视频的第一帧标定一个最初种子点，然后通过本图割算法自动对avi视频逐帧处理分割瞳孔图像，并且对每帧图像自动测量瞳孔直径。自动测量瞳孔直径的原理与手动测量的一致，也是测量瞳孔的左右径，测量时选取整个瞳孔垂直方向正中的位置测量瞳孔直径。

按照上述方法，对瞳孔avi分别手动和自动测量瞳孔直径，手动测量重复次数为5次，取5次的平均值作为手动测量的结果，如图7所示。

图7 瞳孔超声图像获取。(a)超声瞳孔图像获取操作；(b)瞳孔直径手动标定Fig.7 The acquirement of ultrasound pupil images. (a)Operation of acquiring ultrasound pupil images; (b)Manual calibration of pupil diameter

图8 瞳孔直径临床测量自动与人工方法的对比Fig.8 The result of clinical pupil diameter measurement: comparison of the proposed automatic and manual method

从图8中的人工与自动测量的结果对比中可以看出，自动测量的结果与人工结果较接近，误差在0.1 mm内。在瞳孔对光反射测试中，移走光源后瞳孔放大，自动测量准确地检测到了瞳孔的放大，但因为运动伪迹的原因与自动测量还是有一定的误差。总的来说，本算法基本实现了瞳孔直径的准确测量。

表3中列出了5个临床瞳孔自动测量与人工测量的相关系数与绝对误差分析。从整体上看，相关系数为0.83～0.96，平均值为0.93，具有统计显著性。本方法与参考基线的绝对误差是(0.11±0.09)mm,均值偏差较小。另外，在进行人工测量瞳孔时，5次重复测量的标准差在0.02～0.1 mm之间，而自动测量的标准差接近于0，可见本研究自动测量的重复性更高。

表3 人工与自动测量方法的对比

Tab.3 Comparison of manual and automatic measurement

序号图像帧数/帧相关系数绝对误差/mm(平均值±标准差)1600.830.20±0.152600.990.04±0.023600.950.08±0.104600.930.10±0.055600.960.06±0.076600.980.04±0.027600.990.03±0.028600.830.08±0.079600.990.06±0.0410600.970.05±0.02

3 讨论和结论

本研究实现了瞳孔直径的自动测量，解决了医生手动测量瞳孔直径交互量大、测量重复性差且不方便的问题。主要运用图割算法对动态超声瞳孔图像进行分割，然后求取瞳孔直径变化曲线。在图割算法上，用自适应阈值和像素数目约束代替人工种子点的选取，仅需要在第一帧确定好一个瞳孔的中心位置，人工交互少，并且运用了超声序列图像每帧之间的连续性；改进了能量函数的数据项，增加了梯度信息，使得图割算法在超声瞳孔图像的弱边缘处减少溢出，并且减少了分割结果中的小区域。另外，自动直径测量时选取直径的方法是对分割结果选取上下正中部位的左右径作为瞳孔的直径，这与医生选取直径的方法基本一致。最后，本方法应用于临床的10位ICU患者的瞳孔动态图像，从得到的瞳孔直径变化曲线与手动测量结果的对比可知，本方法的鲁棒性更强，并且能够代替手动测量，有效减少了人工交互量。

瞳孔大小的监测是临床监护的重要组成部分，医院监护多采用定时对病人瞳孔直径进行人工测量记录的方法，需要进行人工测量的次数多，并且人工测量受人为因素影响大；另外，医生通常通过肉眼观察瞳孔直径的动态变化，仅能观察到瞳孔大致的变化，定性地评判瞳孔变化无法定量。本研究实现的瞳孔自动测量能有效解决上述问题。运用本方法可自动得到瞳孔直径动态变化曲线，在减少人工交互量的同时增加了诊查信息。医生或是看护人员对于患者的瞳孔监测不再局限于定时的瞳孔直径记录，可以有瞳孔的动态变化信息，简化瞳孔测量并拓展了瞳孔观测，为病情变化提供了信息。

与红外拍摄瞳孔来获取瞳孔直径的动态变化曲线相比，利用超声图像的分割来测量瞳孔直径的方法具有安全、实时、成本低廉、操作简单的特点。但是，超声图像的分辨率低，要想准确地测量瞳孔直径就需要能对超声图像实现较精确的分割，所以本研究对原有的图割算法进行改进，以实现对超声图像的准确分割。

本方法主要处理瞳孔的动态超声图片，即序列超声图片。在对序列超声图片的处理过程中，要充分利用序列图像的连续性，避免人工交互简化后不能保证对后部序列图像的分割准确性；对于传统图割算法，其能量函数的设计与具体图像特征有关，并且关系到分割结果的好坏，所以从能量函数着手改进图割算法是比较常见的切入点。最后，本方法可用于医生根据瞳孔大小的变化来反映生理状态的研究。

[1] Wingate M, Hood B, Shi H,etal. Automated pupil size determination for evaluating fluctuations in physiological arousal [C] //Yue Minzhu, Feng Yang, eds. Signal Processing Proceedings. Beijing: IEEE, 1998:1666-1669.

[2] Parnandi A, Gutierrez-Osuna R. Contactless measurement of heart rate variability from pupillary fluctuations [C] //Suter G, eds. Affective Computing and Intelligent Interaction. Geneva: IEEE, 2013:191-196.

[3] Yoshida H, Mizuta H, Gouhara T,etal. Statistical properties of simultaneously recorded fluctuations in pupil diameter and heart rate [C] //Hood L, eds. Engineering in Medicine and Biology Society. Montreal : IEEE, 1995: 165-166.

[4] Yamaji K, Hiratai Y, Usui S. The pupil as a possible monitor as the autonomic nervous system [C] //Hood L, eds. Engineering in Medicine and Biology Society. Chicago: IEEE, 1997: 2777-2781.

[5] 汪育文，保金华，吕帆. 三种瞳孔直径测量方法在不同测量状态下测量结果的比较和分析 [J]. 眼视光学杂志, 2009, 11(3): 218-220.

[6] Starck T, Liu Yun, Prewett AL,etal. Comparison of scotopic pupil measurement with slitlamp-based cobalt blue light and infrared video-based system [J]. Journal of Cataract and Refractive Surgery, 2002, 28(11): 1952-1956.

[7] 杨亚波, 俞一波, 付清. Ⅱ型糖尿病患者瞳孔相关性自主神经病变的研究 [J]. 中华眼科杂志, 2006, 42(7): 616-619.

[8] Yang H, Lee M, Kim JB,etal. Burst-shot infrared digital photography to determine scotopic pupil diameter [J]. Journal of Cataract and Refractive Surgery, 2006, 32(12): 2113-2117.

[9] Schnitzler EM, Baumeister M, Kohnen T. Scotopic measurement of normal pupils: colvard versus video vision analyzer infrared pupillometer [J]. Journal of Cataract and Refractive Surgery, 2000, 26(6): 859-866.

[10] Arthur C, Dennis L. Comparison of the colvard pupillometer and the zywave for measuring scotopic pupil diameter [J]. Journal of Refractive Surgery, 2004, 20(3): 248-252.

[11] Fotiou F, Fountoulakis KN, Goulas A,etal. Automated standardized pupillometry with optical method for purposes of clinical practice and research [J]. Clinical Physiology, 2000, 20(5): 336-347.

[12] Pascal SCJ, Rosen P. Ocular examination technique for the emergency department [J]. The Journal of Emergency Medicine, 1997, 15(6): 793-810.

[13] Harries A, Shah S, Teismann N,etal. Ultrasound assessment of extraocular movements and pupillary light reflex in ocular trauma [J]. American Journal of Emergency Medicine, 2010, 28(8): 956-959.

[14] Terrence D. Ultrasound physics and instrumentation [J]. Surgical Clinics of North America, 1998, 78(2): 197-217.

[15] 汪一休. 关于交互式图像分割算法的研究 [D]. 合肥: 中国科学技术大学, 2009.

[16] Chang Huali, Chen Zhenping, Huang Qinghua, et al. Graph-based learning for segmentation of 3D ultrasound images [J]. Neurocomputing, 2015, 151(2): 632-644.

[17] Malcolm J, Rathi Y, Tannenbaum A. Graph cut segmentation with nonlinear shape priors [C] //Mahesh P, eds. Image Processing. San Antonio: IEEE, 2007: 365-368.

[18] Slabaugh G, Unal G. Graph cuts segmentation using an elliptical shape prior [C] //Mahesh P, eds. Image Processing. San Antonio: IEEE, 2005: 365-368.

[19] Veksler O. Star shape prior for graph-cut image segmentation [C] //David F, Philip T, Andrew Z, eds. Computer Vision. Marseille: Springer Berlin Heidelberg, 2008:454-467.

[20] 胡焦, 陈辉. 基于像素相似性与Graph-Cut 的图像自动分割 [J]. 现代电子技术, 2004, 37(2): 45-47.

[21] 徐昱琳, 张雷, 陈万米,等. 基于SIFT 特征和图割算法的图像分割方法研究 [J]. 计算机测量与控制, 2009, 17(12): 2491-2493.

[22] Kolmogorov V, Zabih R. What energy functions can be minimized via graph cuts? [J]. IEEE Transactions on Pattern Analysis and Machine Intelligenc, 2004, 26(2): 147-159.

[23] Boykov Y, Veksler O, Zabih R. Fast approxamate energy minimization via graph cuts [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2001, 23(11): 1222-1239.

[24] Boykov Y, Kolmogorov V. An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2004, 26(9): 1124-1137.

Pupil Ultrasound Images Segmentation and Diameter Measurements Based on Improved Graph Cuts

Zhou Yali1Gao Yuzhi2Zhou Hao1Zhou Mi1Zheng Yinfei1*

1(BiomedicalEngineeringandInstrumentScience,ZhejiangUniversity,Hangzhou310027,China)2(EmergencyDepartmentofSecondAffiliatedHospitalofZhejiangUniversity,Hangzhou310027,China)

The traditional measurement of pupil diameter is manually measured, but for patients who have ocular trauma or loss of consciousness this method is inconvenience. Aiming to solve the problem of large manual interaction in pupil diameter measurement and the weak measuring robustness, we used improved graph cut segmentation algorithm to segment pupil ultrasound images and measure pupil diameters. In this paper, we improved the traditional graph cut algorithm in two aspects. One is using adaptive threshold region growing to take place of manually seeds selection, which ensures the segmentation results while reducing the amount of manual interaction. The other one is increasing gradient information of the image into data entry portion of the energy function, which reduces the small area in the segmentation results and enhances the weak edges of segmentation. The proposed method realized automatic segmentation of image and automatic measurement of the pupil diameter. By employing the method we acquired dynamic changes of patient pupil diameter and provide a basis for clinical diagnosis. To verify the validity of the algorithm, we used the method to measure the diameter of ten patients’ dynamic pupil ultrasound images, and compared the results with that obtained by manual measurements. It is shown that the absolute error is less than 0.2 mm, the correlation coefficient is at least 0.83. In conclusion, the modified graph-cut algorithm improved segmentation results, achieved automatical measurements of pupil diameters using dynamic ultrasound pupil images, and can be expected to substitute manual measurements and reduce the amount of human interaction.

graph-cut; pupil ultrasound image; image segmentation; pupil diameter

10.3969/j.issn.0258-8021. 2015. 04.003

2015-04-01， 录用日期:2015-05-30

中央高校基本科研业务费专项资金(2015FZA5019)；国家“十二五”科技支撑计划资助项目(2011BAI12B02)

R318

A

0258-8021(2015) 04-0399-08

*通信作者 (Corresponding author)，E-mail: zyfnjupt@126.com