
% Adam Wang, 17 Mar 2008
% Subsampled 2D S-Transform

function S = st2d_sub(h, subx, suby, subkx, subky)
%
% S = st2d_sub(h, subx, suby, subkx, subky)
%
% h is an image, sub* are subsampling factors
%
% st2d_sub performs a subsampled 2D Stockwell Transform
%
[N M] = size(h);
if (N*M)^2/(subx*suby*subkx*subky)>64^4
    disp('WARNING: Full 2D Stockwell Transform has (n*m)^2 elements! (>16 million)');
end

S = zeros(N/subx,M/suby,N/subkx,M/subky);
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):floor(N/2);
m = -(ceil(M/2)-1):floor(M/2);
nn = -(ceil(N/2/subkx)-1):floor(N/2/subkx);
mm = -(ceil(M/2/subky)-1):floor(M/2/subky);


tic
G1 = zeros(N,N); G2 = zeros(M,M);
for i = 1:N
    if(n(i)~=0)
        G1(i,:) = exp(-2*pi^2*a.^2/n(i)^2);
    else
        G1(i,1) = 1; %kronecker delta function
    end
end
G1 = G1./(sqrt(sum(G1.^2,2))*ones(1,N));

for j = 1:M
    if(m(j)~=0)
        G2(j,:) = exp(-2*pi^2*b.^2/m(j)^2);
    else
        G2(j,1) = 1;
    end
end
G2 = G2./(sqrt(sum(G2.^2,2))*ones(1,M));

for ii = 1:N/subkx
    for jj = 1:M/subky
        tot = zeros(N,M);
        for i = (ii-1)*subkx+1:ii*subkx
            for j = (jj-1)*subky+1:jj*subky
                G = G1(i,:)'*G2(j,:);
                tot = tot + abs(ifft2(G.*circshift(H,[-n(i) -m(j)]))).^2;
            end
        end
        totSum = conv2(tot,ones(subx,suby),'same');
        S(:,:,floor(N/subkx/2)+1-nn(ii),floor(M/2/subky)+1-mm(jj)) = totSum(ceil(subx/2):subx:end,ceil(suby/2):suby:end);
    end
end

S = sqrt(S);
toc


% 
% 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),[])
% 
% 
% testIm2 = 2*ones(64,1)*cos(2*pi*(1:64)/8.*(1:64)/64) + cos(2*pi*(1:64)/4.*(1:64)/64)'*ones(1,64);
% figure, imshow(testIm2,[])
% S2 = st2d(testIm2);
% s2 = squeeze(sum(sum(S2,1),2));
% figure, imshow(abs(s2),[])
% figure, imshow(abs(s2)-abs(fft2c(testIm2)),[])


