
% Adam Wang, 17 Mar 2008.


function S = st2d(h)
%
% S = st2d(h)
%
% h is an image
%
% st2d performs a 2D Stockwell Transform
%
[N M] = size(h);
if N*M>64^2
    disp('WARNING: Full 2D Stockwell Transform has (n*m)^2 elements! (>16 million)');
end

S = zeros(N,M,N,M);
H = fft2(h); % following IEEE TMI 2005, Distributed Vector Processing ..., Brown Zhu and Mitchell

a = [0:floor(N/2) -(ceil(N/2)-1):-1];
b = [0:floor(M/2) -(ceil(M/2)-1):-1];
n = [-(ceil(N/2)-1):-1 1:floor(N/2)];
m = [-(ceil(M/2)-1):-1 1:floor(M/2)];

for i = 1:N-1
    tic
    for j = 1:M-1
        G = exp(-2*pi^2*a.^2/n(i)^2)'*exp(-2*pi^2*b.^2/m(j)^2);
        S(:,:,floor(N/2)+1-n(i),floor(M/2)+1-m(j)) = ifft2(G.*circshift(H,[-n(i) -m(j)]));
    end
    toc
end

S(:,:,floor(N/2)+1,floor(M/2)+1) = sum(h(:))/(N*M); %DC;

% 
% testIm = cos(2*pi*8*(1:64)/64)'*cos(2*pi*16*(1:64)/64);
% figure, imshow(abs(fft2c(testIm)),[])
% S = st2d(testIm);
% s = squeeze(sum(sum(S,1),2));
% figure, imshow(abs(s),[])

