CN101488219A - 一种快速视频图像双边滤波的方法 - Google Patents

Abstract

本发明公开了一种快速视频图像双边滤波的方法，通过将滤波滑动窗分解出多个重叠的常数空间滤波平面，而同时对每个常数空间滤波平面进行像素密度值的直方图统计，在直方图统计结果的基础上把像素值域的滤波处理以常规的加权求和来实现。因为不必再逐像素的计算像素值域上滤波核系数，因此滤波处理的时间复杂度显著降低。在累积直方图计算的基础上，以逐像素扫描的方式进行滑动窗直方图更新，而不是针对每个滤波处理像素都要构建一个独立的直方图，从而避免了大量的直方图统计计算复杂度。

Description

一种快速视频图像双边滤波的方法

技术领域

本发明涉及视频图像处理技术领域，特别是在视频图像增强处理系统中，对视频图像进行快速双边滤波的方法。

背景技术

大多数视频图像处理技术需要相当多的计算功耗，不管图像处理技术是否以互动式或自动式使用，为完成任务都需要计算处理时间，特别时当有大量数据需要处理时。为加快视频图像处理速度一个方法是通过实现更快的电子芯片或增加处理任务的处理器数目，然而该方案特别的昂贵。另一个增加处理速度方法是利用别的更快的处理方法以产生与现存方法相同的结果，典型的来说，通过改进视频图像处理技术所需要的关键算法的执行时间，来完成这一目标。

视频图像处理中常用的一个技术是双边滤波器，是一个能在保持图像边缘细节的同时平滑图像的非迭代方法。双边滤波器实质是一个加权卷积运算，它每个像素相关的权重系数不只依赖于到中心像素的距离，也依赖于它相对的像素密度值距离。经典双边滤波核一般公式是：

f_classic(s，s₀)＝g_spatial(s-s₀)·g_tone(I(s)-I(s₀))       (1)

其中，s₀表示滤波核中心位置，s表示滤波核元素位置，I(s₀)表示滤波位置像素某通道的灰度值，I(s)表示滤波核其他位置上像素某通道的灰度值。

两个权重函数都是高斯函数。

g_spatial(s)＝g(x，σ_s)·g(y，σ_s)                        (2)

g_tone(I)＝g(I，σ_t)                                    (3)

其中：σ_s为双边滤波空间方差参数，σ_t为像素值域控制参数，x，y为像素坐标，

$g(x,{\mathrm{\σ}}_s)=\frac{1}{{\mathrm{\σ}}_s\sqrt{2\mathrm{\π}}}{e}^{-x^2/2{{\mathrm{\σ}}_s}^2}$

$g(y,{\mathrm{\σ}}_s)=\frac{1}{{\mathrm{\σ}}_s\sqrt{2\mathrm{\π}}}{e}^{-y^2/2{{\mathrm{\σ}}_s}^2}$

$g(I,{\mathrm{\σ}}_t)=\frac{1}{{\mathrm{\σ}}_t\sqrt{2\mathrm{\π}}}{e}^{-I^2/2{{\mathrm{\σ}}_t}^2}---\left(4\right)$

我们应该关注的本质属性是每个像素的运行时间，而这个运行时间是双边滤波器核的函数。当调节滤波器核大小γ，即核半径时用户会体验到相应的性能变化，而这是双边滤波器技术主要的特征差异。

现有双边滤波是结合了空间域与像素值域的滤波，所以它计算量要求非常大。一个标准双边滤波器的实现计算每个输出像素的复杂度是O(γ²)，随着滤波器参数调节滤波核半径变大会非常缓慢。

另外有一种快速双边滤波技术，其思路是通过在空间域和像素值域上进行下采样，随后在图像下采样后的版本上进行滤波处理，随后使用线性插值重新组合起滤波处理的结果，以近似逼近标准双边滤波处理结果而达到显著减少计算需求量的目标。换句话说，该方法把密度图像看作是3D空间向量，随后对该向量空间应用高斯光滑滤波，随后再插值重构滤波后的图像。3D向量空间有更多的下采样则滤波核半径越小，所以滤波处理速度越快。然而该方法的缺点是滤波处理的结果依赖于下采样阶段以及上采样的插值重构阶段，也丢失了一部分精确度。

另外还有一种快速双边滤波技术，使用层次式分布直方图方法能获得运行时间收敛到O(logγ)。然而，尽管其运行复杂度减少，因为对特定实现所需要的滤波器大小和直方图记数等，丢失了标准双边滤波技术原来的简洁性。

发明内容

本发明的目的在于克服现有快速视频图像双边滤波方法的不足，提供一种滤波计算结果近似于标准双边滤波计算的结果，计算复杂度显著降低的快速视频图像双边滤波方法。

为达到上述发明目的，本发明的快速视频图像双边滤波的方法，其特征在于，包括以下步骤：

(1)、定义一个大小由双边滤波空间方差参数σ_s决定的矩形区域作为滤波滑动窗，其中心像素P0是当前平滑滤波计算的像素；

(2)、计算图像像素累积直方图；

(3)、计算当前滤波滑动窗内像素直方图，方法为：

h_s＝H(x-，y-)-H(x，y-)-H(x-，y)+H(x，y)

其中，P(x-，y-)，P(x，y-)，P(x-，y)，P(x，y)分别表示滤波滑动窗内左上、右上、左下、右下四个角点像素，其对应的累积直方图为H(x-，y-)，H(x，y-)，H(x-，y)，H(x，y)，根据步骤(2)的计算结果可以得到；

(4)、将当前滤波滑动窗执行常数空间滤波平面分解

当前滤波滑动窗的宽度，即高斯滤波核宽度为kw＝2*r+1，相应的高斯核g_spatial(s-s₀)系数共有r+1个不同值，从内到外系数依次表示为：co0，co1，co2，...，cor，并依据以下公式计算出c0，c1，c2，...，cr等r+1个空间滤波常数系数：

c0＝co0；

c1＝co1-co0；

c2＝co2-c1-c0＝co2-co1；

c3＝co3-c2-c1-c0＝co3-co2；

……

cr＝cor-co(r-1)；

得到c0，c1，c2，...，cr等r+1个空间滤波常数系数后，将r+1个常数空间滤波平面分别表示为S(0)，S(1)，S(2)，…，S(r)；

(5)、对分解后的常数空间滤波平面计算其像素值直方图

按步骤(3)所述的方法计算分解后的r+1个常数空间滤波平面像素值直方图：

h_S(r)＝H(x-，y-)-H(x，y-)-H(x-，y)+H(x，y)；

h_S(r-1)＝H((x-)+1，(y-)+1)-H((x)-1，(y-)+1)-H((x-)+1，(y)-1)+H((x)-1，(y)-1)；

h_S(r-2)＝H((x-)+2*1，(y-)+2*1)-H((x)-2*1，(y-)+2*1)

        -H((x-)+2*1，(y)-2*1)+H((x)-2*1，(y)-2*1)；

……

h_S(0)＝H((x)-r*1，(y)-r*1)；

这里2r+1表示双边滤波核的宽度，h_S(0)，h_S(1)，...，h_S(r-2)，h_S(r-1)，h_S(r)分别表示r+1个常数空间滤波平面所对应的像素值直方图。

(6)、对标准双边滤波的图像平滑滤波近似计算：

通过求和分别计算出S(0)，S(1)，S(2)，…，S(r)等r+1个常数空间滤波的中间双边计算结果依次是I(s₀)₀，I(s₀)₁，I(s₀)₂，...，I(s₀)_r，随着在由这些中间双边计算结果I(s₀)₀，I(s₀)₁，I(s₀)₂，...，I(s₀)_r计算出图像的最终平滑滤波结果I(s₀)，

I(s₀)＝k^-1·(I(s₀)₀+I(s₀)₁+…+I(s₀)_r)

其中中间双边计算结果为：

I(s₀)₀＝c0·∑_i∈S(0)i·h_S(0)(i)·g_tone(i-I(s₀))；

I(s₀)₁＝c1·∑_i∈S(1)i·h_S(1)(i)·g_tone(i-I(s₀))；

……

I(s₀)_r＝cr·∑_i∈S(r)i·h_S(r)(i)·g_tone(i-I(s₀))；

归一化算子为

k＝k₀+k₁+…+k_r

k₀＝c0·∑_i∈S(0)h_S(0)(i)·g_tone(i-I(s₀))；

k₁＝c1·∑_i∈S(1)h_S(1)(i)·g_tone(i-I(s₀))；

……

k_r＝cr·∑_i∈S(r)h_S(r)(i)·g_tone(i-I(s₀))；

本发明滤波滑动窗分解出多个重叠的常数空间滤波平面，而同时对每个常数空间滤波平面进行像素密度值的直方图统计，在直方图统计结果的基础上把像素值域的滤波处理以常规的加权求和来实现。因为不必再逐像素的计算像素值域上滤波核系数，因此滤波处理的时间复杂度显著降低。在累积直方图计算的基础上，以逐像素扫描的方式进行滑动窗直方图更新，而不是针对每个滤波处理像素都要构建一个独立的直方图，从而避免了大量的直方图统计计算复杂度。

附图说明

图1是本发明滤波滑动窗示意图；

图2是当前像素灰度值直方图统计示意图；

图3是当前像素灰度值直方图统计的之字形扫描方式示意图；

图4是本发明滤波滑动窗组成情况示意图；

图5是本发明一种具体实施方式流程图；

图6是滤波滑动窗所有像素的累积直方图扫描存储示意图；

图7是r+1个常数空间滤波平示意图。

具体实施方式

为更好地理解本发明，下面结合附图和具体实施方式对本发明进行更为详细描述。在以下的描述中，当已有的现有技术的详细描述也许会淡化本发明的主题内容时，这些描述在这儿将被忽略。

本发明意欲提供一种创新技术，在视频图像增强处理系统中，以一种基于滑动滤波窗直方图统计的处理技术来近似标准双边滤波器的处理效果，却能获得处理时间复杂度上的显著优化效果。

为了实现这一快速双边滤波器的创新技术，本发明的基本思路是使用一个矩形滑动窗在待滤波处理图像上滑动，该滑动窗覆盖了双边滤波在空间域上滤波核所覆盖到的图像像素区域。在这个矩形滑动窗内分解出多个重叠的常数空间滤波平面，而同时对每个常数空间滤波平面进行像素密度值的直方图统计，在直方图统计结果的基础上把像素值域的滤波处理以常规的加权求和来实现。因为不必再逐像素的计算像素值域上滤波核系数，因此滤波处理的时间复杂度显著降低。在累积直方图计算的基础上，以逐像素扫描的方式进行滑动窗直方图更新，而不是针对每个滤波处理像素都要构建一个独立的直方图，从而避免了大量的直方图统计计算复杂度。

为达到本发明的目标，实现本发明的核心思想，本发明主要内容的一个方面是提供一种近似标准双边滤波效果的图像平滑滤波方法框架和装置，处理流程如图5所示意。进一步具体来说，标准双边滤波在每个像素位置首先计算滤波核如式(1)，随后进行加权积分计算如式(5)。

I(s₀)＝k(s₀)^-1·∫∫_s∈kernelf_classic(s，s₀)·I(s)ds     (5)

k(s₀)＝∫∫_s∈kernelf_classic(s，s₀)ds                    (6)

式(2)空间滤波核g_spatial(s-s₀)是线性滤波器，它只依赖于参数σ_s。在其方差参数σ_s确定后，空间滤波核g_spatial(s-s₀)可被视为恒定常数。标准的双边滤波计算可以变形为：

I(s₀)＝k^-1·∫∫_s∈kernelg_spatial(s-s₀)·g_tone(I(s)-I(s₀))·I(s)ds

                                                    (7)

＝k^-1·∫∫_s∈kernelker_spatial(s-s₀)·g_tone(I(s)-I(s₀))·I(s)ds

k＝∫∫_s∈kernelker_spatial(s-s₀)·g_tone(I(s)-I(s₀))ds     (8)

ker_spatial(x)是根据空间滤波核g_spatial(s-s₀)计算出来后的滤波参数依据当前点到中心点的距离|s-s₀|进行索引的函数。

设若高斯滤波核宽度为kw＝2*r+1，那么相应的高斯核g_spatial(s-s₀)系数共有r+1个，从中心位置依次到外围系数依次表示为：co0，co1，co2，...，cor。据此计算出c0，c1，c2，...，cr等r+1个空间滤波常数系数。具体计算方法为：

c0＝co0；

c1＝co1-co0；

c2＝co2-c1-c0＝co2-co1；

c3＝co3-c2-c1-c0＝co3-co2；

……

cr＝cor-co(r-1)；                           (9)

得到c0，c1，c2，...，cr等r+1个空间滤波常数系数后，对应了如图7所示意的r+1个常数空间滤波平面，分别表示为S(0)，S(1)，S(2)，…，S(r)。如此，便可把(7)，(8)的双边滤波计算分解为r+1个常数空间滤波的双边计算，如下：

I(s₀)₀＝c0·∫∫_s∈S(0)g_tone(I(s)-I(s₀))·I(s)ds

I(s₀)₁＝c1·∫∫_s∈S(1)g_tone(I(s)-I(s₀))·I(s)ds

……

I(s₀)_r＝cr·∫∫_s∈S(r)g_tone(I(s)-I(s₀))·I(s)ds

I(s₀)＝k^-1·(I(s₀)₀+I(s₀)₁+…+I(s₀)_r)       (10)

相应的分解后归一化算子计算为：

k₀＝c0·∫∫_s∈S(0)g_tone(I(s)-I(s₀))ds

k₁＝c1·∫∫_s∈S(1)g_tone(I(s)-I(s₀))ds

……

k_r＝cr·∫∫_s∈S(r)g_tone(I(s)-I(s₀))ds

k＝k₀+k₁+…+k_r                    (11)

本发明主要内容的另一个方面是提供一种数字图像像素累积直方图计算、扫描、存储与读取方法和装置，如图2、图3以及图6所示意。进一步具体来说，如图2所示在数字图像200上，当前像素为P(x，y)，其左上像素P(x-1，y-1)，上像素为P(x，y-1)，左像素为P(x-1，y)。每个像素都将对应一个累积直方图H，将其定义为按图3所示意的之字形扫描方式，从图像最左上像素到当前像素依据如此扫描所能走过最少的像素集合的像素值直方图统计。设若当前像素P(x，y)紧邻的前置像素(左，上，左上)分别对应的累积直方图为H(x-1，y)，H(x，y-1)，H(x-1，y-1)等。那么当前像素的累积直方图计算方法为：

H(x，y)＝H(x-1，y)+H(x，y-1)-H(x-1，y-1)+H(I(x，y))         (12)

H(I(x，y))表示当前像素灰度值的直方图统计。H(x-1，y)，H(x，y-1)，H(x-1，y-1)分别对应图2所示意中H(x-1，y)，H(x，y-1)，H(x-1，y-1)所指示的几个矩形区域所覆盖像素值的直方图统计。

为了计算随后可能的滑动窗内的像素直方图，需要把如图6中601所示意的从滑动窗左上像素P(x-，y-)到滑动窗右下像素P(x，y)的图像扫描线所覆盖所有像素的累积直方图都存储到如602所示意的链表存储区。602是一个先进先出的单链表存储区，存储方法如图6中602所示意。链表尾存储最老的像素P(x-，y-)的累积直方图，一个滤波核宽度kw的偏移存储当前滑动窗右上角像素P(x，y-)，而链表头存储最新的当前像素P(x，y)的累积直方图，向后一个滤波核宽度kw的偏移存储当前滑动窗左下角像素P(x-，y)。

本发明主要内容的另一个方面是提供一种数字图像滤波核滑动窗建立，滑动，以及滑动窗内图像像素灰度值直方图的计算方法和装置，如图1、图4以及图6所示意。进一步具体来说，滑动窗定义为如图4所示401是大小由双边滤波空间方差参数σ_s所决定的矩形区域，其滤波核区域中心像素P0是当前平滑滤波计算的像素，而滑动窗内最右下像素P(x，y)就是滤波滑动窗当前像素。随着平滑滤波核区域的移动，该滑动窗随之移动。如图1所示中100是当前待滤波数字图像，101表示当前滤波滑动窗，而P(x-，y-)，P(x，y-)，P(x-，y)，P(x，y)分别表示滤波滑动窗内左上、右上、左下、右下四个角点像素。当这四个角点累积直方图H(x-，y-)，H(x，y-)，H(x-，y)，H(x，y)都已知时，那么当前滑动窗所覆盖的区域S像素值直方图的计算方法是：

h_s＝H(x-，y-)-H(x，y-)-H(x-，y)+H(x，y)               (13)

滤波滑动窗在数字图像上的移动方法如图6中601所示意。设若完成当前像素的滤波后，滤波滑动窗随之右移一个像素位置，则图6中像素累积直方图存储链表尾P(x-，y-)被移除，链表尾指针指向其右紧邻存储位置，链表头指针指向刚被置空的存储位置，在该空存储位置存入滑动窗滑动后的右下新像素P(x，y)所对应的累积直方图。

本发明主要内容的另一个方面是提供一种分解多个常数空间滤波平面后的双边滤波计算方法和装置，如图7所示意。进一步具体来说，设若当前滤波核滑动窗所对应的最大的常数空间滤波平面是S(r)，其宽度应为当前像素滤波核宽度kw，该常数空间滤波平面四个顶角所对应的像素分别是P(x-，y-)，P(x，y-)，P(x-，y)，P(x，y)如图1所示意。而常数空间滤波平面S(r-1)，其宽度为kw-2，相应的其四个顶角像素分别应该是P((x-)+1，(y-)+1)，P((x)-1，(y-)+1)，P((x-)+1，(y)-1)，P((x)-1，(y)-1)。同理常数空间滤波平面S(r-2)，其宽度为kw-2*2，相应的其四个顶角像素分别应该是P((x-)+2*1，(y-)+2*1)，P((x)-2*1，(y-)+2*1)，P((x-)+2*1，(y)-2*1)，P((x)-2*1，(y)-2*1)。依次类推，直至常数空间滤波平面S(0)，此时宽度为1，四个顶角像素完全重合，也即P((x-)+r*1，(y-)+r*1)＝P((x)-r*1，(y-)+r*1)＝P((x-)+r*1，(y)-r*1)＝P((x)-r*1，(y)-r*1)。依据本发明前述的矩形区域内图像像素灰度值直方图的计算方法，分解后的r+1个常数空间滤波平面所对应的像素值直方图均能计算出来，如下：

h_S(r)＝H(x-，y-)-H(x，y-)-H(x-，y)+H(x，y)；

h_s(r-1)＝H((x-)+1，(y-)+1)-H((x)-1，(y-)+1)-H((x-)+1，(y)-1)+H((x)-1，(y)-1)；

h_S(r-2)＝H((x-)+2*1，(y-)+2*1)-H((x)-2*1，(y-)+2*1)-H((x-)+2*1，(y)-2*1)+H((x)-2*1，(y)-2*1)；

……

h_S(0)＝H((x)-r*1，(y)-r*1)；               (14)

这里2r+1表示双边滤波核的宽度。随后就可以把r+1个常数空间滤波平面上的双边滤波运算(10)和(11)应用(14)式所计算的r+1个相应的直方图统计进行改写如下：

I(s₀)₀＝c0·∑_i∈S(0)i·h_S(0)(i)·g_tone(i-I(s₀))；

I(s₀)₁＝c1·∑_i∈S(1)i·h_S(1)(i)·g_tone(i-I(s₀))；

I(s₀)_r＝cr·∑_i∈S(r)i·h_S(r)(i)·g_tone(i-I(s₀))；

I(s₀)＝k^-1·(I(s₀)₀+I(s₀)₁+…+I(s₀)_r)            (15)

相应的分解后归一化算子计算为：

k₀＝c0·∑_i∈S(0)h_S(0)(i)·g_tone(i-I(s₀))；

k₁＝c1·∑_i∈S(1)h_S(1)(i)·g_tone(i-I(s₀))；

……

k_r＝cr·∑_i∈S(r)h_S(r)(i)·g_tone(i-I(s₀))；

k＝k₀+k₁+…+k_r                          (16)

这样就实现了标准双边滤波的快速近似。

实施例

本发明的具体实施过程以下步骤进行：

(1)、如图5所示意执行滤波滑动窗定义以及扫描501过程。进一步具体来说，滑动窗定义为如图4所示401是大小由双边滤波空间方差参数σ_s所决定的矩形区域，其滤波核区域中心像素P0是当前平滑滤波计算的像素，而滑动窗内最右下像素P(x，y)就是滤波滑动窗当前像素。随着平滑滤波核区域的移动，该滑动窗随之移动。如图1所示中100是当前待滤波数字图像，101表示当前滤波滑动窗，而P(x-，y-)，P(x，y-)，P(x-，y)，P(x，y)分别表示滤波滑动窗内左上、右上、左下、右下四个角点像素。

(2)、如图5所示意执行图像像素累积直方图统计502过程。进一步具体来说，如图2所示在数字图像200上，当前像素为P(x，y)，其左上像素P(x-1，y-1)，上像素为P(x，y-1)，左像素为P(x-1，y)。每个像素都将对应一个累积直方图H，将其定义为按图3所示意的之字形扫描方式，从图像最左上像素到当前像素依据如此扫描所能走过最少的像素集合的像素值直方图统计。设若当前像素P(x，y)紧邻的前置像素(左，上，左上)分别对应的累积直方图为H(x-1，y)，H(x，y-1)，H(x-1，y-1)等。那么当前像素的累积直方图计算方法为：

H(x，y)＝H(x-1，y)+H(x，y-1)-H(x-1，y-1)+H(I(x，y))       (12)

H(I(x，y))表示当前像素灰度值的直方图统计。H(x-1，y)，H(x，y-1)，H(x-1，y-1)分别对应图2所示意中(x-1，y)，H(x，y-1)，H(x-1，y-1)所指示的几个矩形区域所覆盖像素值的直方图统计。

为了计算随后的滑动窗内的像素直方图，需要把如图6中601所示意的从滑动窗左上像素P(x-，y-)到滑动窗右下像素P(x，y)的图像扫描线所覆盖所有像素的累积直方图都存储到如602所示意的链表存储区。602是一个先进先出的单链表存储区，存储方法如图6中602所示意。链表尾存储最老的像素P(x-，y-)的累积直方图，一个滤波核宽度kw的偏移存储当前滑动窗右上角像素P(x，y-)，而链表头存储最新的当前像素P(x，y)的累积直方图，向后一个滤波核宽度kw的偏移存储当前滑动窗左下角像素P(x-，y)。

(3)、如图5所示意执行图像滤波核滑动窗内像素直方图统计计算503过程。如图1所示中100是当前待滤波数字图像，101表示当前滤波滑动窗，而P(x-，y-)，P(x，y-)，P(x-，y)，P(x，y)分别表示滤波滑动窗内左上、右上、左下、右下四个角点像素。当这四个角点像素的累积直方图H(x-，y-)，H(x，y-)，H(x-，y)，H(x，y)都已知时，那么当前滑动窗所覆盖的区域S像素值直方图的计算方法是：

h_s＝H(x-，y-)-H(x，y-)-H(x-，y)+H(x，y)        (13)

滤波滑动窗在数字图像上的移动方法如图6中601所示意。设若完成当前像素的滤波后，滤波滑动窗随之右移一个像素位置，则图6中像素累积直方图存储链表尾P(x-，y-)被移除，链表尾指针指向其右紧邻存储位置，链表头指针指向刚被置空的存储位置，在该空存储位置存入滑动窗滑动后的右下新像素P(x，y)所对应的累积直方图。

(4)、如图5所示意执行常数空间滤波平面分解504过程，以对标准双边滤波效果的图像平滑滤波进行近似计算。设若高斯滤波核宽度为kw＝2*r+1，那么相应的高斯核g_spatial(s-s₀)系数共有r+1个不同值，从中心位置依次到外围系数依次表示为：co0，co1，co2，...，cor。据此计算出c0，c1，c2，...，cr等r+1个空间滤波常数系数。具体计算方法为：

c0＝co0；

c1＝co1-co0；

c2＝co2-c1-c0＝co2-co1；

c3＝co3-c2-c1-c0＝co3-co2；

……

cr＝cor-co(r-1)；                                  (9)

得到c0，c1，c2，...，cr等r+1个空间滤波常数系数后，对应了如图7所示意的r+1个常数空间滤波平面，分别表示为S(0)，S(1)，S(2)，…，S(r)。

(5)、如图5所示意执行对分解后常数空间滤波平面计算其像素直方图统计505过程，设若当前滤波核滑动窗所对应的最大的常数空间滤波平面是S(r)，其宽度应为当前像素滤波核宽度kw，该常数空间滤波平面四个顶角所对应的像素分别是P(x-，y-)，P(x，y-)，P(x-，y)，P(x，y)如图1所示意。而常数空间滤波平面S(r-1)，其宽度为kw-2，相应的其四个顶角像素分别应该是P((x-)+1，(y-)+1)，P((x)-1，(y-)+1)，P((x-)+1，(y)-1)，P((x)-1，(y)-1)。同理常数空间滤波平面S(r-2)，其宽度为kw-2*2，相应的其四个顶角像素分别应该是P((x-)+2*1，(y-)+2*1)，P((x)-2*1，(y-)+2*1)，P((x-)+2*1，(y)-2*1)，P((x)-2*1，(y)-2*1)。依次类推，直至常数空间滤波平面S(0)，此时宽度为1，四个顶角像素完全重合，也即P((x-)+r*1，(y-)+r*1)＝P((x)-r*1，(y-)+r*1)＝P((x-)+r*1，(y)-r*1)＝P((x)-r*1，(y)-r*1)。依据本发明前述的矩形区域内图像像素灰度值直方图的计算方法，分解后的r+1个常数空间滤波平面所对应的像素值直方图均能计算出来，如下：

h_S(r)＝H(x-，y-)-H(x，y-)-H(x-，y)+H(x，y)；

h_S(r-1)＝H((x-)+1，(y-)+1)-H((x)-1，(y-)+1)-H((x-)+1，(y)-1)+H((x)-1，(y)-1)；

h_S(r-2)＝H((x-)+2*1，(y-)+2*1)-H((x)-2*1，(y-)+2*1)

-H((x-)+2*1，(y)-2*1)+H((x)-2*1，(y)-2*1)；

……

h_S(0)＝H((x)-r*1，(y)-r*1)；                (14)

这里2r+1表示双边滤波核的宽度，h_S(0)，h_S(1)，...，h_S(r-2)，h_S(r-1)，h_S(r)分别表示r+1个常数空间滤波平面所对应的像素值直方图。

(6)、如图5所示意执行对标准双边滤波效果的图像平滑滤波近似计算的506过程，随后就可以把r+1个常数空间滤波平面上的双边滤波运算(10)和(11)应用(14)式所计算的r+1个相应的直方图统计进行改写，用如下计算进行快速近似实现：

I(s₀)₀＝c0·∑_i∈S(0)i·h_S(0)(i)·g_tone(i-I(s₀))；

I(s₀)₁＝c1·∑_i∈S(1)i·h_S(1)(i)·g_tone(i-I(s₀))；

……

I(s₀)_r＝cr·∑_i∈S(r)i·h_S(r)(i)·g_tone(i-I(s₀))；

I(s₀)＝k^-1·(I(s₀)₀+I(s₀)₁+…+I(s₀)_r)            (15)

也就是说先通过求和分别计算出S(0)，S(1)，S(2)，…，S(r)等r+1个常数空间滤波的中间双边计算结果依次是I(s₀)₀，I(s₀)₁，I(s₀)₂，...，I(s₀)_r，随着在由这些中间双边计算结果I(s₀)₀，I(s₀)₁，I(s₀)₂，...，I(s₀)_r计算出图像的最终平滑滤波结果I(s₀)。相应的分解后归一化算子k计算为：

k₀＝c0·∑_i∈S(0)h_S(0)(i)·g_tone(i-I(s₀))；

k₁＝c1·∑_i∈S(1)h_S(1)(i)·g_tone(i-I(s₀))；

……

k_r＝cr·∑_i∈S(r)h_S(r)(i)·g_tone(i-I(s₀))；

k＝k₀+k₁+…+k_r                            (16)

也就是说先通过求和分别计算出S(0)，S(1)，S(2)，…，S(r)等r+1个常数空间滤波所对应的归一化算子k₀，k₁，…，k_r，然后由之求和得到归一化算子k。

这样就实现了标准双边滤波的快速近似实现。

在本实施例中，本发明的显著特点是提供了一个近似标准双边滤波效果的图像平滑滤波方法框架和系统装置。进一步来说使用一个矩形滑动窗在待滤波处理图像上滑动，该滑动窗覆盖了双边滤波在空间域上滤波核所覆盖到的图像像素区域。在这个矩形滑动窗内分解出多个重叠的常数空间滤波平面，而同时对每个常数空间滤波平面进行像素密度值的直方图统计，在直方图统计结果的基础上把像素值域的滤波处理以常规的加权求和来实现。

在本实施例中，本发明的另一个特点是提供一种数字图像像素累积直方图计算、扫描、存储与读取的方法和装置。该方法与装置能利用图像滤波的扫描顺序方便的依次计算出扫描线上像素点的累积直方图，而且使用一个固定长度的链表存储必要的一段扫描线上像素的累积直方图数据，就能满足随后计算的需要。

在本实施例中，本发明的再一个特点是提供一种分解多个常数空间滤波平面后的双边滤波计算方法和装置。在分解后的多个常数空间滤波平面上可以便利计算出相应的像素值直方图，随后可以方便的计算出近似标准双边滤波的滤波后像素值。

因此，尽管上面对本发明说明性的具体实施方式进行了描述，但应当清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (4)

1、一种快速视频图像双边滤波的方法，其特征在于，包括以下步骤：

(1)、定义一个大小由双边滤波空间方差参数σ_s决定的矩形区域作为滤波滑动窗，其中心像素P0是当前平滑滤波计算的像素；

(2)、计算图像像素累积直方图；

(3)、计算当前滤波滑动窗内像素直方图，方法为：

h_s＝H(x-，y-)-H(x，y-)-H(x-，y)+H(x，y)

其中，P(x-，y-)，P(x，y-)，P(x-，y)，P(x，y)分别表示滤波滑动窗内左上、右上、左下、右下四个角点像素，其对应的累积直方图为H(x-，y-)，H(x，y-)，H(x-，y)，H(x，y)，根据步骤(2)的计算结果可以得到；

(4)、将当前滤波滑动窗执行常数空间滤波平面分解

当前滤波滑动窗的宽度，即高斯滤波核宽度为kw＝2*r+1，相应的高斯核g_spatial(s-s₀)系数共有r+1个不同值，从内到外系数依次表示为：co0，co1，co2，...，cor，并依据以下公式计算出c0，c1，c2，...，cr等r+1个空间滤波常数系数：

c0＝co0；

c1＝co1-co0；

c2＝co2-c1-c0＝co2-co1；

c3＝co3-c2-c1-c0＝co3-co2；

......

cr＝cor-co(r-1)；

得到c0，c1，c2，...，cr等r+1个空间滤波常数系数后，将r+1个常数空间滤波平面分别表示为S(0)，S(1)，S(2)，…，S(r)；

(5)、对分解后的常数空间滤波平面计算其像素值直方图

按步骤(3)所述的方法计算分解后的r+1个常数空间滤波平面像素值直方图：

h_S(r)＝H(x-，y-)-H(x，y-)-H(x-，y)+H(x，y)；

h_S(r-1)＝H((x-)+1，(y-)+1)-H((x)-1，(y-)+1)-H((x-)+1，(y)-1)+H((x)-1，(y)-1)；

h_S(r-2)＝H((x-)+2*1，(y-)+2*1)-H((x)-2*1，(y-)+2*1)

         -H((x-)+2*1，(y)-2*1)+H((x)-2*1，(y)-2*1)；

.....

h_S(0)＝H((x)-r*1，(y)-r*1)；

这里2r+1表示双边滤波核的宽度，h_S(0)，h_S(1)，...，h_S(r-2)，h_S(r-1)，h_S(r)分别表示r+1个常数空间滤波平面所对应的像素值直方图。

(6)、对标准双边滤波的图像平滑滤波近似计算：

通过求和分别计算出S(0)，S(1)，S(2)，…，S(r)等r+1个常数空间滤波的中间双边计算结果依次是I(s₀)₀，I(s₀)₁，I(s₀)₂，...，I(s₀)_r，随着在由这些中间双边计算结果I(s₀)₀，I(s₀)₁，I(s₀)₂，...，I(s₀)_r计算出图像的最终平滑滤波结果I(s₀)，

I(s₀)＝k^-1·(I(s₀)₀+I(s₀)₁+…+I(s₀)_r)

其中中间双边计算结果为：

I(s₀)₀＝c0·Σ_i∈S(0)i·h_S(0)(i)·g_tone(i-I(s₀))；

I(s₀)₁＝c1·Σ_i∈S(1)i·h_S(1)(i)·gtone(i-I(s₀))；

.....

I(s₀)_r＝cr·Σ_i∈S(r)i·h_S(r)(i)·g_tone(i-I(s₀))；

归一化算子为

k＝k₀+k₁+…+k_r

k₀＝c0·Σ_i∈S(0)h_S(0)(i)·g_tone(i-I(s₀))；

k₁＝c1·Σ_i∈S(1)h_S(1)(i)·g_tone(i-I(s₀))；

.....

k_r＝cr·Σ_i∈S(r)h_S(r)(i)·g_tone(i-I(s₀))；。

2、根据权利要求1所述的快速视频图像双边滤波的方法，其特征在于，所述的图像像素累积直方图计算为：

H(x，y)＝H(x-1，y)+H(x，y-1)-H(x-1，y-1)+H(I(x，y))

其中，H(I(x，y))表示当前像素灰度值的直方图统计，H(x-1，y)、H(x，y-1)、H(x-1，y-1)之字形扫描方式，从图像最左上像素到当前像素依据如此扫描所能走过最少的像素集合的像素值直方图统计。

3、根据权利要求1所述的快速视频图像双边滤波的方法，其特征在于，所述的步骤(3)计算当前滤波滑动窗内像素直方图时，须需要从滑动窗左上像素P(x-，y-)到滑动窗右下像素P(x，y)的图像扫描线所覆盖所有像素的累积直方图都存储到链表存储区，链表存储区是一个先进先出的单链表存储区，链表尾存储最老的像素P(x-，y-)的累积直方图，一个滤波核宽度kw的偏移存储当前滑动窗右上角像素P(x，y-)，而链表头存储最新的当前像素P(x，y)的累积直方图，向后一个滤波核宽度kw的偏移存储当前滑动窗左下角像素P(x-，y)。

Images