demo3: group-sparse thresholding, 3D data
This example illustrates 3D sparse signal denoising using non-convex overlapping group sparsity (OGS).
This demo illustrates the non-convex form of OGS, with two different non-convex penalties: atan and MCP.
Ivan Selesnick, July 2014. Revised, May 2016.
Reference: P. Y. Chen and I. W. Selesnick, "Group-Sparse Signal Denoising: Non-Convex Regularization, Convex Optimization," IEEE Transactions on Signal Processing, vol. 62, no. 13, pp. 3464-3478, July 1, 2014. doi: 10.1109/TSP.2014.2329274
Contents
Set up
clear close all % Set random number generator for reproducibility rng('default') rng(1) rmse = @(err) sqrt(mean(abs(err(:)).^2));
Make data
% N : size of 3D data N = [30, 40, 20]; % s : noise-free data s = zeros(N); s(10+(1:4), 10+(1:4), 10+(1:4)) = rand(4,4,4)-0.5; s(15+(1:4), 25+(1:4), 8+(1:4)) = rand(4,4,4)-0.5; % y : noisy data sigma = 0.1; y = s + sigma * randn(N);
Denoising using OGS
K = [3 3 3]; % K : group size % Note: the group size does not need to be the same as % the true groups in the data. Usually, the specified % group size should be smaller. (Also, in many real data, % groups are of different sizes anyway, there is not % single group size to identify.) % OGS (atan penalty (default)) lam = 0.027; x1 = ogs3(y, K, lam); rmse(x1 - s) % OGS (mcp penalty) pen = 'mcp'; lam = 0.027; x2 = ogs3(y, K, lam, pen); rmse(x2 - s)
ans = 0.01 ans = 0.01
Display 2D slices
Show 2D slices through the data as images
Clim = [-1 1]/2; % color limits figure(1) clf subplot(2, 2, 1) imagesc(s(:,:,11), Clim) title('Noise-free data'); colormap(gray(256)) axis image off subplot(2, 2, 2) imagesc(y(:,:,11), Clim) title('Noisy data'); axis image off subplot(2, 2, 3) imagesc(x1(:,:,11), Clim) title('Denoised using OGS (atan penalty (default))'); axis image off subplot(2, 2, 4) imagesc(x2(:,:,11), Clim) title('Denoised using OGS (mcp penalty) '); axis image off orient landscape print -dpdf -fillpage figures/demo3_fig1
Display 1D lines
Display 1D lines through the data as plots. The 'mcp' penalty does a little better than the default ('atan') penalty. Both are non-convex penalties.
figure(2) clf subplot(2, 2, 1) plot(s(:,26,11)) title('Noise-free data'); ylim(Clim) subplot(2, 2, 2) plot(y(:,26,11)) title('Noisy data'); ylim(Clim) subplot(2, 2, 3) plot(x1(:,26,11)) title('Denoised using OGS (atan penalty (default))'); ylim(Clim) subplot(2, 2, 4) plot(x2(:,26,11)) title('Denoised using OGS (mcp penalty) '); ylim(Clim) orient landscape print -dpdf -fillpage figures/demo3_fig2