compare_nla.m

为stanford大学donoho教授所编写的redgit变换源代码。是用c编写的

function [wsnr, fsnr, risnr, rwsnr] = compare_nla(im1)
% Comparison of "non-linear approximation" curves
% image is of size 257 by 257

% Wavelet setting
wname = 'sym4';
dwtmode('per');
n = 8;

% Dyadic length
im = im1(1:end-1, 1:end-1);

% Remove DC
im = im - mean2(im);
im1 = im1 - mean2(im1);

% Take the maximum wavelet decomposition level
% n = wmaxlev(size(im), wname);
% n = log2(im_size(1));

% Fourier (DCT) transform
f = dct2_iv(im);
fsort = sort(abs(f(:))); 
ferr = flipud(cumsum(fsort.^2));

% Wavelet transform
[w, s] = wavedec2(im, n, wname);
wsort = sort(abs(w(:)));
werr = flipud(cumsum(wsort.^2));

% Radon transform
ra = frat(im1);
rasort = sort(abs(ra(:)));
raerr = flipud(cumsum(rasort.^2));

% Ridgelet transform
ri = frito(im1, wname);
% ri = frit(im1, n, wname);
% ri = wavedecc(ra1, n, wname);
% Take out one FRAT coefficient on each projection
% ra1 = ra(1:end-1, :);

% Wavelet and DCT on each projection
rc = dct(ra);

% Remove DC component
rc(1, :) = [];

% Compute entropies (l^1 norm) of each method
ei = col_entropy(ri, 'norm', 1);
ec = col_entropy(rc, 'norm', 1);

% Replace wavelet by DCT if DCT has smaller entropy
ind = find(ei > ec)
ri(:, ind) = rc(:, ind);

risort = sort(abs(ri(:)));
rierr = flipud(cumsum(risort.^2));

% Compute SNR's
wsnr = 10 * log10(werr(1) ./ werr);
fsnr = 10 * log10(ferr(1) ./ ferr);
rasnr = 10 * log10(raerr(1) ./ raerr);
risnr = 10 * log10(rierr(1) ./ rierr);


% Plot
total = prod(size(im));
ncoefs = 1:total;
rate = ncoefs / total * 100;

plot(rate, fsnr(ncoefs), ':', ...
     rate, wsnr(ncoefs), '-.',...
     rate, rasnr(ncoefs), '--', ...
     rate, risnr(ncoefs), '-');

% title('Comparison of "Rate-Distortion"', 'FontSize', 14);
xlabel('Percentage of Retained Coefficients (%)', 'FontSize', 12);
ylabel('SNR (dB)', 'FontSize', 12)
legend('DCT', 'Wavelet', 'Radon', 'Ridgelet', 2);

% NOTE:
% To "zoom in" a particular part of the plot, use the AXIS command
% e.g., axis([0 5 0 20]);