一种自适应的图像细节增强方法

注意：这篇文章上次更新于1247天前，文章内容可能已经过时。

算法简介

算法的核心思想为：使用 DoG 模型提取图像边缘信息并与原图相加来获得细节增强的效果^[1]。

$output(x,y) = input(x,y) + (input * DoG_{\sigma_c,\sigma_s,k})(x,y)$

其中：

$DoG_{\sigma_c,\sigma_s,k} = max(0,\frac{1}{2\pi\sigma_c^2}exp(-\frac{x^2+y^2}{2\sigma_c^2}) - k\times \frac{1}{2\pi\sigma_s^2}exp(-\frac{x^2+y^2}{2\sigma_s^2}))$

为了能够在增强细节的同时不放大噪声，依据图像的梯度信息对 k 的取值进行自适应化。图像的梯度越大时，k 的取值越小。理论上可以实现对图像的每一个像素都使用不同 k 值的 DoG 模型进行运算，但这会使得运算的时间代价特别大。因此，我们对梯度图量化到 N ，这样只需要 N 个DoG 模型作用到原图像上，可有效减少计算的时间花费^[2]。

首先，梯度图可以使用如下公式进行计算：

$G(x) = \sqrt{G_h(x)^2 + G_v(x)^2}$

其中，$ 和 $ 分别为输入图像使用 $3 \times 3$ 的 Sobel 算子计算的结果。我们假设在一个局部区域内图像的细节丰富程度是相同的，因此，在对梯度图进行量化操作之前，首先对梯度图进行扩张（dilate）操作^[3]。

我们通过对扩张后的梯度图按其强度所在的区间进行标记来实现量化操作，例如，当梯度图的强度在 [0 , 0.2] 范围内时，将该点标记为 1 ；当梯度图的强度在 [0.2 , 0.4] 范围内时，将其标记为 2 …… 以此类推。接下来对输入图像进行 DoG 滤波，当该像素点的标记为 1 时，此处 DoG 模型的 k 值取 1 ；当该像素点的标记为 2 时，此处 DoG 模型的 k 值取 0.95 …… 当该像素点的标记为 5 时，此处的 DoG 模型的 k 值取 0.8。

最后，将自适应 DoG 模型提取的细节图与原图相加，在不放大噪声的情况下实现细节增强。

PS：对于像素值超过 1 的情况，可以使用线性拉伸，非线性拉伸，裁剪等操作将其缩放到 [0 , 1]

Matlab实现

二维非对称高斯函数

function gauss = GaussFun2D(xsigma,ysigma)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 二维非对称高斯函数
% 参数列表：
%    xsigma:水平方向标准差 
%    ysigma:竖直方向标准差（缺省值为xsigma）
% 高斯核的直径为标准差的 6 倍
% GaussFun2D(0.5) 等效于 fspecial('gaussian',[3,3],0.5)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin == 1
    ysigma = xsigma;
end

xsize = xsigma * 6;
ysize = ysigma * 6;

if mod(round(xsize),2) == 0
    if xsize < round(xsize)
        xsize = floor(xsize);
    elseif xsize > round(xsize)
        xsize = ceil(xsize);
    else
        xsize = xsize + 1;
    end
else
    xsize = round(xsize);
end
if mod(round(ysize),2) == 0
    if ysize < round(ysize)
        ysize = floor(ysize);
    elseif ysize > round(ysize)
        ysize = ceil(ysize);
    else
        ysize = ysize + 1;
    end
else
    ysize = round(ysize);
end

gauss = zeros(xsize,ysize);

x_center = (xsize + 1)/2;
y_center = (ysize + 1)/2;
A = 1/(2*pi*xsigma*ysigma);
for i = 1:xsize
    for j = 1:ysize
        gauss(i,j) = A*exp(-0.5*( ((i-x_center)/xsigma).^2 + ((j-y_center)/ysigma).^2  ));
    end
end

gauss = gauss ./ sum(gauss(:));
end

DoG 模型

function output = MyDoG(input,sigma_c,sigma_s,k)
    gauss_c = MyGaussFun2D(sigma_c);
    gauss_s = MyGaussFun2D(sigma_s);
    c = imfilter(input,gauss_c);
    s = imfilter(input,gauss_s);
    output = max(0,c - k .* s);
end

梯度图计算

function [Grad] = getGrad(I, win)
r0 = 0;
r1 = 1;
Strch = @(x) (x-min(x(:))).*((r1-r0)/(max(x(:))-min(x(:)))) + r0;
I = im2double(I);
gray = rgb2gray(I);
[Gmag, ~] = imgradient(gray, 'sobel');

radius = win;
se = strel('square', radius); % ???

imG = imdilate(Gmag, se); % ???
Grad = Strch(imfill(imG, 8, 'holes'));
end

量化操作

function label = Quantisation(input,N)

[m,n] = size(input);
label = ones(m,n);
levels = [0:1/N:1];

for i = 1:N
    for j = 1:m*n
        if input(j) > levels(i) && input(j) <= levels(i+1)
            label(j) = i;
        end
    end
end
label = int32(label);
end

数据归一化

function [outputs] = Normalize(arg,str)
if nargin == 1
    str = 'Line';
end
if strcmp(str,'Line')
MAX = max(arg(:));
MIN = min(arg(:));
outputs = (arg - MIN)./ (MAX - MIN);
end
if strcmp(str,'sigmoid')
    u = mean(arg(:));
    s = std(arg(:));
    tmp = (arg - u)./ s;
    outputs = 1 ./ (1 + exp(-tmp));  
end
if strcmp(str,'log')
    alpha = max(arg(:)) + 1;
    outputs = log(arg + 1) ./ log(alpha);
end
end

细节增强

function output = DetailEnhancement(input)
%% 参数
N = 5;
k = 1:-0.05:0.8;
win = 15;

%%
I = im2double(input);
grad = GetGrad(I,win);
label = Quantisation(grad,N);

for i = 1:N
    dog{i} = MyDoG(I,0.5,1.5,k(i));
end
dog_out = zeros(size(I));
for i = 1:size(I,1)
    for j = 1:size(I,2)
        dog_out(i,j,:) = dog{label(i,j)}(i,j,:);
    end
end
output = I + dog_out;

output = Normalize(output,'sigmoid');
% output = Normalize(output,'Line');
% output = Normalize(output,'log');
% output = min(1,output);
end