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');