一种基于点集匹配的图像配准方法

摘要

本发明公开了一种基于点集匹配的图像配准方法，通过对浮动图像和参考图像进行特征点提取，得到浮动图像特征点集和参考图像特征点集，采用了一种构造中间变换的方法，使得待匹配的二个点集之间的旋转角度小于90度。本发明通过构造由12个四元数表征的空间映射信息的集合，对于浮动图像特征点集作用这12个初始变换会得到12个中间图像特征点集，构造的这12个中间特征点集中至少会存在一个点集与参考图像特征点集的绕轴旋转角度小于等于90度。本发明可以解决任意角度的旋转点集配准，从而实现图像刚性配准。本发明在医学图像配准领域有重要的应用价值。

说明书

技术领域

本发明属于图像配准技术领域，具体涉及一种基于点集匹配的图像配 准的方法，应用于医学图像配准领域。

背景技术

图像配准是医学图像处理和分析的一个关键问题，它可以帮助人们进 行三维重建、图像融合、图像分割以及病理分析，在医学图像领域有重要 的应用价值。

图像配准一般可以分为二种，基于图像灰度信息的配准(Wells et al. 1996.Multi-modal volume registration by maximization of mutual information) 和基于图像特征的配准(Wu et al.2007.Learning best features and  deformation statistics for hierarchical registration of MR brain images)。基于 图像特征配准的方法有四个主要步骤：一是对于图像进行特征提取，二是 对于图像特征进行匹配，三是构造映射函数，四是对图像进行映射变换。 基于图像特征配准的方法，可以不需要为特征点进行属性构造，仅仅保留 特征点的位置信息，这将节省大量的时间，现有的点集匹配方法包括最近 点迭代算法(Besl et al.1992.Method for registration of3-D shapes)和高斯 混合模型点集匹配方法(Jian et al.2011.Robust point set registration using  gaussian mixture models)等，但是现有的点集匹配方法都无法解决点集旋 转大于90度的问题。

发明内容

(一)要解决的技术问题

本发明所要解决的技术问题是现有的图像配准方法对较大角度旋转 的三维图像的配准效果不佳的问题。

(二)技术方案

(三)有益效果

本发明通过构造由12个四元数表征的空间映射信息的集合，对于用 NSIFT方法从浮动图像中提取的浮动图像特征点集作用这12个初始变换 会得到12个中间图像特征点集，构造的这12个中间特征点集中至少会存 在一个点集与用NSIFT方法从参考图像中提取的特征点集的绕轴旋转角 度小于等于90度。然后利用高斯混合模型点集匹配方法对于这12个中间 图像点集分别与参考图像点集进行匹配，并求出几何变换。然后从这12 个几何变换中挑选出最优的空间映射使得浮动图像经此映射变换后得到 的配准图像与参考图像的差异最小。解决了现有的图像配准方法对较大角 度旋转的三维图像的配准效果不佳的问题。

附图说明

图1示意了分别采用ICP、GMM和本发明的方法对于浮动图像与参 考图像旋转角度为10°到180°的配准误差。

图2示意了分别采用ICP、GMM和本发明的方法对于浮动图像与参 考图像旋转角度为60°时的配准直观视觉结果，(a)是参考图像切片，(b)是 浮动图像切片，(c)是参考图像切片和浮动图像切片的棋盘差异图像，(d)～ (f)分别表示用ICP、GMM、和本发明的方法得到的配准图像的切片，(g)～ (i)分别表示用ICP、GMM和本发明的方法得到的配准图像的切片与参考 图像切片的棋盘差异图像。

图3示意了分别采用ICP、GMM和本发明的方法对于浮动图像与参 考图像旋转角度为90°时的配准直观视觉结果，(a)是参考图像切片，(b)是 浮动图像切片，(c)是参考图像切片和浮动图像切片的棋盘差异图像，(d)～ (f)分别表示用ICP、GMM和本发明的方法得到的配准图像的切片，(g)～ (i)分别表示用ICP、GMM和本发明的方法得到的配准图像的切片与参考 图像切片的棋盘差异图像。

图4示意了分别采用ICP、GMM和本发明的方法对于浮动图像与参 考图像旋转角度为120°时的配准直观视觉结果，(a)是参考图像切片，(b) 是浮动图像切片，(c)是参考图像切片和浮动图像切片的棋盘差异图像， (d)～(f)分别表示用ICP、GMM和本发明的方法得到的配准图像的切片， (g)～(i)分别表示用ICP、GMM和本发明的方法得到的配准图像的切片与 参考图像切片的棋盘差异图像。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实 施例，并参照附图，对本发明作进一步的详细说明。

本发明可在具有可编程能力的GPU硬件平台上实现。在一个实施例 中，在NVIDIA GTX480显卡上通过CUDA接口实现。在实施该实施例 之前，首先实现NSIFT在高斯差分空间对于特征点的提取，以及根据现有 技术实现一个基本的高斯混合模型点集匹配方法。

步骤S1：对于浮动图像和参考图像进行特征点的提取，浮动图像和参 考图像就可以表示为浮动图像特征点集和参考图像特征点集。

利用GPU加速(CUDA)根据文献(Cheung et al.2007.Scale invariant  feature transform for n-dimensional images)中的方法构造高斯差分金字塔 对浮动图像和参考图像快速检测关键点，在获取关键点时，我们仅保留其 位置信息，而不对关键点进行描述子构造这一步。得到浮动图像特征点集 参考图像特征点集

步骤S2：构造中间图像特征点集，根据构造的12个四元数表征的空 间映射信息的集合，对于浮动点集分别作用这12个空间映射，我们可以 得到12个中间图像特征点集。

$\left(\begin{array}{c}\Omega =\{{\left(\mathrm{1,0,0,0}\right)}^T,{\left(\mathrm{0,1,0,0}\right)}^T,{\left(\mathrm{0,0,1,0}\right)}^T,{\left(\mathrm{0,0,0,1}\right)}^T,\\ {(\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2})}^T,{(-\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2})}^T,{(\frac{1}{2},-\frac{1}{2},\frac{1}{2},\frac{1}{2})}^T,{(\frac{1}{2},\frac{1}{2},-\frac{1}{2},\frac{1}{2})}^T,\\ {(-\frac{1}{2},\frac{1}{2},-\frac{1}{2},\frac{1}{2})}^T,{(\frac{1}{2},-\frac{1}{2},-\frac{1}{2},\frac{1}{2})}^T,{(-\frac{1}{2},-\frac{1}{2},-\frac{1}{2},\frac{1}{2})}^T,{(-\frac{1}{2},-\frac{1}{2},\frac{1}{2},\frac{1}{2})}^T\}\end{array}\right)$

假定四元数q₁，q₂….q₁₂分别代表上述集合Ω的第一个到第12个元素， 四元数和旋转矩阵满足一下关系式：

q＝(r₁，r₂，r₃，r4)^T的空间旋转信息可以用旋转矩阵R表示为：

$R=\left(\begin{array}{ccc}{r_4}^2+{r_1}^2-{r_2}^2-{r_3}^2& 2(r_1{r}_2+r_3{r}_4)& 2(r_3{r}_1-r_4{r}_2)\\ 2(r_1{r}_2-r_3{r}_4)& {r_4}^2-{r_1}^2+{r_2}^2-{r_3}^2& 2(r_2{r}_3+r_1{r}_4)\\ 2(r_1{r}_3+r_2{r}_4)& 2(r_2{r}_3-r_1{r}_4)& {r_4}^2-{r_1}^2-{r_2}^2+{r_3}^2\end{array}\right)$

可以求出q₁，q₂……q₁₂所对应的旋转矩阵R₁，R₂……R₁₂。对浮动点集分别作 用这12个旋转矩阵，我们得到12个中间图像的特征点集为

步骤S3：利用高斯混合模型点集匹配方法对于这12个中间图像点集 分别与参考图像点集进行匹配。求出几何变换。这12个几何变换分别与 步骤2所对应的空间映射进行复合运算，得到从浮动图像特征点集到参考 图像特征点集的最终空间映射12个。

因为在P₁，P₂……P₁₂这12个点集中至少存在着一个点集，它与参考 图像特征点集之间的绕轴旋转角度小于或等于90度，而高斯混合点集匹 配方法能完美解决旋转角度小于90的点集匹配方法。对这P₁，P₂……P₁₂与 F分别利用高斯混合模型点集匹配方法算出几何变换分别为 R_PF1，t₁；R_PF2，t₂，……R_PF12，t₁₂(R_PF1；R_PF2，……R_PF12是旋转矩阵，t₁；t₂，……t₁₂是平移向量)。

这12个几何变换分别与步骤2所对应的空间映射进行复合运算，得 到的最终的12个空间映射为：R_PF1R₁，t₁；R_PF2R₂，t₂，……R_PF12R₁₂，t₁₂(R_PF1R₁；R_PF2R₂，……R_PF12R₁₂是复合后旋转矩阵，t₁；t₂，……t₁₂是平移向量)。

步骤S4我们对浮动图像用以R_PF1R₁为旋转矩阵，t₁为平移向量进行几 何变换得到配准图像registration1，以R_PF2R₂为旋转矩阵，t₂为平移向量进 行几何变换得到配准图像registration2……以R_PF12R₁₂为旋转矩阵，t₁₂为平移 向量进行几何变换得到配准图像registration12。然后在这12幅配准图像 registration1，registration2……registration12中挑选与参考图像最相近的图 片。

至此，我们就完成了图像配准。

图1示意了ICP，GMM，和本发明的方法对于浮动图像与参考图像旋 转角度为10°到180°的配准误差(当误差较大时，认为没有完成配准，不 记录此次误差)。

从图1中可以看出，当旋转角度大于60度时，ICP的配准方法会失效。 当旋转角度大于90度时，GMM的配准方法会失效，本发明的方法在任意 旋转角度都能对图像做很好的配准。

图2示意了ICP，GMM，和本发明的方法对于浮动图像与参考图像旋 转角度为60°时的配准直观视觉结果，(a)是参考图像切片，(b)是浮动图像 切片，(c)是参考图像切片和浮动图像切片的棋盘差异图像，(d)-(f)分别表 示用ICP，GMM，和本发明的方法得到的配准图像的切片，(g)-(i)分别表 示用ICP，GMM，和本发明的方法得到的配准图像的切片与参考图像切片 的棋盘差异图像。

从图2中可以看出，当旋转角度为60度的时候，ICP，GMM，和本发明 的方法都很很好的对图像进行配准。

图3示意了ICP，GMM，和本发明的方法对于浮动图像与参考图像旋 转角度为90°时的配准直观视觉结果，(a)是参考图像切片，(b)是浮动图像 切片，(c)是参考图像切片和浮动图像切片的棋盘差异图像，(d)-(f)分别表 示用ICP，GMM，和本发明的方法得到的配准图像的切片，(g)-(i)分别表 示用ICP，GMM，和本发明的方法得到的配准图像的切片与参考图像切片 的棋盘差异图像。

从图3中可以看出，当旋转角度为90度的时候，ICP的方法无法对图 像进行有效的配准，而GMM和本发明的方法均能对图像进行配准。

图4示意了ICP，GMM，和本发明的方法对于浮动图像与参考图像旋 转角度为120°时的配准直观视觉结果，(a)是参考图像切片，(b)是浮动图像 切片，(c)是参考图像切片和浮动图像切片的棋盘差异图像，(d)-(f)分别表 示用ICP，GMM，和本发明的方法得到的配准图像的切片，(g)-(i)分别表 示用ICP，GMM，和本发明的方法得到的配准图像的切片与参考图像切片 的棋盘差异图像。

从图4中可以看出，当旋转角度为120度的时候，ICP和GMM的方 法无法对图像进行有效的配准，仅有本发明的方法能够对图像进行配准。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行 了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已， 并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、 等同替换、改进等，均应包含在本发明的保护范围之内。