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

摘要

本发明公开了一种基于血管的眼底图像与频域光学相干断层（SD-OCT）投影图像自动配准方法，属于图像处理技术领域。该方法首先对输入的SD-OCT图像进行组织层分割，定位光感受器内节/外节层（IS-OS）上边界和视网膜色素上皮（RPE）下边界，根据IS-OS层上边界和RPE层下边界生成SD-OCT投影图像，利用改进的水平集方法分割SD-OCT投影图像，然后采用形态学多尺度顶帽（top-hat）方法分割彩色眼底图像，最后使用迭代分块方法搜索最大相似函数的局部最大值来提取控制点，并根据控制点计算仿射变换系数实现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，则用当前点的局部邻域内最大梯度位置修正该点 的边界值，即

$\mathrm{ILMposition}\left(x\right)=\left(\begin{array}{cc}\underset{z\in S}{\max }d2l(z,x)& \left(\begin{array}{cc}\mathrm{if}& \mathrm{Diff}>T\end{array}\right)\\ \mathrm{ILMposition}\left(x\right)& \mathrm{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投影图像；所 用公式为：

$\mathrm{SVP}(x,y)=\frac{\underset{i=\mathrm{ONL}/\mathrm{IS}(x,y)}{\overset{\mathrm{RPE}/\mathrm{Choroid}(x,y)}{\Sigma }}\mathrm{Img}(x,i,y)}{\mathrm{RPE}/\mathrm{Choroid}(x,y)-\mathrm{ONL}/\mathrm{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₁)+μΡ(φ)

其中$P\left(\phi \right)=\int \frac{1}{2}{\left(\right|▿\phi \left(x\right)|-1)}^2\mathrm{dx},f_i\left(x\right)=\frac{K_{\sigma }\left(x\right)*\left[M_i\left(\phi \left(x\right)\right)I\left(x\right)\right]}{K_{\sigma }\left(x\right)*M_i\left(\phi \left(x\right)\right)}$

$ϵ(\phi ,f_1,f_2)=\underset{i=1}{\overset{2}{\Sigma }}{\lambda }_i\int (\int K_{\sigma }(x-y){|I\left(y\right)-f_i\left(x\right)|}^2{M}_i\left(\phi \left(y\right)\right)\mathrm{dy})\mathrm{dx}+v\int |▿H\left(\phi \left(x\right)\right)|\mathrm{dx}$

$M_1=H\left(\phi \right),M_2=1-H\left(\phi \right),H\left(x\right)=\frac{1}{2}[1+\frac{2}{\pi }\arctan \left(\frac{x}{\sigma }\right)]$

$K_{\sigma }\left(u\right)=\frac{1}{{\left(2\pi \right)}^{n/2}{\sigma }^n}{e}^{{-\left|u\right|}^2/{2\sigma }^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]Θ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的彩色眼 底血管图像块，之后计算与模板的重叠率，具体公式为：

$\mathrm{overlaprat}(m,n)=\frac{\mathrm{overlapCount}}{\mathrm{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，则用当前点的局部邻域内最大梯度位置修正该点 的边界值，即

$\mathrm{ILMposition}\left(x\right)=\left(\begin{array}{cc}\underset{z\in S}{\max }d2l(z,x)& \left(\begin{array}{cc}\mathrm{if}& \mathrm{Diff}>T\end{array}\right)\\ \mathrm{ILMposition}\left(x\right)& \mathrm{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投影图像，所 用公式为：

$\mathrm{SVP}(x,y)=\frac{\underset{i=\mathrm{ONL}/\mathrm{IS}(x,y)}{\overset{\mathrm{RPE}/\mathrm{Choroid}(x,y)}{\Sigma }}\mathrm{Img}(x,i,y)}{\mathrm{RPE}/\mathrm{Choroid}(x,y)-\mathrm{ONL}/\mathrm{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₁)+μΡ(φ)

其中$P\left(\phi \right)=\int \frac{1}{2}{\left(\right|▿\phi \left(x\right)|-1)}^2\mathrm{dx},f_i\left(x\right)=\frac{K_{\sigma }\left(x\right)*\left[M_i\left(\phi \left(x\right)\right)I\left(x\right)\right]}{K_{\sigma }\left(x\right)*M_i\left(\phi \left(x\right)\right)}$

$ϵ(\phi ,f_1,f_2)=\underset{i=1}{\overset{2}{\Sigma }}{\lambda }_i\int (\int K_{\sigma }(x-y){|I\left(y\right)-f_i\left(x\right)|}^2{M}_i\left(\phi \left(y\right)\right)\mathrm{dy})\mathrm{dx}+v\int |▿H\left(\phi \left(x\right)\right)|\mathrm{dx}$

$M_1=H\left(\phi \right),M_2=1-H\left(\phi \right),H\left(x\right)=\frac{1}{2}[1+\frac{2}{\pi }\arctan \left(\frac{x}{\sigma }\right)]$

$K_{\sigma }\left(u\right)=\frac{1}{{\left(2\pi \right)}^{n/2}{\sigma }^n}{e}^{{-\left|u\right|}^2/{2\sigma }^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]Θ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的彩色眼 底血管图像块，之后计算与模板的重叠率，具体公式为：

$\mathrm{overlaprat}(m,n)=\frac{\mathrm{overlapCount}}{\mathrm{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可知：本发明得到的配准结果精度高， 能够满足临床医生的判断依据；而且当两幅图像融合在一起后，眼底图像的信息 增加了，能更好的分析图像的信息特征。