CN111798500B - 一种基于层次邻域谱特征的微分同胚非刚性配准算法 - Google Patents

Abstract

本发明属于图像处理技术领域，具体公开了一种基于层次邻域谱特征的微分同胚非刚性配准算法。由于层次邻域谱特征具有旋转不变性、亮度线性变化不变性和抗噪声能力强的优点，首先提取参考图像和浮动图像的层次邻域谱特征；其次，利用层次邻域谱匹配算法将这种新的谱特征融合到微分同胚Demons配准框架的能量函数中，从而提高微分同胚Demons算法捕获复杂形变的能力，并保证微分同胚变换和配准精度；最后，在整个配准过程中采用基于小波分解的多分辨率策略，进一步提高配准精确性和效率。本发明不仅能有效地生成光滑和可逆的形变场，而且比经典配准算法具有更好的性能。

Description

一种基于层次邻域谱特征的微分同胚非刚性配准算法

技术领域

本发明属于图像处理技术领域，具体涉及一种基于层次邻域谱特征的微分同胚非刚性配准算法。

背景技术

医学图像配准是指对浮动图像进行合理的空间变换，使其与参考图像上相应的解剖点或至少诊断点达到相同空间关系的过程。根据空间变换的不同，医学图像配准一般分为刚性配准和非刚性配准两类。刚性配准只描述仅限于全局旋转和平移的运动，而非刚性配准通常包含非常复杂的局部和全局弹性形变。实际上，在图像采集过程中，同一个或不同患者往往会受到肺部运动、膀胱充盈等因素影响，导致复杂的非刚性形变。非刚性配准在自适应放疗、图像引导手术、疾病诊断、病变跟踪和治疗效果评价等医学临床应用中起着重要的作用。同时，非刚性配准仍然是医学图像分析中最具挑战性的课题之一。

目前，医学图像非刚性配准有很多经典模型，如基于B样条的自由形变(FFD)模型、有限元模型(FEM)、粘性流体模型、Demons模型等。Demons配准是一种非常经典和流行的基于光流场理论的非刚性配准算法，由于其完整的数学理论，在医学图像非刚性配准中得到了广泛地应用。然而，Demons算法存在以下三个缺点：(1)由于只利用参考图像的梯度信息来驱动形变，算法的收敛速度比较慢；(2)当参考图像无纹理区域的梯度接近于零时，很容易得到错误的配准结果；(3)Demons算法很难估计复杂形变。为了加快Demons算法的收敛性，将浮动图像的梯度信息引入到扩散方程中，同时利用均匀化系数来调整驱动力的强度，提出了Active Demons算法。然而，Active Demons算法只通过均匀化系数来调整驱动力的强度，无法实现大形变和小形变之间的平衡。为了平衡大形变和小形变，在Active Demons中引入了一个新的称为平衡系数的参数来调节大形变和小形变之间的驱动力，提出了改进的Active Demons算法。为了保持形变场的拓扑结构，避免物理上的不合理形变，将配准过程看作是能量函数的优化过程，提出了基于李群理论的微分同胚Demons非刚性配准算法。尽管微分同胚Demons是一个很好的基于图像灰度的配准框架，但它也不能捕获复杂的大尺度形变。为了解决这一问题，同时提高配准精度，许多研究者针对微分同胚Demons算法提出了很多改进算法，大致分为两种策略：

(1)使用新的相似性测度。基于微分同胚Demons，通过使用新的相似性测度，提出了很多改进算法。例如：互相关测度被集成到对称的微分同胚Demons配准框架中，改进了医学图像非刚性配准性能。局部相关系数LCC被引入到Log-Demons配准中，提出了LCC-Demons非刚性配准算法。

(2)引入新特征到驱动力中。这些特征主要包括：灰度梯度场的相似性、几何约束特征、鲁棒性特征描述子SIFT、SURF和ASIFT、基于SAE的无监督深度特征、模态独立邻域描述子和谱特征等。例如，将谱特征引入到Demons配准算法的驱动力中，提出了SpectralDemons非刚性配准算法。Spectral Demons是一种新的捕获大形变的配准框架，与传统的Demons算法相比，在配准精度和对大形变的鲁棒性方面都有明显提高。然而，SpectralDemons得到的最终变换只适用于拓扑结构相似的图像配准，不适用于肿瘤切除前后的图像配准，不利于术后康复的监测和评价。此外，高维矩阵的谱分解导致Spectral Demons非常耗时，限制了其在许多临床应用中的实用性。

然而，由于引入了新的相似性或新的特征，导致利用上述两种策略改进的许多非刚性配准算法通常比较耗时。在实际应用中，加速这些改进算法的方法主要有两种：多分辨率策略和并行计算技术。研究人员将多分辨率策略应用于配准过程中，以加快配准收敛速度。在采用多分辨率策略进行配准时，首先对低分辨率图像进行较为复杂的相似性度量，既快速又能抑制局部极值；然后将粗配准阶段的配准参数设置为精细配准阶段的初始参数，可以缩短优化过程的时间，有效提高配准效率。此外，基于GPU等多核硬件的并行计算技术可以加速许多非刚性配准算法，是提高配准性能的重要途径。

综上所述，上述算法要么不保留细节，要么不具有微分同胚性，要么非常耗时，对于复杂而严重的器官形变，它们不能提供令人满意的结果。因此，如何对具有复杂形变的医学图像实现精确地非刚性配准还需进一步深入研究。

发明内容

本发明的目的在于提供了一种基于层次邻域谱特征的微分同胚非刚性配准算法，针对具有复杂形变的图像，能够实现精确地非刚性配准。

为了实现以上目的，本发明提供了一种基于层次邻域谱特征的微分同胚非刚性配准算法，包括以下步骤：

(1)对参考图像R和浮动图像M分别生成参考图像金字塔R_i(i＝1，2，…，N)和浮动图像金字塔M_i(i＝1，2，…，N)；

(2)初始化应用于第一层浮动图像M₁上的速度场v₁；

(3)设置当前待配准的图像层数i＝1，从第一层开始进行配准；

(4)判断i≤N是否成立，如果是则继续步骤(5)，否则转至步骤(14)；

(5)如果i＞1，初始化当前层的速度场v_i；

(6)判断配准条件是否收敛，如果是则转至步骤(13)，否则继续步骤(7)；

(7)利用层次邻域谱匹配算法计算从R到

的更新场u_R→M；

(8)利用层次邻域谱匹配算法计算从M到

的更新场u_M→R；

(9)计算平均更新场u←(u_R→M-u_M→R)/2；

(10)平滑平均更新场u←K_fluid*u；

(11)利用平滑后的更新场u来计算速度场

；

(12)平滑速度场v_i←K_diffuse*v_i，转至步骤(6)执行；

(13)输出当前层的最优配准变换φ_i ^*，然后i＝i+1并转至步骤(4)执行；

(14)输出最终的最优变换φ^*。

进一步地，所述步骤(1)中采用基于小波分解的多分辨率策略分别对参考图像R和浮动图像M进行多分辨率分解，分别生成参考图像金字塔R_i(i＝1，2，…，N)和浮动图像金字塔M_i(i＝1，2，…，N)。

进一步地，所述步骤(7)和步骤(8)中层次邻域谱匹配算法，包括：

(a)对参考图像R和浮动图像M分别逐像素遍历，并判断是否遍历结束，如果没有结束，设定当前中心像素p并顺序执行步骤(b)，如果当前图像已遍历结束，则跳转至步骤(i)；

(b)初始化层次邻域层数k＝10，当前层数l＝1；

(c)判断l≤k是否成立，如果成立，顺序执行步骤(d)，针对当前中心像素p开始进行层次邻域谱特征λ_l计算，否则转至步骤(h)；

(d)针对当前中心像素p，构造当前层对应的线图L(S)；

(e)生成线图L(S)对应的邻接矩阵A；

(f)对邻接矩阵A进行谱分解，从而计算当前第l层的层次邻域谱特征

；

(g)从当前第l层的层次邻域谱特征λ_l中选择该层的主谱特征

；

(h)连接当前中心像素p的每层对应的主谱特征

来生成当前中心像素p的最终层次邻域谱特征χ_R(p)或χ_M(p)，转至步骤(a)更新中心像素p；

(i)根据参考图像R对浮动图像M的层次邻域谱特征χ_R(p)或χ_M(p)进行排序；

(j)将参考图像R和浮动图像M的层次邻域谱特征χ_R(p)和χ_M(p)分别嵌入到各自的特征空间R＝(I_R,X_R,χ_R)和M＝(I_M,X_M,χ_M)中；

(k)寻找映射c实现从参考图像R到浮动图像M的映射

进一步地，所述步骤(d)中构造当前层对应的线图L(S)，包括：

首先构建星图S；然后在星图S中，用边pq_j的权重r_j相减来构造线图L(S)。

更进一步地，所述星图S构建时，中心像素p的第l(1≤l≤k,l∈R⁺)层的8l邻域像素为q_j(j＝1,2,…,8l)，利用中心像素p和8l邻域像素q_j(j＝1,2,…,8l)构建星图S，在星图S中，中心像素p与所有邻域像素q_j(j＝1,2,…,8l)之间的边为pq_j，每条边的权重r_j为：

其中，I(p)和X(p)分别是中心像素点p的灰度值和空间坐标，I(q_j)和X(q_j)分别是像素点q_j的灰度值和空间坐标。

进一步地，所述步骤(e)生成线图L(S)对应的邻接矩阵A，其邻接矩阵A的元素为：

其中，1≤i,j≤8l,i,j∈R⁺。

更进一步地，所述邻接矩阵A的8l个实数特征通过谱分解来获得，且λ₁(p)≥λ₂(p)≥…≥λ_8l(p)，使用λ_i(p)来代替像素点p的灰度值并正则化到[0，255]，计算大小为(m-1)×(n-1)的特征图I_i(i＝1,2,…,8l)。

进一步地，所述步骤(f)对邻接矩阵A进行谱分解，从而计算当前第l层的层次邻域谱特征

包括：

S1：初始化采样列数m，秩n和过采样参数o；

S2：从矩阵A中随机采样m列和m行生成矩阵W；

S3：生成高斯随机矩阵Ω_m×(n+o)；

S4：Y←WΩ；

S5：通过QR分解来计算Q：Y＝QQ^TY；

S6：B←Q^TWQ；

S7：通过SVD分解来计算V：B＝VΛV ^T；

S8：

S9：输出U和Λ，并将λ_l←Λ。

进一步地，所述步骤(g)从当前第l层的层次邻域谱特征λ_l中选择该层的主谱特征

包括：设置r_l＝2(1≤l≤k)，从当前第l层的层次邻域谱特征λ_l中选择前r_l项构成该层的主谱特征/>

进一步地，所述步骤(h)包括：顺序组合每层的主谱特征λ'_l(l＝1,2,…,k)，从而生成像素点p的最终层次邻域谱特征χ(p)＝(λ'₁/||λ'₁||,λ'₂/||λ'₂||,…,λ'_k/||λ'_k||)。

相比于现有技术，本发明具有以下有益效果：

本发明方法是基于层次邻域谱特征的微分同胚非刚性配准算法，由于层次邻域谱特征具有旋转不变性、亮度线性变化不变性和抗噪声能力强的优点，因此利用层次邻域谱匹配算法提取参考图像和浮动图像的层次邻域谱特征，并对这两幅图像的层次邻域谱特征进行相似性度量，从而将这种新的谱特征融合到微分同胚Demons配准框架的能量函数中提高捕获复杂形变的能力，不仅解决了复杂的形变问题，而且保证了微分同胚变换和配准精度。

进一步地，在层次邻域谱匹配算法中采用基于随机奇异值分解(SVD)的

近似方法，有效地减少了高维矩阵的谱分解时间，提高了本发明的配准效率。

进一步地，本发明中采用基于小波分解的多分辨率策略对参考图像和浮动图像进行多分辨率分解，在配准过程中，采用基于小波分解的多分辨率策略，能够进一步提高配准精确性和效率。

本发明能有效地生成光滑和可逆的形变场，比经典配准算法具有更好的性能。

附图说明

图1为本发明的流程原理图；

图2为层次邻域谱特征的构造过程图；

图3为使用七种算法对六组不同模态和位置的医学图像进行配准的实验结果图，其中，(a)为参考图像；(b)为浮动图像；(c)-(i)分别为采用基于B样条的FFD、Demons、Active Demons、改进Active Demons、Log-Demons、Spectral Demons和本发明的方法对待配准图像进行配准的结果；

图4为使用七种算法对四组真实图像进行配准的实验结果图，其中，(a)为参考图像；(b)为浮动图像；(c)-(i)分别为采用基于B样条的FFD、Demons、Active Demons、改进Active Demons、Log-Demons、Spectral Demons和本发明方法对待配准图像进行配准的结果；

图5为针对两组待配准图像使用三种算法进行配准时的能量函数值随迭代次数变化的曲线图，其中，(a)为一组医学图像配准过程中的能量函数随迭代次数变化曲线图；(b)为一组真实图像配准过程中的能量函数随迭代次数变化曲线图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细的描述。

参见图1和图2，本发明提出的基于层次邻域谱特征的微分同胚非刚性配准算法，具体步骤包括：

步骤(1)：采用基于小波分解的多分辨率策略分别对参考图像R和浮动图像M进行多分辨率分解，分别生成参考图像金字塔R_i(i＝1，2，…，N)和浮动图像金字塔M_i(i＝1，2，…，N)；

步骤(2)：初始化应用于第一层浮动图像M₁上的速度场v₁；

步骤(3)：设置当前待配准的图像层数i＝1，从第一层开始进行配准；

步骤(4)：判断i＜＝N是否成立，如果是则继续步骤(5)，否则转至步骤(14)；

步骤(5)：如果i＞1，初始化当前层的速度场v_i；

步骤(6)：判断配准条件是否收敛，如果是则转至步骤(13)，否则继续步骤(7)；

步骤(7)：利用层次邻域谱匹配算法计算从R到

的更新场u_R→M，其中，层次邻域谱匹配算法的具体步骤主要包括：

(a)对参考图像R和浮动图像

分别进行逐像素遍历，并判断是否遍历结束，如果没有结束，设定当前中心像素p并顺序执行步骤(b)，如果当前图像已遍历结束，则跳转至步骤(i)；

(b)初始化层次邻域层数k＝10，当前层数l＝1；

(c)判断l≤k是否成立，如果成立，顺序执行步骤(d)，针对当前中心像素p开始进行层次邻域谱特征λ_l计算，否则转至步骤(h)；

(d)针对当前中心像素p，构造当前层对应的线图L(S)，具体步骤主要包括：

首先，构建星图S，针对大小为m×n的图像I，任何非边界像素p的特征可以用它与邻域像素的关系来定义，因此，像素p的特征用图来表示，假设p的第l(1≤l≤k,l∈R⁺)层的8l邻域像素为q_j(j＝1,2,…,8l)，用像素p和8l邻域像素q_j(j＝1,2,…,8l)构建星图S。在星图S中，中心像素p与所有邻域像素q_j(j＝1,2,…,8l)之间的边为pq_j，每条边的权重r_j为：

其中，I(p)和X(p)分别是像素点p的灰度值和空间坐标，I(q_j)和X(q_j)分别是像素点q_j的灰度值和空间坐标；

其次，为进一步反映像素点p的特征，在星图S中，用边pq_j的权重r_j相减来构造线图L(S)；

(e)生成线图L(S)对应的邻接矩阵A，其中邻接矩阵A的元素为：

其中，1≤i,j≤8l,i,j∈R⁺。由于邻接矩阵A是一个实对称矩阵，8l个实数特征可以通过谱分解来获得，并且λ₁(p)≥λ₂(p)≥…≥λ_8l(p)，因此使用λ_i(p)来代替像素点p的灰度值并正则化到[0，255]，可以计算大小为(m-1)×(n-1)的特征图I_i(i＝1,2,…,8l)；

(f)对邻接矩阵A进行谱分解，从而计算当前第l层的层次邻域谱特征

具体步骤包括：

S1：初始化采样列数m，秩n和过采样参数o；

S2：从矩阵A中随机采样m列和m行生成矩阵W；

S3：生成高斯随机矩阵Ω_m×(n+o)；

S4：Y←WΩ；

S5：通过QR分解来计算Q：Y＝QQ^TY；

S6：B←Q^TWQ；

S7：通过SVD分解来计算V：B＝VΛV^T；

S8：

S9：输出U和Λ，并将λ_l←Λ；

(g)从当前第l层的层次邻域谱特征λ_l中选择该层的主谱特征

具体步骤主要包括：设置r_l＝2(1≤l≤k)，从当前第l层的层次邻域谱特征λ_l中选择选择前r_l项构成该层的主谱特征λ'_l，通过选择主谱特征λ'_l，既可以保留层次邻域谱特征λ_l的核心特征信息，同时可以减少特征匹配时的时间代价；

(h)连接当前中心像素p的每层对应的主谱特征

来生成当前中心像素p的最终层次邻域谱特征χ_R(p)或χ_M(p)，具体步骤主要包括：顺序组合每层的主谱特征λ'_l(l＝1,2,…,k)，从而生成中心像素点p的最终层次邻域谱特征χ(p)＝(λ'₁/||λ'₁||,λ'₂/||λ'₂||,…,λ'_k/||λ'_k||)，然后转至步骤(a)更新中心像素p；

(i)根据参考图像R对浮动图像

的层次邻域谱特征χ_R(p)或χ_M(p)进行排序；

(j)将参考图像R和浮动图像

的层次邻域谱特征χ_R(p)和χ_M(p)分别嵌入到各自的特征空间R＝(I_R,X_R,χ_R)和M＝(I_M,X_M,χ_M)中；

(k)寻找映射u_R→M实现从参考图像R到浮动图像

的映射；

步骤(8)：利用层次邻域谱匹配算法计算从M到

的更新场u_M→R；其中，层次邻域谱匹配算法的具体步骤主要包括：

(a)对参考图像M和浮动图像

分别逐像素遍历，并判断是否遍历结束，如果没有结束，设定当前中心像素p并顺序执行步骤(b)，如果当前图像已遍历结束，则跳转至步骤(i)；

(b)初始化层次邻域层数k＝10，当前层数l＝1；

(c)判断l≤k是否成立，如果成立，顺序执行步骤(d)，针对当前中心像素p开始进行层次邻域谱特征λ_l计算，否则转至步骤(h)；

(d)针对当前中心像素p，构造当前层对应的线图L(S)，具体步骤主要包括：

首先，构建星图S，针对大小为m×n的图像I，任何非边界的中心像素p的特征可以用它与邻域像素的关系来定义，因此，中心像素p的特征用图来表示，假设p的第l(1≤l≤k,i∈R⁺)层的8l邻域像素为q_j(j＝1,2,…,8l)，用中心像素p和8l邻域像素q_j(j＝1,2,…,8l)构建星图S。在星图S中，中心像素p与所有邻域像素q_j(j＝1,2,…,8l)之间的边为pq_j，每条边的权重r_j为：

其中，I(p)和X(p)分别是中心像素点p的灰度值和空间坐标，I(q_j)和X(q_j)分别是像素点q_j的灰度值和空间坐标；

其次，为进一步反映中心像素点p的特征，在星图S中，用边pq_j的权重r_j相减来构造线图L(S)；

(e)生成线图L(S)对应的邻接矩阵A，其中邻接矩阵A的元素为：

其中，1≤i,j≤8l,i,j∈R⁺。由于邻接矩阵A是一个实对称矩阵，8l个实数特征可以通过谱分解来获得，并且λ₁(p)≥λ₂(p)≥…≥λ_8l(p)，因此使用λ_i(p)来代替像素点p的灰度值并正则化到[0，255]，可以计算大小为(m-1)×(n-1)的特征图I_i(i＝1,2,…,8l)；

(f)对邻接矩阵A进行谱分解，从而计算当前第l层的层次邻域谱特征

具体步骤包括：

S1：初始化采样列数m，秩n和过采样参数o；

S2：从矩阵A中随机采样m列和m行生成矩阵W；

S3：生成高斯随机矩阵Ω_m×(n+o)；

S4：Y←WΩ；

S5：通过QR分解来计算Q：Y＝QQ^TY；

S6：B←Q^TWQ；

S7：通过SVD分解来计算V：B＝VΛV^T；

S8：

S9：输出U和Λ，并将λ_l←Λ；

(g)从当前第l层的层次邻域谱特征λ_l中选择该层的主谱特征

具体步骤主要包括：设置r_l＝2(1≤l≤k)，从当前第l层的层次邻域谱特征λ_t中选择选择前r_k项构成该层主谱特征λ'_t，通过选择主谱特征λ'_t，既可以保留层次邻域谱特征λ_l的核心特征信息，同时可以减少特征匹配时的时间代价；

(h)连接当前中心像素p的每层对应的主谱特征来生成当前中心像素p的最终层次邻域谱特征χ_R(p)或χ_M(p)，具体步骤主要包括：顺序组合每层的主谱特征λ'_l(l＝1,2,…,k)，从而生成中心像素点p的最终层次邻域谱特征χ(p)＝(λ'₁/||λ'₁||,λ'₂/||λ'₂||,…,λ'_k/||λ'_k||)，然后转至步骤(a)重新确定中心像素；

(i)根据参考图像M对浮动图像

的层次邻域谱特征进行排序；

(j)将参考图像M和浮动图像

的层次邻域谱特征分别嵌入到各自的特征空间R＝(I_R,X_R,χ_R)和M＝(I_M,X_M,χ_M)中；

(k)寻找映射u_M→R实现从参考图像M到浮动图像

的映射；

步骤(9)：计算平均更新场u←(u_R→M-u_M→R)/2；

步骤(10)：平滑平均更新场u←K_fluid*u；

步骤(11)：利用平滑后的更新场u来计算速度场

；

步骤(12)：平滑速度场v_i←K_diffuse*v_i，转至步骤(6)执行；

步骤(13)：输出当前层的最优配准变换φ_i ^*，然后i＝i+1并转至步骤(4)执行；

步骤(14)：输出最终的最优变换φ^*。

为了验证本发明的图像配准效果，在两类图像上完成实验：医学图像和真实图像，其中医学图像分别从Vanderbilt大学的RIRE数据库中选择和从医院收集。在这些图像上人工生成复杂形变作为待配准图像。用于实验的测试图像如下表所示。为了展示本发明在不同模态和不同身体部位的泛化能力，测试数据包括不同图像模态：CT、MRI-T1、MRI-T2和MRI-T2-retified，并包括不同身体部位：头部、胸部和腹部。所有算法和实验评估在配置Intel(R)Core(TM)CPU i7-9700,4.2GHz×8和32GB RAM的PC机上执行。

为了评估本发明方法的有效性，将本发明与六个非刚性配准算法进行比较，六个非刚性配准算法分别是：基于B样条的FFD、Demons、Active Demons、改进Active Demons、Log-Demons和Spectral Demons。为了公平比较，基于B样条的FFD算法选择使其性能最好的配准参数。所有算法都采用三层分辨率，从高分辨率到低分辨率每层的最大迭代次数分别为100、50和20。为了平衡本发明的效率和精确性，层次邻域谱特征仅仅在粗分辨率层上进行计算，并且设置层次邻域层数k＝10来捕获复杂形变。另外，在选择主谱特征时设置主谱特征个数r_k＝2来降低相似度计算的复杂度。所有Demons算法采用相同的参数：σ_fluid＝1，σ_diffuse＝1和α_x＝1。由于空域一致性权重α_s惩罚大形变，α_s较大时会降低配准精确性。因此，α_s应该比较小但是非0，参数的最佳实验值应该满足条件1/10＜α_s/α_g＜1/3。在实际中，Spectral Demons常常选择参数α_i＝0.8，α_s＝0.05和α_g＝0.15来支持几何和灰度一致性，防止每次迭代中发生不合理形变。为了保证公平，本发明选择相同的权重参数α_i＝0.8，α_s＝0.05和α_g＝0.15。

通过主观可视化和客观量化两种方式评估本发明的配准结果。使用五种客观评价指标：归一化互信息NMI、差值平方和SSD、差值绝对值和SAD、归一化互相关NCC和Dice相似性系数DSC。归一化互信息NMI描述参考图像和浮动图像的统计相关性，定义如下所示：

SSD可以在灰度匹配和灰度形变之间获得较好的平衡，定义如下所示：

SAD在评估配准性能方面与SSD相似，定义如下所示：

NCC描述参考图像和浮动图像的相关程度，定义如下所示：

DSC描述参考图像和浮动图像的重叠程度，定义如下所示：

针对这些测度，好的配准结果对应较大的NMI、NCC和DSC值，对应较小的SSD和SAD值。

使用以上七种算法对医学图像进行非刚性配准实验。在实验之前，首先对医学图像进行刚性配准，然后进行归一化处理来避免由于不同设备定位产生的数据不一致性的影响。配准实验使用基于B样条的FFD算法、Demons、Active Demons、改进Active Demons、Log-Demons、Spectral Demons和本发明。图3显示了使用以上七种算法对六组不同模态和身体部位的医学图像进行配准的实验结果图，其中，(a)为参考图像；(b)为浮动图像；(c)-(i)分别为采用基于B样条的FFD、Demons、Active Demons、改进Active Demons、Log-Demons、Spectral Demons和本发明对待配准图像进行配准的结果。

在图3中，由于在浮动图像中存在复杂形变，基于B样条的FFD算法、Demons、ActiveDemons、改进Active Demons和Log-Demons容易使配准过程陷入局部最优，影响配准精确性。因此，这些算法产生的配准结果具有较低的精确性。与基于B样条的FFD算法、Demons、Active Demons、改进Active Demons和Log-Demons相比，Spectral Demons产生了更好的配准结果。然而，这些结果仍旧包含明显的误配准区域。与其它配准算法相比，使用本发明进行配准得到的差分图像没有明显的配准误差，可以获得最好的平滑性和最小的配准误差。

另外，量化评价本发明在医学图像配准时的配准精确性。以上七种算法分别配准六组待配准的医学图像，针对配准结果统计五种客观评价指标如下表所示。从下表可以看出，根据五种测度，本发明可以获得最好的配准结果，这与主观评价方式得到的结果一致。因此，本实验证明针对不同模态和不同身体部位的医学图像进行配准，本发明比其它六种算法可以获得更好的配准精确性。

/>

为了进一步评估本发明的配准精确性和鲁棒性，在从公开数据集上选择的四组真实图像上进行配准实验(浮动图像通过在参考图像上人工应用形变来生成)，实验结果如图4所示。图4显示了使用以上七种算法对四组真实图像进行配准的实验结果图，其中，(a)为参考图像；(b)为浮动图像；(c)-(i)分别为采用基于B样条的FFD、Demons、Active Demons、改进Active Demons、Log-Demons、Spectral Demons和本发明对待配准图像进行配准的结果。

从图4中看出，与Spectral Demons和本发明相比，Demons、Active Demons、改进Active Demons、Log-Demons和基于B样条的FFD算法产生了较差的配准结果。一个重要原因是这些算法仅使用图像的灰度和梯度信息来产生驱动力，并使用高斯滤波器对位移形变场进行更新。因此，由于不充分的驱动力导致一些误配准区域出现在图4中(c)-(g)中。Spectral Demons使用谱特征可以捕获浮动图像的全局几何相似性，能够处理全局大的非刚性形变，但是不能很好地改善局部小形变。本发明可以更好地处理浮动图像中的复杂形变。同时，由于层次邻域谱特征的使用，本发明能有效处理局部非刚性形变并改善配准精确性。

更进一步，下表中的客观量化评价结果与图4中的主观可视化评价结果一致。从下表中五种评价指标可以看出，由于集成了层次邻域谱特征到微分同胚Demons配准框架中，本发明针对真实图像进行配准时优于其它六种算法。

通过比较本发明和Log-Demons、Spectral Demons来分析和验证本发明的收敛性。分别针对一组医学图像和一组真实图像执行非刚性配准。为了显示能量函数随迭代次数变化的趋势，归一化每次迭代的能量函数值。由于需要分析收敛性，这三种算法都使用了相同的三层分辨率策略。在第三层分辨率图像上，计算能量函数随迭代次数的变化值。能量函数值的变化曲线如图5所示。

图5展示了本发明在12-13次迭代后开始接近最优配准结果，因此，本发明具有更快的收敛性。本发明提出一种新的谱匹配相似性测度并把它集成到微分同胚框架的能量函数中，可以产生有先验信息的形变场，这些先验信息具有旋转和亮度不变性和对噪声的鲁棒性。这些先验信息用来引导整个配准过程，可以克服驱动力不足的问题，加速能量函数的收敛并改善配准性能。

进一步验证本发明的配准效率。在配置为Intel(R)Core(TM)CPU i7-9700，4.2GHz×8的PC机上执行实验。在以上10组图像上执行七种配准算法。所有算法的参数设置和上面实验一致。为了消除系统运行误差，所有算法针对每组待配准图像执行10次，然后比较平均运行时间。例如：Demons、Active Demons和改进Active Demons在尺寸为256×256的Lena图像上进行配准，仅需要不到30s，具有较高的执行效率。但是，基于B样条的FFD算法和Log-Demons分别需要大概100s。然而，Spectral Demons需要大概330s，本发明需要大概270s。

由于在Demons、Active Demons和改进Active Demons的能量函数中仅仅使用了灰度信息和梯度信息，这三个算法需要相对较短的执行时间。Log-Demons为了获得可逆的形变场，因此需要比Demons算法更多的时间。基于B样条的FFD算法为了处理复杂形变需要用到很多参数，所以需要比Demons、Active Demons和改进Active Demons消耗更多的执行时间。Spectral Demons和本发明引入谱特征到配准过程中，在每次迭代中计算谱相似度。因此，Spectral Demons和本发明比其它配准算法需要更多的执行时间。然而，SpectralDemons和本发明引入谱信息到能量函数中，可以使配准过程有更强的驱动力并获得更快的收敛性。根据收敛性分析，本发明在12-13次迭代之后近似最好的配准结果，因此本发明可以通过减少迭代次数获得和Log-Demons相同的结果，并且本发明能够提供更好的配准精确性。与Spectral Demons相比，本发明通过利用基于随机SVD的

近似方法来近似高维邻接矩阵，可以有效降低谱分解的计算复杂度。因此，本发明比Spectral Demons获得更好的配准性能。/>

Claims (9)

1.一种基于层次邻域谱特征的微分同胚非刚性配准算法，其特征在于，包括：

(1)对参考图像R和浮动图像M分别生成参考图像金字塔R_i(i＝1，2，...，N)和浮动图像金字塔M_i(i＝1，2，...，N)；

(2)初始化应用于第一层浮动图像M₁上的速度场v₁；

(3)设置当前待配准的图像层数i＝1，从第一层开始进行配准；

(4)判断i≤N是否成立，如果是则继续步骤(5)，否则转至步骤(14)；

(5)如果i>1，初始化当前层的速度场v_i；

(6)判断配准条件是否收敛，如果是则转至步骤(13)，否则继续步骤(7)；

(7)利用层次邻域谱匹配算法计算从R到

的更新场u_R→M；

(8)利用层次邻域谱匹配算法计算从M到

的更新场u_M→R；

(9)计算平均更新场u←(u_R→M-u_M→R)/2；

(10)平滑平均更新场u←K_fluid*u；

(11)利用平滑后的更新场u来计算速度场

(12)平滑速度场v_i←K_diffuse*v_i，转至步骤(6)执行；

(13)输出当前层的最优配准变换

然后i＝i+1并转至步骤(4)执行；

(14)输出最终的最优变换φ^*；

所述步骤(7)和步骤(8)中层次邻域谱匹配算法，包括：

(a)对参考图像R和浮动图像M分别逐像素遍历，并判断是否遍历结束，如果没有结束，设定当前中心像素p并顺序执行步骤(b)，如果当前图像已遍历结束，则跳转至步骤(i)；

(b)初始化层次邻域层数k＝10，当前层数l＝1；

(c)判断l≤k是否成立，如果成立，顺序执行步骤(d)，针对当前中心像素p开始进行层次邻域谱特征λ_l计算，否则转至步骤(h)；

(d)针对当前中心像素p，构造当前层对应的线图L(S)；

(e)生成线图L(S)对应的邻接矩阵A；

(f)对邻接矩阵A进行谱分解，从而计算当前第l层的层次邻域谱特征

(g)从当前第l层的层次邻域谱特征λ_l中选择该层的主谱特征

(h)连接当前中心像素p的每层对应的主谱特征

来生成当前中心像素p的最终层次邻域谱特征χ_R(p)或χ_M(p)，转至步骤(a)更新中心像素p；

(i)根据参考图像R对浮动图像M的层次邻域谱特征χ_R(p)或χ_M(p)进行排序；

(j)将参考图像R和浮动图像M的层次邻域谱特征χ_R(p)和χ_M(p)分别嵌入到各自的特征空间R＝(I_R,X_R,χ_R)和M＝(I_M,X_M,χ_M)中；

(k)寻找映射c实现从参考图像R到浮动图像M的映射

2.根据权利要求1所述的一种基于层次邻域谱特征的微分同胚非刚性配准算法，其特征在于，所述步骤(1)中采用基于小波分解的多分辨率策略分别对参考图像R和浮动图像M进行多分辨率分解，分别生成参考图像金字塔R_i(i＝1，2，…，N)和浮动图像金字塔M_i(i＝1，2，…，N)。

3.根据权利要求1所述的一种基于层次邻域谱特征的微分同胚非刚性配准算法，其特征在于，所述步骤(d)中构造当前层对应的线图L(S)，包括：

首先构建星图S；然后在星图S中，用边pq_j的权重r_j相减来构造线图L(S)。

4.根据权利要求3所述的一种基于层次邻域谱特征的微分同胚非刚性配准算法，其特征在于，所述星图S构建时，中心像素p的第l(1≤l≤k,l∈R⁺)层的8l邻域像素为q_j(j＝1,2,…,8l)，利用中心像素p和8l邻域像素q_j(j＝1,2,…,8l)构建星图S，在星图S中，中心像素p与所有邻域像素q_j(j＝1,2,…,8l)之间的边为pq_j，每条边的权重r_j为：

其中，I(p)和X(p)分别是中心像素点p的灰度值和空间坐标，I(q_j)和X(q_j)分别是像素点q_j的灰度值和空间坐标。

5.根据权利要求1所述的一种基于层次邻域谱特征的微分同胚非刚性配准算法，其特征在于，所述步骤(e)生成线图L(S)对应的邻接矩阵A，其邻接矩阵A的元素为：

其中，1≤i,j≤8l,i,j∈R⁺。

6.根据权利要求5所述的一种基于层次邻域谱特征的微分同胚非刚性配准算法，其特征在于，所述邻接矩阵A的8l个实数特征通过谱分解来获得，且λ₁(p)≥λ₂(p)≥…≥λ_8l(p)，使用λ_i(p)来代替像素点p的灰度值并正则化到[0，255]，计算大小为(m-1)×(n-1)的特征图I_i(i＝1,2,…,8l)。

7.根据权利要求1所述的一种基于层次邻域谱特征的微分同胚非刚性配准算法，其特征在于，所述步骤(f)对邻接矩阵A进行谱分解，从而计算当前第l层的层次邻域谱特征

包括：

S1：初始化采样列数m，秩n和过采样参数o；

S2：从矩阵A中随机采样m列和m行生成矩阵W；

S3：生成高斯随机矩阵Ω_m×(n+o)；

S4：Y←WΩ；

S5：通过QR分解来计算Q：Y＝QQ^TY；

S6：B←Q^TWQ；

S7：通过SVD分解来计算V：B＝VΛV ^T；

S8：

S9：输出U和Λ，并将λ_l←Λ。

8.根据权利要求1所述的一种基于层次邻域谱特征的微分同胚非刚性配准算法，其特征在于，所述步骤(g)从当前第l层的层次邻域谱特征λ_l中选择当前层的主谱特征

包括：设置r_l＝2(1≤l≤k)，从当前第l层的层次邻域谱特征λ_l中选择前r_l项构成该层的主谱特征/>

9.根据权利要求1所述的一种基于层次邻域谱特征的微分同胚非刚性配准算法，其特征在于，所述步骤(h)包括：顺序组合每层的主谱特征λ'_l(l＝1,2,…,k)，从而生成中心像素点p的最终层次邻域谱特征χ(p)＝(λ'₁/||λ'₁||,λ'₂/||λ'₂||,…,λ'_k/||λ'_k||)。

Images