lpq函数matlab,LPQ算法总结

本文为原创，转载请说明出处哈~

LPQ(Local Phase Quantization)算法的原理是假设平滑函数h(x)中心对称，则其傅里叶变换为H(u)，对所有H(u)≥0有∠G(u)= ∠F(u),其中F(u)和G(u)分别为原图像和平滑后图像的傅里叶变换，因而在H(u)≥0条件下，图像对平滑有不变性。

为使H(u)≥0，a取为不超过第一个过零点的频率点，其值为a=1/winSize(winSize为输入参数)。分别用f(x)对u1=(a,0),u2=(0,a),u3=(a,a),u4=(a,-a)四个点作STFT，然后分别把四个点的实部和虚部分开,形成一个向量W=[Re{F(u1,x)},Re{F(u2,x)},Re{F(u3,x)},Re{F(u4,x)},Im{F(u1,x)},Im{F(u2,x)},Im{F(u3,x)},Im{F(u4,x)}]T。

最终得到LPQ变换即为Fx=W*fx。然后对参数进行统计分析，若参数是相关的，则用奇异值分解去相关并且量化。

LPQ算法的matlab实现如下(代码引用自UNIVERSITY OF OULU,DEPARTMENT OF COMPUTER SCIENCE AND ANGINEERING)：

function LPQdesc = lpq(img,winSize,decorr,freqestim,mode)

% Funtion LPQdesc=lpq(img,winSize,decorr,freqestim,mode) computes the Local Phase Quantization (LPQ) descriptor

% for the input image img. Descriptors are calculated using only valid pixels i.e. size(img)-(winSize-1).

%

% Inputs: (All empty or undefined inputs will be set to default values)

% img = N*N uint8 or double, format gray scale image to be analyzed.

% winSize = 1*1 double, size of the local window. winSize must be odd number and greater or equal to 3 (default winSize=3).

% decorr = 1*1 double, indicates whether decorrelation is used or not. Possible values are:

% 0 -> no decorrelation,

% (default) 1 -> decorrelation

% freqestim = 1*1 double, indicates which method is used for local frequency estimation. Possible values are:

% (default) 1 -> STFT with uniform window (corresponds to basic version of LPQ)

% 2 -> STFT with Gaussian window (equals also to Gaussian quadrature filter pair)

% 3 -> Gaussian derivative quadrature filter pair.

% mode = 1*n char, defines the desired output type. Possible choices are:

% (default) 'nh' -> normalized histogram of LPQ codewords (1*256 double vector, for which sum(result)==1)

% 'h' -> un-normalized histogram of LPQ codewords (1*256 double vector)

% 'im' -> LPQ codeword image ([size(img,1)-r,size(img,2)-r] double matrix)

%

% Output:

% LPQdesc = 1*256 double or size(img)-(winSize-1) uint8, LPQ descriptors histogram or LPQ code image (see "mode" above)

%

% Example usage:

% img=imread('cameraman.tif');

% LPQhist = lpq(img,3);

% figure; bar(LPQhist);

%

% Version published in 2010 by Janne Heikkil? Esa Rahtu, and Ville Ojansivu

% Machine Vision Group, University of Oulu, Finland

%% Defaul parameters

% Local window size

if nargin<2 || isempty(winSize)

winSize=3; % default window size 3

end

% Decorrelation

if nargin<3 || isempty(decorr)

decorr=1; % use decorrelation by default

end

rho=0.90; % Use correlation coefficient rho=0.9 as default

% Local frequency estimation (Frequency points used [alpha,0], [0,alpha], [alpha,alpha], and [alpha,-alpha])

if nargin<4 || isempty(freqestim)

freqestim=1; %use Short-Term Fourier Transform (STFT) with uniform window by default

end

STFTalpha=1/winSize; % alpha in STFT approaches (for Gaussian derivative alpha=1)

sigmaS=(winSize-1)/4; % Sigma for STFT Gaussian window (applied if freqestim==2)

sigmaA=8/(winSize-1); % Sigma for Gaussian derivative quadrature filters (applied if freqestim==3)

% Output mode

if nargin<5 || isempty(mode)

mode='nh'; % return normalized histogram as default

end

% Other

convmode='valid'; % Compute descriptor responses only on part that have full neigborhood. Use 'same' if all pixels are included (extrapolates image with zeros).

%% Check inputs

if size(img,3)~=1

error('Only gray scale image can be used as input');

end

if winSize<3 || rem(winSize,2)~=1

error('Window size winSize must be odd number and greater than equal to 3');

end

if sum(decorr==[0 1])==0

error('decorr parameter must be set to 0->no decorrelation or 1->decorrelation. See help for details.');

end

if sum(freqestim==[1 2 3])==0

error('freqestim parameter must be 1, 2, or 3. See help for details.');

end

if sum(strcmp(mode,{'nh','h','im'}))==0

error('mode must be nh, h, or im. See help for details.');

end

%% Initialize

img=double(img); % Convert image to double

r=(winSize-1)/2; % Get radius from window size

x=-r:r; % Form spatial coordinates in window

u=1:r; % Form coordinates of positive half of the Frequency domain (Needed for Gaussian derivative)

%% Form 1-D filters

if freqestim==1 % STFT uniform window

% Basic STFT filters

w0=(x*0+1);

w1=exp(complex(0,-2*pi*x*STFTalpha));

w2=conj(w1);

elseif freqestim==2 % STFT Gaussian window (equals to Gaussian quadrature filter pair)

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

elseif freqestim==3 % Gaussian derivative quadrature filter pair

% Frequency domain definition of filters

G0=exp(-x.^2*(sqrt(2)*sigmaA)^2);

G1=[zeros(1,length(u)),0,u.*exp(-u.^2*sigmaA^2)];

% Normalize to avoid small numerical values (do not change the phase response we use)

G0=G0/max(abs(G0));

G1=G1/max(abs(G1));

% Compute spatial domain correspondences of the filters

w0=real(fftshift(ifft(ifftshift(G0))));

w1=fftshift(ifft(ifftshift(G1)));

w2=conj(w1);

% Normalize to avoid small numerical values (do not change the phase response we use)

w0=w0/max(abs([real(max(w0)),imag(max(w0))]));

w1=w1/max(abs([real(max(w1)),imag(max(w1))]));

w2=w2/max(abs([real(max(w2)),imag(max(w2))]));

end

%% Run filters to compute the frequency response in the four points. Store real and imaginary parts separately

% Run first filter

filterResp=conv2(conv2(img,w0.',convmode),w1,convmode);%卷积

% Initilize frequency domain matrix for four frequency coordinates (real and imaginary parts for each frequency).

freqResp=zeros(size(filterResp,1),size(filterResp,2),8);

% Store filter outputs

freqResp(:,:,1)=real(filterResp);

freqResp(:,:,2)=imag(filterResp);

% Repeat the procedure for other frequencies

filterResp=conv2(conv2(img,w1.',convmode),w0,convmode);

freqResp(:,:,3)=real(filterResp);

freqResp(:,:,4)=imag(filterResp);

filterResp=conv2(conv2(img,w1.',convmode),w1,convmode);

freqResp(:,:,5)=real(filterResp);

freqResp(:,:,6)=imag(filterResp);

filterResp=conv2(conv2(img,w1.',convmode),w2,convmode);

freqResp(:,:,7)=real(filterResp);

freqResp(:,:,8)=imag(filterResp);

% Read the size of frequency matrix

[freqRow,freqCol,freqNum]=size(freqResp);

%% If decorrelation is used, compute covariance matrix and corresponding whitening transform

if decorr == 1

% Compute covariance matrix (covariance between pixel positions x_i and x_j is rho^||x_i-x_j||)

[xp,yp]=meshgrid(1:winSize,1:winSize);%直角坐标系

pp=[xp(:) yp(:)];

dd=dist(pp,pp');%欧式几何距离

C=rho.^dd;

% Form 2-D filters q1, q2, q3, q4 and corresponding 2-D matrix operator M (separating real and imaginary parts)

q1=w0.'*w1;

q2=w1.'*w0;

q3=w1.'*w1;

q4=w1.'*w2;

u1=real(q1); u2=imag(q1);

u3=real(q2); u4=imag(q2);

u5=real(q3); u6=imag(q3);

u7=real(q4); u8=imag(q4);

M=[u1(:)';u2(:)';u3(:)';u4(:)';u5(:)';u6(:)';u7(:)';u8(:)'];

% Compute whitening transformation matrix V

D=M*C*M';

A=diag([1.000007 1.000006 1.000005 1.000004 1.000003 1.000002 1.000001 1]); % Use "random" (almost unit) diagonal matrix to avoid multiple eigenvalues.

[U,S,V]=svd(A*D*A);%奇异值分解

% In order to avoid any sign problems in SVM, force the sign of the largest magnitude element in each singular vector to be positive

[~,ii]=max(abs(V),[],1);

V=V*diag(ones(1,size(V,2))-2*double(V((ii+(0:(length(ii)-1))*size(V,1)))

% Reshape frequency response

freqResp=reshape(freqResp,[freqRow*freqCol,freqNum]);%改变输入数组的维数

% Perform whitening transform

freqResp=(V.'*freqResp.').';

% Undo reshape

freqResp=reshape(freqResp,[freqRow,freqCol,freqNum]);

end

%% Perform quantization and compute LPQ codewords

LPQdesc=zeros(freqRow,freqCol); % Initialize LPQ code word image (size depends whether valid or same area is used)

for i=1:freqNum

LPQdesc=LPQdesc+(double(freqResp(:,:,i))>0)*(2^(i-1));%编码

end

%% Switch format to uint8 if LPQ code image is required as output

if strcmp(mode,'im')

LPQdesc=uint8(LPQdesc);

end

%% Histogram if needed

if strcmp(mode,'nh') || strcmp(mode,'h')

LPQdesc=hist(LPQdesc(:),0:255);

end

%% Normalize histogram if needed

if strcmp(mode,'nh')

LPQdesc=LPQdesc/sum(LPQdesc);

end

在利用LPQ进行图像匹配的时候，若效果不明显，可以尝试使用分块图像，就是先把图像划分成子图像，再对每个子图像计算LPQ直方图，把这些子图像的LPQ直方图拼接起来作为图像的LPQ特征值再进行匹配，效果会明显改善。