图像分析之梯度L0范数平滑

本文是Image Smoothing via L0 Gradient Minimization一文的笔记。L0 Gradient Smoothing的formulation与TV和WLS等基于变分的模型很相似，所以本文重在推导。读者需注意，本文采用的符号标记与原论文不同，笔者觉得本文采用的符号标记表达力更强些，且不容易产生歧义。本文重写了原论文中的问题描述，推导了原论文中的公式（8），笔者还推导了一个新的向量形式的Solver，并编码验证了该Solver的正确性（遗憾的是，效率不及原作者用FFT实现的代码）。

更新记录

本文持续更新！如文中有错误，或你对本文有疑问或建议，欢迎留言或发邮件至quarrying#qq.com！

2016年01月12日，发布博文。

2016年01月13日，修改代码二，效果与代码一一致。感谢网友guide（QQ昵称）指出错误。

参考

[1] L. Xu, C. Lu, Y. Xu, and J. Jia, “Image smoothing via L0 gradient minimization,” Proc. 2011 SIGGRAPH Asia Conf. - SA ’11, vol. 30, no. 6, p. 1, 2011.

相关代码

代码一，这是原作者采用FFT实现的L0 Gradient Smoothing

%   Distribution code Version 1.0 -- 09/23/2011 by Jiaya Jia Copyright 2011, The Chinese University of Hong Kong.
%
%   The Code is created based on the method described in the following paper 
%   [1] "Image Smoothing via L0 Gradient Minimization", Li Xu, Cewu Lu, Yi Xu, Jiaya Jia, ACM Transactions on Graphics, 
%   (SIGGRAPH Asia 2011), 2011. 
%  
%   The code and the algorithm are for non-commercial use only.


function S = L0Smoothing(Im, lambda, kappa)
%L0Smooth - Image Smoothing via L0 Gradient Minimization
%   S = L0Smoothing(Im, lambda, kappa) performs L0 graidient smoothing of input
%   image Im, with smoothness weight lambda and rate kappa.
%
%   Paras: 
%   @Im    : Input UINT8 image, both grayscale and color images are acceptable.
%   @lambda: Smoothing parameter controlling the degree of smooth. (See [1]) 
%            Typically it is within the range [1e-3, 1e-1], 2e-2 by default.
%   @kappa : Parameter that controls the rate. (See [1])
%            Small kappa results in more iterations and with sharper edges.   
%            We select kappa in (1, 2].    
%            kappa = 2 is suggested for natural images.  
%
%   Example
%   ==========
%   Im  = imread('pflower.jpg');
%   S  = L0Smoothing(Im); % Default Parameters (lambda = 2e-2, kappa = 2)
%   figure, imshow(Im), figure, imshow(S);


if ~exist('kappa','var')
    kappa = 2.0;
end
if ~exist('lambda','var')
    lambda = 2e-2;
end
S = im2double(Im);
betamax = 1e5;
fx = [-1, 1];
fy = [-1; 1];
[N,M,D] = size(Im);
sizeI2D = [N,M];
otfFx = psf2otf(fx,sizeI2D);
otfFy = psf2otf(fy,sizeI2D);
Normin1 = fft2(S);
Denormin2 = abs(otfFx).^2 + abs(otfFy ).^2;
if D>1
    Denormin2 = repmat(Denormin2,[1,1,D]);
end
beta = 2*lambda;
while beta < betamax
    Denormin   = 1 + beta*Denormin2;
    % h-v subproblem
    h = [diff(S,1,2), S(:,1,:) - S(:,end,:)];
    v = [diff(S,1,1); S(1,:,:) - S(end,:,:)];
    if D==1
        t = (h.^2+v.^2)<lambda/beta;
    else
        t = sum((h.^2+v.^2),3)<lambda/beta;
        t = repmat(t,[1,1,D]);
    end
    h(t)=0; v(t)=0;
    % S subproblem
    Normin2 = [h(:,end,:) - h(:, 1,:), -diff(h,1,2)];
    Normin2 = Normin2 + [v(end,:,:) - v(1, :,:); -diff(v,1,1)];
    FS = (Normin1 + beta*fft2(Normin2))./Denormin;
    S = real(ifft2(FS));
    beta = beta*kappa;
    fprintf('.');
end
fprintf('
');
end

代码一的效果图

代码二，本代码对应于向量化显式求解，是实验代码，只能处理单通道，效率不及代码一，仅为示例。

% Author: Kang Kai( Nickname: quarryman)
% Update: 2016-01-13
% References: 
% [1] "Image Smoothing via L0 Gradient Minimization", Li Xu, 
% Cewu Lu, Yi Xu, Jiaya Jia, ACM Transactions on Graphics, 
% (SIGGRAPH Asia 2011), 2011. 
%  
% This code is only for non-commercial use .

function U = kcvL0Smooth(U0, lambda, kappa)
% kcvL0Smooth - Image Smoothing via L0 Gradient Minimization
%   U = kcvL0Smooth(U0, lambda, kappa) performs L0 gradient smoothing of input
%   image U0, with smoothness weight lambda and rate kappa.
%
%   Paras: 
%   @U0    : Input UINT8 image, only accept grayscale images.
%   @lambda: Smoothing parameter controlling the degree of smooth. (See [1]) 
%            Typically it is within the range [1e-3, 1e-1], 2e-2 by default.
%   @kappa : Parameter that controls the rate. (See [1])
%            Small kappa results in more iterations and with sharper edges.   
%            We select kappa in (1, 2].    
%            kappa = 2 is suggested for natural images.  
%
%   Example
%   ==========
%   U0  = imread('pflower.jpg');
%   U  = kcvL0Smooth(U0);
%   figure, imshow(U0), figure, imshow(U);

if ~exist('U0','var')
    U0 = imread('lena.jpg');
    U0 = rgb2gray(U0);
end
if ~exist('lambda','var')
    lambda = 0.005;
end
if ~exist('kappa','var')
    kappa = 2.0;
end

betaMax = 1e5;
beta = 2 * lambda;

U = im2double(U0);

while beta < betaMax
    % v subproblem
    Vx = padarray(diff(U, 1, 2), [0 1], 'post');
    Vy = padarray(diff(U, 1, 1), [1 0], 'post');
    t = (Vx.^2 + Vy.^2) < lambda / beta;
    Vx(t) = 0; Vy(t) = 0;
    % U subproblem
    U = updateU(U, Vx, Vy, beta);
    beta = beta * kappa;
    imshow(U); pause(1)
    fprintf('.');
end

end

function U = updateU(U0, Vx, Vy, beta)
    [m, n] = size(U0); k = m * n;
    
    dVx = padarray(-diff(Vx, 1, 2), [0 1], 'pre');
    dVy = padarray(-diff(Vy, 1, 1), [1 0], 'pre');
    B = U0 + beta * (dVx + dVy); b = B(:);
    
    dx = [ones(m, n - 1), zeros(m, 1)];
    dy = [ones(m - 1, n); zeros(1, n)];
    dx = beta * dx(:); dy = beta * dy(:);
    T1 = spdiags([dx, dy], [-m, -1], k, k);
    
    W = padarray(dx, m, 'pre'); W = W(1 : end - m);
    N = padarray(dy, 1, 'pre'); N = N(1 : end - 1);
    % T2 = 1 + (dx + dy + W + N);
    T2 = 1 + (2 * beta + W + N);
    A = spdiags(T2, 0, k, k) - T1 - T1';
    % deprecated, out of memory
    % Sm = diag(ones(m - 1, 1), 1) - eye(m);
    % Sn = diag(ones(n - 1, 1), 1) - eye(n);
    % A = eye(m * n) + beta * (kron(Sn'*Sn, eye(m)) +  ...
    %     kron(eye(n), Sm'*Sm));
    U = reshape(A  b, m, n);
end

代码二的效果图：

正文