CN104392452A - 一种基于应用的dicom医学图像处理方法 - Google Patents

Abstract

本发明提供了一种基于应用的DICOM医学图像处理方法，该方法包括：步骤一，输入两幅待配准的DICOM图像，依次对其进行格式转换、图像去噪处理和图像分割，从分割后的图像中抽取出特征信息组成特征空间；步骤二，给定初始变换，采用线性插值方法对特征空间进行线性变换；步骤三，选择优化策略，在特征空间中寻找最优变换；步骤四，根据预定义相似度函数来测量两幅待配准的DICOM图像的相似度，如果当前变换满足预设相似度阈值，则输出该变换，否则重复步骤二和步骤三。本发明提出了改进的图像配准方法，有效提高了DICOM医学图像配准的速度和精度。

Description

一种基于应用的DICOM医学图像处理方法

技术领域

本发明涉及医学图像处理，特别涉及一种基于应用的DICOM医学图像处理方法。

背景技术

随着计算机科学的发展，医学影像信息学已成为国际上发展迅速的一个领域。DICOM作为医学图像数据标准的使用使不同厂家，不同设备之间医学图像数据传输，交换和存储成为可能。DICOM标准强大的兼容性使它越来越广泛地应用到医学领域，逐步成为一种标准的医学图像数据格式。DICOM医学图像为临床诊断医学、计算机辅助手术、病灶监控等提供了高质量的信息支持，但是，DICOM图像也带来了一些问题，如时序差异、位置差异、获得设备差异等导致影像版本信息的不一致性。怎样消除各种差异通过图像配准获取最大有效信息成为亟待解决的问题。医学图像的配准关键在于界定医学图像之间的相似性。在临床方面，这种方法常用于匹配一段时间内来自不同或同一患者的医学图像。如果将普通图像配准直接应用于结构复杂的DICOM医学图像，则不仅计算量非常大，效率低；而且容易产生误差导致结果不准确。

因此，针对相关技术中所存在的上述问题，目前尚未提出有效的解决方案。

发明内容

为解决上述现有技术所存在的问题，本发明提出了一种基于应用的DICOM医学图像处理方法，包括：

步骤一，输入两幅待配准的DICOM图像，依次对其进行格式转换、图像去噪处理和图像分割，从分割后的图像中抽取出特征信息组成特征空间；

步骤二，给定初始变换，采用线性插值方法对特征空间进行线性变换；

步骤三，选择优化策略，在特征空间中寻找最优变换；

步骤四，根据预定义相似度函数来测量两幅待配准的DICOM图像的相似度，如果当前变换满足预设相似度阈值，则输出该变换，否则重复步骤二和步骤三。

优选地，所述线性插值进一步包括：

利用最近邻域点灰度值的加权平均值来表示该点灰度值，通过线性插值计算直接修改直方图的项，将每对像素对联合直方图的贡献分散到联合直方图上相邻的各个像素对，联合直方图上各点的值以小于1的小数增加，从而得到一个平滑的配准函数。

优选地，所述图像去噪进一步包括：

对于每个种子像素点I(k，l，n)，选定最大搜索区域为M×M，其中k和l代表像素坐标，n代表帧数；

种子像素点的8邻域像素表示为I(i,j，n)，i,j∈Ω，根据以下准则为每个像素确定一个权值：

$|I(k,l,n)-I(i,j,n)|<\sqrt{2}{\mathrm{\&sigma;}}_N$

其中的σ_N为图像噪声标准差的估计值；

通过加权平均得到该种子像素点的估计灰度值：

$I_e(k,l,n)=\frac{\underset{i,j\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}p(i,j)I(i,j,n)}{\underset{i,j\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}p(i,j)}$

其中 $p(i,j)=\left\{\begin{array}{cc}1& |I(k,l,n)-I(i,j,n)|<\sqrt{2}{\mathrm{\&sigma;}}_N\\ 0& I(k,l,n)-I(i,j,n)|\&GreaterEqual;\sqrt{2}{\mathrm{\&sigma;}}_N\end{array}\right.,$ Ω表示8邻域像素坐标空间；

然后在选定的搜索区域M×M内，以I_e(k，l，n)为种子像素点构建邻域：将邻域内的像素组成自适应邻域，得到一个像素点集合；

设置空域的高斯核函数为：

$g_s(i,j)=\exp (-\frac{D(I(k,l,n),I(i,j,n))}{2{\mathrm{bn}}^2})$

其中D代表种子像素点I(k,l,n)和自适应邻域像素点集合中任意一点I(i,j,n)之间的欧氏距离，bn的值取最大搜索区域尺寸M的一半；

设置灰度域的高斯核函数为：

$g_g(i,j)=\exp (-\frac{{(I_e(k,l,n),I(i,j,n))}^2}{2{\mathrm{bn}}^2})$

其中σ_N仍为图像噪声标准差的估计值，通过空间域和灰度域的高斯核函数，按下式进行滤波：

$\hat{I}(k,l,n)=\frac{\underset{i,j\&Element;V}{\mathrm{\&Sigma;}}\left[g_s(i,j)g_g(i,j)I(i,j,n)\right]}{\underset{i,j\&Element;V}{\mathrm{\&Sigma;}}\left[g_s(i,j)g_g(i,j)\right]}$

其中V表示自适应邻域像素点集合的坐标空间。

优选地，对于待处理的第n帧，建立一个图像缓冲器，用来暂存第(n-1)、n、(n+1)三帧图像，首先计算以两个对应像素点为中心，大小为3×3的前后两帧图像块之间灰度差的绝对值和SAD_prev和SAD_next：

${\mathrm{SAD}}_{\mathrm{prev}}=\underset{i,j\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}|\hat{I}(i,j,n)-\hat{I}(i,j,n-1)|$

${\mathrm{SAD}}_{\mathrm{next}}=\underset{i,j\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}|\hat{I}(i,j,n)-\hat{I}(i,j,n+1)|$

其中表示待处理的帧，它与其前后两帧图像都经过的空域和灰度域滤波，Ω表示3×3图像块的像素坐标空间；

设置时域的高斯核函数为：

$g_t(\mathrm{SAD},T)=\exp \left(\frac{{\mathrm{SAD}}^2}{2T^2}\right)$

其中T是由噪声引起的误差决定的常量参数，最终的时域滤波按下式进行：

${\hat{I}}^{\&prime;}(k,l,n)=\frac{\hat{I}(k,l,n)+g_t({\mathrm{SAD}}_{\mathrm{prev}},T)\hat{I}(k,l,n-1)+g_t({\mathrm{SAD}}_{\mathrm{next}},T)\hat{I}(k,l,n+1)}{1+g_t({\mathrm{SAD}}_{\mathrm{prev}},T)+g_t({\mathrm{SAD}}_{\mathrm{next}},T)}.$

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

提出了改进的图像配准方法，有效提高了DICOM医学图像配准的速度和精度。

附图说明

图1是根据本发明实施例的基于应用的DICOM医学图像处理方法的流程图。

具体实施方式

下文与图示本发明原理的附图一起提供对本发明一个或者多个实施例的详细描述。结合这样的实施例描述本发明，但是本发明不限于任何实施例。本发明的范围仅由权利要求书限定，并且本发明涵盖诸多替代、修改和等同物。在下文描述中阐述诸多具体细节以便提供对本发明的透彻理解。出于示例的目的而提供这些细节，并且无这些具体细节中的一些或者所有细节也可以根据权利要求书实现本发明。

本发明提出一种DICOM医学图像配准方法。先对DICOM医学图像进行预处理，然后测量待配准图像间的相似度，最后优化图像。图1是根据本发明实施例的基于应用的DICOM医学图像处理方法流程图。首先抽取出图像的特征信息组成特征空间；然后根据所抽取的特征空间确定出一个空间变换，使一幅图像经过该变换后能够达到所定义的相似性测度。在确定变换的过程中若采取一定的优化措施或搜索策略，则可使相似性测度更好、更快地达到最优值。

本发明寻找最佳空间和灰度变换函数，便于图像使用相似性测度进行匹配，其过程就是寻找最优变换，其中该最优变换将图像B进行变换，以使其尽可能多的包含图像A的信息。用互信息表示相似性测度对两图像进行精确配准，通过对完整的医学图像寻找一个仿射变换即旋转参数和平移参数，该仿射变换对应的两个图像的互信息最大，此时图像实现精确配准互信息用熵来表示：

对于概率分布函数为p(a)的随机变量集A，其熵H(A)定义如下：

H(A)＝-∑p(a)logp(a)a∈A

对两个图像，它们的离散的随机变量集分别是A和B，它们的边缘概率分布函数分别是p(a)和p(b)，联合概率分布函数是P_AB(a，b)，则随机变量A和B的联合熵定义如下：

H(A)＝-∑p(a)logp(a)logp(b)a∈Ab∈B

如果H(A/B)表示已知B是A的条件熵，那么H(A)与H(A/B)的差值，就代表了在B中包含的A的信息，即互信息；因此两个随机变量集间的互信息定义为：

I(A，B)＝H(A)+H(B)-H(A，B)

＝H(A)-H(A/B)

＝H(B)-H(B/A)

在图像配准中，当两幅图像的空间位置达到完全一致时，其中一幅图像A表达的关于另一幅图像B的信息，也就是对应像素灰度的互信息应为最大；用联合概率分布和完全独立时的概率分布间的广义距离来估计互信息：

$I(A,B)=\mathrm{\&Sigma;p}(a,b)\log \frac{p(a,b)}{p\left(a\right)p\left(b\right)}$

对于待配准的两幅图像，可以认为它们是关于图像灰度的两个随机变量集参考图像和浮动图像，a和b是两幅图像中相关的像素灰度值，a和b通过图像之间的坐标变换相联系，对于离散的数字图像，联合概率分布p_AB(a，b)可以用归一化的联合直方图表示，则有：

$I(A,B)=\underset{i,j}{\mathrm{\&Sigma;}}{p}_{\mathrm{AB}}(i,j)\log \frac{p_{\mathrm{AB}}(i,j)}{p_A\left(i\right)p_B\left(j\right)}$

接下来寻找一个合适的线性变换，使得一幅试验图像经过此变换后和参考图像的互信息最大。

本发明的配准方法总体步骤为：

1、输入两幅待配准的DICOM图像，即参考图像和试验图像，对其进行格式转换；然后进行图像去噪和图像分割，从处理后的图像中抽取出特征信息组成特征空间。

2、给定初始变换，使用线性插值方法进行线性变换。

3、选择优化策略，在特征空间中寻找最优变换。

4、测量参考图像和试验图像的相似度，如果当前最优变换满足预设相似度阈值，则输出该变换，否则重复步骤2和3。

在图像预处理的图像格式转换过程中，首先在DIOCM文件中找到数据域的起始位置，用pdata表示，根据DICOM的组标号，单位标号，获取数据的帧数和数据长度。判断当前字节排序方式，如果排序方式是BIG_ENDIAN，则分别交换组标号和单位标号的2个单元。若组标号为0x0002，单元标号是0x0010，则获取UID，并根据UID判断数据的字节顺序；若组标号为0x0002，获取单元标号，根据所获取的单元标号的值确定数据帧数，每帧中的行数和列数；若组标号＝0x7fe0，则设置位图的文件头、信息头。判断是否需要设定颜色表，如果需要，则获取单元标号，根据所获取的单元标号的值确定数据的存储长度，并获取数据起始位置。设置颜色表并按照逆序方式存储pdata数据到位图文件，如果不需要，则按照行逆序和RGB逆序的方式存储pdata数据到位图文件，最后对位图文件进行显示。

DICOM文件可能封装了多幅图像，对于这样的DICOM文件，需要找到每幅图像的偏移量，根据偏移量分次转换。在DICOM图像文件转换为位图文件只需把DICOM中的像素数据输出到BMP文件即可，但在转换前需考虑窗宽和窗位。窗宽是指图像数据显示的范围，窗位是指图像数据显示的中心值，两指标的值根据不同的生理组织结构设定不同的值，当把DICOM文件的全部像素读入内存后，就要考虑如何把DICOM显示到屏幕上，由于DICOM图像像素的灰阶可达到2¹⁶，而系统的调色板只能显示2⁸种颜色，所以为了较逼真的显示医学图像的全部信息，必须对图像数据进行窗宽、窗位的调整，把窗口区域图像数据线形的转换到显示器的最大显示范围内，这样就可以动态调整窗宽和窗位，观察医学图像的全部信息。可对图像的每个数据进行转换，把数据压缩到256个灰度等级。

其中：y为显示的位图灰度值，x为图像的数据，w为窗宽，c为窗位

$y\left(x\right)=\left\{\begin{array}{cc}0& c-w/2>x\\ [x-(c-w/2)]255/w& c-w/2<x<c+w/2\\ 255& c+w/2<x\end{array}\right.$

由于DICOM图像和位图所规定的图像数据的存储顺序不同，所以在存放位图数据时，其顺序是先从最后一行开始，依次向上，直到第一行。

位图文件有4部分数据项组成：位图文件头、位图信息头、颜色表和字节阵列。初始化位图文件实现了对位图文件数据项的赋值和对颜色表的设置。设置彩色表是将数据的灰度值映射为较少级数的*’+图像下的灰度值。位图文件头主要包括文件类型，文件大小，字节偏移量等。每个数据项的值有DICOM格式的行数、列数、帧数等确定。

要实现位图数据的写入，应首先确定DICOM文件中数据域的起始地址值，该起始地址是通过DICOM数据预处理得到的，用pdata来记录该值。位图数据的写入分两种情况，一种是DICOM格式下真彩色存储的情况，另一种是DICOM格式下非真彩色存储的情况。两种情况在由DICOM格式转存为位图格式时都将采取逆序方式存储原来数据域的值，不同点在于非真彩色存储时只需将DICOM数据按行逆序存储即可转换为相应位图，而真彩色除了按行逆序存储外还需对RGB值也逆序存储。将DICOM格式的图像转换成非真彩色存储时，每个数据单元的值不会发生变化，但是存储位置都将发生变化，主要表现在逆序存储每一行的数据值，即把DICOM数据中最后一行的值存储在位图非真彩色存储的第一行上，第一行的数据存储在最后一行上。真彩色格式像素的存储比非真彩色要复杂得多，除了采取行逆序存储外，还要对每个像素中RGB值逆序存储。

在图像预处理的图像去噪过程中，医学图像通常存在大量的斑点噪声，使DICOM图像的质量明显下降，甚至掩盖了某些组织细节，这对图像特征的提取、识别与分析带来了极大困难。通常采用滤波的方法来抑制噪声频率成分，从而提高图像的信噪比，改善图像质量。本发明在现有技术的基础上，相继在空域、灰度域、时域中进行滤波，有效抑制图像中的噪声，并保留边缘细节信息。提高了算法的鲁棒性。有利于特征提取，识别和分析。具体的去噪过程为：

对于一个种子像素点，I(k，l，n)，选定最大搜索区域为M×M，n代表帧数。先考虑种子像素点的8邻域像素，I(i,j，n)，并根据以下准则为每个像素确定一个权值： $|I(k,l,n)-I(i,j,n)|<\sqrt{2}{\mathrm{\&sigma;}}_N$

其中的σ_N为图像噪声标准差的估计值。然后，通过加权平均得到一个种子像素点的估计灰度值：

$I_e(k,l,n)=\frac{\underset{i,j\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}p(i,j)I(i,j,n)}{\underset{i,j\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}p(i,j)}$

其中 $p(i,j)=\left\{\begin{array}{cc}1& |I(k,l,n)-I(i,j,n)|<\sqrt{2}{\mathrm{\&sigma;}}_N\\ 0& I(k,l,n)-I(i,j,n)|\&GreaterEqual;\sqrt{2}{\mathrm{\&sigma;}}_N\end{array}\right.,$ Ω表示8邻域像素坐标空间。这一步是为了减小种子像素点本身噪声较大时对自适应邻域构建的影响。

然后在选定的搜索区域M×M内，以I_e(k，l，n)为种子像素点构建邻域：将满足上式的像素，将其纳入最终的自适应邻域，从而得到一个像素点集合。在这一步中，阈值放大到了配合高斯加权滤波，使自适应邻域的构建不会对这个阈值过分依赖。

设置空域的高斯核函数为：

$g_s(i,j)=\exp (-\frac{D(I(k,l,n),I(i,j,n))}{2{\mathrm{bn}}^2})$

其中D代表种子像素点I(k,l,n)和自适应邻域像素点集合中任意一点I(i,j,n)之间的欧氏距离，bn的值取最大搜索区域尺寸M的一半。

灰度域的高斯核函数为：

$g_g(i,j)=\exp (-\frac{{(I_e(k,l,n),I(i,j,n))}^2}{2{\mathrm{bn}}^2})$

其中σ_N仍为图像噪声标准差的估计值。结合空间域和灰度域的高斯核函数，最终的滤波按下式进行：

$\hat{I}(k,l,n)=\frac{\underset{i,j\&Element;V}{\mathrm{\&Sigma;}}\left[g_s(i,j)g_g(i,j)I(i,j,n)\right]}{\underset{i,j\&Element;V}{\mathrm{\&Sigma;}}\left[g_s(i,j)g_g(i,j)\right]}$

其中V表示自适应邻域像素点集合的坐标空间。

为了避免ghost和提高噪声抑制能力，对于待处理的第n帧，其前后两帧与它的相关性最大。为此本发明建立一个图像缓冲器，用来暂存第(n-1)、n、(n+1)三帧图像。首先，计算以两个对应像素点为中心，大小为3×3的图像块之间灰度差的绝对值和SAD_prev和SAD_next：

${\mathrm{SAD}}_{\mathrm{prev}}=\underset{i,j\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}|\hat{I}(i,j,n)-\hat{I}(i,j,n-1)|$

${\mathrm{SAD}}_{\mathrm{next}}=\underset{i,j\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}|\hat{I}(i,j,n)-\hat{I}(i,j,n+1)|$

其中表示待处理的帧，它与其前后两帧图像都已经过了上一步的空域和灰度域滤波，Ω表示3×3图像块的像素坐标空间。时域的高斯核函数为：

$g_t(\mathrm{SAD},T)=\exp \left(\frac{{\mathrm{SAD}}^2}{2T^2}\right)$

其中T是由噪声引起的误差决定的常量参数。最终的时域滤波按下式进行：

${\hat{I}}^{\&prime;}(k,l,n)=\frac{\hat{I}(k,l,n)+g_t({\mathrm{SAD}}_{\mathrm{prev}},T)\hat{I}(k,l,n-1)+g_t({\mathrm{SAD}}_{\mathrm{next}},T)\hat{I}(k,l,n+1)}{1+g_t({\mathrm{SAD}}_{\mathrm{prev}},T)+g_t({\mathrm{SAD}}_{\mathrm{next}},T)}.$

进一步地，在步骤2中，本发明采用线性变换和刚性配准。在医学图像配准中，称刚体变换和仿射变换为刚性配准。当试验图像变换后坐标位置不在网格点，通常采用插值技术对其灰度值进行估计。插值方法影响图像的平滑性、优化性能和总体计算时间。为了提高优化效率和节省总体时间，本发明采用了线性插值和为建立联合灰度直方输入DICOM图像而进行部分体积插值。利用最近邻域点灰度值的加权平均值来表示该点灰度值。部分体积插值使用线性插值权值，但并不通过邻域点加权来计算该点灰度值和更新一个简单直方图条目，而是通过插值计算直接修改几个直方图的项，将每对像素对联合直方图的贡献分散到联合直方图上与之相邻的各个像素对，联合直方图上各点的值以小于1的小数增加，从而得到一个平滑的配准函数，有利于优化搜索。

步骤3中，配准方法的核心是一个含有能搜索到最小化测度函数的变换模型参数的优化框架。在基于特征配准方法中，例如一系列点相关的点匹配，使用直接优化技术，如最小范数的解决方案。在较为复杂的情况下，使用迭代优化方法：

步骤A1、给定允许误差ε>0，选取初始点X₀(0，0，0)和一组线性无关的搜索方向向量e⁽ⁱ⁾(i＝1，2，3)，其中e⁽¹⁾＝(1，0，0)，e⁽²⁾＝(0，1，0)，e⁽³⁾＝(0，0，1)，即令e⁽ⁱ⁾分别为沿坐标轴的方向；设置迭代计数器K＝1；

步骤A2、从X₀出发沿方向e⁽ⁱ⁾(i＝1，2，3)依次进行一维搜索，得到

X⁽ⁱ⁾＝X_(i-1)+λ_i×e⁽ⁱ⁾

f(X_(i))＝min(f(X_(i-1)+λ_i×e⁽ⁱ⁾))

其中λ_i为预设常数权值，在完成了这3次一维搜索后，得到X₍₃₎；

步骤A3、检验是否已满足终止准则：若|f(X_(i))-f(X₍₀₎)|≤ε或者迭代次数大于等于U，则停止迭代，退出；其中，ε为预先设定的阈值，U是预设的最大迭代次数；

步骤A4、按照下式计算最速下降方向上的函数的变化量：

Δf＝max|f(X(i))-f(X(0))|(1<i<N)，并记最快下降方向为e_max；

步骤A5、取e＝X₍₃₎-X₍₀₎，Pcs＝2X₍₃₎-X₍₀₎，计算f(e)＝f(X₍₃₎-X₍₀₎)，f(Pcs)＝f(2X₍₃₎-X₍₀₎)，若f(e)≥f(x₍₀₎)或f(Pcs)<f(x₍₀₎)，且2(f(x₍₀₎)-2f(x_(n))+f(e))[(f(x₍₀₎)-f(x_(n)))-Δf]²≥(f(x₍₀₎)-f(e))²Δf，则将X₍₃₎赋值给X₍₀₎作为新的起点，沿上面的一组旧方向向量e⁽ⁱ⁾(i＝1，2，3)重复步骤A2至A5；若上述条件均不满足，则沿方向e＝X₍₃₎-X₍₀₎，以X₍₃₎作为起始点进行搜索得到目标函数在此方向上的极小值点P_min；将原来的下降方向最大的e_max去掉，而保留其它的2个方向，再加上方向e，仍得到3个方向：e⁽¹⁾，e⁽²⁾，e⁽³⁾，把此时的P_min赋给X₍₀₎作为新起点，重复步骤A2至A5。

本发明采用的迭代优化方法进行配准参数的寻优，不必对目标函数求导，具有二次收敛性，收敛速度快，可加快搜索最大值的速度，且便于编程实现。

综上所述，本发明提出了改进的图像配准方法，有效提高了DICOM医学图像配准的速度和精度。

显然，本领域的技术人员应该理解，上述的本发明的各模块或各步骤可以用通用的计算系统来实现，它们可以集中在单个的计算系统上，或者分布在多个计算系统所组成的网络上，可选地，它们可以用计算系统可执行的程序代码来实现，从而，可以将它们存储在存储系统中由计算系统来执行。这样，本发明不限制于任何特定的硬件和软件结合。

应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和边界、或者这种范围和边界的等同形式内的全部变化和修改例。

Claims (4)

1.一种基于应用的DICOM医学图像处理方法，其特征在于，包括：

步骤一，输入两幅待配准的DICOM图像，依次对其进行格式转换、图像去噪处理和图像分割，从分割后的图像中抽取出特征信息组成特征空间；

步骤二，给定初始变换，采用线性插值方法对特征空间进行线性变换；

步骤三，选择优化策略，在特征空间中寻找最优变换；

步骤四，根据预定义相似度函数来测量两幅待配准的DICOM图像的相似度，如果当前变换满足预设相似度阈值，则输出该变换，否则重复步骤二和步骤三。

2.根据权利要求1所述的方法，其特征在于，所述线性插值进一步包括：

利用最近邻域点灰度值的加权平均值来表示该点灰度值，通过线性插值计算直接修改直方图的项，将每对像素对联合直方图的贡献分散到联合直方图上相邻的各个像素对，联合直方图上各点的值以小于1的小数增加，从而得到一个平滑的配准函数。

3.根据权利要求2所述的方法，其特征在于，所述图像去噪进一步包括：

对于每个种子像素点I(k,l,n)，选定最大搜索区域为M×M，其中k和l代表像素坐标，n代表帧数；

种子像素点的8邻域像素表示为I(i,j,n)，i,j∈Ω，根据以下准则为每个像素确定一个权值：

$|I(k,1,n)-I(i,j,n)|<\sqrt{2}{\mathrm{\&sigma;}}_N$

其中的σ_N为图像噪声标准差的估计值；

通过加权平均得到该种子像素点的估计灰度值：

$I_e(k,l,n)=\frac{\underset{i,j\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}p(i,j)I(i,j,n)}{\underset{i,j\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}p(i,j)}$

其中 $p(i,j)=\left\{\begin{array}{cc}1& |I(k,l,n)-I(i,j,n)|<\sqrt{2}{\mathrm{\&sigma;}}_N\\ 0& I(k,l,n)-I(i,j,n)|\&GreaterEqual;\sqrt{2}{\mathrm{\&sigma;}}_N\end{array}\right.,$ Ω表示8邻域像素坐标空间；

然后在选定的搜索区域M×M内，以I_e(k,l,n)为种子像素点构建邻域：将邻域内的像素组成自适应邻域，得到一个像素点集合；

设置空域的高斯核函数为：

$g_s(i,j)=\exp (-\frac{D(I(k,l,n),I(i,j,n))}{2{\mathrm{bn}}^2})$

其中D代表种子像素点I(k,l,n)和自适应邻域像素点集合中任意一点I(i,j,n)之间的欧氏距离，bn的值取最大搜索区域尺寸M的一半；

设置灰度域的高斯核函数为：

$g_g(i,j)=\exp (-\frac{{(I_e(k,l,n),I(i,j,n))}^2}{2{\mathrm{bn}}^2})$

其中σ_N仍为图像噪声标准差的估计值，通过空间域和灰度域的高斯核函数，按下式进行滤波：

$\hat{I}(k,l,n)=\frac{\underset{i,j\&Element;V}{\mathrm{\&Sigma;}}\left[g_s(i,j)g_g(i,j)I(i,j,n)\right]}{\underset{i,j\&Element;V}{\mathrm{\&Sigma;}}\left[g_s(i,j)g_g(i,j)\right]}$

其中V表示自适应邻域像素点集合的坐标空间。

4.根据权利要求3所述的方法，其特征在于，还包括：

对于待处理的第n帧，建立一个图像缓冲器，用来暂存第(n-1)、n、(n+1)三帧图像，首先计算以两个对应像素点为中心，大小为3×3的前后两帧图像块之间灰度差的绝对值和SAD_prev和SAD_next：

${\mathrm{SAD}}_{\mathrm{prev}}=\underset{i,j\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}|\hat{I}(i,j,n)-\hat{I}(i,j,n-1)|$

${\mathrm{SAD}}_{\mathrm{next}}=\underset{i,j\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}|\hat{I}(i,j,n)-\hat{I}(i,j,n+1)|$

其中表示待处理的帧，它与其前后两帧图像都经过的空域和灰度域滤波，Ω表示3×3图像块的像素坐标空间；

设置时域的高斯核函数为：

$g_t(\mathrm{SAD},T)=\exp (-\frac{{\mathrm{SAD}}^2}{2T^2})$

其中T是由噪声引起的误差决定的常量参数，最终的时域滤波按下式进行：

${\hat{I}}^{\&prime;}(k,l,n)=\frac{\hat{I}(k,l,n)+g_t({\mathrm{SAD}}_{\mathrm{prev}},T)\hat{I}(k,l,n-1)+g_t({\mathrm{SAD}}_{\mathrm{next}},T)\hat{I}(k,l,n+1)}{1+g_t({\mathrm{SAD}}_{\mathrm{prev}},T)+g_t({\mathrm{SAD}}_{\mathrm{next}},T)}.$