CN103927725A - 基于分数阶微分的电影核磁共振图像序列运动场估计方法 - Google Patents

Abstract

基于分数阶微分的电影核磁共振图像序列运动场估计方法，属于核磁共振成像数据处理领域，本发明的目的是为了解决应用现有技术利用整数阶微分图像增强和直接建立光流方程的方法估计电影核磁共振图像序列运动场，存在丢失图像的纹理细节，估计结果受光照变化影响，不具有旋转不变性，抗噪声性能差，导致基于电影核磁共振图像运动估计的精确度低的问题。实现该方法的主要步骤为：一、对电影核磁共振图像利用分数阶微分进行纹理增强；二、通过Riesz变换提取图像的单演信号，即单演相位、单演方位、单演振幅；三、利用单演信号的相位向量建立光流方程；四、通过光流方程估计电影核磁共振图像序列的运动场。本发明用于电影核磁共振图像估计成像对象的运动。

Description

基于分数阶微分的电影核磁共振图像序列运动场估计方法

技术领域

本发明属于核磁共振成像数据处理领域。

背景技术

核磁共振成像(Magnetic Resonance Imaging,MRI)技术已经成为临床诊断的重要辅助手段。基于电影核磁共振图像(Cine-MRI)运动场估计是一个不可忽视的研究方向。由于Cine-MRI的灰度非常相近，难以找到密集的对应特征点，因而对于Cine-MRI的运动估计比对于加标记核磁共振图像(Tagged-MRI)的运动估计更具有难度，目前国内外研究较少。针对Cine-MRI的灰度相近，特征稀少的挑战，微分运算增强图像已得到广泛应用，应用整数阶微分图像增强方法(如Sobel、Prewitt、Laplacian算子)在增强图像的同时丢失了图像的纹理细节，尤其是平滑区域的纹理在整数阶微分后损失严重；而分数阶微分可以在提高图像信号高频信息的同时非线性地保留低频信息(纹理细节)。近年来，分数阶微积分逐渐应用于图像处理的各个方面。

单演信号是一种二维解析信号，是由Felsberg提出的，它是一维解析信号的一个扩展，原始的二维图像信息通过具有旋转不变性的Riesz变换，以一种非线性的方式被映射到了虚平面。与一维解析信号类似，该模型涵盖了原始图像在实平面和虚平面上的不同信息，可以据此提取出相应的单演幅度，单演相位和单演方向特征信息。与其他提取图像相位信息的方法如基于Gabor滤波器的相位计算方法相比，该模型的相位计算不需要对方向进行采样。相对于边缘轮廓信息，其中的单演相位特征包含了该图像的绝大部分重要特征信息，只需要这个单演相位向量信息就能建立光流方程并估计运动场。单演相位特征最大的优点是独立于亮度信息，不受光照变化的影响，而且该特征量具有旋转不变性同时也具有很好的抗噪声性能。

应用现有技术整数阶微分图像增强方法和直接建立光流方程的方法，估计电影核磁共振图像序列运动场，存在丢失图像的纹理细节，结果受光照变化影响，不具有旋转不变性，抗噪声性能差，基于电影核磁共振图像运动估计的精确度低。

发明内容

发明目的：为了解决应用现有技术利用整数阶微分图像增强和直接建立光流方程的方法估计电影核磁共振图像序列运动场，存在丢失图像的纹理细节，估计结果受光照变化影响，不具有旋转不变性，抗噪声性能差，导致基于电影核磁共振图像运动估计的精确度低的问题，本发明提出了一种基于分数阶微分的电影核磁共振图像序列运动场估计方法。

技术方案：本发明基于分数阶微分的电影核磁共振图像序列运动场估计方法，实现该方法的步骤如下：

一、利用分数阶微分对电影核磁共振图像序列进行纹理增强；

二、通过Riesz变换提取图像的单演信号，即单演相位、单演方位、单演振幅；

三、利用单演信号的相位向量建立光流方程；

四、通过光流方程估计电影核磁共振图像序列的运动场。

有益效果：

(1)利用分数阶微分对电影核磁共振图像进行纹理增强，克服了现有技术应用整数阶微分图像增强方法在增强图像的同时产生的丢失图像纹理细节的缺陷，可以在提高图像信号高频信息的同时非线性地保留低频信息，低频信息即图像纹理细节。

(2)利用单演信息构造光流方程，单演相位信息独立于亮度信息，不受光照变化的影响。

(3)本方法利用的单演信息特征量具有旋转不变性，具有很好的抗噪声性能，有效提高电影核磁共振图像序列的运动场的估计精度，比不利用单演信息构造光流方程估计精度高10％以上。

本发明用于电影核磁共振图像估计成像对象的运动。

附图说明

图1基于分数阶微分的电影核磁共振图像序列运动场估计方法的流程图；

图2任意一点(x,y)的8个对称方向，分别为x轴正方向，x轴负方向，y轴正方向，y轴负方向,左下方向，左上方向，右上方向，右下方向；其中水平向右为x轴正方向，水平向左为为x轴负方向，竖直向上为y轴正方向，竖直向下为y轴负方向；

图3分数阶微分增强掩模构造示意图，以兴趣点为中心，分别向8个对称方向构造掩模

图4单演信息提取示意图，p，q₁，q₂构成正交三维向量空间，Α＝(p,q₁,q₂)^T，q＝(q₁,q₂)^T；A与p的夹角为q与q₁的夹角为θ。

具体实施方式

具体实施方式一：结合附图1，本实施方式基于分数阶微分的电影核磁共振图像序列运动场估计方法的具体实施步骤如下：

一、利用分数阶微分对电影核磁共振图像序列进行纹理增强；

二、通过Riesz变换提取图像的单演相位，单演方位，单演振幅，构造单演信号；

三、利用单演信号的相位向量建立光流方程；

四、通过光流方程估计电影核磁共振图像序列的运动场。

具体实施方式二：本实施方式是对具体实施方式一中的基于分数阶微分的电影核磁共振图像序列运动场估计方法的具体实施步骤一的进一步限定，结合附图2、附图3，步骤一所述利用分数阶微分对电影核磁共振图像序列进行纹理增强的过程为：

从连续函数整数阶导数的出发，将微分的阶数由整数扩展到分数，构建基本的v阶Grümwald–Letnikov分数阶微分方程，

$D_{G-L}^{v}s\left(u\right)=\frac{d^v}{[d(u-a){]}^v}s\left(u\right){|}_{G-L}=\underset{N\→\∞}{\lim }\{\frac{{\left(\frac{u-a}{N}\right)}^{-v}}{\mathrm{\Γ}(-v)}\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}\frac{\mathrm{\Γ}(k-v)}{\mathrm{\Γ}(k+1)}\×s(u-k\frac{u-a}{N})\}---\left(1\right)$

这里连续函数s(u)表示一维图像信号，s(u)∈[a,u],a<u,a∈R,u∈R，u表示一维信号变化量，R表示实数集，[v]表示v的整数部分，表示整数集合；当v>0时，k不小于[v]，表示信号的长度；是Gamma函数，表示Grümwald–Letnikov-分数阶微分操作符；当N值足够大时，推导出一维图像信号分数阶微分表达式，

$\frac{d^v}{d{x}^v}s\left(u\right)\&cong;\frac{x^{-v}{N}^v}{\mathrm{\&Gamma;}(-v)}\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\frac{\mathrm{\&Gamma;}(k-v)}{\mathrm{\&Gamma;}(k+1)}s(u-\frac{\mathrm{ku}}{N})=\frac{u^{-v}{N}^v}{\mathrm{\&Gamma;}(-v)}\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\frac{\mathrm{\&Gamma;}(k-v)}{\mathrm{\&Gamma;}(k+1)}\&times;s_k---\left(2\right)$

针对二维图像s(x,y)，有以下两个表达式，(x,y)表示图像中一点的坐标，

$\begin{array}{c}\frac{{\&PartialD;}^{v}s(x,y)}{\&PartialD;x^v}\&cong;s(x,y)+(-v)s(x-1,y)+\frac{(-v)(-v+1)}{2}s(x-2,y)+\\ \frac{(-v)(-v+1)(-v+2)}{g}s(x-3,y)+\&CenterDot;\&CenterDot;\&CenterDot;\frac{\mathrm{\&Gamma;}(n-v-1)}{(n-1)!\mathrm{\&Gamma;}(-v)}s(x-n+1,y)\end{array}---\left(3\right)$

$\begin{array}{c}\frac{{\&PartialD;}^{v}s(x,y)}{\&PartialD;x^v}\&cong;s(x,y)+(-v)s(x,y-1)+\frac{(-v)(-v+1)}{2}s(x,y-2)+\\ \frac{(-v)(-v+1)(-v+2)}{g}s(x,y-3)+\&CenterDot;\&CenterDot;\&CenterDot;\frac{\mathrm{\&Gamma;}(n-v-1)}{(n-1)!\mathrm{\&Gamma;}(-v)}s(x,y-n+1)\end{array}---\left(4\right)$

为了使分数阶微分掩模具有45°旋转不变性，针对二维图像中的任意一点(x,y)设计在8个对称方向下的分数阶微分，分别是X轴正方向，X轴负方向，Y轴正方向，Y轴负方向,左下方向，左上方向，右上方向，右下方向；

根据公式(2)和公式(3)构造分数阶微分增强掩模公式如下：

$\left\{\begin{array}{c}{C}_{s_0}=1\\ C_{s_1}=-v\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ C_{s_k}=\frac{\mathrm{\&Gamma;}(k-v)}{k!\mathrm{\&Gamma;}(-v)}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ C_{s_{n-1}}=\frac{\mathrm{\&Gamma;}(n-v-1)}{(n-1)!\mathrm{\&Gamma;}(-v)}\\ C_{s_n}=\frac{\mathrm{\&Gamma;}(n-v)}{n!\mathrm{\&Gamma;}(-v)}\end{array}\right.---\left(5\right)$

这里公式(5)中各式是兴趣点掩模公式，s₀＝s(x₀,y₀)，(x₀,y₀)表示兴趣点坐标，用掩模依次卷积需增强图像的每个像素点。

具体实施方式三：本实施方式是对具体实施方式一中的基于分数阶微分的电影核磁共振图像序列运动场估计方法的具体实施步骤二的进一步限定，结合附图4，步骤二所述通过Riesz变换提取图像的单演相位，单演方位，单演振幅，构造单演信号的过程为：

a、通过3个二维空间正交滤波器将图像对应点的局部振幅、局部相位和局部方向分离出来；所采用的二维空间正交滤波器为差分泊松滤波器；空间正交滤波器由1个旋转不变的偶数带通滤波器b_e(Z)和2个奇数带通滤波器b_o1(Z)、b_o2(Z)组成；

b、求单演相位单演方向θ(Z)和单演振幅A(Z)：

$\mathrm{\&theta;}\left(Z\right)=\arctan \left(\frac{q_2\left(Z\right)}{q_1\left(Z\right)}\right)$

$A\left(Z\right)=\sqrt{p^2\left(Z\right)+|q\left(Z\right){|}^2}$

这里Z＝(x,y),p(Z)＝(I*b_e)(Z)，q₁(Z)＝(I*b_o1)(Z)，q₂(Z)＝(I*b_o2)(Z)，q(Z)＝[q₁(Z),q₂(Z)]^T，“*”，符号表示二维卷积，I表示对应Z点灰度。

具体实施方式四：本实施方式是对具体实施方式一中的基于分数阶微分的电影核磁共振图像序列运动场估计方法的具体实施步骤三的进一步限定，步骤三所述利用单演信号的相位向量建立光流方程的方法为：

基于局部约束，假定以p点为中心的一个小区域Ω内各点的光流相同，对区域内不同的点给予不同的权重，离p点越近，权重越高；

列出光流的计算方程： $\underset{x\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{W}^2\left(p\right)[\&dtri;I(p,t)\&CenterDot;v+I_t(p,t){]}^2=0---\left(7\right)$

上式中，Ω代表2维图像中以p点为中心的一个小的区域，W(p)为窗函数，代表区域中各点的权重，离p点越近，权重越高，表示t时刻p点灰度的梯度，v[U V]^T表示运动向量，U表示在X轴方向位移，V表示在Y轴方向位移，I_t(p,t)表示图像中p点在t时刻的灰度；

对于邻域Ω内的n点，其中，

$A={(\&dtri;I\left(p_1\right),\&dtri;I\left(p_2\right),...,\&dtri;I\left(p_n\right))}^T,$

W＝diag(W(p₁),W(p₂),...,W(p_n))，diag表示主对角矩阵是一个主对角线之外的元素皆为0的矩阵，

b＝-(I_t(p₁),I_t(p₂),...,I_t(p_n))^T。

具体实施方式五：本实施方式是对具体实施方式四中的基于分数阶微分的电影核磁共振图像序列运动场估计方法的具体实施步骤四的进一步限定，步骤四通过光流方程估计电影核磁共振图像序列的运动场的方法为：

由光流计算方程： $\underset{x\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{W}^2\left(p\right)[\&dtri;I(p,t)\&CenterDot;v+I_t(p,t){]}^2=0,$

解得运动向量v＝(A^TW2A)^-1A^TW²b；

其中A^TW²A为2×2矩阵：

$A^T{W}^{2}A=[\left.\begin{array}{cc}\mathrm{\&Sigma;}{W}^2\left(p\right)I_x^2& \mathrm{\&Sigma;}{W}^2\left(p\right)I_x{I}_y\\ \mathrm{\&Sigma;}{W}^2\left(p\right)I_y{I}_x& \mathrm{\&Sigma;}{W}^2\left(p\right)I_y^2\end{array}\right],$

上式中所有的求和都是在邻域Ω内所有点上进行的；

Ix和I_y表示灰度在X轴和Y轴方向的变化率，本方法用单演信息相位沿方向n的变化率代替Ix和I_y；

代入单演信息求出单演相位信息下的运动向量v，对电影核磁共振图像进行区域划分，依次求出各个区域内单演相位信息下的运动向量v，根据各个运动向量构造出电影核磁共振图像的运动场。

Claims (5)

1.基于分数阶微分的电影核磁共振图像序列运动场估计方法，其特征在于，实现该方法的步骤如下：

一、利用分数阶微分对电影核磁共振图像序列进行纹理增强；

二、通过Riesz变换提取图像的单演信号，即单演相位、单演方位、单演振幅；

三、利用单演信号的相位向量建立光流方程；

四、通过光流方程估计电影核磁共振图像序列的运动场。

2.根据权利要求1所述的基于分数阶微分的电影核磁共振图像序列运动场估计方法，其特征在于，

步骤一所述利用分数阶微分对电影核磁共振图像进行纹理增强的过程为：

从连续函数整数阶导数的出发，将微积分的阶数由整数扩展到分数，构建基本的v阶Grümwald–Letnikov分数阶微分方程，

$D_{G-L}^{v}s\left(u\right)=\frac{d^v}{[d(u-a){]}^v}s\left(u\right){|}_{G-L}=\underset{N\&RightArrow;\&infin;}{\lim }\{\frac{{\left(\frac{u-a}{N}\right)}^{-v}}{\mathrm{\&Gamma;}(-v)}\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\frac{\mathrm{\&Gamma;}(k-v)}{\mathrm{\&Gamma;}(k+1)}\&times;s(u-k\frac{u-a}{N})\}---\left(1\right)$

这里连续函数s(u)表示一维图像信号，s(u)∈[a,u],a<u,a∈R,u∈R，u表示一维信号变化量，R表示实数集，[v]表示v的整数部分，表示整数集合；当v>0时，k不小于[v]，表示信号的长度；是Gamma函数，表示Grümwald–Letnikov-分数阶微分操作符；当N值足够大时，推导出一维图像信号分数阶微分表达式，

$\frac{d^v}{d{x}^v}s\left(u\right)\&cong;\frac{x^{-v}{N}^v}{\mathrm{\&Gamma;}(-v)}\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\frac{\mathrm{\&Gamma;}(k-v)}{\mathrm{\&Gamma;}(k+1)}s(u-\frac{\mathrm{ku}}{N})=\frac{u^{-v}{N}^v}{\mathrm{\&Gamma;}(-v)}\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\frac{\mathrm{\&Gamma;}(k-v)}{\mathrm{\&Gamma;}(k+1)}\&times;s_k---\left(2\right)$

针对二维图像s(x,y)，有以下两个表达式，(x,y)表示图像中一点的坐标，

$\begin{array}{c}\frac{{\&PartialD;}^{v}s(x,y)}{\&PartialD;x^v}\&cong;s(x,y)+(-v)s(x-1,y)+\frac{(-v)(-v+1)}{2}s(x-2,y)+\\ \frac{(-v)(-v+1)(-v+2)}{g}s(x-3,y)+\&CenterDot;\&CenterDot;\&CenterDot;\frac{\mathrm{\&Gamma;}(n-v-1)}{(n-1)!\mathrm{\&Gamma;}(-v)}s(x-n+1,y)\end{array}---\left(3\right)$

$\begin{array}{c}\frac{{\&PartialD;}^{v}s(x,y)}{\&PartialD;x^v}\&cong;s(x,y)+(-v)s(x,y-1)+\frac{(-v)(-v+1)}{2}s(x,y-2)+\\ \frac{(-v)(-v+1)(-v+2)}{g}s(x,y-3)+\&CenterDot;\&CenterDot;\&CenterDot;\frac{\mathrm{\&Gamma;}(n-v-1)}{(n-1)!\mathrm{\&Gamma;}(-v)}s(x,y-n+1)\end{array}---\left(4\right)$

针对二维图像中的任意一点(x,y)设计在8个对称方向下的分数阶微分，分别是X轴正方向，X轴负方向，Y轴正方向，Y轴负方向,左下方向，左上方向，右上方向，右下方向；

根据公式(2)和公式(3)构造分数阶微分增强掩模公式如下：

$\left\{\begin{array}{c}{C}_{s_0}=1\\ C_{s_1}=-v\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ C_{s_k}=\frac{\mathrm{\&Gamma;}(k-v)}{k!\mathrm{\&Gamma;}(-v)}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ C_{s_{n-1}}=\frac{\mathrm{\&Gamma;}(n-v-1)}{(n-1)!\mathrm{\&Gamma;}(-v)}\\ C_{s_n}=\frac{\mathrm{\&Gamma;}(n-v)}{n!\mathrm{\&Gamma;}(-v)}\end{array}\right.---\left(5\right)$

这里公式(5)中各式是兴趣点掩模公式，s₀＝s(x₀,y₀)，(x₀,y₀)表示兴趣点坐标，用掩模依次卷积需增强图像的每个像素点。

3.根据权利要求1或2所述的基于分数阶微分的电影核磁共振图像序列运动场估计方法，其特征在于，步骤二通过Riesz变换提取图像的单演信号，即单演相位、单演方位、单演振幅的过程为：

a、通过3个二维空间正交滤波器将图像对应点的局部振幅、局部相位和局部方向分离出来；所采用的二维空间正交滤波器为差分泊松滤波器；空间正交滤波器由1个旋转不变的偶数带通滤波器b_e(Z)和2个奇数带通滤波器b_o1(Z)、b_o2(Z)组成；

b、求单演相位单演方向θ(Z)和单演振幅A(Z)：

$\mathrm{\&theta;}\left(Z\right)=\arctan \left(\frac{q_2\left(Z\right)}{q_1\left(Z\right)}\right)$

$A\left(Z\right)=\sqrt{p^2\left(Z\right)+|q\left(Z\right){|}^2}$

这里Z＝(x,y),p(Z)＝(I*b_e)(Z)，q₁(Z)＝(I*b_o1)(Z)，q₂(Z)＝(I*b_o2)(Z)，q(Z)＝[q₁(Z),q₂(Z)]^T，“*”，符号表示二维卷积，I表示对应Z点灰度。

4.根据权利要求3所述的基于分数阶微分的电影核磁共振图像序列运动场估计方法，其特征在于，步骤三所述利用单演信号的相位向量建立光流方程的过程为：

基于局部约束，假定以p点为中心的一个小区域Ω内各点的光流相同，对区域内不同的点给予不同的权重，离p点越近，权重越高；

列出光流的计算方程： $\underset{x\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{W}^2\left(p\right)[\&dtri;I(p,t)\&CenterDot;v+I_t(p,t){]}^2=0---\left(7\right)$

上式中，Ω代表2维图像中以p点为中心的一个小的区域，W(p)为窗函数，代表区域中各点的权重，离p点越近，权重越高，表示t时刻p点灰度的梯度，v＝[U V]^T表示运动向量，U表示在X轴方向位移，V表示在Y轴方向位移，I_t(p,t)表示图像中p点在t时刻的灰度；

对于邻域Ω内的n点，其中，

$A={(\&dtri;I\left(p_1\right),\&dtri;I\left(p_2\right),...,\&dtri;I\left(p_n\right))}^T,$

W＝diag(W(p₁),W(p₂),...,W(p_n))，diag表示主对角矩阵是一个主对角线之外的元素皆为0的矩阵，

b＝-(I_t(p₁),I_t(p₂),...,I_t(p_n))^T。

5.根据权利要求4所述的基于分数阶微分的电影核磁共振图像序列运动场估计方法，其特征在于，所述步骤四通过光流方程估计电影核磁共振图像序列的运动场的过程为：

由光流计算方程： $\underset{x\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{W}^2\left(p\right)[\&dtri;I(p,t)\&CenterDot;v+I_t(p,t){]}^2=0,$

解得运动向量v＝(A^TW²A)^-1A^TW²b；

其中A^TW²A为2×2矩阵：

$A^T{W}^{2}A=[\left.\begin{array}{cc}\mathrm{\&Sigma;}{W}^2\left(p\right)I_x^2& \mathrm{\&Sigma;}{W}^2\left(p\right)I_x{I}_y\\ \mathrm{\&Sigma;}{W}^2\left(p\right)I_y{I}_x& \mathrm{\&Sigma;}{W}^2\left(p\right)I_y^2\end{array}\right],$

上式中所有的求和都是在邻域Ω内所有点上进行的；

Ix和Iy表示灰度在X轴和Y轴方向的变化率，本方法用单演信息相位沿方向n的变化率代替I_x和I_y；

代入单演信息和求出单演相位信息下的运动向量v；

对电影核磁共振图像进行区域划分，依次求出各个区域内单演相位信息下的运动向量v，根据各个运动向量构造出电影核磁共振图像的运动场。