Group-Sparse Total-Variation Denoising (GS-TVD) : Example 2

Comparison of TV denoising and group-sparse TV denoising. The test signal is a single row extracted from an image.

Ivan Selesnick, selesi@poly.edu, 2012

Contents

Start

clear
close all

printme = @(filename) print('-dpdf', sprintf('figures/Example2_%s', filename));

set(0,'DefaultAxesColorOrder', [1 1 1]*0);   % make all plot lines black

Load data

a = imread('data/lena.png');    % a : 512 by 512 image
a = double(a);
N = size(a,2);                  % N : 512
n = 0:N-1;

s = a(256,:)';                  % s : noise-free signal

sigma = 10.0;                   % sigma : noise standard deviation
randn('state',0);               % Set state so example is reproducible
noise = sigma * randn(size(s));
y = s + noise;                  % y : noisy signal

% Display noise-free signal
figure(1)
clf
subplot(2,1,1)
plot(s)
xlim([0 N])
ax1 = [0 N -30 300];
axis(ax1)
title('Image scanline')

% Display noisy signal
subplot(2,1,2)
plot(y)
axis(ax1)
title('Signal plus noise');

TV Denoising (TVD)

Peform total variation denoising

lam1 = 7.6;                         % lam1 : regularization parameter
Nit = 30;                           % Nit : number of iterations
[x1, cost] = tvd(y, lam1, Nit);     % TVD algorithm

txt1 = sprintf('Total Variation Denoising,  \\lambda = %.2f', lam1);
rmse_tv = sqrt( mean( (s - x1).^2 ) )
txt_rmse_tv = sprintf('RMSE = %.3f', rmse_tv);

% Display cost function history to verify convergence of TVD algorithm
figure(2)
clf
subplot(2,1,1)
plot(1:Nit, cost, '.-')
title('Cost function history');
xlabel('Iteration')

% Display output of TVD
subplot(2,1,2)
plot(x1)
axis(ax1)
title(txt1);
axnote(txt_rmse_tv)
rmse_tv =

    6.8525

Group-Sparse TV Denoising (GS-TVD)

Peform group-sparse total variation denoising with group size 6

K = 6;                                  % K : group size
lam = 2.6;                              % lam : regularization parameter
[xg, costg] = gstvd(y, K, lam, Nit);    % GS-TVD algorithm

txtg = sprintf('Group-Sparse TV Denoising. Group size K = %d,  \\lambda = %.2f', K, lam);
rmse_tvgs = sqrt( mean( (s - xg).^2 ) )
txt_rmse_tvgs = sprintf('RMSE = %.3f', rmse_tvgs);

% Display cost function history to verify convergence of GS-TVD algorithm
figure(2)
clf
subplot(2,1,1)
plot(costg,'.-')
title('Cost function history');
xlabel('Iteration')

% Display output of GS-TVD
subplot(2,1,2)
plot(xg)
axis(ax1)
title(txtg);
axnote(txt_rmse_tvgs)
rmse_tvgs =

    6.4102

Vary group size and lambda

Calculate RMSE as a function of group size K and regularization parameter lambda

lam_vec = 1.0:0.1:14;           % lam_vec : range of lambda

if exist('Example2_rmse_data.mat', 'file')
    % If RMSE is already computed, then load from previously saved file
    load Example2_rmse_data
else
    % Otherwise, compute RMSE for each (K,lam) pair
    L = length(lam_vec);
    rmse = zeros(10, L);
    for K = 1:10
        for l = 1:L
            [xgv, costg] = gstvd(y, K, lam_vec(l), Nit);
            rmse(K,l) = sqrt( mean( (s - xgv).^2 ) );
        end
    end
    save Example2_rmse_data rmse
end

Display RMSE versus lambda

Display RMSE as a function of lambda and group size K

figure(1)
clf
plot(lam_vec, rmse')
xlabel('\lambda')
title('RMSE : Group size 1 through 10')
xlim([1 10])
ylim([6.2 7.5])

[rmse_min, i] = min(rmse(:));
[K, l] = ind2sub(size(rmse), i);

hold on
l1 = 25;
l2 = 16;
l10 = 6;
p = plot(...
    lam_vec(l1), rmse(1,l1), 'o-', ...
    lam_vec(l2), rmse(2,l2), 's-', ...
    lam_vec(l10), rmse(10,l10), 'v-', ...
    'markerfacecolor','white');
hold off
legend(p, 'K = 1', 'K = 2', 'K = 10', 'location', 'se')

printme('RMSE')

Display signals

figure(10)
clf

subplot(4,1,1)
plot(s)
xlim([0 N])
axis(ax1)
title('Signal');

subplot(4,1,2)
plot(y)
axis(ax1)
title('Signal plus noise');
axnote(sprintf('\\sigma = %.1f', sigma))

subplot(4,1,3)
plot(x1)                    % Display TVD output
axis(ax1)
title(txt1);
axnote(txt_rmse_tv)

subplot(4,1,4)
plot(xg)                    % Display GS-TVD output
axis(ax1)
title(txtg);
axnote(txt_rmse_tvgs)

orient tall
printme('all')

Display detail

xlim1 = [260 350];

figure(1)
clf

FS = get(gca, 'fontsize');
subplot(2,2,1)
plot(x1)
set(gca,'fontsize', FS-2)
xlim(xlim1)
ax2 = axis;
axis(ax2)
title('TV Denoising (detail)')
subplot(2,2,2)
plot(xg)
set(gca,'fontsize', FS-2)
title('Group-sparse TVD (detail)')
axis(ax2)

printme('compare')
set(0,'DefaultAxesColorOrder', 'remove');