function [h]= lmd(im,winSize,patchSize,bins)
% a=imread('D:\face.jpg');
% im=rgb2gray(a);
% winSize=7;
% bins=16;
% patchSize=[4,4];
load(strcat('lmd_winSize',num2str(winSize),'_',num2str(bins),'bins.mat'));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
im = double(im);
im = padarray(im,[floor(winSize/2),floor(winSize/2)],'symmetric','both');
%padarray的目的是扩展图像尺寸，因为要做卷积
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% initialize parameters
r=(winSize-1)/2;
x=-r:r;
%% Form 1-D filters
STFTalpha=1/winSize; % alpha in STFT approaches
sigmaS=(winSize-1)/4; % Sigma for STFT Gaussian window
% Basic STFT filters
w0=(x*0+1);
w1=exp(complex(0,-2*pi*x*STFTalpha));
w2=conj(w1);
% Gaussian window
gs=exp(-0.5*(x./sigmaS).^2)./(sqrt(2*pi).*sigmaS);
% Windowed filters
w0=gs.*w0;
w1=gs.*w1;
w2=gs.*w2;
% Normalize to zero mean
w1=w1-mean(w1);
w2=w2-mean(w2);
%12~33 在设置STFT的窗口
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
convmode='valid';%设置卷积模式，help covn2了解相关问题
%% 通过两次一维卷积实现四个频率的STFT
Qa=conv2(conv2(im,w0.',convmode),w1,convmode);% 频率u1的STFT,其中，wo.'表示非共轭转置
Qb=conv2(conv2(im,w1.',convmode),w0,convmode);% 频率u2的STFT
Qc=conv2(conv2(im,w1.',convmode),w1,convmode);% 频率u3的STFT
Qd=conv2(conv2(im,w1.',convmode),w2,convmode);% 频率u4的STFT
%38~41这四行是通过卷积运算实现了STFT，至于为什么STFT可通过卷积实现，可参看我发给你的文章及其参考文献5，6，11
%% 求不同频率STFT的幅度
Ma=abs(Qa);
Mb=abs(Qb);
Mc=abs(Qc);
Md=abs(Qd);
clear Qa Qb Qc Qd;
%用abs命令求出幅值
%% LMD：Local Magnitude Descriptor
ww = [1,2,4,8,16,32,64,128];
% 频率u1
LMD_u1=zeros(size(Ma,1)-2,size(Ma,2)-2);
for ii=2:size(Ma,1)-1
for jj=2:size(Ma,2)-1
tmp=Ma(ii,jj);
LMD_u1(ii-1,jj-1)=(tmp<=Ma(ii-1,jj-1))+(tmp<=Ma(ii-1,jj))*ww(2)+(tmp<=Ma(ii-1,jj+1))*ww(3)+(tmp<=Ma(ii,jj+1))*ww(4) ...
+(tmp<=Ma(ii+1,jj+1))*ww(5)+(tmp<=Ma(ii+1,jj))*ww(6)+(tmp<=Ma(ii+1,jj-1))*ww(7)+(tmp<=Ma(ii,jj-1))*ww(8);
end
end
% 频率u2
LMD_u2=zeros(size(Mb,1)-2,size(Mb,2)-2);
for ii=2:size(Mb,1)-1
for jj=2:size(Mb,2)-1
tmp=Mb(ii,jj);
LMD_u2(ii-1,jj-1)=(tmp<=Mb(ii-1,jj-1))+(tmp<=Mb(ii-1,jj))*ww(2)+(tmp<=Mb(ii-1,jj+1))*ww(3)+(tmp<=Mb(ii,jj+1))*ww(4) ...
+(tmp<=Mb(ii+1,jj+1))*ww(5)+(tmp<=Mb(ii+1,jj))*ww(6)+(tmp<=Mb(ii+1,jj-1))*ww(7)+(tmp<=Mb(ii,jj-1))*ww(8);
end
end
% 频率u3
LMD_u3=zeros(size(Mc,1)-2,size(Mc,2)-2);
for ii=2:size(Mc,1)-1
for jj=2:size(Mc,2)-1
tmp=Mc(ii,jj);
LMD_u3(ii-1,jj-1)=(tmp<=Mc(ii-1,jj-1))+(tmp<=Mc(ii-1,jj))*ww(2)+(tmp<=Mc(ii-1,jj+1))*ww(3)+(tmp<=Mc(ii,jj+1))*ww(4) ...
+(tmp<=Mc(ii+1,jj+1))*ww(5)+(tmp<=Mc(ii+1,jj))*ww(6)+(tmp<=Mc(ii+1,jj-1))*ww(7)+(tmp<=Mc(ii,jj-1))*ww(8);
end
end
% 频率u4
LMD_u4=zeros(size(Md,1)-2,size(Md,2)-2);
for ii=2:size(Md,1)-1
for jj=2:size(Md,2)-1
tmp=Md(ii,jj);
LMD_u4(ii-1,jj-1)=(tmp<=Md(ii-1,jj-1))+(tmp<=Md(ii-1,jj))*ww(2)+(tmp<=Md(ii-1,jj+1))*ww(3)+(tmp<=Md(ii,jj+1))*ww(4) ...
+(tmp<=Md(ii+1,jj+1))*ww(5)+(tmp<=Md(ii+1,jj))*ww(6)+(tmp<=Md(ii+1,jj-1))*ww(7)+(tmp<=Md(ii,jj-1))*ww(8);
end
end
%%%%以LBP方式从4个幅值图中提取LBP描述子
%% uniform pattern of histogram
clear Ma Mb Mc Md;
LMD=struct;
LMD(1).u=LMD_u1;LMD(2).u=LMD_u2;LMD(3).u=LMD_u3;LMD(4).u=LMD_u4;%
lmdhist1=[];
for ii=1:length(LMD)
tmp_LMD = imPart(LMD(ii).u,patchSize);%
for jj=1:numel(tmp_LMD)
tmp=hist(tmp_LMD{jj}(:),0:255);
histgram=tmp/sum(tmp);
tmp=tmp/sum(tmp);
[histgram]=scale_bins(tmp,hist_idx_mag,bins);
lmdhist1=[lmdhist1,histgram];
end
end
h=lmdhist1';
% 使用直方图信息，各个频率串接使用
%97~111行在压缩直方图
