[发明专利]能够减少物体运动造成的成像伪影的生物医学图像配准方法和计算机程序产品

权利要求书

1.用于减少物体运动造成的成像伪影的生物医学图像配准方法，其包括以下步骤：

a)提供同一物体的至少第一幅和第二幅数字或数字化图像，或者至少第一组和第二组该物体的截面图像，所述图像由像素或体素的二维或三维阵列形成；

b)通过选取一定数量的被设为标志或特征的像素或体素，在第一幅图像或第一组图像中定义一定数量的标志，即所谓的特征，并生成要跟踪的所述特征的列表；

c)通过对被选定为特征的每个像素或体素确定从第一幅到第二幅图像或从第一组到第二组图像的光流矢量，从第一幅到第二幅图像或从第一组到第二组图像跟踪选定为特征的每个像素或体素的位置；

d)通过对第二幅或第二组图像的像素或体素施加逆光流，来对第一幅和第二幅图像、或对第一组和第二组图像进行配准；

其特征在于

执行自动特征选择步骤，该步骤包括：

B1)在第一幅图像或第一组截面图像的每个像素或体素周围定义像素或体素邻域，所述的像素或体素邻域包括有限数量的像素或体素；对于单幅图像，选用二维邻域；对于截面图像组，则选用三维邻域；

B2)对于每个像素或体素，确定所述像素或体素邻域的所有像素或体素的平均信号强度值；

B3)定义平均信号强度值阈值；

B4)比较在步骤B2中为每个像素或体素邻域确定的平均信号强度，并将所述平均信号强度值与平均信号强度值阈值进行比较；

B5)如果步骤B4中所述邻域的平均信号强度值大于阈值，则将所述像素或体素定义为要跟踪特征，并将其添加到要跟踪特征列表中。

2.根据权利要求1所述的方法，其特征在于，借助实验，通过在一组测试或样本图像上执行所述方法，并比较取不同阈值时获得的结果，来确定所述平均信号强度阈值。

3.根据权利要求1或2所述的方法，其特征在于，通过所谓的Lucas & Kanade自动特征跟踪方法或算法，来自动地跟踪被选作在第二幅图像中要跟踪的特征或标志的图像的像素或体素。

4.根据权利要求3所述的方法，其特征在于，将发明出来仅作用于二维图像的Lucas & Kanade算法的原始实现推广到作用于三维图像体积而非二维图像，以适用于生物医学截面图像体积，其包含如下步骤：

a)提供由三维体素阵列构成的第一三维图像体积I；

b)对于每个目标体素定义体素邻域，该邻域包围所述目标体素，并具有中心在所述目标体素的体素阵列，该体素阵列具有可调节的大小，例如3×3×3；

c)定义笛卡尔坐标系统，以及相对于所述笛卡尔坐标系统的每个轴，计算所谓的梯度体积，由第一图像体积的每个体素相对于其邻域体素的强度梯度来形成所述梯度体积，且定义为：

a)X方向的梯度体积：$I_{\mathrm{\Δx}}(x,y,z)=\frac{I(x+1,y,z)-I(x+1,y,z)}{2}$

b)Y方向的梯度体积：$I_{\mathrm{\Δy}}(x,y,z)=\frac{I(x,y+1,z)-I(x,y-1,z)}{2}$

c)Z方向的梯度体积：$I_{\mathrm{\Δz}}(x,y,z)=\frac{I(x,y,z+1)-I(x,y,z-1)}{2}$

d)对于每个目标体素，计算所谓的梯度矩阵，所述梯度矩阵定义为：

$G=\left[\begin{array}{ccc}{m}_{200}& m_{110}& m_{101}\\ m_{110}& m_{020}& m_{011}\\ m_{101}& m_{011}& m_{002}\end{array}\right],$其中$\begin{array}{c}{m}_{200}=I_{\mathrm{\Δx}}^2(x,y,z)\\ m_{020}=I_{\mathrm{\Δy}}^2(x,y,z),m_{002}=I_{\mathrm{\Δz}}^2(x,y,z),\\ m_{110}=I_{\mathrm{\Δx}}(x,y,z)\·I_{\mathrm{\Δy}}(x,y,z),\\ m_{101}=I_{\mathrm{\Δx}}(x,y,z)\·I_{\mathrm{\Δz}}(x,y,z),\\ m_{011}=I_{\mathrm{\Δy}}(x,y,z)\·I_{\mathrm{\Δz}}(x,y,z)\end{array}$

e)对于三维图像体积的每个目标体素的每个梯度矩阵，应用下式来计算最小特征值λ_m：

定义：

c＝m₂₀₀·m₀₂₀$d=m_{011}^2$e＝m₁₁₀·m₁₀₁ f＝m₁₀₁·m₁₀₁ p＝-m₂₀₀-m₀₂₀-m₀₀₂

q＝c+(m₂₀₀+m₀₂₀)m₀₀₂-d-e-f

r＝(e-c)m₀₀₂+dm₂₀₀-2(m₁₁₀·m₀₁₁·m₁₀₁)+f·m₀₂₀

$a=q-\frac{p^2}{3}$$b=\frac{2p^3}{27}-\frac{\mathrm{pq}}{3}+r$

$\mathrm{\θ}=\frac{1}{3}{\cos }^{-1}\left(\frac{3b}{\mathrm{an}}\right),$其中$n=2\sqrt{\frac{-a}{3}}$

注意：由于G是中心矩矩阵，所以$\left|3\frac{b}{\mathrm{an}}\right|\≤1^a\≤0$

如下，计算所述梯度矩阵的最小特征值：

${\mathrm{\λ}}_1=n\·\cos \mathrm{\θ}-\frac{p}{3}$

${\mathrm{\λ}}_2=n\left[\frac{\cos \mathrm{\θ}+\sqrt{3}\sin \mathrm{\θ}}{2}\right]-\frac{p}{3}$λ_m＝minimum(λ₁，λ₂，λ₃)

${\mathrm{\λ}}_3=n\left[\frac{\cos \mathrm{\θ}-\sqrt{3}\sin \mathrm{\θ}}{2}\right]-\frac{p}{3}$

f)定义最小特征值λ_m的阈值；

g)将对应梯度矩阵的最小特征值λ_m满足一定标准的每个体素选择作为图像I中要跟踪的特征，这些标准如下：

i)λ_m大于所述阈值；

ii)λ_m不小于所有最小特征值λ_m中的最大值的某个百分比；

iii)如果在保持为选定特征的目标体素周围的已定义体素邻域内存在另一个被选作特征的体素，则仅将其梯度矩阵具有较大的最小特征值的一个体素选为特征，而将另一个体素从可跟踪特征列表中删除；

iv)选作特征的体素周围的三维邻域的平均信号强度值大于可调节平均信号强度值阈值；