frito.m

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

function [r, l, m] = frito(a, wname)
% FRITO	   Orthonormal finite ridgelet transform
%	[r, l, m] = frito(a, wname)
%
% Input:
%	a:	image matrix of size P by P, P is a prime number
%	wname:	wavelet name
%
% 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 FRITO is 257 by 257.
%
%      	The wavelet transform takes the maximum posible	number of 
%	decomposition levels.	      
%
% See also:	IFRITO, FRIT

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

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

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
    [r, l] = wavedecc(ra(1:end-1, :), n, wname);
    
    % Incooperate the Radon coefficient to the waveklet approx. coeff
    r(1, :) = (r(1, :) - sqrt(p-1) * ra(end, :)) / sqrt(p);    

else
    error('Have not support this size of image yet!');
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