fritm.m

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

function [r, l, m] = fritm(a, wname, maxR)
% FRITM	   Orthonormal finite ridgelet transform mixed with DCT
%	[r, l, m] = fritm(a, wname, [maxR])
%
% Input:
%	a:	image matrix of size P by P, P is a prime number
%	wname:	wavelet name
%	maxR:	[optional] maximum radius of directions that used by wavelet
%		Default is 3.
%
% Output:
%       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
%
% Note:
%	This implement the orthonormal version of finite ridgelet
%	transform.  The result (P-1)x(P+1) coefficients in r together
%	with m makes up total of exactly PxP coefficients.
%
%	However there is a restriction on the size P: P+1 must be dyadic
%	number (the method can easily be extended for P = 2^n + m, though).
%	The typical size of the input image for FRITM is 257 by 257, or
%	17 by 17 image blocks.
%
%      	The wavelet transform takes the maximum posible	number of 
%	decomposition levels.	      
%
% See also:	IFRITM,

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

if ndims(a) ~= 2
    error('Input must be a matrix of 2 dimensions');
end

[p, q] = size(a);
if (p ~= q) | ~isprime(p)
    error('Input must be a P by P matrix, P is a prime number')
end

% Subtract the DC component
m = mean(a(:));
a = a - m;

% Normalize for unit norm
m = p * m;

% Finite Radon transform
ra = frat(a);

% 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);

% 1D wavelet transform on the selected columns of the Radon transform 
% -> "Ridgelets".  Care is taken for the non-dyadic size to ensure
% orthonormal condition

r = zeros(p-1, p+1);

if isdyadic(p - 1)
    % Number of wavelet decomposition levels
    n = log2(p - 1);

    % Make sure using the periodic extension mode
    st = dwtmode('status', 'nodisp');
    if ~strcmp(st, 'per')
	dwtmode('per');
    end    
    
    % Take wavelet transform at each direction except the last coefficients
    [ri, l] = wavedecc(ra(1:end-1, wind), n, wname);
    
    % Incooperate the Radon coefficient to the wavelet approx. coeff
    ri(1, :) = (ri(1, :) - sqrt(p-1) * ra(end, wind)) / sqrt(p);    

else
    error('Have not support this size of image yet!');
end

% Save "true" ridgelet coefficients into a mixed matrix
r = zeros(p-1, p+1);
r(:, wind) = ri;

% The remaining projections are transformed by the DCT
% The DC component is removed since projections have zero mean
if ~isempty(cind)
    rc = dct(ra(:, cind));
    r(:, cind) = rc(2:end, :);
end

%----------------------------------------------------------------------------%
% 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