compare_denoise.m

脊波工具和相关实验以及有关实验数据

function [im, wim, riim] = compare_denoise(org_im, rho, th, sorh)
% COMPARE_DENOISE   Comparison of Denosing using Wavelet and Ridgelet
%
%	[im, wim, riim] = compare_denoise(org_im, rho, th, sorh)
%
% Input:
%	org_im:	original square image of size 257 x 257
%	rho:	signal to noise ratio
%	th:	threshold for denoising, T = sigma * th
%
% Output:
%	im:	additive Gaussian white noise image
%	wim:	denoised image using wavelet
%	riim:	denoised image using ridgelet

%%% Favourite run:
% >> im = imread('obj.png'); im = double(im) / 256;
% >> th = thselect(im, 'sqtwolog')                 
% >> [nim, wim, riim] = compare_denoise(im, 3, th, 'h');
              

% Wavelet setting
wname = 'sym4';
edge = 'per';
st = dwtmode('status', 'nodisp');
if ~strcmp(st, edge)
    dwtmode(edge);
end

% Signal to noise ratio
if ~exist('rho', 'var')
    rho = 5;
end

im_size = size(org_im);
if any(im_size ~= [257, 257])
    error('org_im is supposed to be of size 257 x 257');
end

% Take the maximum wavelet decomposition level
nlevels = log2(im_size(1) - 1);

% Generate noisy image. 
init = 2055615866; 
randn('seed', init); 
sig = std(org_im(:));
sigma = sig / rho;
im = org_im + sigma * randn(size(org_im));

% Image of dyadic size for wavelet transform
im1 = im(1:end-1, 1:end-1);
org_im1 = org_im(1:end-1, 1:end-1);

% Find default thresholding values (see ddencmp). 
% [thr, sorh, keepapp] = ddencmp('den', 'wv', im);
thr = th * sigma;
keepapp = 0;

if ~exist('sorh', 'var')
    sorh = 'h';
end


%%%%% WAVELET DENOSING %%%%%
% Global denoising
wim = wdencmp('gbl', im1, wname, nlevels, thr, sorh, keepapp);


%%%%% RIDGELET DENOSING %%%%%
%%% Ridgelet transform
[ri, s, m] = frito(im, wname);

% Thresholding
rit = wthresh(ri, sorh, thr);

% Inverse ridgelet transform
riim = ifrito(rit, s, m, wname);

% Optional step: Wiener denoising
window = [11, 11];      % Window size
wn_riim = wiener2(riim, window);


%%%%% Plot
clf, colormap('gray')
range = [-.5, 1.25];

subplot(2,2,1), imagesc(im, range); axis image
set(gca, 'FontSize', 8);
title(sprintf('Noisy Image (SNR = %.2f dB)', ...
 	      SNR(org_im, im)), 'FontSize', 10);

subplot(2,2,2), imagesc(wim, range); axis image
set(gca, 'FontSize', 8);
title(sprintf('Denoise using 2D DWT (SNR = %.2f dB)', ...
	      SNR(org_im1, wim)), 'FontSize', 10);

subplot(2,2,3), imagesc(riim, range); axis image
set(gca, 'FontSize', 8);
title(sprintf('Denoise using FRIT (SNR = %.2f dB)', ...
	      SNR(org_im, riim)), 'FontSize', 10);

subplot(2,2,4), imagesc(wn_riim, range); axis image
set(gca, 'FontSize', 8);
title(sprintf('Denoise using FRIT + Wiener (SNR = %.2f dB)', ...
	      SNR(org_im, wn_riim)), 'FontSize', 10);