ifritm.m

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

function a = ifritm(r, l, m, wname, maxR)
% IFRITM	Inverse orthogonal finite ridgelet transform mixed with DCT
%	a = ifritm(r, l, m, wname, [maxR])
%
% Input:
%       r:      ridgelet coefficients in a (P-1) by (P+1) matrix,
%		one column for each direction
%	l:	structure of the wavelet decomposition that is
%		needed for reconstruction
%	m:	normalized mean value of the image
%	wname:	wavelet name
%	maxR:	[optional] maximum radius of directions that used by wavelet
%		Default is 2.
%
% Output:
%	a:	reconstructed image
%
% See also:	FRITM

if ~exist('maxR', 'var')
    maxR = 3;
end

p = size(r, 2) - 1;
if (size(r, 1) ~= (p - 1)) | ~isprime(p)
    error('Ridgelet coefficients must be in a (P-1) by (P+1) matrix.');
end

% Directions that has normal vector up to maxR are used by DWT,
% the rest are used by DCT instead.
M = bestdir(p);
L = max(abs(M(1,:)), abs(M(2,:)));
wind = find(L <= maxR);
cind = setdiff(1:(p+1), wind);

% Back to Radon domain by inverting the wavelet transform.
% By definition, Radon coefficients should have zero mean
% at each direction

ra = zeros(p, p + 1);

if isdyadic(p - 1)
    % Take out the "true" ridgelet colums
    ri = r(:, wind);
    
    % Number of wavelet decomposition levels
    n = log2(p - 1);

    % Recorrect the wavelet approx. cofficients
    ri(1, :) = sqrt(p) * ri(1, :) / p;

    % Make sure using the periodic extension mode
    st = dwtmode('status', 'nodisp');
    if ~strcmp(st, 'per')
	dwtmode('per');
    end    
        
    % Inverse the wavelet transform
    ra(1:end-1, wind) = waverecc(ri, l, wname);
    
    % Compute the last Radon coefficients
    ra(end, wind) = -sqrt(p-1) * ri(1, :);
    
else
    error('Have not support this size of image yet!');
end

% The remaining projections were transformed by the DCT
if ~isempty(cind)
    rc = zeros(p, length(cind));
    rc(2:end, :) = r(:, cind);
    ra(:, cind) = idct(rc);
end

% Inverse the finite Radon transform with the mean corrected
a = ifrat(ra, 0);

% Add back the DC component
a = a + m / (size(r, 2) - 1);


%----------------------------------------------------------------------------%
% Internal Function(s)
%----------------------------------------------------------------------------%
function id = isdyadic(n)
% True for dyadic (power of two) number

l = log2(n);

if (l == round(l)) & (l > 0)
    id = 1;
else
    id = 0;
end