本发明属于图像识别技术领域，涉及一种基于主动轮廓和能量约束的筛板前表面深度测量方法。
背景技术：
青光眼是世界第二大致盲疾病，它通过对人视乳头周围的视神经节细胞轴突的破坏，导致人的视野缺失。由于青光眼的不可逆性，对青光眼的早检测、早发现和早治疗能够减缓疾病发展的进程。但因为青光眼的发病机制尚未完全明确，因此对青光眼的风险因素的研究仍然是目前的热点问题。
视神经筛板(laminacribrosa)位于巩膜，是视网膜神经节细胞轴突穿出眼睛的位置，呈筛状结构，负责对视网膜神经节细胞轴突提供营养和结构支持，并有助于保持眼内和眼外之间的压力梯度。筛板的解剖学结构如下：大约1.2-1.5万个视网膜神经节细胞轴突聚集在视乳头，通过视神经管的内部(即bruch's膜开口)和外部(即巩膜)部分离开眼睛。视神经是中枢神经系统的一部分，外面包有由三层脑膜延续而来的三层被膜(硬脑膜、蛛网膜、软脑膜)，硬脑膜下腔和蛛网膜下腔也随之延续到视神经周围，蛛网膜下腔充满脑脊液。在巩膜后孔处，巩膜组织前1/3形成薄层纤维隔，隔板内有很多筛孔允许视网膜神经节细胞轴突穿行，后2/3与硬脑膜相接。研究表明筛板与青光眼最大的风险因素——眼内压存在一定关系，在青光眼的发病过程中通常会伴随着眼内压的升高，而眼内压的升高同时会造成筛板的向后移位和弯曲变形等。王宁利等人提出视神经节细胞轴突损伤与跨筛板压力差(眼内压与脑内压之差)有关，而在这种压力差之下筛板也会表现出向后移位的特征。
光学相干断层扫描技术(opticcoherencetomography,oct)技术在眼科的临床应用不过二十年，其技术创新迅速，现已经成为眼科临床最重要的检查之一。该技术通过对组织发射相干光，回收组织的反射光与散射光，并与时间的延迟相结合，可以获取到组织的二维断层信息或整体视网膜三维图像的信息。除了实时监测，无创等优点外，oct最重要的特点是其分辨率高，并且观察到的细微结构为横截面结构，符合病理学的常规观察习惯，为眼科研究者的活体视网膜形态研究提供技术支持。oct的工作过程如下：由低相干光源发出的低相干光线，由干涉仪分成两束光线，一束进入探测光路，直接入射眼内，被眼内不同组织的界面发射回来，用以提供眼内各种组织的厚度与距离信息；另一束进入参照光路，由已知空间距离的参照镜反射回来。两束光在光纤偶联器中整合为一，由于发射到参照镜的来回距离和发射到眼内给定结构的距离精确匹配，故产生干涉现象，被光敏探测器探测到。调解后的信号输入计算机进行运算，得到被测物体的光学相干断层图像。由于眼内组织具有不同的深度和空间结构，所以两束光线之间会产生时间差，这种时间差叫做光学延迟时间。把探测器所计算出的时间差利用低相干光涉度量学原理，可获得组织反射的信息。获取到这些信息后，由计算机计算出一维的扫描图像信息，通常为二维图像中的一列。
近年来，oct技术的进步逐渐解决了深层组织无法成像的问题。海德堡公司oct的加强深度成像(enhanceddepthimaging,edi)模式能够将光线汇聚在眼底的更深层位置，如巩膜、脉络膜和筛板等。这种成像技术的出现有利于对活体筛板的病理性变化如筛板向后移位等进行进一步研究，提高青光眼早期诊断的精确度。我们通常使用筛板前表面深度来表示筛板的向后移位，这种筛板深度的测量方法为：以bruch膜开口点(bruch’smembraneopening,bmo)组成的bmo参考平面作为基准面，测量基准面与筛板前表面(即筛板组织与筛板前组织之间的分界面)之间的距离。因此，筛板前表面深度的测量分为两步：bmo测定和筛板前表面分割。得到bmo参考平面和筛板前表面的位置后，筛板前表面深度可以通过测量两个面之间的距离得出。
由于对筛板的研究才刚刚起步，一些其他现有的方法也无法很好地应用于筛板图像处理中，所以相关的研究成果相对较少。
a)bmo测定
bmo点是oct图像中最重要的生物标记之一，它在oct断层图像中代表视盘边缘的位置，因此被广泛应用于视盘分割，盘沿测量等问题中。
2013年，fu等人提出了一种基于字典学习和低秩矩阵重建的bmo测定方法，该方法首先分割出视网膜的内界膜(ilm)层和视网膜色素上皮(rpe)层并把图像分为训练区域和候选区域，然后用图像的训练区域训练出低秩字典，并用候选区域的数据重建出整个图像，得到重建误差曲线，包括亮度误差，局部二值模式(lbp)误差和距离偏置。最后通过三个误差曲线得出视盘的位置即bmo点的位置。该方法的缺点在于无法有效测定青光眼病变的bmo位置。
2014年，akram等人提出了一种基于反卷积的模型来测定bmo的位置。该模型认为oct图像中的每一层可以由两个成分建模，分别是：用一条曲线来建模层的骨架和用一个或一组滤波器来建模层的厚度。只要得到了这两个成分的参数，就可以唯一确定oct图像上每一层的位置。该模型使用蒙特卡洛马尔科夫链(montecarlomarkovchain,mcmc)方法进行参数测定，达到了较高的准确率，但是同时也带来了效率较低的问题。
2015年，wang等人提出了一种bmo的测定方法，该方法使用了五种方法分别测定bmo的位置，第一种方法基于ilm层的形状特征，其余四种方法基于bmo附近位置的纹理特征。得到五种方法的结果后，选择其中三个最接近的值并求它们的平均值作为最终的bmo测定结果。该方法存在效率较低，鲁棒性较差等问题。
hussian等人提出了一种基于图论的bmo测定方法，该方法首先测出三个参照层的位置，分别是视网膜神经纤维(rnfl)层，外丛状层(onl)层和rpe层。然后使用层的位置和特征作为图的权重，使用一种图搜索算法得出bmo的位置，达到了较高的准确率。
b)筛板前表面分割
筛板前表面是筛板组织与筛板前组织的分界面，由于人的筛板形态各异，没有固定的形状特征，对比度不明显，有时还会因为青光眼等疾病导致缺失，所以传统的基于梯度特征、颜色和形状的分割方法效果并不好，一些现有的层次分割方法也不适用于筛板前表面分割中。
2015年，akram等人提出了第一种筛板前表面的自动分割方法，该方法基于马尔科夫随机场和mcmc采样算法，将分割问题转化为分类问题，只需要知道每一个点是否属于筛板前表面即可。以每一个点为单位，提取出该点及其邻居点的特征来建立一个马尔科夫随机场，然后利用mcmc算法来测定马尔科夫随机场的参数，并通过迭代采样得到最终的表面。
然而，上述现有方法通常存在效率较低，对青光眼病变图像不鲁棒等问题，导致它们在青光眼的筛查和临床诊断中的实用性较差。
技术实现要素：
本发明提供了一种基于主动轮廓和能量约束的筛板前表面深度测量方法，其目的在于，克服现有技术中效率低以及鲁棒性不强的问题。
一种基于主动轮廓和能量约束的筛板前表面深度测量方法，包括以下步骤：
步骤1：将待检测的灰度图像进行反色处理后，对反色处理图像中的所有像素点采用k均值聚类方法进行聚类，得到聚类图；
步骤2：构造c-v主动轮廓模型的能量函数，利用聚类图的轮廓作为c-v主动轮廓模型的初始曲线，以c-v主动轮廓模型的能量函数最小时的曲线作为目标轮廓，以目标轮廓的下边界的终止点作为bmo点，连接两个bmo点得到bmo参考投影线；
步骤3：从反射处理图像中选取位于bmo参考投影线以下的区域，以所选区域对应的最小矩形区域作为筛板前表面分割的感兴趣区域，其中，所选区域的上边界为bmo参考投影线的在水平方向上的投影线；
步骤4：基于能量函数从感兴趣区域中各列中提取一个筛板前表面候选点，并依据完整性约束条件，对所有的筛板前表面候选点进行剔除处理，得到控制点集；将控制点集进行曲线拟合得到筛板表面轮廓线；
步骤5：基于bmo参考投影线和筛板表面轮廓线，计算筛板前表面深度；
以bmo参考投影线的中点及其鼻侧和颞侧各100微米的两个点作为测量点，在三个测量点上对bmo参考投影线作垂线并与筛板前表面轮廓线相交，三条垂线长度的平均值为待测量的筛板前表面深度。
进一步地，所述聚类图为m(i,j,k)，是将反色处理图像中属于像素值最小的聚类中心代表的类的像素点保留，其余像素点的灰度值置0：
其中，cm和ck分别表示第m个和第k个聚类中心，且cm代表像素点(i,j)的当前聚类中心，i(i,j)为图像中的像素点(i,j)的像素值，i(cm)和i(ck)分别为聚类中心cm和ck处的像素值。
进一步地，所述c-v主动轮廓模型的能量函数为：
e(c)＝μl(c)+v*area(inside(c))+λ∫inside(c)(i(x,y)-ii)2dxdy
+ω∫outside(c)(i(x,y)-io)2dxdy
其中，i(x,y)为聚类图中的像素点(x,y)的像素值，l(c)代表曲线的长度，area(inside(c))表示轮廓曲线所围成的面积，ii和io分别为轮廓外和轮廓内的像素点的像素值；inside(c)表示轮廓内的像素点，outside(c)表示聚类图中轮廓外的像素点；μ为长度参数，取值范围为(0,1)；v为面积参数，取值为0；λ和ω分别为内能量和外能量的权重系数，λ＝ω＝1。
进一步地，所述基于筛板前表面分割的能量函数从感兴趣区域中各列中按照以下公式提取一个筛板前表面候选点：
其中，sj表示第j列中的筛板前表面候选点，n表示感兴趣区域中像素的总列数；ij(i)表示感兴趣区域中第j列像素的第i个像素点的亮度特征值，gradient(ij(i))表示ij(i)的梯度特征值，α和β分别表示像素点的亮度特征和梯度特征的权重参数，α的取值范围为[0.3，0.4]，β的取值范围为[0.5，0.6]，argmax()表示寻找使能量函数最大的点的位置。
感兴趣区域中像素的总列数即为感兴趣区域图像的宽度；
argmax()即为使得α*ij(i)+β*gradient(ij(i))取得最大值时，对应的像素点作为sj；
进一步地，所述完整性约束条件为：
其中，control(j)表示控制点集，d表示偏移距离控制参数，取值范围为(10,15)；s为整个候选点纵坐标集合，yj代表第j个a-scan上的候选点的纵坐标，average(s)为整个候选点集合纵坐标的平均值。
进一步地，采用b样条拟合算法将控制点集进行曲线拟合得到筛板表面轮廓线。
有益效果
本发明提出了一种基于主动轮廓和能量约束的筛板前表面深度测量方法，该方法基于主动轮廓模型，提取出oct图像中的bmo点和bmobmo参考投影线，再使用能量约束方法分割出筛板前表面。最后利用bmo测定和筛板前表面分割的结果，测量出筛板前表面深度。本方法在bmo点测定中使用了基于区域的主动轮廓模型，避免使用传统的梯度特征使测定结果受到血管阴影，病变等异常的影响；在筛板前表面分割中，本方法采用了能量函数的思想，解决了筛板前表面没有固定的形状，大小等传统特征的问题，并采用曲线拟合来拟合筛板的缺失部分；实验结果表明，该方法在bmo测定的精确度和筛板前表面分割的准确率和平滑性上均优于现有方法，更为接近专家手动标定的结果。本发明能够解决专家在临床诊断时需要手动标定测量筛板前表面深度的费时费力的问题，对青光眼的早期筛查和临床诊断具有积极意义。
附图说明
图1为本发明的流程图，其中，(a)为其逻辑流程图，(b)为其执行步骤的流程图；
图2为筛板前表面深度测量的示意图；
图3为本方法的bmo测定结果，其中，(a)为本方法的测定结果图1；(b)为与(a)对应的专家手工标定的结果；(c)为本方法的测定结果图2；(d)为与(c)对应的专家手工标定的结果；(e)为本方法的测定结果图3；(f)为与(e)对应的专家手工标定的结果；
图4为本方法的筛板前表面分割结果；
图5为本方法测量的筛板前表面深度与专家手工标定的对比示意图，其中，(a)为本方法的测量结果图1；(b)为与(a)对应的专家手工标定的结果；(c)为本方法的测量结果图2；(d)为与(c)对应的专家手工标定的结果；(e)为本方法的测量结果图3；(f)为与(e)对应的专家手工标定的结果。
具体实施方式
以下结合附图对本发明做进一步说明。
本发明提出了一种基于主动轮廓和能量约束的筛板前表面深度测量方法，包括四个步骤：首先利用k均值算法的聚类图作为c-v主动轮廓模型的初始轮廓，再取轮廓下边界的终止点，即为bmo点，然后提取出图像中包含bmo参考投影线以下区域的最小矩形作为筛板前表面分割的感兴趣区域，对图像感兴趣区域的每一列根据能量函数提取出候选点，将不符合约束条件的候选点删去，并用曲线拟合方法分割出整个表面，最后计算出筛板前表面深度。
具体流程如图1(a)所示，该方法的执行过程如图1(b)所示，其具体实施步骤如下：
步骤1：将待检测的灰度图像进行反色处理后，对反色处理图像中的所有像素点采用k均值聚类方法进行聚类，得到聚类图；
给定反色后的待检测图像i(m,n)，k均值算法把图像中的m×n个像素分成预定义的k个类，根据每个类类内像素的差异性较小，而类与类之间的差异性较大的原理，我们可以把问题转化为使价值函数(代表像素之间的差异性)最小的问题，然后得到每一类的聚类中心，使得每个点与最近的聚类中心的差异性最小。所使用的价值函数如下：
其中c为预定义的类别数，d(i,j,k)＝||i(i,j)-i(ck)||2为描述图像上两个点之间的差异度的度量函数，因为oct图像具有层次性特征，故使用图像的像素值特征作为差异性度量，i(i,j)为图像中的像素点，ck为聚类中心。那么使价值函数最小的聚类图m(i,j,k)根据下公式计算得出。
其中，cm和ck分别表示第m个和第k个聚类中心，且cm代表像素点(i,j)的当前聚类中心，i(i,j)为图像中的像素点(i,j)的像素值，i(cm)和i(ck)分别为聚类中心cm和ck处的像素值。将图中属于像素值最小的聚类中心代表的类的像素点保留，其余的置0，作为待测定聚类图。得到聚类图后，每一类的所有像素的中心点就是新的聚类中心。
步骤2：构造c-v主动轮廓模型的能量函数，利用聚类图的轮廓作为c-v主动轮廓模型的初始曲线，以c-v主动轮廓模型的能量函数最小时的曲线作为目标轮廓，以目标轮廓的下边界的终止点作为bmo点，连接两个bmo点得到bmo参考投影线；
由于主动轮廓模型的结果很大程度上受到所给定的初始曲线的影响，因此将上一步所得到的k均值聚类图作为主动轮廓模型的初始曲线c，c-v主动轮廓模型的能量函数如下式：
e(c)＝μl(c)+v*area(inside(c))+einside+eoutside
＝μl(c)+v*area(inside(c))+λ∫inside(c)(i(x,y)-ii)2dxdy
+ω∫outside(c)(i(x,y)-io)2dxdy
其中，i(x,y)为聚类图中的像素点(x,y)的像素值，l(c)代表曲线的长度，area(inside(c))表示轮廓曲线所围成的面积，ii和io分别为轮廓外和轮廓内的像素点的像素值；inside(c)表示轮廓内的像素点，outside(c)表示聚类图中轮廓外的像素点；μ为长度参数，它的取值由图像中目标的大小决定，如果目标的尺寸较大，则μ的取值也较大，反之μ的取值就会较小，取值范围为0-1，通常取μ＝0.5。v为面积参数，λ和ω分别为内能量和外能量的权重系数。通常取λ＝ω＝1,v＝0。c-v模型通过构造能量函数得到其能量的分布，以能量函数不断减小的方法来推动曲线逼近目标的边缘，最终将目标从背景中分割出来。它摆脱了经典边缘分割算法如log和传统snake模型中对图像梯度的依赖，对oct图像具有很好的分割能力。
步骤3：从反射处理图像中选取位于bmo参考投影线以下的区域，以所选区域对应的最小矩形区域作为筛板前表面分割的感兴趣区域，其中，所选区域的上边界为bmo参考投影线的在水平方向上的投影线；
获得bmo参考投影线的位置后，以bmo参考投影线的水平投影线作为感兴趣区域的上边界，即矩形区域的宽度；区域下边界为图像底部，左右边界分别为两个bmo点的位置。
步骤4：基于能量函数从感兴趣区域中各列中提取一个筛板前表面候选点，并依据完整性约束条件，对所有的筛板前表面候选点进行剔除处理，得到控制点集；将控制点集进行曲线拟合得到筛板表面轮廓线；
围绕筛板前表面的亮度特征和梯度特征来构建筛板前表面分割的能量函数，并在每一个a-scan(即每一列，下同)中根据能量函数选择一个筛板前表面候选点，如下式：
其中，sj表示第j列中的筛板前表面候选点，n表示感兴趣区域中像素的总列数；ij(i)表示感兴趣区域中第j列像素的第i个像素点的亮度特征值，gradient(ij(i))表示ij(i)的梯度特征值，α和β分别表示像素点的亮度特征和梯度特征的权重参数，α的取值范围为[0.3，0.4]，β的取值范围为[0.5，0.6]，argmax()表示表示寻找使能量函数最大的点的位置。得到候选点集合之后，由于筛板的形态特征因人而异，以及血管阴影的影响，有一些候选点事实上不在筛板前表面的位置上，需要将这部分候选点的共同特征找出来并将其删去。完整性约束条件如下：
其中，control(j)表示控制点集，d表示偏移距离控制参数，用这个参数控制需要移除的偏移值过大的候选点，取值在10-15之间。s为整个候选点纵坐标集合，yj代表第j个a-scan上的候选点的纵坐标，average(s)为整个候选点集合纵坐标的平均值。为了保证得到的筛板前表面平滑，我们在得到的控制点集control(j)上进行曲线拟合，使用的拟合算法为b样条拟合算法。
步骤5：基于bmo参考投影线和筛板表面轮廓线，计算筛板前表面深度；
得到bmo测定结果和筛板前表面分割结果后，以bmo参考投影线的中点及其鼻侧和颞侧各100微米的两个点作为测量点，在三个测量点上对bmo参考投影线作垂线并与筛板前表面相交。如图2，三条垂线长度的平均值即为待测量的筛板前表面深度。
为验证本发明提出的方法的性能，采集18个样本的30幅图像作为数据集，并以专家手工标定的bmo点和筛板前表面作为对应的groundtruth图。基于步骤3所述，α和β参数很大程度上决定了方法是否有效，通过实验发现，当α增大时，曲线会接近筛板区域；当β增大时，曲线会接近内界膜区域。为了保证方法的鲁棒性，设置α＝0.3-0.4,β＝0.5-0.6。
本发明将提出的计算方法与2015年hussian等人提出的bmo测定方法及nlsc筛板前表面分割方法进行比较。在bmo测定中，我们以测定点和groundtruth的平均误差和方差作为评价标准；在筛板前表面分割中，我们则采用错误率作为评价标准。
表1bmo测定结果
从表1中可以看出，本发明提出的方法在精确度上(平均误差48.17微米)要优于hussian等人所提出的bmo测定方法(平均误差54.18微米)，而在稳定性上，本方法的标准差是51.32微米，也略优于它的53.74微米。图3列出了本方法得到分割结果与专家手动分割的bmo点及bmo参考平面，从图3中可以看出由于避开了梯度特征，本发明可以很好的避免由大血管阴影带来的图像模糊问题，在血管阴影多的edi-oct图像上也能得到非常准确的分割结果，并与专家手动标定的结果相一致。
在筛板前表面分割中采用筛板前表面样点的错误率(与groundtruth存在误差的点个数与样本总数之比)作为评价标准，错误率等级从0到2，定义为：等级0代表平均误差小于3个像素的点占样本总数之比，等级1代表平均误差在3个像素和5个像素之间的点占样本总数之比，等级2代表平均误差大于5个像素的点占样本总数之比。将本发明提出方法的分割结果与一种基于马尔科夫链蒙特卡洛的bfps边缘分割方法以及目前唯一的筛板前表面自动分割方法——nlsc方法进行对比，结果如表2所示。
表2筛板前表面分割结果
在每张图片上采样100-150个评价点。从表2中可以看出与groundtruth相比，在错误率等级0和等级1中，本发明所提出方法分别达到76.5％和17.6％，高于nlsc方法的73.7％和16.1％和bfps方法的40.9％和22.8％，而在误差较高的错误率等级2中，本发明所提出方法则只有5.9％，要优于nlsc方法的10.2％和bfps方法的36.3％。由此可见本发明所提出的方法在整体上优于nlsc方法和bfps方法。图4列出了实验结果与groundtruth的对比，可以看出实验结果与专家手动标定的结果在图像上非常接近，而且曲线也非常平滑，从组织学上来看也具有合理性。
最后本发明方法所测得的筛板深度与groundtruth进行比较，仍然使用平均误差和标准差作为评价标准，其结果如表3所示。
表3筛板前表面自动测量结果
从表3中可以看到，本发明测得的筛板深度和专家标定的结果误差较小，基本一致，图5给出了较为直观的结果，从图中也能够看出本方法与专家的手工标定基本一致，且准确率较高，平均有符号误差为22.1微米，无符号误差为28.7微米。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。