一种基于亮度分层的快速三边滤波器色调映射方法

摘要

本发明涉及一种基于亮度分层的快速三边滤波器色调映射方法，属于数字图像处理领域。本发明在现有三边滤波器色调映射方法的基础上，仅对获得平滑后的水平方向梯度图像L′x(p)和经平滑后的垂直方向梯度图像L′y(p)的方法进行了改进，即：通过对水平方向梯度图像Lx(p)的像素值进行分级，然后根据Wk(p)＝s(Lx(p)-Ik)和Jk(p)＝Wk(p)·Ix(p)分别获得图像Wk(p)和图像Jk(p)；再利用图像Wk(p)和图像Jk(p)对水平方向梯度图像Lx(p)进行平滑处理，得到平滑后的水平方向梯度图像L′x(p)；使用相同方法，得到平滑后的垂直方向梯度图像L′y(p)。本发明方法与已有的三边滤波器色调映射方法相比较，具有处理速度快，但不损失图像细节的优点。

说明书

技术领域

本发明涉及一种高动态范围图像(HDRI-high dynamic range image)的色调映射技术，特别涉及一种基于亮度分层的快速三边滤波器色调映射方法，属于数字图像处理领域。

背景技术

动态范围是指图像中像素最大亮度与最小亮度的比值。传统的RGB图像格式使用8位二进制整数表示每个颜色分量，其动态范围为100∶1，而真实场景的动态范围远远超出这一范围。传统图像格式无法真实地表达动态范围很高的场景。高动态范围(High Dynamic Range，HDR)图像是针对这一问题所设计的新型图像格式，通常使用浮点数来表示每个颜色分量，这样能够表达更高的动态范围。

色调映射是指将一幅高动态范围图像转换为低动态范围图像的方法。它提供了一种将现实场景的亮度值映射到显示设备能显示的范围的方法，除了压缩亮度范围之外，它还必须充分保留原始图像的感观质量，比如保留对比度、图像细节等对于表现原始场景非常重要的信息。色调映射算法是高动态范围图像研究领域中最为关键的部分。其主要原因在于目前占主导地位的有限动态范围媒介，譬如CRT和LCD显示器或者投影仪等，只能显示有限动态范围的图像，使人们不能观察高动态范围图像的全部信息。

现有的色调映射技术可分为两大类：全局方法和局部方法。全局方法是将一个映射函数应用于图像中的所有像素，而忽略像素的在图像中的位置，例如直方图调整法、自适应指数映射法等。全局方法的优点在于整体架构简洁和运算效率高，缺点是会导致图像结果在细节、颜色和对比度损失。局部方法中某个像素的映射与其在图像中的位置相关，或者说是与其自身亮度及周围像素的亮度相关，例如双边滤波器法、三边滤波器法等。局部方法的特点是能够较好地保留图像细节，但速度慢。其中三边滤波器是一种效果非常好的方法。

P.Choudhury等人在文献《The trilateral filter for high contrast images and meshes》(P.Choudhury and J.Tumblin，Proceedings ofEurographics Symposium on Rendering，pp.1-11，2003.)中介绍了三边滤波器色调映射方法，其具体步骤为：

第1步：输入一幅高动态范围图像。

第2步：根据该图像的红色通道R、绿色通道G、蓝色通道B的通道值，获得该图像中每个像素p的对数域亮度L(p)。

第3步：根据前向差分法计算图像中每个像素p的对数域亮度L(p)的水平方向梯度图像L_x(p)和垂直方向梯度图像L_y(p)。

第4步：根据公式1和公式2计算经平滑后的水平方向梯度图像L′_x(p)。

$L_x^{\prime }\left(p\right)=\frac{1}{k_x\left(p\right)}\underset{\zeta \in N\left(p\right)}{\Sigma }c\left(\zeta \right)s(L_x(\zeta +p)-L\left(p\right))L_x\left(p\right)---\left(1\right)$

$k_x\left(p\right)=\underset{\zeta \in N\left(p\right)}{\Sigma }c\left(\zeta \right)s(L_x(\zeta +p)-L\left(p\right))---\left(2\right)$

其中，c为一高斯函数，s为一高斯函数；高斯函数c和高斯函数s的区别在于其均值和方差根据具体情况取不同值；N(p)表示以像素p为中心的n×n大小的区域，其中n≥3且n为正整数；ζ表示N(p)区域内的每一个像素。

第5步：根据公式3和公式4计算经平滑后的垂直方向梯度图像L′_y(p)，其方法是：

$L_y^{\prime }\left(p\right)=\frac{1}{k_y\left(p\right)}\underset{\zeta \in N\left(p\right)}{\Sigma }c\left(\zeta \right)s(L_y(\zeta +p)-L\left(\zeta \right))L_y\left(p\right)---\left(3\right)$

$k_y\left(p\right)=\underset{\zeta \in N\left(p\right)}{\Sigma }c\left(\zeta \right)s(L_y(\zeta +p)-L\left(p\right))---\left(4\right)$

第6步：对于每一个像素p及其周围n×n大小区域内的每一个像素ζ，根据公式5计算局部细节值L_Δ(p，ζ)。

L_Δ(p，ζ)＝L(ζ)-L(p)-(L′_x(p)，L′_y(p))·(ζ-p)               (5)

第7步：根据公式6计算像素p经三边滤波的结果L_out(p)，L_out(p)即为图像中每个像素p的对数域亮度L(p)的基本层图像。

$I_{\mathrm{out}}\left(p\right)=L\left(p\right)+\frac{1}{k_{\Delta }\left(x\right)}\underset{\zeta \in N\left(p\right)}{\Sigma }{L}_{\Delta }(p,\zeta )c(\zeta -p)s\left(L_{\Delta }(p,\zeta )\right)f(p,\zeta )---\left(6\right)$

$k_{\Delta }\left(x\right)=\underset{\zeta \in N\left(p\right)}{\Sigma }c(\zeta -p)s\left(L_{\Delta }(p,\zeta )\right)f(p,\zeta )---\left(7\right)$

其中，f(p，ζ)为取值为0或1的二值函数，其定义如公式8所示。

$f(p,\zeta )=\left(\begin{array}{cc}1& (\mathrm{if}\sqrt{{(L_x^{\prime }\left(\zeta \right)-L_x^{\prime }\left(p\right))}^2+{(L_y^{\prime }\left(\zeta \right)-L_y^{\prime }\left(p\right))}^2}<R)\\ 0& \left(\mathrm{otherwise}\right)\end{array}\right)---\left(8\right)$

R为一个阈值，其通过公式9来确定。

R＝0.15(MaxG(L)-MinG(L))                                    (9)

其中，MaxG(L)和MinG(L)分别为图像中每个像素p的对数域亮度L(p)中的最大梯度幅值和最小梯度幅值。

第8步：通过公式10计算L(p)的细节层图像L_detail(p)。

L_detail(p)＝L(p)-L_out(p)                                    (10)

第9步：通过公式11计算经动态范围压缩后的基本层图像L_base(p)。

L_base(p)＝L_out(p)/factor                                    (11)

其中，factor为压缩比例因子，取值为大于1的正实数。

第10步：通过公式12合成经动态范围压缩后的对数域亮度图像L_mapped(p)。

L_mapped(p)＝L_base(p)+L_detail(p)                             (12)

第11步：使用公式13～15从L_mapped(p)还原出彩色图像I_out(p)的三个颜色通道值。

I_{out_R}(p)＝exp(L_mapped(p))·I_{in_R}(p)/exp(L(p))              (13)

I_{out_G}(p)＝exp(L_mapped(p))·I_{in_G}(p)/exp(L(p))              (14)

I_{out_B}(p)＝exp(L_mapped(p))·I_{in_B}(p)/exp(L(p))              (15)

其中，exp为指数函数，其底数与步骤2中所取对数的底数相同，I_{out_R}(p)是指I_out(p)的红色通道值，I_{in_R}(p)是指输入图像的红色通道值，I_{out_G}(p)是指I_out(p)的绿色通道值，I_{in_G}(p)是指输入图像的绿色通道值，I_{out_B}(p)是指I_out(p)的蓝色通道值，I_{in_B}(p)是指输入图像的蓝色通道值。

第12步：将第11步中的颜色通道值线性平移并缩放至0到255的区间，并四舍五入取整，再合并成一幅新的图像，即为输入的高动态范围图像经色调映射后的结果。

该方法的优点是：处理效果非常好；但其缺点也显而易见：处理速度非常慢；在处理一幅图像时，往往需要好几分钟或者好几十分钟的计算时间。

发明内容

本发明的目的是克服已有三边滤波器色调映射方法存在的处理速度慢的缺陷，提出一种基于亮度分层的快速三边滤波器色调映射方法。

本发明的目的是通过下述技术方案实现的。

一种基于亮度分层的快速三边滤波器色调映射方法，其具体操作步骤为：

第1步：输入一幅高动态范围图像。

第2步：根据该图像的红色通道R、绿色通道G、蓝色通道B的通道值，获得该图像中每个像素p的对数域亮度L(p)。

第3步：利用前向差分法获得每个像素p的对数域亮度L(p)的水平方向梯度图像L_x(p)和垂直方向梯度图像L_y(p)。

第4步：获得经平滑后的水平方向梯度图像L′_x(p)；具体为：

步骤①：找出第3步得到的水平方向梯度图像L_x(p)中的最大像素值I_max和最小像素值I_min，把从I_min到I_max的区间等分为m个亮度等级，分别记为I₀、I₁、......、I_m-1；其中，m≥4且m为正整数。

步骤②：依次对第k(k＝1，2，......，m且k为正整数)个亮度等级，利用第3步得到的水平方向梯度图像L_x(p)，分别根据公式16和公式17得到图像W_k(p)和图像J_k(p)；

W_k(p)＝s(L_x(p)-I_k)                                             (16)

J_k(p)＝W_k(p)·L_x(p)                                            (17)

其中：s为一高斯函数。

步骤③：对第3步得到的水平方向梯度图像L_x(p)进行平滑处理，得到平滑后的水平方向梯度图像L′_x(p)；具体为：

如果L_x(p)∈[I_t，I_t+1]，1≤t≤m-2且t为正整数；则：

$L_x^{\prime }\left(p\right)=\frac{I_{t+1}-L_x\left(p\right)}{k_t\left(p\right)(I_{t+1}-I_t)}\underset{\zeta \in N\left(p\right)}{\Sigma }c(\zeta -p)J_t\left(\zeta \right)+\frac{L_x\left(p\right)-I_t}{k_{t+1}\left(x\right)(I_{t+1}-I_t)}\underset{\zeta \in N\left(p\right)}{\Sigma }c(\zeta -p)J_{t+1}\left(\zeta \right)---\left(18\right)$

$k_t\left(p\right)=\underset{p\in N\left(p\right)}{\Sigma }c(\zeta -p)W_t\left(p\right)---\left(19\right)$

$k_{t+1}\left(p\right)=\underset{p\in N\left(p\right)}{\Sigma }c(\zeta -p)W_{t+1}\left(p\right)---\left(20\right)$

其中，c为一高斯函数，N(p)表示以像素p为中心的n×n大小的区域，其中n≥3且n为正整数；ζ表示N(p)区域内的每一个像素。

经过第4步的操作，即可获得经平滑后的水平方向梯度图像L′_x(p)。

第5步：获得经平滑后的垂直方向梯度图像L′_y(p)；具体为：

步骤①：找出第3步得到的垂直方向梯度图像L_y(p)中的最大像素值I′_max和最小像素值I′_min，把从I′_min到I′_max的区间等分为m′个亮度等级，分别记为I′₀、I′₁、......、I′_m′-1；其中，m′≥4且m′为正整数。

步骤②：依次对于第k′(k′＝1，2，......，m′且k′为正整数)个亮度等级，利用第3步得到的垂直方向梯度图像L_y(p)，分别根据公式21和公式22得到图像W′_k′(p)和图像J′_k′(p)；

W′_k′(p)＝s(L_y(p)-I′_k′)                                        (21)

J′_k(p)＝W′_k′(p)·L_y(p)                                         (22)

步骤③：对第3步得到的垂直方向梯度图像L_y(p)进行平滑处理，得到平滑后的垂直方向梯度图像L′_y(p)；具体为：

如果L_y(p)∈[I′_t′，I′_t′+1]，其中1≤t′≤m′-2且t′为正整数；则：

$L_y^{\prime }\left(p\right)=\frac{{I^{\prime }}_{t^{\prime }+1}-L_y\left(p\right)}{{k^{\prime }}_{t^{\prime }}\left(p\right)({I^{\prime }}_{t^{\prime }+1}-{I^{\prime }}_{t^{\prime }})}\underset{\zeta \in N\left(p\right)}{\Sigma }c(\zeta -p){J^{\prime }}_{t^{\prime }}\left(\zeta \right)+\frac{{L^{\prime }}_{x^{\prime }}\left(p\right)-{I^{\prime }}_{t^{\prime }}}{{k^{\prime }}_{t^{\prime }+1}\left(x\right)({I^{\prime }}_{t^{\prime }+1}-{I^{\prime }}_{t^{\prime }})}\underset{\zeta \in N\left(p\right)}{\Sigma }c(\zeta -p)J_{t+1}\left(\zeta \right)---\left(23\right)$

${k^{\prime }}_{t^{\prime }}\left(p\right)=\underset{p\in N\left(p\right)}{\Sigma }c(\zeta -p){W^{\prime }}_{t^{\prime }}\left(p\right)---\left(24\right)$

${k^{\prime }}_{t^{\prime }+1}\left(p\right)=\underset{p\in N\left(p\right)}{\Sigma }c(\zeta -p){W^{\prime }}_{t^{\prime }+1}\left(p\right)---\left(25\right)$

经过第5步的操作，即可获得经平滑后的垂直方向梯度图像L′_y(p)。

第6步：对于每一个像素p及其周围n×n大小区域内的每一个像素ζ，计算局部细节值L_Δ(p，ζ)；

第7步：获得每一个像素p经三边滤波的结果L_out(p)，L_out(p)即为L(p)的基本层图像；

第8步：获得L(p)的细节层图像L_detail(p)；

第9步：获得经动态范围压缩后的基本层图像I_base(p)；

第10步：合成经动态范围压缩后的对数域亮度图像L_mapped(p)；

第11步：从L_mapped(p)还原出彩色图像的三个颜色通道值；

第12步：将第11步中的颜色通道值平移并缩放至0到255的区间，并四舍五入取整，再合并成一幅新的图像，即为输入的高动态范围图像经色调映射后的结果。

有益效果

本发明方法与已有的三边滤波器色调映射方法相比较，具有处理速度快，但不损失图像细节的优点。

附图说明

图1为本发明具体实施方式中采用的原始图像；

图2为本发明具体实施方式中采用本发明方法对图1处理过程中得到的基本层图像；

图3为本发明具体实施方式中采用本发明方法对图1处理过程中得到的细节层图像；

图4为本发明具体实施方式中采用本发明方法对图1处理后的输出结果；

图5为采用已有的三边滤波器色调映射方法对图1处理后的输出结果；

图6为本发明具体实施方式中采用6幅多种场景类型输入图像的原始图像a～f。

具体实施方式

下面结合具体实施例对本发明进行详细说明。实施例中，采用本发明方法处理一幅如1所示的高动态范围图像，其处理过程如下：

第1步：输入一幅高动态范围图像，如图1所示。

第2步：根据该图像的红色通道R、绿色通道G、蓝色通道B的通道值，获得该图像中每个像素p的对数域亮度L(p)。

第3步：利用前向差分法获得每个像素p的对数域亮度L(p)的水平方向梯度图像L_x(p)和垂直方向梯度图像L_y(p)。

L_x(p)＝L(p_right)-L(p)                                        (26)

L_y(p)＝L(p_below)-L(p)                                        (27)

其中p_right表示p右侧的像素，p_below表示p下方的像素。

第4步：获得经平滑后的水平方向梯度图像L′_x(p)；具体为：

步骤①：找出第3步得到的水平方向梯度图像L_x(p)中的最大像素值I_max和最小像素值I_min，把从I_min到I_max的区间等分为10个亮度等级，分别记为I₀、I₁、......、I₉；。

步骤②：依次对于第k个亮度等级，利用第3步得到的水平方向梯度图像L_x(p)，分别根据公式16和公式17得到图像W_k(p)和图像J_k(p)；其中1≤k≤10，公式16中的高斯函数s的中值为0，标准差为公式9得到的R值。

步骤③：根据公式18、公式19和公式20对第3步得到的水平方向梯度图像L_x(p)进行平滑处理，得到平滑后的水平方向梯度图像L′_x(p)；具体为：

如果L_x(p)∈[I_t，I_t+1]，其中1≤t≤8且t为正整数；则：

$L_x^{\prime }\left(p\right)=\frac{I_{t+1}-L_x\left(p\right)}{k_t\left(p\right)(I_{t+1}-I_t)}\underset{\zeta \in N\left(p\right)}{\Sigma }c(\zeta -p)J_t\left(\zeta \right)+\frac{L_x\left(p\right)-I_t}{k_{t+1}\left(x\right)(I_{t+1}-I_t)}\underset{\zeta \in N\left(p\right)}{\Sigma }c(\zeta -p)J_{t+1}\left(\zeta \right)---\left(18\right)$

$k_t\left(p\right)=\underset{p\in N\left(p\right)}{\Sigma }c(\zeta -p)W_t\left(p\right)---\left(19\right)$

$k_{t+1}\left(p\right)=\underset{p\in N\left(p\right)}{\Sigma }c(\zeta -p)W_{t+1}\left(p\right)---\left(20\right)$

其中，高斯函数c的中值为0，标准差为21，N(p)表示以像素p为中心的5×5大小的区域，ζ表示N(p)区域内的每一个像素。

经过第4步的操作，即可获得经平滑后的水平方向梯度图像L′_x(p)。

第5步：获得经平滑后的垂直方向梯度图像L′_y(p)；具体为：

步骤①：找出第3步得到的垂直方向梯度图像L_y(p)中的最大像素值I′_max和最小像素值I′_min，把从I′_min到I′_max的区间等分为10个亮度等级，分别记为I′₀、I′₁、......、I′₉。

步骤②：依次对于第k′个亮度等级，利用第3步得到的垂直方向梯度图像L_y(p)，分别根据公式21和公式22得到图像W′_k′(p)和图像J′_k′(p)；其中1≤k′≤10，公式21中的高斯函数s的中值为0，标准差为公式9得到的中的R值。

步骤③：根据公式23、公式24和公式25对第3步得到的垂直方向梯度图像L_y(p)进行平滑处理，得到平滑后的垂直方向梯度图像L′_y(p)；具体为：

如果L_y(p)∈[I′_t′，I′_t′+1]，其中1≤t′≤8且t′为正整数；则：

$L_y^{\prime }\left(p\right)=\frac{{I^{\prime }}_{t^{\prime }+1}-L_y\left(p\right)}{{k^{\prime }}_{t^{\prime }}\left(p\right)({I^{\prime }}_{t^{\prime }+1}-{I^{\prime }}_{t^{\prime }})}\underset{\zeta \in N\left(p\right)}{\Sigma }c(\zeta -p){J^{\prime }}_{t^{\prime }}\left(\zeta \right)+\frac{{L^{\prime }}_{x^{\prime }}\left(p\right)-{I^{\prime }}_{t^{\prime }}}{{k^{\prime }}_{t^{\prime }+1}\left(x\right)({I^{\prime }}_{t^{\prime }+1}-{I^{\prime }}_{t^{\prime }})}\underset{\zeta \in N\left(p\right)}{\Sigma }c(\zeta -p)J_{t+1}\left(\zeta \right)---\left(23\right)$

${k^{\prime }}_{t^{\prime }}\left(p\right)=\underset{p\in N\left(p\right)}{\Sigma }c(\zeta -p){W^{\prime }}_{t^{\prime }}\left(p\right)---\left(24\right)$

${k^{\prime }}_{t^{\prime }+1}\left(p\right)=\underset{p\in N\left(p\right)}{\Sigma }c(\zeta -p){W^{\prime }}_{t^{\prime }+1}\left(p\right)---\left(25\right)$

其中，高斯函数c的中值为0，标准差为21，N(p)表示以像素p为中心的5×5大小的区域，ζ表示N(p)区域内的每一个像素。

经过第5步的操作，即可获得经平滑后的垂直方向梯度图像L′_y(p)。

第6步：根据公式5对于每一个像素p及其周围5×5大小区域内的每一个像素ζ，计算局部细节值L_Δ(p，ζ)；

第7步：根据公式6获得每一个像素p经三边滤波的结果L_out(p)，L_out(p)即为L(p)的基本层图像，如图2所示；

第8步：根据公式10获得L(p)的细节层图像L_detail(p)，如图3所示；

第9步：根据公式11获得经动态范围压缩后的基本层图像L_base(p)，factor取值为5；

第10步：根据公式12合成经动态范围压缩后的对数域亮度图像L_mapped(p)；

第11步：根据公式13、公式14、公式15从L_mapped(p)还原出彩色图像的三个颜色通道值；

第12步：将第11步中的颜色通道值平移并缩放至0到255的区间，并四舍五入取整，再合并成一幅新的图像，即为输入的高动态范围图像经色调映射后的结果，如图4所示。

采用已有的三边滤波器色调映射方法对图1进行处理，其结果如图5所示，图4与图5的比较说明本发明方法与已有的三边滤波器色调映射方法的结果非常接近，但本发明方法速度有明显提高。

为了更好地说明本发明效果，分别采用对图6中的6幅场景本发明方法和已有三边滤波器色调映射方法类型输入图像进行处理，得到本发明方法与已有三边滤波器色调映射方法对图6中6幅图片进行处理所花费时间的对比，如表1所示。

表1本发明方法和已有方法对图6处理所花费时间的对比表

从表1中可以看出，本发明方法的处理速度比已有的三边滤波器色调映射方法快很多。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进，或者对其中部分技术特征进行等同替换，这些改进和替换也应视为本发明的保护范围。