CN103810709B - 基于血管的眼底图像与sd-oct投影图像配准方法 - Google Patents

Abstract

本发明公开了一种基于血管的眼底图像与频域光学相干断层投影图像自动配准方法，属于图像处理技术领域。该方法首先对输入的SD‑OCT图像进行组织层分割，定位光感受器内节/外节层上边界和视网膜色素上皮下边界，根据IS‑OS层上边界和RPE层下边界生成SD‑OCT投影图像，利用改进的水平集方法分割SD‑OCT投影图像，然后采用形态学多尺度顶帽方法分割彩色眼底图像，最后使用迭代分块方法搜索最大相似函数的局部最大值来提取控制点，并根据控制点计算仿射变换系数实现SD‑OCT投影血管图像与彩色眼底血管图像的配准。实验结果表明，本发明中所得到的不同模态下的图像配准精度高。

Description

基于血管的眼底图像与SD-OCT投影图像配准方法

技术领域

本发明涉及一种图像配准处理的方法，特别是一种基于血管的眼底图像与SD-OCT投影图像的配准方法。

背景技术

眼底照相可以快速获取不同视野范围的彩色眼底图像，而SD-OCT能有效地呈现视网膜组织层的细微变化(厚度变化、组织层结构变化)，不同的成像设备可以提供不同的图像信息。图像配准技术能够将来自不同设备的信息进行结合。图像配准方法大体上可分为两类：

(1)基于灰度的方法通常利用图像灰度信息的不同实现图像的配准。利用图像分割或根据相关性、互信息测量等测度作为目标函数。然而基于灰度的方法容易陷入一个局部极小值，而且变换模型的搜索空间非常大。

(2)基于特征的配准是以待配准图像对之间存在的匹配点为基础的，通过匹配点之间的相似度来完成图像配准。

由于彩色眼底图像与SD-OCT投影图像的灰度不在同一个灰度级内，很难用基于灰度的方法实现图像的配准。而血管是这两种模态下图像共有的特征，所以采用基于特征的方法完成图像配准。

目前国际上只有一种方法能够实现彩色眼底图像与SD-OCT投影图像的配准，该方法将血管脊作为配准的特征，然后定义一种距离函数，在函数中加入血管脊的法方向信息，计算相似匹配点的距离，ICP搜索相似函数的最大值，从而完成图像的配准。但是该方法的计算复杂度较高。

发明内容

本发明的目的在于提供一种基于血管的眼底图像与SD-OCT投影图像配准方法。

实现本发明的目的的技术解决方案为：一种基于血管的眼底图像与SD-OCT投影图像配准方法，包括以下步骤：

步骤1、采集彩色眼底图像与三维SD-OCT视网膜图像；

步骤2、用基于双梯度信息和空间约束信息的视网膜组织层分割算法对SD-OCT视网膜图像进行IS-OS层和RPE层分割；具体为：

步骤2-1、采用基于非局部均值思想的三维块匹配去噪算法对SD-OCT图像进行降噪处理，然后用基于边缘流的方法提取视网膜图像的边缘，所述边缘的像素值为255，其它为黑色背景；

步骤2-2、对于步骤2-1结果的每个A-scan，从上向下搜索第一个值为255的像素点，提取ILM边界；所述A-scan表示图像的列，ILM为视网膜图像内界膜；之后在垂直方向上使用确定图像梯度d2l来修正ILM边界，I为视网膜图像，具体方法如下：

定义ILM边界位置向量ILMposition(x)，x∈[1,M]，M为A-scan方向的宽度；之后判断相邻两列间的位置差：Diff＝|ILMposition(x)-ILMposition(x-1)|，如果位置差大于给定的阈值T，则用当前点的局部邻域内最大梯度位置修正该点的边界值，即

$ILMposition\left(x\right)=\left\{\begin{array}{cc}\underset{z\∈S}{max}d2l(z,x)& \begin{array}{cc}if& Diff>T\end{array}\\ ILMposition\left(x\right)& Others\end{array}\right.$

其中S是以ILMposition(x)为中心在x列窗口大小为l区域内的像素点集合，d2l(z,x)为由暗到亮梯度图像，z表示图像的高度，x表示图像的宽度，这里所述的阈值T＝20；

将ILM边界位置向量构成ILM位置矩阵ILMmatrix(y,x)＝ILMposition_y(x)，y∈[1,N]，N为B-scan方向的序列数；利用帧间相关性，在B-scan方向上使用一维中值滤波平滑约束ILM位置矩阵，去除误差点，模板大小为5；

步骤2-3、由于RPE位于SD-OCT视网膜图像的超反射层，因此搜索每个A-scan中像素最亮值来定位RPE初始位置，具体方法如下：

由于采用灰度最大值定位的RPE初始位置中含有大量RNFL层的像素点，因此根据ILM边界与RPE层的距离限制RPE搜索区域以分割RPE层；定义RPE位置矩阵RPEmatrix(x,y)，x∈[1,M]，y∈[1,N]，用I_y(z,x)表示B-scan方向上每一帧图像的灰度，搜索图像I_y(z,x)中每一列的亮度最大值{z|maxI_y(z,x),x∈M}，构造RPE的初始位置矩阵：

RPEmatrix₀(x,y)＝{z|maxI_y(z,x),x∈[1,M]}，y∈[1,N]

采用模板为7×7的开运算和闭运算对RPE初始位置矩阵RPEmatrix₀滤波，使用OTSU阈值分割滤波后的RPE位置矩阵，得到二值图像，二值图像中为零的像素点表示定位错误的RPE边界点；然后采用8邻域填充灰度值为零的像素点，对填充后的RPE位置矩阵进行中值滤波以消除误点；

之后根据RPE厚度构造第二次搜索的范围为40pixels，即[RPEmatrix₀(x,y)-10,RPEmatrix₀(x,y)+30]，然后采用上述方法重新搜索RPE，记为RPEmatrix₁，所述RPE厚度是10±1pixels，即20±2μm；此时得到的RPE位置位于最亮层，利用RPEmatrix₁+5估计RPE-Choroid边界，即RPE层的下边界；构造OS-RPE搜索范围为5pixels，即[RPEmatrix₀(x,y)-3,RPEmatrix₀(x,y)-2]，然后采用从下向上搜索第一个亮点，若不存在该亮点，则用[RPEmatrix₀(x,y)-5]估计该位置，从而定位出OS-RPE边界；

步骤2-4、定位出OS-RPE边界后，采用从下向上的策略检测ONL-IS边界；限制OS-RPE边界面的搜索区域为40pixels，在[OS/RPE-30，OS/RPE-10]内搜索边缘图像中第一个峰值点，从而定位出ONL-IS边界；利用帧间相关性，在B-scan方向上使用模板为5×1的一维均值滤波平滑ONL-IS位置矩阵，校正ONL-IS边界深度误差点。

步骤3、根据IS-OS层上边界和RPE层下边界确定SD-OCT投影图像；所用公式为：

$SVP(x,y)=\frac{\underset{i=ONL/IS(x,y)}{\overset{RPE/Choroid(x,y)}{\Σ}}\mathrm{Im}g(x,i,y)}{RPE/Choroid(x,y)-ONL/IS(x,y)}$

式中x表示A-scan方向的宽度，y表示B-scan方向的序列数，Img为三维SD-OCT图像数据，RPE/Choroid表示RPE层下边界位置矩阵，ONL/IS表示IS-OS层上边界位置矩阵。

步骤4、应用改进的水平集方法分割SD-OCT投影图像得到血管图像；具体为：

步骤4-1、采用模板大小为3的均值滤波对SD-OCT投影图像进行噪声预处理；

步骤4-2、用改进的水平集方法分割投影图像，所用公式为：

F(φ,f₁,f₂)＝ε(φ,f₁,f₁)+μΡ(φ)

其中

$\mathrm{\ϵ}(\mathrm{\φ},f_1,f_2)=\underset{i=1}{\overset{2}{\Σ}}{\mathrm{\λ}}_i\∫(\∫K_{\mathrm{\σ}}(x-y\left){\left|I\left(y\right)-f_i\left(x\right)\right|}^2{M}_i\right(\mathrm{\φ}\left(y\right)\left)dy\right)dx+\mathrm{\ν}\∫|\▿H\left(\mathrm{\φ}\right(x\left)\right)|dx$

$M_1=H\left(\mathrm{\φ}\right),M_2=1-H\left(\mathrm{\φ}\right),H\left(x\right)=\frac{1}{2}\[1+\frac{2}{\mathrm{\π}}arctan\left(\frac{x}{\mathrm{\σ}}\right)\]$

$K_{\mathrm{\σ}}\left(u\right)=\frac{1}{{\left(2\mathrm{\π}\right)}^{n/2}{\mathrm{\σ}}^n}{e}^{-{\left|u\right|}^2/2{\mathrm{\σ}}^2}$

其中φ表示水平集函数，I(x)表示输入图像，σ表示标准差，该方法参数设置为迭代次数为200，λ₁＝λ₂＝1.0，ν＝0.005*255*255，μ＝1，σ＝5.0。

步骤5、提取彩色眼底图像的绿色通道，并对绿色通道取反；具体为：

步骤5-1、提取彩色眼底图像的绿色通道：Img(x,y)＝im(x,y,2，)其中im彩色图像，Img为灰度图像；

步骤5-2、对灰度图像取反，Img＝255-Img。

步骤6、采用形态学多尺度top-hat方法提取基于步骤5的图像的血管图像；具体步骤为：

步骤6-1、对灰度图像进行形态学闭操作，公式为：

$A\·S=\[A\⊕S\]\mathrm{\Θ}S$

其中·表示闭运算符,Θ表示腐蚀运算符，A为输入图像，S表示半径为2个像素的圆盘形结构元；

步骤6-2、对步骤6-1结果进行顶帽变换，每次顶帽变换使用一组不同尺度的圆形卷积核，圆形卷积核的半径从1个像素变化到8个像素，于是得到8个结果，公式为：

tophat(Img)＝Img-min(open(closeImg,s),Img)

式中closeimg是步骤6-1的闭运算结果，s为卷积核半径，从1变化到8；

步骤6-3、对基于步骤6-2的连续图像对进行取均值，得到四个结果，并根据不同的权重取四组图像中的最大值，具体公式为：

mimg＝max(mimg,5*(8-2*i+1)avgimg_i)

式中avgimg_i表示第i组均值图像；

步骤6-4、阈值分割基于步骤6-3的结果，其阈值为整幅的均值；

步骤6-5、对基于步骤6-4的分割结果进行中值滤波，其模板大小为7×7；

步骤6-6、用形态学方法对基于步骤6-5的结果进行腐蚀运算，卷积核半径为1个像素；

步骤6-7、根据连通区域特性去除孤立的像素点，其阈值为200。

步骤7、利用快速并行细化算法骨骼化血管图像得到细化图像；具体是假设已知目标点标记为1，背景点标记为0，边界点定义为本身标记为1而其8-连通邻域中至少有1个标记点为0的点，处理方法如下：

步骤7-1、对以边界点为中心的8-邻域进行处理，记中心点为p1，其邻域的8个点顺时针绕中心点分别记为p2,p3,……p9，其中p2在p1的正上方；首先标记同时满足下列条件的边界点：

2≤N(p1)≤6；S(p1)＝1；p2×p4×p6＝0；p4×p6×p8＝0；

其中N(p1)是p1的非零邻点的个数，S(p1)是以p2,p3,……p9为序时这些点的值从0到1变化的次数，当对所有边界点都检验完毕后，除去所有标记的点；

步骤7-2、对以边界点为中心的8-邻域进行处理，记中心点为p1，其邻域的8个点顺时针绕中心点分别记为p2,p3,……p9，其中p2在p1的正上方；首先标记同时满足下列条件的边界点：

2≤N(p1)≤6；S(p)＝1；p2×p4×p8＝0；p2×p6×p8＝0；

其中N(p1)是p1的非零邻点的个数，S(p1)是以p2,p3,……p9为序时这些点的值从0到1变化的次数，当对所有边界点都检验完毕后，除去所有标记的点；

步骤7-3、以上两步操作构成一次迭代，直至没有点再满足标记条件，这时剩下的点组成区域的骨架。

步骤8、对SD-OCT血管图像按照4×4进行均匀分块，然后对每一局部图像块作为模板遍历彩色眼底血管图像并确定血管重叠像素数，最后根据重叠率提取控制点；具体方法为：

步骤8-1、将每一个局部图像块作为模板遍历彩色眼底血管图像，然后统计每个像素点的血管重叠像素数，具体公式为：overlapCount＝sum(CImg*M)，式中*是卷积算子，CImg是彩色眼底血管图像，M是局部图像块，图像中的值由0和1组成，sum是求和运算符；之后求取最大重叠像素数的像素位置；

步骤8-2、根据步骤8-1的像素位置为中心，提取大小为128×128的彩色眼底血管图像块，之后计算与模板的重叠率，具体公式为：

$overlaprat(m,n)=\frac{overlapCount}{Count}$

式中overlapCount表示两个图像块血管重叠的像素数，Count表示两个图像块血管像素总数；m,n分别表示4×4矩阵的横坐标和纵坐标；

步骤8-3、将重叠率与给定阈值进行比较，如果小于给定阈值，则不记录控制点，否则记录控制点，其阈值为0.33，具体为：首先提取血管重叠的部分，并采用步骤7的方法进行骨骼化，然后根据区域连通性对骨骼化图像进行去除孤立点，其阈值为55，最后提取SD-OCT投影图像和彩色眼底图像每个重叠区域的3个不相邻的控制点，该三个控制点分别为正数第10个，中间位置，倒数第10个；

步骤8-4、判断SD-OCT图像四个局部块是否存在控制点，如果不存在控制点，则根据四个局部块与提取的第一个控制点的相对位置计算出SD-OCT投影图像和彩色眼底图像的控制点，所述SD-OCT图像四个局部块是指位于SD-OCT投影图像左上角、右上角、右下角、左下角的图像块。

步骤9、根据提取的控制点计算仿射变换系数，然后利用仿射变换系数计算SD-OCT投影图像在彩色眼底图像的位置，得到粗配准结果；

步骤10、将基于步骤9的粗配准结果作为控制点搜索区域，重新提取四个局部图像块的控制点，之后确定仿射变换系数，再次进行图像匹配，得到精配准结果；所述重新提取的四个局部图像块是指位于SD-OCT投影图像左上角、右上角、右下角、左下角的图像块。具体方法为：

步骤10-1、提取粗配准结果的彩色眼底图像，并对其按照2×2进行分块；

步骤10-2、对这四个局部图像块分别遍历对应步骤10-1的四个图像块，求取最大重叠像素数的位置，并根据位置提取大小为128×128的图像块；

步骤10-3、利用步骤7的方法细化重叠的图像块；

步骤10-4、根据面积去除步骤10-3细化图像的孤立点，其面积阈值为50；

步骤10-5、基于步骤10-4结果获取重叠像素的位置，并提取SD-OCT投影图像和彩色眼底图像每个重叠区域的2个控制点，该两个控制点分别是第一个控制点和中间位置的控制点。

本发明与现有技术相比，其显著优点为：(1)本发明给出了一种基于血管的眼底图像与SD-OCT投影图像配准方法，该方法采用分块的思想遍历彩色眼底图像，根据相似函数计算像素的重叠率，之后利用重叠率提取图像匹配点，该方法能够准确地提取出图像匹配点；(2)本发明利用逐步迭代的思想优化图像控制点，在粗匹配的基础上进一步提取不存在血管区域的控制点，这样能够保证配准的精度。

下面结合附图对本发明作进一步详细描述。

附图说明

图1是本发明基于血管的眼底图像与SD-OCT投影图像配准方法的流程图。

图2是基于双梯度信息和空间约束信息的SD-OCT视网膜图像IS-OS层和RPE层分割的流程图。

图3是SD-OCT投影图像的血管分割流程图。

图4是彩色眼底图像的血管分割流程图。

图5是提取图像的控制点的流程图。

图6是重提取四个局部图像块的控制点的流程图。

图7是眼底图像。

图8是三维SD-OCT视网膜图像。

图9是SD-OCT投影图像。

图10是SD-OCT投影图像的血管分割示意图。

图11是彩色眼底图像的血管分割示意图。

图12是配准结果示意图。

具体实施方式

结合图1，本发明的基于血管的眼底图像与SD-OCT投影图像配准方法包括以下步骤：

步骤1、采集彩色眼底图像与三维SD-OCT视网膜图像；

步骤2、用基于双梯度信息和空间约束信息的视网膜组织层分割算法对SD-OCT视网膜图像进行IS-OS层和RPE层分割；具体为：

步骤2-1、采用基于非局部均值思想的三维块匹配去噪算法对SD-OCT图像进行降噪处理，然后用基于边缘流的方法提取视网膜图像的边缘，所述边缘的像素值为255，其它为黑色背景；

步骤2-2、对于步骤2-1结果的每个A-scan，从上向下搜索第一个值为255的像素点，提取ILM边界；所述A-scan表示图像的列，ILM为视网膜图内界膜；之后在垂直方向上使用确定图像梯度d2l来修正ILM边界，I为视网膜图像，具体方法如下：

定义ILM边界位置向量ILMposition(x)，x∈[1,M]，M为A-scan方向的宽度；之后判断相邻两列间的位置差：Diff＝|ILMposition(x)-ILMposition(x-1)|，如果位置差大于给定的阈值T，则用当前点的局部邻域内最大梯度位置修正该点的边界值，即

$ILMposition\left(x\right)=\left\{\begin{array}{cc}\underset{z\∈S}{max}d2l(z,x)& \begin{array}{cc}if& Diff>T\end{array}\\ ILMposition\left(x\right)& Others\end{array}\right.$

其中S是以ILMposition(x)为中心在x列窗口大小为l区域内像素点集合，d2l(z,x)为由暗到亮梯度图像，z表示图像的高度，x表示图像的宽度，阈值T＝20；

将ILM边界位置向量构成ILM位置矩阵ILMmatrix(y,x)＝ILMposition_y(x)，y∈[1,N]，N为B-scan方向的序列数；利用帧间相关性，在B-scan方向上使用一维中值滤波平滑约束ILM位置矩阵，去除误差点，模板大小为5；

步骤2-3、由于RPE位于SD-OCT视网膜图像的超反射层，因此利用图像每个A-scan中像素最亮值定位RPE初始位置，具体方法如下：

由于采用灰度最大值定位的RPE初始位置中含有大量的RNFL层像素点，因此根据ILM边界与RPE层的距离限制RPE搜索区域以分割RPE层；定义RPE位置矩阵RPEmatrix(x,y)，x∈[1,M]，y∈[1,N]，用I_y(z,x)表示B-scan方向上每一帧图像的灰度，搜索图像I_y(z,x)中每一列的亮度最大值{z|maxI_y(z,x),x∈M}，构造RPE的初始位置矩阵：

RPEmatrix₀(x,y)＝{z|maxI_y(z,x),x∈[1,M]},y∈[1,N]

采用模板为7×7的开运算和闭运算对RPE初始位置矩阵RPEmatrix0滤波，使用OTSU阈值分割滤波后的RPE位置矩阵，得到二值图像，二值图像中为零的像素点表示定位错误的RPE边界点；然后采用8邻域填充灰度值为零的像素点，对填充后的RPE位置矩阵进行中值滤波以消除误点；

之后根据RPE厚度构造第二次搜索的范围为40pixels，即[RPEmatrix₀(x,y)-10,RPEmatrix₀(x,y)+30]，然后采用上述方法重新搜索RPE，记为RPEmatrix₁，所述RPE厚度是10±1pixels，20±2μm；此时得到的RPE位置位于最亮层，利用RPEmatrix₁+5估计RPE-Choroid边界，即RPE层的下边界；构造OS-RPE搜索范围为5pixels，即[RPEmatrix₀(x,y)-3,RPEmatrix₀(x,y)-2]，然后采用从下向上搜索第一个亮点，若不存在该亮点，则用[RPEmatrix₀(x,y)-5]估计该位置，从而估计出OS-RPE边界；

步骤2-4、检测出OS-RPE边界后，采用从下向上的策略检测ONL-IS边界；限制OS-RPE边界面的搜索区域为40pixels，在[OS/RPE-30，OS/RPE-10]内搜索边缘图像中第一个峰值点，从而定位出ONL-IS边界；利用帧间相关性，在B-scan方向上使用模板为5×1的一维均值滤波平滑ONL-IS位置矩阵，校正ONL-IS边界深度误差点。

步骤3、根据IS-OS层上边界和RPE层下边界确定SD-OCT投影图像，所用公式为：

$SVP(x,y)=\frac{\underset{i=ONL/IS(x,y)}{\overset{RPE/Choroid(x,y)}{\Σ}}\mathrm{Im}g(x,i,y)}{RPE/Choroid(x,y)-ONL/IS(x,y)}$

式中x表示A-scan方向的宽度，y表示B-scan方向的序列数，Img为三维SD-OCT图像数据，RPE/Choroid表示RPE层下边界位置矩阵，ONL/IS表示IS-OS层上边界位置矩阵。

步骤4、应用改进的水平集方法分割SD-OCT投影图像得到血管图像，具体为：

步骤4-1、采用模板大小为3的均值滤波对SD-OCT投影图像进行噪声预处理；

步骤4-2、用改进的水平集方法分割投影图像，所用公式为：

F(φ,f₁,f₂)＝ε(φ,f₁,f₁)+μΡ(φ)

其中

$\mathrm{\ϵ}(\mathrm{\φ},f_1,f_2)=\underset{i=1}{\overset{2}{\Σ}}{\mathrm{\λ}}_i\∫(\∫K_{\mathrm{\σ}}(x-y\left){\left|I\left(y\right)-f_i\left(x\right)\right|}^2{M}_i\right(\mathrm{\φ}\left(y\right)\left)dy\right)dx+\mathrm{\ν}\∫|\▿H\left(\mathrm{\φ}\right(x\left)\right)|dx$

$M_1=H\left(\mathrm{\φ}\right),M_2=1-H\left(\mathrm{\φ}\right),H\left(x\right)=\frac{1}{2}\[1+\frac{2}{\mathrm{\π}}arctan\left(\frac{x}{\mathrm{\σ}}\right)\]$

$K_{\mathrm{\σ}}\left(u\right)=\frac{1}{{\left(2\mathrm{\π}\right)}^{n/2}{\mathrm{\σ}}^n}{e}^{-{\left|u\right|}^2/2{\mathrm{\σ}}^2}$

其中φ表示水平集函数，I(x)表示输入图像，σ表示标准差，该方法参数设置为迭代次数为200，λ₁＝λ₂＝1.0，ν＝0.005*255*255，μ＝1，σ＝5.0。

步骤5、提取彩色眼底图像的绿色通道，并对绿色通道取反，具体为：

步骤5-1、提取彩色眼底图像的绿色通道：Img(x,y)＝im(x,y,2，)其中im彩色图像，Img为灰度图像；

步骤5-2、对灰度图像取反，Img＝255-Img。

步骤6、采用形态学多尺度top-hat方法提取基于步骤5结果的血管图像，具体步骤为：

步骤6-1、对灰度图像进行形态学闭操作，公式为：

$A\·S=\[A\⊕S\]\mathrm{\Θ}S$

其中·表示闭运算符,Θ表示腐蚀运算符，A为输入图像，S表示半径为2个像素的圆盘形结构元；

步骤6-2、对步骤6-1结果进行顶帽变换，每次顶帽变换使用一组不同尺度的圆形卷积核，圆形卷积核的半径从1个像素变化到8个像素，于是得到8个结果，公式为：

tophat(Img)＝Img-min(open(closeImg,s),Img)

式中closeimg是步骤6-1的闭运算结果，s为卷积核半径，从1变化到8；

步骤6-3、对基于步骤6-2的连续图像对进行取均值，得到四个结果，并根据不同的权重取四组图像中的最大值，具体公式为：

mimg＝max(mimg,5*(8-2*i+1)avgimg_i)

式中avgimg_i表示第i组均值图像；

步骤6-4、阈值分割基于步骤6-3的结果，其阈值为整幅的均值；

步骤6-5、对基于步骤6-4的分割结果进行中值滤波，其模板大小为7×7；

步骤6-6、用形态学方法对基于步骤6-5的结果进行腐蚀运算，卷积核半径为1个像素；

步骤6-7、根据连通区域特性去除孤立的像素点，其阈值为200。

步骤7、利用快速并行细化算法骨骼化血管图像得到细化图像，具体是假设已知目标点标记为1，背景点标记为0，边界点定义为本身标记为1而其8-连通邻域中至少有1个标记点为0的点，处理方法如下：

步骤7-1、对以边界点为中心的8-邻域进行处理，记中心点为p1，其邻域的8个点顺时针绕中心点分别记为p2,p3,……p9，其中p2在p1的正上方；首先标记同时满足下列条件的边界点：

2≤N(p1)≤6；S(p1)＝1；p2×p4×p6＝0；p4×p6×p8＝0；

其中N(p1)是p1的非零邻点的个数，S(p1)是以p2,p3,……p9为序时这些点的值从0到1变化的次数，当对所有边界点都检验完毕后，除去所有标记的点；

步骤7-2、对以边界点为中心的8-邻域进行处理，记中心点为p1，其邻域的8个点顺时针绕中心点分别记为p2,p3,……p9，其中p2在p1的正上方；首先标记同时满足下列条件的边界点：

2≤N(p1)≤6；S(p)＝1；p2×p4×p8＝0；p2×p6×p8＝0；

其中N(p1)是p1的非零邻点的个数，S(p1)是以p2,p3,……p9为序时这些点的值从0到1变化的次数，当对所有边界点都检验完毕后，除去所有标记的点；

步骤7-3、以上两步操作构成一次迭代，直至没有点再满足标记条件，这时剩下的点组成区域的骨架。

步骤8、对SD-OCT血管图像按照4×4进行均匀分块，每个局部图像块的大小为128×128，然后以每个局部图像块为模板遍历彩色眼底血管图像并计算血管重叠率，根据重叠率提取控制点，具体方法为：

步骤8-1、以每个局部图像块为模板遍历彩色眼底血管图像，然后统计每个像素点的血管重叠像素数，具体公式为：overlapCount＝sum(CImg*M)，式中*是卷积算子，CImg是彩色眼底血管图像，M是局部图像块，图像中的值由0和1组成，sum是求和运算符；之后求取最大重叠像素数的像素位置；

步骤8-2、根据步骤8-1的像素位置为中心，提取大小为128×128的彩色眼底血管图像块，之后计算与模板的重叠率，具体公式为：

$overlaprat(m,n)=\frac{overlapCount}{Count}$

式中overlapCount表示两个图像块血管重叠的像素数，Count表示两个图像块血管像素总数；m,n分别表示4×4矩阵的横坐标和纵坐标；

步骤8-3、将重叠率与给定阈值进行比较，如果小于给定阈值，则不记录控制点，否则记录控制点，其阈值为0.33，具体为：首先提取血管重叠的部分，并采用步骤7的方法进行骨骼化，然后根据区域连通性对骨骼化图像进行去除孤立点，其阈值为55，最后提取SD-OCT投影图像和彩色眼底图像每个重叠区域的3个不相邻的控制点，该三个控制点分别为正数第10个，中间位置，倒数第10个；

步骤8-4、判断SD-OCT图像四个局部块是否存在控制点，如果不存在控制点，则根据第一控制点与四个局部块的相对位置计算出SD-OCT投影图像和彩色眼底图像的控制点，所述SD-OCT图像四个局部块是指位于SD-OCT投影图像左上角、右上角、右下角、左下角的图像块。

步骤9、根据提取的控制点计算仿射变换系数，然后利用仿射变换系数计算SD-OCT投影图像在彩色眼底图像的位置，得到粗图像配准结果；

步骤10、将基于步骤9的粗配准结果限制搜索范围，重新提取四个局部图像块的控制点，之后确定仿射变换系数，再次进行图像匹配，具体方法为：

步骤10-1、提取粗配准结果的彩色眼底图像，并对其按照2×2进行分块；

步骤10-2、对这四个局部图像块分别遍历对应步骤10-1的四个图像块，求取最大重叠像素数的位置，并根据位置提取大小为128×128的图像块；

步骤10-3、利用步骤7的方法细化重叠的图像块；

步骤10-4、根据面积去除步骤10-3细化图像的孤立点，其面积阈值为50；

步骤10-5、基于步骤10-4结果获取重叠像素的位置，并提取SD-OCT投影图像和彩色眼底图像每个重叠区域的2个控制点，该两个控制点分别是第一个控制点和中间位置的控制点。

下面结合实施例对本发明做进一步详细的说明：

本发明以彩色眼底图像和三维SD-OCT视网膜图像作为输入，采用图像处理手段对不同模态下的图像进行自动配准。

本实施例的流程如图1所示，通过OCT成像设备采集到的彩色眼底图像大小为2392×2048，及三维SD-OCT视网膜图像大小为1024×512×128，图7和图8分别给出了一幅彩色眼底图像和一组三维SD-OCT视网膜图像。首先对三维SD-OCT视网膜图像进行IS-OS层和RPE层分割，然后根据IS-OS层上边界和RPE层下边界，对两边界之间的灰度累加取平均，从而生成图9所示的SD-OCT投影图像。由于生成的SD-OCT投影图像含有噪声，首先采用均值滤波对投影图像进行去噪，然后应用改进的水平集方法对SD-OCT投影图像进行血管分割，于是得到SD-OCT投影的血管图像，如图10所示。为了便于彩色眼底图像血管的提取，首先取彩色眼底图像的绿色通道，并对绿色通道取反，然后用形态学方法对其进行闭操作，并对结果采用8尺度的顶帽滤波，得到8个结果，将连续图像对进行取均值，得到四个结果，并根据不同的权重取四组图像的最大值，得到血管图像，最后对血管图像进行后续处理，如中值滤波、腐蚀操作、根据面积去除孤立的像素点，得到最终的血管图像(图11)。

将图10均匀分成4×4的图像块，每个图像块的大小为128×128。将每个图像块作为模板去遍历彩色眼底的血管图像(图11)，并计算以每个像素为中心重叠血管的像素数，找出最大的重叠数，利用公式计算出重叠率，将重叠率与给定阈值进行比较(阈值为0.33)，如果大于给定阈值，则取模板大小的血管重叠部分，对其进行细化，为了提取更加稳定的控制点，根据面积对细化图像进行去除孤立点，面积阈值为55，然后从重叠的血管中选择三个不相邻的控制点。

由于四个局部图像块(左上角、左下角、右上角、右下角)的控制点对配准结果影响较大，所以需要判断这四个局部图像块是否存在控制点，如果不存在控制点，根据第一个控制点与四个局部图像块的相对位置计算出SD-OCT投影图像和彩色眼底图像的控制点。

根据提取的控制点计算仿射变换系数，然后利用变换系数计算SD-OCT投影图像在彩色眼底图像的位置，将投影图像叠加到彩色眼底图像上，得到粗略的图像配准结果。由于四个局部图像块控制点造成仿射变换系数误差大，导致配准精度不高，所以需要根据粗配准结果进一步提取四个局部图像块的控制点，以达到精确配准。首先提取粗匹配部分的彩色眼底图像，并对其按照2×2进行分块，对这四个局部图像块分别遍历对应粗匹配的局部图像块，求取最大重叠血管像素数的位置，并以该位置为中心提取大小为128×128的图像块，对其进行细化，根据面积去除孤立的像素点，然后根据重叠部分的血管选择2个不相邻的控制点，最后用提取的控制点替代上述估计的控制点，再次利用仿射变换进行配准，得到精确配准，配准结果如图12所示。从图12可知：本发明得到的配准结果精度高，能够满足临床医生的判断依据；而且当两幅图像融合在一起后，眼底图像的信息增加了，能更好的分析图像的信息特征。

Claims (9)

1.一种基于血管的眼底图像与SD-OCT投影图像配准方法，其特征在于，包括以下步骤：

步骤1、采集彩色眼底图像与三维SD-OCT视网膜图像；

步骤2、用基于双梯度信息和空间约束信息的视网膜组织层分割算法对SD-OCT视网膜图像进行ONL.IS边界和RPE下边界分割；

步骤3、根据ONL.IS边界和RPE层下边界确定SD-OCT投影图像；

步骤4、应用改进的水平集方法分割SD-OCT投影图像得到血管图像；

步骤5、提取彩色眼底图像的绿色通道，并对绿色通道取反；

步骤6、采用形态学多尺度top-hat方法提取基于步骤5图像的血管图像；

步骤7、利用快速并行细化算法骨骼化血管图像得到细化图像；

步骤8、对SD-OCT血管图像按照4×4进行均匀分块，然后将每一局部图像块作为模板遍历彩色眼底血管图像并确定血管重叠像素数，最后根据重叠率提取控制点；

步骤9、根据提取的控制点确定仿射变换系数，然后利用仿射变换系数确定SD-OCT投影图像在彩色眼底图像的位置，得到粗配准结果；

步骤10、将基于步骤9的粗配准结果作为控制点搜索区域，重新提取四个局部图像块的控制点，之后确定仿射变换系数，再次进行图像匹配，得到精配准结果；所述重新提取的四个局部图像块是指位于SD-OCT投影图像左上角、右上角、右下角、左下角的图像块。

2.根据权利要求1所述的基于血管的眼底图像与SD-OCT投影图像配准方法，其特征在于，步骤2用基于双梯度信息和空间约束信息的视网膜组织层分割算法对SD-OCT视网膜图像进行ONL.IS边界和RPE下边界分割，具体为：

步骤2-1、采用基于非局部均值思想的三维块匹配去噪算法对SD-OCT图像进行降噪处理，然后用基于边缘流的方法提取视网膜图像的边缘，所述边缘的像素值为255，其它为黑色背景；

步骤2-2、对于步骤2-1结果的每个A-scan，从上向下搜索第一个值为255的像素点，提取ILM边界；所述A-scan表示图像的列，ILM为视网膜图像内界膜；之后在垂直方向上使用确定图像梯度d2l来修正ILM边界，I为视网膜图像，具体方法如下：

定义ILM边界位置向量ILMposition(x)，x∈[1,M]，M为A-scan方向的宽度；之后判断相邻两列间的位置差：Diff＝|ILMposition(x)-ILMposition(x-1)|，如果位置差大于给定的阈值T，则用当前点的局部邻域内最大梯度位置修正该点的边界值，即

$ILMposition\left(x\right)=\left\{\begin{array}{cc}\underset{z\∈S}{\max }d2l(z,x)& ifDiff>T\\ ILMposition\left(x\right)& Others\end{array}\right.$

其中S是以ILMposition(x)为中心在x列窗口大小为l区域内的像素点集合，d2l(z,x)为由暗到亮梯度图像，z表示图像的高度，x表示图像的宽度，阈值T＝20；

将ILM边界位置向量构成ILM位置矩阵ILMmatrix(y,x)＝ILMposition_y(x)，y∈[1,N]，N为B-scan方向的序列数；利用帧间相关性，在B-scan方向上使用一维中值滤波平滑约束ILM位置矩阵，去除误差点，模板大小为5；

步骤2-3、由于RPE位于SD-OCT视网膜图像的超反射层，因此搜索每个A-scan中像素最亮值来定位RPE初始位置，具体方法如下：

由于采用灰度最大值定位的RPE初始位置中含有大量RNFL层的像素点，因此根据ILM边界与RPE层的距离限制RPE搜索区域以分割RPE层；定义RPE位置矩阵RPEmatrix(x,y)，x∈[1,M]，y∈[1,N]，用I_y(z,x)表示B-scan方向上每一帧图像的灰度，搜索图像I_y(z,x)中每一列的亮度最大值{z|max I_y(z,x),x∈M}，构造RPE的初始位置矩阵：

RPEmatrix₀(x,y)＝{z|max I_y(z,x),x∈[1,M]}，y∈[1,N]

采用模板为7×7的开运算和闭运算对RPE初始位置矩阵RPEmatrix₀滤波，使用OTSU阈值分割滤波后的RPE位置矩阵，得到二值图像，二值图像中为零的像素点表示定位错误的RPE边界点；然后采用8邻域填充灰度值为零的像素点，对填充后的RPE位置矩阵进行中值滤波以消除误点；

之后根据RPE厚度构造第二次搜索的范围为40pixels，即[RPEmatrix₀(x,y)-10,RPEmatrix₀(x,y)+30]，然后采用上述方法重新搜索RPE，记为RPEmatrix₁，所述RPE厚度是10±1pixels，即20±2μm；此时得到的RPE位置位于最亮层，利用RPEmatrix₁+5估计RPE-Choroid边界，即RPE层的下边界；构造OS.RPE搜索范围为5pixels，即[RPEmatrix₀(x,y)-3,RPEmatrix₀(x,y)-2]，然后采用从下向上搜索第一个亮点，若不存在该亮点，则用[RPEmatrix₀(x,y)-5]估计该位置，从而定位出OS.RPE边界；

步骤2-4、定位出OS.RPE边界后，采用从下向上的策略检测ONL-IS边界；限制OS.RPE边界面的搜索区域为40pixels，在[OS.RPE-30，OS.RPE-10]内搜索边缘图像中第一个峰值点，从而定位出ONL-IS边界；利用帧间相关性，在B-scan方向上使用模板为5×1的一维均值滤波平滑ONL-IS位置矩阵，校正ONL-IS边界深度误差点。

3.根据权利要求1所述的基于血管的眼底图像与SD-OCT投影图像配准方法，其特征在于，步骤3中根据ONL.IS边界和RPE层下边界确定SD-OCT投影图像，所用公式为：

$SVP(x,y)=\frac{\underset{i=ONL/IS(x,y)}{\overset{RPE/Choroid(x,y)}{\Σ}}\mathrm{Im}g(x,i,y)}{RPE/Choroid(x,y)-ONL/IS(x,y)}$

式中x表示A-scan方向的宽度，y表示B-scan方向的序列数，Img为三维SD-OCT图像数据，RPE/Choroid表示RPE层下边界位置矩阵，ONL/IS表示ONL.IS边界位置矩阵。

4.根据权利要求1所述的基于血管的眼底图像与SD-OCT投影图像配准方法，其特征在于，步骤4应用改进的水平集方法分割SD-OCT投影图像得到血管图像，具体为：

步骤4-1、采用模板大小为3的均值滤波对SD-OCT投影图像进行噪声预处理；

步骤4-2、用改进的水平集方法分割投影图像，所用公式为：

F(φ,f₁,f₂)＝ε(φ,f₁,f₁)+μP(φ)

其中

$\mathrm{\ϵ}(\mathrm{\φ},f_1,f_2)=\underset{i=1}{\overset{2}{\Σ}}{\mathrm{\λ}}_i\∫(\∫K_{\mathrm{\σ}}(x-y\left)\right|I\left(y\right)-f\left(x\right){|}^2{M}_i\left(\mathrm{\φ}\left(y\right)\right)dy)dx+v\∫|\▿H\left(\mathrm{\φ}\right(x\left)\right)|dx$

$M_1=H\left(\mathrm{\φ}\right),M_2=1-H\left(\mathrm{\φ}\right),H\left(x\right)=\frac{1}{2}\[1+\frac{2}{\mathrm{\π}}arctan\left(\frac{x}{\mathrm{\σ}}\right)\]$

$K_{\mathrm{\σ}}\left(u\right)=\frac{1}{{\left(2\mathrm{\π}\right)}^{n/2}{\mathrm{\σ}}^n}{e}^{-|u{|}^2/2{\mathrm{\σ}}^2}$

其中φ表示水平集函数，I(x)表示输入图像，σ表示标准差，该方法参数设置为迭代次数为200，λ₁＝λ₂＝1.0，v＝0.005*255*255，μ＝1，σ＝5.0。

5.根据权利要求1所述的基于血管的眼底图像与SD-OCT投影图像配准方法，其特征在于，步骤5提取彩色眼底图像的绿色通道，并对绿色通道取反，具体为：

步骤5-1、提取彩色眼底图像的绿色通道：Img(x,y)＝im(x,y,2，其中im彩色图像，Img为灰度图像；

步骤5-2、对灰度图像取反，Img＝255-Img。

6.根据权利要求1所述的基于血管的眼底图像与SD-OCT投影图像配准方法，其特征在于，步骤6采用形态学多尺度top-hat方法提取基于步骤5结果的血管图像，具体步骤为：

步骤6-1、对灰度图像进行形态学闭操作，公式为：

$A\·S=\[A\⊕S\]\mathrm{\Θ}S$

其中·表示闭运算符,Θ表示腐蚀运算符，A为输入图像，S表示半径为2个像素的圆盘形结构元；

步骤6-2、对步骤6-1结果进行顶帽变换，每次顶帽变换使用一组不同尺度的圆形卷积核，圆形卷积核的半径从1个像素变化到8个像素，于是得到8个结果，公式为：

tophat(Img)＝Img-min(open(closeImg,s),Img)

式中closeimg是步骤6-1的闭运算结果，s为卷积核半径，从1变化到8；

步骤6-3、对基于步骤6-2的连续图像对进行取均值，得到四个结果，并根据不同的权重取四组图像中的最大值，具体公式为：

mimg＝max(mimg,5*(8-2*i+1)avgimg_i)

式中avgimg_i表示第i组均值图像；

步骤6-4、阈值分割基于步骤6-3的结果，其阈值为整幅的均值；

步骤6-5、对基于步骤6-4的分割结果进行中值滤波，其模板大小为7×7；

步骤6-6、用形态学方法对基于步骤6-5的结果进行腐蚀运算，卷积核半径为1个像素；

步骤6-7、根据连通区域特性去除孤立的像素点，其阈值为200。

7.根据权利要求1所述的基于血管的眼底图像与SD-OCT投影图像配准方法，其特征在于，步骤7利用快速并行细化算法骨骼化血管图像得到细化图像，具体是假设已知目标点标记为1，背景点标记为0，边界点定义为本身标记为1而其8-连通邻域中至少有1个标记点为0的点，处理方法如下：

步骤7-1、对以边界点为中心的8-邻域进行处理，记中心点为p1，其邻域的8个点顺时针绕中心点分别记为p2,p3,……p9，其中p2在p1的正上方；首先标记同时满足下列条件的边界点：

2≤N(p1)≤6；S(p1)＝1；p2×p4×p6＝0；p4×p6×p8＝0；

其中N(p1)是p1的非零邻点的个数，S(p1)是以p2,p3,……p9为序时这些点的值从0到1变化的次数，当对所有边界点都检验完毕后，除去所有标记的点；

步骤7-2、对以边界点为中心的8-邻域进行处理，记中心点为p1，其邻域的8个点顺时针绕中心点分别记为p2,p3,……p9，其中p2在p1的正上方；首先标记同时满足下列条件的边界点：

2≤N(p1)≤6；S(p)＝1；p2×p4×p8＝0；p2×p6×p8＝0；

其中N(p1)是p1的非零邻点的个数，S(p1)是以p2,p3,……p9为序时这些点的值从0到1变化的次数，当对所有边界点都检验完毕后，除去所有标记的点；

步骤7-3、以上两步操作构成一次迭代，直至没有点再满足标记条件，这时剩下的点组成区域的骨架。

8.根据权利要求1所述的基于血管的眼底图像与SD-OCT投影图像配准方法，其特征在于，步骤8对SD-OCT血管图像按照4×4进行均匀分块，然后将每一个局部图像块作为模板遍历彩色眼底血管图像并确定血管重叠像素数，根据重叠率提取控制点，具体方法为：

步骤8-1、将每一个局部图像块作为模板遍历彩色眼底血管图像，然后统计每个像素点的血管重叠像素数，具体公式为：overlapCount＝sum(CImg*M)，式中*是卷积算子，CImg是彩色眼底血管图像，M是局部图像块，图像中的值由0和1组成，sum是求和运算符；之后求取最大重叠像素数的像素位置；

步骤8-2、根据步骤8-1的像素位置为中心，提取大小为128×128的彩色眼底血管图像块，之后计算与模板的重叠率，具体公式为：

$overlaprat(m,n)=\frac{overlapCount}{Count}$

式中overlapCount表示两个图像块血管重叠的像素数，Count表示两个图像块血管像素总数；m,n分别表示4×4矩阵的横坐标和纵坐标；

步骤8-3、将重叠率与给定阈值进行比较，如果小于给定阈值，则不记录控制点，否则记录控制点，其阈值为0.33，具体为：首先提取血管重叠的部分，并采用步骤7的方法进行骨骼化，然后根据区域连通性对骨骼化图像进行去除孤立点，其阈值为55，最后提取SD-OCT投影图像和彩色眼底图像每个重叠区域的3个不相邻的控制点，该三个控制点分别为正数第10个，中间位置，倒数第10个；

步骤8-4、判断SD-OCT图像四个局部块是否存在控制点，如果不存在控制点，则根据四个局部块与提取的第一个控制点的相对位置计算出SD-OCT投影图像和彩色眼底图像的控制点，所述SD-OCT图像四个局部块是指位于SD-OCT投影图像左上角、右上角、右下角、左下角的图像块。

9.根据权利要求1所述的基于血管的眼底图像与SD-OCT投影图像配准方法，其特征在于，步骤10将基于步骤9的粗配准结果作为控制点搜索区域，重新提取四个局部图像块的控制点，之后确定仿射变换系数，再次进行图像匹配，得到精配准结果；所述重新提取的四个局部图像块是指位于SD-OCT投影图像左上角、右上角、右下角、左下角的图像块，具体方法为：

步骤10-1、提取粗配准结果的彩色眼底图像，并对其按照2×2进行分块；

步骤10-2、对这四个局部图像块分别遍历对应步骤10-1的四个图像块，求取最大重叠像素数的位置，并根据位置提取大小为128×128的图像块；

步骤10-3、利用步骤7的方法细化重叠的图像块；

步骤10-4、根据面积去除步骤10-3细化图像的孤立点，其面积阈值为50；

步骤10-5、基于步骤10-4结果获取重叠像素的位置，并提取SD-OCT投影图像和彩色眼底图像每个重叠区域的2个控制点，该两个控制点分别是第一个控制点和中间位置的控制点。