%   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('\n');

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