function [tdst_at_certain_point,tdst_at_certain_freq]=tdst(image)
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%输入参数：image：待变换图像，不能是彩色图像，也就是说必须是一个二维矩阵
%输出参数：tdst_at_certain_point：各点处的S变换结果。是一个元胞矩阵，行列数与图像矩阵相同
% tdst_at_certain_freq：各频率处对整幅图像进行S变换的矩阵，行列数比图像矩阵多1
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
tic;
TRUE=1;
FALSE=0;
factor_x=1; %x方向上分辨率控制因子，通常取1，其值增加可以提高分辨率
factor_y=1; %y方向上分辨率控制因子，通常取1，其值增加可以提高分辨率
suppress_kx0ky0=TRUE; %抑制kx=0，ky=0处的值
[N,M]=size(image);
fft2_image=fft2(image); %对图像进行傅里叶变换，为后面移位做准备
tdst_at_certain_freq=cell(N+1,M+1); %初始化元胞矩阵，每个元素是一个矩阵，代表着在某频率点分析整幅图像的结果
kx0ky0=mean(mean(image)); %kx=0,ky=0处二维S变换
%计算二维S变换矩阵
for kx=0:N
for ky=0:M
tdst_at_certain_freq{kx+1,ky+1}=ones(N,M); %初始化该频率点处的变换矩阵
if (kx==0&&ky==0) %零频处的值为图像特点均值
if suppress_kx0ky0==1 %检查是否选择抑制kx=0，ky=0处的值
tdst_at_certain_freq{kx+1,ky+1}=kx0ky0*zeros(N,M); %抑制该点的值
elseif suppress_kx0ky0==0
tdst_at_certain_freq{kx+1,ky+1}=kx0ky0*ones(N,M); %正确生成该点的值
else disp(sprintf('Error.Please choose whether to suppress the value at kx=0,ky=0'));
end
else
tdst_at_certain_freq{kx+1,ky+1}=ifft2(circshift(fft2_image,[-kx-1,-ky-1]).*gauss(N,M,kx+1,ky+1,factor_x,factor_y));%移位是用circshift完成的，x是竖直向下的方向，y是水平向右的方向。移位后与高斯函数窗相乘，再进行傅里叶反变换
end
end
end
%上面得到的变换结果是在各频率点处的整幅图像的变换值，下面将其转换为在各左边点处整个频率平面上的变换值
tdst_at_certain_point=cell(N,M); %初始化矩阵
for x=1:N
for y=1:M
tdst_at_certain_point{x,y}=extract(tdst_at_certain_freq,x,y); %这里无法避免使用循环，对某个坐标点，从tdst_at_certain_freq中抽取相应元素，得到该点在各频率处的变换值
end
end
time=toc;
disp(sprintf('Estimated time is %f',time));
function gauss_w=gauss(length_x,length_y,freq_x,freq_y,factor_x,factor_y)
%二维高斯函数窗
%输入参数：length_x：x方向上长度，即图像矩阵的行数
% length_y：y方向上长度，即图像矩阵的列数
% freq_x：需要计算高斯函数点处x方向上的频率值
% freq_y：需要计算高斯函数点处y方向上的频率值
% factor_x:x方向上分辨率控制因子
% factor_y:y方向上分辨率控制因子
%输出参数：gauss_w：二维高斯函数，是一个矩阵
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
if freq_x~=0 %频率为0时，高斯函数需要单独处理
vector_x=[0:length_x-1;-length_x:-1]; %定义行,这里使用两行来计算高斯函数，是参照Stockwell给出的一维S变换程序，具体原因不详，但是这样可以获得更好的分辨率
vector_x=factor_x*vector_x.^2/freq_x^2;
vector_x=sum(vector_x); %两行相加
vector_x=vector_x'; %得到列
matrix_x=vector_x*ones(1,length_y); %把列变为矩阵
else matrix_x=zeros(length_x,length_y); %频率为0时，没有该部分
end
if freq_y~=0 %频率为0时，高斯函数需要单独处理
vector_y=[0:length_y-1;-length_y:-1]; %定义行
vector_y=factor_y*vector_y.^2/freq_y^2;
vector_y=sum(vector_y); %两行相加
matrix_y=ones(length_x,1)*vector_y; %把行变成矩阵
else matrix_y=zeros(length_x,length_y); %频率为0时，没有该部分
end
matrix=matrix_x+matrix_y; %得到高斯函数系数雏形
gauss_w=matrix*(-2)*pi^2; %整个二维高斯函数窗的指数
gauss_w=exp(gauss_w); %二维高斯函数窗
function matrix=extract(matrix_to_be_extracted,x,y)
%将tdst_at_certain_freq变成tdst_at_certain_point用到的抽取函数，从tdst_at_certain_freq各矩阵元素的同一位置抽取元素按次序组成tdst_at_certain_point的某一矩阵元素
%输入参数：matrix_to_be_extracted:待抽取的元胞矩阵
% x：需要抽取的位置的坐标x
% y：需要抽取的位置的坐标y
%输出参数：matrix：由输入矩阵抽取相应元素得到的新矩阵
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
[N M]=size(matrix_to_be_extracted); %确定输入矩阵的行列数，即确定在x，y方向上共有多少频率点
matrix=ones(N,M); %初始化矩阵
for j=1:N
for k=1:M
temp=matrix_to_be_extracted{j,k}; %找出输入矩阵某一位置的矩阵
matrix(j,k)=temp(x,y); %将该矩阵的元素按输入的要求取出并赋给matrix
end
end
