发明名称 基于Schwalbe线的眼前房角参数自动测量方法
摘要 本发明涉及一种基于Schwalbe线的眼前房角参数自动测量方法。本发明首先对图像进行增强与分割操作;根据房角图像特性分别检测角膜前边缘,并对图像进行折射校正、进而检测角膜后边缘和后弹力膜前边缘,以及虹膜前边缘;利用角膜后弹力膜的形状和梯度特征,及其与角膜后边缘的位置关系实现对Schwalbe线位置的准确检测;最后结合所检测到的角膜后边缘与虹膜前边缘可计算出最终的前房角相关解剖结构参数的测量结果。该发明能实现Schwalbe线准确位置的自动检测和前房角相关解剖结构参数的精确测量,并减少手动测量各参数所引入的误差。
申请公布号 CN103971348A 申请公布日期 2014.08.06
申请号 CN201410138890.2 申请日期 2014.04.08
申请人 杭州电子科技大学 发明人 武薇;陆晓娟;王梦蕾;范影乐
分类号 G06T5/20(2006.01)I;G06T5/00(2006.01)I;G06T7/00(2006.01)I;A61B3/117(2006.01)I 主分类号 G06T5/20(2006.01)I
代理机构 杭州求是专利事务所有限公司 33200 代理人 杜军
主权项 基于Schwalbe线的眼前房角参数自动测量方法,其特征在于该方法包括如下步骤:步骤(1)对获得的前房角OCT图像通过配准后进行帧平均操作;步骤(2)对帧平均图像进行增强和图像分割操作,具体是:2‑1、对帧平均图像I采用7×3的自适应Wiener滤波器滤波,得到图像P;2‑2、对图像P做形态学灰度重建得到I<sub>R</sub>,其计算范围为各像素点的四邻域;2‑3、将I<sub>R</sub>按比例缩小至原图的四分之一大小得到I'<sub>R</sub>;2‑4、对I'<sub>R</sub>图像进行幂次变换得到图像F,其中将输入灰度值范围为[0.06,0.9]的像素映射至输出灰度值范围[0,1],幂次参数γ=0.18;2‑5、对图像F采用Otsu自动分割算法进行图像分割,得到分割后的图像I'<sub>BW</sub>;步骤(3)根据前房角结构特性,对分割后的图像和原始图像的梯度图分别做边缘粗检测和细检测;其中在检测到角膜前边缘的基础上对图像进行折射校正,而后再检测出虹膜前边缘、角膜后边缘以及角膜后弹力膜前边缘,具体是:定义边缘为垂直方向上连续的6个取值为0的像素与连续的6个取值为1的像素交界处;基于上述的边缘定义,分割后的图像的边缘粗检测具体是:3‑1、从图像I'<sub>BW</sub>顶端到底部逐行扫描,检测得到第一条初始边缘,即为角膜前边缘,根据第一条边缘对图像I'<sub>BW</sub>做折射校正,得到校正后的图像I'<sub>BWd</sub>;3‑2、从I'<sub>BWd</sub>底部向上至第一条边缘逐行扫描检测第三条边缘,即为虹膜前边缘;3‑3、从第三条边缘向上逐行扫描至第一条边缘,检测得到第二条边缘,即为角膜后边缘;3‑4、对检测到的三条边缘分别用局部加权线性最小二乘的二次多项式回归模型进行平滑,得到初始边缘;使用双三次插值法将得到的初始边缘在横坐标方向按比例放大两倍,得到参考边缘;基于放大后的初始边缘,在原始图像的梯度图中重新进行边缘细检测,得到最终的精确边缘,其步骤如下:3‑5、搜索范围以参考边缘为中心,向上和向下分别以n<sub>1</sub>和n<sub>2</sub>的距离搜索相应梯度最大值;对于搜索到的新边缘在垂直方向上的坐标值应加上一个偏移量Δd;按照边缘粗检测的步骤先检测出第一条边缘,进而对梯度图进行折射校正;3‑6、继续检测第二和第三条边缘,分别对检测出来的三条边缘采用局部加权线性最小二乘的二次多项式回归模型进行平滑处理;在梯度图中,基于检测到的第二条边缘,向其上方距离为[n'<sub>1</sub>,n'<sub>2</sub>]的区域内搜索梯度最大值位置,对搜索到的边缘采用局部加权线性最小二乘的二次多项式回归模型做平滑处理,得到最终的角膜后弹力膜前边缘;步骤(4)利用角膜后弹力膜的形状和灰度特征对Schwalbe线的检测区间做出限定,并在限定区间内进一步检测得到准确的Schwalbe线位置,具体是:4‑1、分别计算后弹力膜前边缘上每一像素点E<sub>t</sub>(x,y)与其上方n<sub>t</sub>个像素之间的灰度差值g(x,y‑i),i=1,2,…,n<sub>t</sub>;x=1,2,…,M;y=1,2,…,N;M和N分别为图像宽度和高度;4‑2、计算后弹力膜前边缘与角膜后边缘间距离D<sub>E</sub>(x,y);4‑3、按式(1),计算后弹力膜前边缘上每一个像素点所对应的特征值f(x,y)<maths num="0001" id="cmaths0001"><math><![CDATA[<mrow><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><mi>A</mi><mo>&CenterDot;</mo><msub><mi>D</mi><mi>E</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>&CenterDot;</mo><mfrac><mn>1</mn><msub><mi>n</mi><mi>t</mi></msub></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>n</mi><mi>t</mi></msub></munderover><mo>[</mo><mi>g</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>-</mo><mi>i</mi><mo>)</mo></mrow><mo>&CenterDot;</mo><mi>w</mi><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow><mo>]</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000488252260000021.GIF" wi="1308" he="132" /></maths>其中A为常数;w为对应每一个g(x,y‑i)值的权值;4‑4、对所有特征值f(x,y)进行8阶多项式拟合,得到方程F(x,y);4‑5、当后弹力膜的位置接近Schwalbe线时,f(x,y)取值为最大;通过计算F(x,y)的一阶和二阶导数得到F(x,y)取最大值的坐标点,该点的横坐标x即为s<sub>begin</sub>;在s<sub>begin</sub>右侧的第一个极小值即为s<sub>end</sub>;倘若无法检测到极小值,则取s<sub>begin</sub>右侧F(x,y)第一个取值小于阈值th<sub>s</sub>的点的横坐标x为s<sub>end</sub>;若仍无法得到s<sub>end</sub>,则取s<sub>end</sub>=c<sub>p</sub>,其中c<sub>p</sub>为第三条边缘的拟合直线经最近邻插值获得的两倍放大直线代替虹膜前边缘,并求其与角膜后边缘的交点;在区间[s<sub>begin</sub>,s<sub>end</sub>]中精确定位Schwalbe线位置,其步骤如下:在区间[s<sub>begin</sub>,s<sub>end</sub>]中查找D<sub>E</sub>(x,y)内取值最大的像素,并得到其中连续性最长的一段像素中最后一个像素坐标值;即搜索位于两条边缘相距距离最大,且离Schwalbe线最近的弹力膜上边缘点的横坐标位置x<sub>begin</sub>;从x<sub>begin</sub>向右继续搜索D<sub>E</sub>(x,y)内取值最小的元素,该点即为Schwalbe线的横坐标位置;由此便可计算得到Schwalbe线在图像中的精确位置;步骤(5)结合所检测到的Schwalbe线位置、角膜后边缘和虹膜前边缘位置计算最终的前房角结构测量参数,具体是:5‑1、计算基于Schwalbe线的房角开放距离AOD_SL、AOD_SL200、AOD_SL550;沿角膜后边缘由Schwalbe线向小梁网方向分别延伸200μm和550μm得到点S<sub>200</sub>和S<sub>550</sub>;分别在Schwalbe线、S<sub>200</sub>和S<sub>550</sub>左右两侧距其各30个像素点的位置对其附近的角膜后边缘进行直线拟合,以此为距离垂线的基准线,分别做垂线得到AOD_SL、AOD_SL200、AOD_SL550这三个距离参数在虹膜前边缘的落点R<sub>SL_i</sub>的坐标(x<sub>r</sub>,y<sub>r</sub>),其中i=1,2,3,开放距离按式(2)计算:<maths num="0002" id="cmaths0002"><math><![CDATA[<mrow><mi>AOD</mi><mo>=</mo><msqrt><msup><mrow><mo>(</mo><msub><mi>x</mi><mi>s</mi></msub><mo>-</mo><msub><mi>x</mi><mi>r</mi></msub><mo>)</mo></mrow><mn>2</mn></msup><mo>+</mo><msup><mrow><mo>[</mo><mrow><mo>(</mo><msub><mi>y</mi><mi>s</mi></msub><mo>-</mo><msub><mi>y</mi><mi>r</mi></msub><mo>)</mo></mrow><mo>&CenterDot;</mo><msub><mi>n</mi><mi>c</mi></msub><mo>]</mo></mrow><mn>2</mn></msup></msqrt><mo>&CenterDot;</mo><msub><mi>d</mi><mi>p</mi></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000488252260000031.GIF" wi="1137" he="101" /></maths>其中(x<sub>s</sub>,y<sub>s</sub>)为S(S<sub>200</sub>或S<sub>550</sub>)的坐标,(x<sub>r</sub>,y<sub>r</sub>)为R<sub>SL_i</sub>坐标,n<sub>c</sub>为角膜折射系数,d<sub>p</sub>为图像分辨率;5‑2、计算基于Schwalbe线的小梁虹膜空间面积TISA_SL500和TISA_SL700;从Schwalbe线向小梁网方向分别延伸500μm和700μm得到点S<sub>500</sub>和S<sub>700</sub>;分别在Schwalbe线S、S<sub>500</sub>和S<sub>700</sub>左右两侧距其各30个像素点的位置对其附近的角膜后边缘进行直线拟合,以此为距离垂线的基准线,分别做垂线与虹膜相交,得到三条垂线段SI、S<sub>500</sub>I<sub>500</sub>和S<sub>700</sub>I<sub>700</sub>;分别令S和S<sub>500</sub>,S和S<sub>700</sub>间的角膜后边缘表示为SS<sub>500</sub>和SS<sub>700</sub>;再将I和I<sub>500</sub>,I和I<sub>700</sub>间的虹膜前边缘表示为II<sub>500</sub>和II<sub>700</sub>;最后,分别计算由线段SS<sub>500</sub>、S<sub>500</sub>I<sub>500</sub>、II<sub>500</sub>和SI,以及由线段SS<sub>700</sub>、S<sub>700</sub>I<sub>700</sub>、II<sub>700</sub>和SI所包围的区域面积即为TISA_SL500和TISA_SL700;5‑3、计算从Schwalbe线开始的房角隐窝面积ARA_SL;取参数SI,并计算由SI和靠隐窝方向一侧的角膜后边缘和虹膜前边缘之间区域的面积,即为所求ARA_SL;其中,若角膜后边缘与虹膜前边缘不相交,则取两者最靠近隐窝一端各50个像素分别做一阶曲线拟合,并令拟合曲线为相应坐标位置处的边缘;5‑4、Schwalbe线下虹膜厚度IT;过虹膜后边缘点I<sub>B</sub>做垂线,使该垂线与点I相交,该垂线段长度即为IT;该垂线需与在I<sub>B</sub>左右两侧虹膜后边缘上各取30个像素点拟合得到的直线相垂直。
地址 310018 浙江省杭州市下沙高教园区2号大街