afritfd.m

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

function [r, m, ind] = afritfd(x, l, wpar, tent, pent)
% AFRITFD  Adaptive Folded Finite Ridgelet Transform
%	[r, m, ind] = afritfd(x, n, wpar, [tent, pent])
%
% Input:
%	x:	image matrix of size N by N where (2N-1) is a prime number
%	l:	number of wavelet decomposition levels
%	wpar:	parameters for spline biorthogonal filters, e.g. [1,3]
%	tent, pent:	[optional] type and parameter for entropy
%
% Output:
%	r:	ridgelet coefficients in a matrix, one column 
%		for each direction (there is N+1 directions)
%	m:	normalized mean value of the folded image
%	ind:	list of projections in which DWT is replaced by DCT
%
% Note:
%	This uses folded FRAT with symmetric extention and biorthogonal
%	symmetric wavelet transform (WaveLab function)
%
%	It adaptively select the best transform for each projections
%
% See also:	FRATFD, FRITFD

% Folded Finite Radon transform
[ra, m] = fratfd(x);

% Columnwise DCT (DCT of the FRAT projections)
rc = dct(ra);

% Columnwise wavelet transform of the folded FRAT

% Retrieve the filters from input parameters
[qmf, dqmf] = MakeBSFilter('CDF', wpar);

%%% UN-OPTIMIZED Part (modified from WaveLab's FWT_SBS)
n = size(ra, 1);
J = ceil(log(n)/log(2));
L = J - l;		% coarset level
dp = dyadpartition(n);

r = zeros(n, n+1);

% Temporary variables
wcoef = zeros(1, n);

% Transpose ra since WaveLab operates on row vector
ra = ra';

for i = 1:n+1
    % FWT_SBS of the i-th "column" of ra
    beta = ra(i, :);  % take samples at finest scale as beta-coeffts

    for j = J-1:-1:L,
	[beta, alfa] = DownDyad_SBS(beta, qmf, dqmf);
	wcoef((dp(j+1)+1):dp(j+2)) = alfa;
    end
    wcoef(1:length(beta)) = beta;
    
    r(:, i) = wcoef';
end


%%% Adaptive transform selection
% Compute entropies of projections for each method
if ~exist('tent', 'var')
    tent = 'norm';
end

if ~exist('pent', 'var')
    pent = 1;
end

ei = col_entropy(r, tent, pent);
ec = col_entropy(rc, tent, pent);

% Replace DWT by DCT in the projections that DCT has smaller entropy
ind = find(ei > ec);
r(:, ind) = rc(:, ind);