一种三维非刚性光学相干断层扫描图像的配准方法与流程

本发明涉及图像配准技术领域，特别是对三维oct(光学相干断层成像)的扫描图像进行配准的三维非刚性光学相干断层扫描图像的配准方法。

背景技术：

视网膜是脑部神经组织的延伸，具有复杂的多层次组织结构。sd-oct技术能提供快速的、高分辨率的、显示视网膜内部分层的三维图像，可为临床眼科医生对疾病的诊断和治疗提供参考和帮助。视网膜oct图像的配准对临床实践具有重要意义。

目前的视网膜oct图像配准算法存在以下的缺陷：(1)大部分算法都不是完全意义上的三维配准算法，都是先对某一维或是某两维进行配准，然后再对剩下的维数进行配准，这样做会使得之前得到的配准结果在另外的维度上发生形变。(2)大部分已有的视网膜配准算法都是针对正常视网膜设计的，当视网膜组织由于病变产生较大的形变时，这些算法将失效。

脉络膜新生血管又称视网膜下新生血管，是来自脉络膜毛细血管的增殖血管，通过bruch膜的裂口而扩展，位于bruch膜与视网膜色素上皮之间、或神经视网膜与视网膜色素上皮之间、或位于视网膜色素上皮与脉络膜之间增殖形成。目前为止，还没有针对脉络膜新生血管的视网膜oct图像的三维非刚性配准方法的相关报道。

技术实现要素：

本发明要解决的技术问题为：提供一种三维非刚性光学相干断层扫描图像的配准方法，将基于灰度和基于特征的配准方法进行结合，对三维oct(光学相干断层成像)的扫描图像进行配准，提高配准结果的精确度和可靠性。

本发明采取的技术方案具体为：一种三维非刚性光学相干断层扫描图像的配准方法，包括步骤：

s1，对oct扫描图像进行预处理，包括依次对oct扫描图像进行oct去噪、图像分层和图像投影，得到视网膜血管的二维投影图像；

s2，从视网膜血管的二维投影图像中提取血管图像；

s3，对提取的血管图像进行二维特征点的提取，根据得到的二维特征点返回到三维空间得到三维特征点；

s4，对固定图像以及待配准图像进行基于特征点的仿射变换，得到粗配准的视网膜图像；

s5，利用基于灰度的非刚性配准方法对粗配准的视网膜图像进行进一步配准，得到精确配准的结果。

进一步的，本发明步骤s1中，对oct扫描图像的预处理还包括病变区域移除。对于有新生血管疾病的病人，对视网膜图像进行投影，其病变区域会造成非血管的阴影区域，这片区域的灰度值与血管的灰度值非常的接近，在检测血管时，会有误检，所以利用图像处理技术中的数学形态学的方法将这片区域移除，更有利于提取血管点

优选的，本发明步骤s1中，oct去噪为采用曲线异向扩散滤波方法对视网膜图像进行滤波，以去除散斑噪声，并保留图像的边界。优选的，曲线异向扩散滤波时，滤波参数包括时间步长、迭代次数和传导率，分别选用经验值0.0625秒、5次和3。

优选的，本发明步骤s1中，采用图搜索算法对oct去燥后的图像进行分层，分为5层，依次为神经纤维层的上边界层、神经纤维层的下边界层、内核层的边界层、内外视网膜的边界层和布鲁赫莫层。现有技术对视网膜图像进行分层时，已经可以分至10层以上，本发明仅根据算法需要从视网膜图像中分出5层。

优选的，本发明步骤s1中，图像投影时，对外层视网膜进行投影，即计算内外视网膜的边界层与布鲁赫莫层的在各位置处的平均灰度值，作为投影图像在相应位置的灰度值，即得到视网膜血管的投影图像。

例如：定义沿三维空间的y轴对视网膜图像进行分层；即分层后的各层图像沿y轴排列，同时各层图像之间视网膜图像相应位置的x轴和z轴坐标不变，仅y轴坐标变化。在计算坐标为(x,z)的投影图像位置的灰度值i(x,z)时，公式如下：

为方便描述，定义内外视网膜的边界层为第四层，布鲁赫莫层为第五层；上式中，x、y、z表示图像的坐标，y1表示第四层的y轴坐标，y2表示第五层的y轴坐标，i表示像素值(灰度值)，n表示视网膜图像的(x,z)坐标处第四层到第五层之间的像素个数。当x、z固定时，即可以通过求第四层与第五层之间的平均灰度值得到投影图像。

进一步的，步骤s1还包括对视网膜血管的投影图像进行上采样至512像素×512像素，以便于血管的提取。病变区域的移除设置于图像投影及上采样之后。

本发明步骤s3中，对提取的血管图像进行二维特征点的提取，即提取血管图像中的二维的血管点作为二维特征点。本发明对二维血管图像进行二维特征点的提取为现有技术。

优选的，本发明步骤s3中，根据得到的二维的特征点返回到三维空间，在三维空间里寻找灰度值最大的点即为血管点；同时获取二维特征点返回三维空间后与视网膜各分层的交接点，包括二维特征点分别与上边界层和布鲁赫莫层的交点，所述血管点以及二维特征点与分层的交点组成三维特征点。当已知三维坐标中的两维坐标，在第三维中找取灰度值最大的那个点，找到该点对应的坐标，即得到准确的三维坐标。本发明在寻找灰度值最大点时，采用中值法从多个灰度值较大点中选取坐标在中间的一个。在选择二维特征点与各分层交点进行应用时，实质上，二维特征点与每个分层皆存在交点，但全部选取会影响算法效率，在算法结果上却没有明显进步，故本发明方法仅选取二维特征点返回三维空间后与上边界层和布鲁赫莫层的交点。

有益效果

与现有技术相比，本发明具有以下优点：

(1)对视网膜图像进行三维的配准，以使得配准结果能够适应三维各维度上的形变；

(2)可应用于病变视网膜的图像配准，即使视网膜组织产生较大形变，本发明的配准方法依然具备高有效性，特别适用于脉络膜新生血管的视网膜oct图像的三维配准。

综上，本发明提供了一种三维oct图像的配准方，结合了基于的特征的仿射变换和基于灰度的非刚性变换两种技术。在采用本发明之后，使用仿射变换得到的dice系数是0.3左右，使用非刚性变换的dice系数是0.7左右，配准有效性得到了较大的提升，能够大大提高配准图像在后续应用中的数据可靠性。

附图说明

图1所示为本发明方法流程示意图；

图2所示为曲线异向扩散滤波去噪示意图；图2-(a)为原图像，图2-(b)为去噪后图像；

图3所示为图像预处理中除oct去噪步骤外的其它步骤示意图；图3-(a)为分层后的图像，图3-(b)为投影图像，图3-(c)为移除病变区域之后的图像，图3-(d)为提取血管后的图像，图3-(e)为下采样之后的图像；

图4所示为部分测试数据层手动分割示意图；rnfl-神经纤维层(retinalnervefiberlayer)；gcl+ipl-节细胞层和内丛状层(ganglioncelllayerandinnerplexiformlayer)；inl+opl-内核层和外丛状层(innernuclearlayerandouterplexiformlayer)；onl-外核层(outernuclearlayer)；is/os+rpe-内/外段层和视网膜色素上皮层(innersegment/outersegmentandretinalpigmentepithelium)；

图5所示为利用本发明进行图像配准的实施例图像过程变化示意图；图中第1列为参考图像，第2列为待配准图像，第3列为仿射变换的粗配准结果，第4列为非刚性的精确配准结果。

具体实施方式

以下结合附图和具体实施例进一步描述。

参考图1所示，本发明的三维非刚性光学相干断层扫描图像的配准方法，包括步骤：

s1，对oct扫描图像进行预处理，包括依次对oct扫描图像进行oct去噪、图像分层和图像投影，得到视网膜血管的二维投影图像；

s2，从视网膜血管的二维投影图像中提取血管图像；

s3，对提取的血管图像进行二维特征点的提取，根据得到的二维特征点返回到三维空间得到三维特征点；

s4，对固定图像以及待配准图像进行基于特征点的仿射变换，得到粗配准的视网膜图像；

s5，利用基于灰度的非刚性配准方法对粗配准的视网膜图像进行进一步配准，得到精确配准的结果。

实施例

以下分别对各个步骤进行详细说明。

s1，图像预处理

图像预处理主要包括以下四个步骤：oct去噪、图像的分层、投影图像的获取以及病变区域的移除。

s11，oct图像去噪：

oct眼部成像仪获取的三维图像含有较多的散斑噪声，为保证后续分割的效果，必须在有效去除噪声的同时尽可能保留图像中的边缘信息。本发明采用曲线异向扩散滤波方法对视网膜图像进行滤波，可以去掉散斑噪声，并且图像的边界能清晰地保留。

曲线异向扩散滤波方法为现有算法，如孙助力等人提出的曲线异向扩散滤波方法。其中曲线扩散方程为：

其中，f是输入图像即原始图像，ft是去噪后的输出图像，▽是梯度算子，c()是传导率变换函数。该滤波方法在用c++代码实现时涉及三个参数：时间步长，迭代次数和传导率，本发明分别对应取值为0.0625秒，5次，3。去噪后结果如图2所示。

s12，图像分层：

使用图搜索算法将视网膜oct图像分为5层，第一层是神经纤维层的上边界，结果如图3(a)中红线所示；第二层是神经纤维层的下边界，结果如图3(a)中绿线所示；第三层是内核层的边界，结果如图3(a)中蓝线所示；第四层是内外视网膜的边界，结果如图3(a)中黄线所示；第五层是布鲁赫莫，结果如图3(a)中紫线所示。图搜索算法可采用现有的石霏等人提出的图搜索算法。

s13，投影图像的获取：

对外层视网膜进行投影，即计算第四层与第五层之间的平均灰度值，就可以得到视网膜血管的投影图像。

定义沿三维空间的y轴对视网膜图像进行分层；即分层后的各层图像沿y轴排列，同时各层图像之间视网膜图像相应位置的x轴和z轴坐标不变，仅y轴坐标变化。在计算坐标为(x,z)的投影图像位置的灰度值i(x,z)时，公式如下：

上式中，x、y、z表示图像的坐标(即图像的大小)，y1表示第四层的y轴坐标，y2表示第五层的y轴坐标，i(.)函数表示像素值(灰度值)，n表示视网膜图像的(x,z)坐标处第四层到第五层之间的像素个数。当x、z固定时，即可以通过求第四层与第五层之间的平均灰度值得到投影图像在坐标为(x,z)位置处的像素灰度值，进而得到全部投影图像。

同理，只要x、y、z中任意两维坐标确定，即可通过另一维变化的坐标，求取两层之间的平均灰度值得到投影图像。

在获得投影图像后，将得到的图像上采样至512×512个像素，可便于血管的提取，结果如图3(b)中所示。

s14，病变区域的移除

对于有新生血管疾病的病人，对视网膜图像进行投影，其病变区域会造成非血管的阴影区域，这片区域的灰度值与血管的灰度值非常的接近，在检测血管时，会有误检，所以本发明利用图像处理技术中的数学形态学方法将这片区域移除，更有利于提取血管点，结果如图3(c)中所示。

s2，血管提取

利用frangi等人提出的血管增强滤波器进行血管的提取。该算法的原理是利用血管的管状结构来提取血管，将血管增强视为一个滤波的过程，寻找三维数据集中类似管状的图形结构。需要分析一个像素是不是血管，就需要分析他的局部区域的特性，通过分析泰勒展开式中的hessian矩阵的特征向量和特征值来对血管进行判断，最终实现血管的提取。

s3，特征点提取

因为要进行三维的视网膜图像的配准，所以要提取三维的特征点。通常情况下，因为血液对于光线的吸收，在oct图像中，血管会出现在高反射区。本发明首先对提取的血管图像提取出二维特征点，即提取血管图像中的二维的血管点作为二维特征点，然后根据得到的二维的特征点返回到三维空间得到三维的特征点。二维特征点的提取为现有技术。

根据得到的二维的特征点返回到三维空间，在三维空间里寻找灰度值最大的点即为血管点；同时获取二维特征点返回三维空间后与视网膜各分层的交接点。本发明仅关注二维特征点分别与第一层(上边界层)和第五层(布鲁赫莫层)的交点。所述血管点以及二维特征点与分层的交点组成三维特征点。

当已知三维坐标中的两维坐标，在第三维中找取灰度值最大的那个点，找到该点对应的坐标，即得到准确的三维坐标。本发明在寻找灰度值最大点时，采用中值法从多个灰度值较大点中选取坐标在中间的一个。在选择二维特征点与各分层交点进行应用时，实质上，二维特征点与每个分层皆存在交点，但全部选取会影响算法效率，在算法结果上却没有明显进步，故本发明方法仅选取二维特征点返回三维空间后与上边界层和布鲁赫莫层的交点。

s4，对固定图像以及待配准图像进行基于特征点的仿射变换，得到粗配准的视网膜图像

本实施例的仿射变换步骤使用的是myronenko等人提出的相干点漂移算法。在相干点漂移算法中，相干点的其中一个点集表示高斯混合模型的质心，另一个点集表示的是数据集，仿射变换的原理即，采用确定性的模拟退火的em(最大期望算法，expectationmaximizationalgorithm)算法对最大似然估计进行最优化，从而找到两个点集的对应关系与仿射变换。仿射变换定义为：

t(ym；r,t,s)＝bym+t(3)

其中bd×d表示仿射变换矩阵，td×1是平移矢量，ym表示高斯混合模型的质心，r、s为矩阵b中的参数，r为旋转矩阵，s为缩放倍数。

仿射变换目标方程为：

其中中间变量：

上述公式中，d表示两个点集的维数，n表示数据集xn中点的个数，m分别表示点集ym中点的个数，xn表示输入数据集，ym表示高斯混合模型的质心，θ表示一系列变换的参数，在此处即表示b和t。参数k为1至m的整数，参数σ在迭代过程中是变化的，初始值可根据需要设定。ω为权重参数，本发明中ω设置为0。上标old表示由迭代计算得到。

仿射变换算法的实现为：步骤s3中提取出的三维特征点作为仿射变换算法的输入，即ym和xn，其中固定图像的三维特征点的点集作为数据集xn，待配准图像的三维特征点的点集作为高斯混合模型的质心ym。通过em算法迭代计算出θ和σ²。首先，估计初始参数，利用贝叶斯公式计算后验概率分布，即p^old(m|xn)；然后，通过最小化前述仿射变换目标方程来更新参数θ和σ²，得到最终的最小的变换参数bd×d和td×1。最后根据参数bd×d，td×1对待配准图像进行变换，得到仿射变换后的粗配准图像。

s5，利用基于灰度的非刚性配准方法对粗配准的视网膜图像进行进一步配准，得到精确配准的结果

由于视网膜的一些变形生长及病变区域的变化，所以仅仅使用仿射变换是不能解决问题的，还需要运用非刚性配准，本实施例运用的是rueckert等人提出的基于b样条的自由变换模型，它是基于灰度的非刚性配准的重要手段。

基于b样条的自由变换模型在应用时，只要是使用样条函数和控制点来描述非线性几何变换域，将位于样条函数栅格上的控制点赋值，当这些控制点产生空间位移的时候，样条函数被弯曲，那么控制点周围的图像空间也会随着样条函数的弯曲而产生不同程度的扭曲。由于限制了节点位移的大小，一层固定网格间距的b样条模型就不能完成大规模形变，所以使用了多层次的b样条模型。

在三维空间域中，图像区域表示为ω＝{(x,y,z)|0≤x＜x,0≤y＜y,0≤z＜z}，(x,y,z)代表组成图像的各点的坐标，x、y、z代表坐标大小；φ是图像域ω中的一组点集，由等间距为δ的控制点φi,j,k组成，等间距为δ的控制点φi,j,k构成nx×ny×nz的网格，控制点φi,j,k即网格中的交点，nx、ny、nz代表网格的尺寸；

则基于b样条函数的自由形变可描述为一阶三次b样条函数的张量积t：

其中中间变量：

βi(.)为第i次b样条的基函数，i即式(7)中的l、m、n，取值为0,1,2,3：

这样就将配准问题转化为三次b样条的局部控制问题，这里三次b样条基函数每个控制点φi,j,k的作用范围在周围的两个像素点之间，只要计算出固定图像以及待配准图像的对应点周围领域的相似性测度即可得到最终的精确配准图像。

方法验证

利用本发明在3个有脉络膜新生血管疾病的病人的12个数据上进行实验。在这三个病人在接受药物治疗的过程中，每隔一个月采集一次数据，一共采集了5次，所以每个病人有五个时间点，即每个病人有5个oct数据，一个oct数据的大小为512像素×1024像素×128像素。其中第一次采集的数据作为固定图像，其它4次作为待配准图像，每人利用本发明进行4次配准。在该方法中，每个病人第一次采集的图像作为固定图像，其余剩下的数据作为待配准图像。部分实验结果如图5所示，图5中第一列是固定图像，第二列是待配准图像，第三列是仿射变换之后的粗配准结果，第四列是非刚性变换之后的精确配准的结果。

采用视网膜各层的dice系数dk来评判配准的结果。dice系数的公式表示为：

其中tk是待配准图像的第k层视网膜层的所有像素点，sk是固定图像的第k层视网膜层的所有像素点。

dice系数越接近1，说明两个图像的重合度越高，即配准更有效，0就表示完全没有重合。

采用本发明之后，使用仿射变换得到的dice系数是0.3左右，进一步使用非刚性变换后的dice系数是0.7左右，得到了较大的提升。本发明不仅在dice系数上游较大的提升，而且图像的一些局部区域也得到很好的配准。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。