将特定的Matlab脚本转换为Python

Question

我需要使用numpy和scipy将以下Matlab脚本转换为Python.

此脚本计算称为LPQ(局部相位定量器)的功能,该功能通常用于面部识别.

记录LPQ特征提取的文件可以在这里找到:http : //www.ee.oulu.fi/mvg/files/pdf/ICISP08.pdf 2

该脚本的Matlab版本如下:

function LPQdesc = lpq(img,winSize,mode)

%% Defaul parameters
% Local window size
if nargin<2 || isempty(winSize)
    winSize=3; % default window size 3
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(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
 % STFT uniform window
    % Basic STFT filters
    w0=(x*0+1);
    w1=exp(complex(0,-2*pi*x*STFTalpha));
    w2=conj(w1);


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

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

试图在python中转换Matlab LPQ函数我生成了以下代码:

# -*- coding: cp1253 -*-
import numpy as np
from scipy.signal import convolve2d

def lpq(img,winSize=3,decorr=0,freqestim=1,mode='nh'):
    rho=0.90;

    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)

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

    img=np.float64(img); # Convert np.image to double
    r=(winSize-1)/2; # Get radius from window size
    x=np.arange(-r,r) # Form spatial coordinates in window
    u=np.arange(1,r) # Form coordinates of positive half of the Frequency domain (Needed for Gaussian derivative)

    if freqestim==1:  #  STFT uniform window
         #  Basic STFT filters
        w0=(x*0+1);
        w1=np.exp(-2*np.pi*x*STFTalpha)
        w2=np.conj(w1);

    ## Run filters to compute the frequency response in the four points. Store np.real and np.imaginary parts separately
    # Run first filter
    filterResp=convolve2d(convolve2d(img,w0,convmode),w1,convmode);
    # Initilize frequency domain matrix for four frequency coordinates (np.real and np.imaginary parts for each frequency).
    shape0, shape1 = filterResp.shape
    freqResp=np.zeros((shape0,shape1,8)); 
    # Store filter outputs
    freqResp[:,:,0]=np.real(filterResp);
    freqResp[:,:,1]=np.imag(filterResp);
    # Repeat the procedure for other frequencies
    filterResp=convolve2d(convolve2d(img,w1.transpose(),convmode),w0,convmode);
    freqResp[:,:,2]=np.real(filterResp);
    freqResp[:,:,3]=np.imag(filterResp);
    filterResp=convolve2d(convolve2d(img,w1.transpose(),convmode),w1,convmode);
    freqResp[:,:,4]=np.real(filterResp);
    freqResp[:,:,5]=np.imag(filterResp);
    filterResp=convolve2d(convolve2d(img,w1.transpose(),convmode),w2,convmode);
    freqResp[:,:,6]=np.real(filterResp);
    freqResp[:,:,7]=np.imag(filterResp);

    # Read the size of frequency matrix
    freqRow,freqCol,freqNum=freqResp.shape;

    ## Perform quantization and compute LPQ codewords
    LPQdesc=np.zeros((freqRow,freqCol)); # Initialize LPQ code word np.image (size depends whether valid or same area is used)
    for i in range(0, freqNum):
        LPQdesc=LPQdesc+((freqResp[:,:,i])>0)*(2^(i-1));

    ## Switch format to uint8 if LPQ code np.image is required as output
    if mode=='im':
        LPQdesc=uint8(LPQdesc);

    ## Histogram if needed
    if mode=='nh' or mode=='h':
        LPQdesc=np.histogram(LPQdesc[:],range(256));

    ## Normalize histogram if needed
    if mode=='nh':
        LPQdesc[0]=LPQdesc[0]/sum(LPQdesc[0]);

    print LPQdesc[0]
    return LPQdesc[0]

但是,在执行相同图像的python脚本后,我收到以下错误:

Traceback (most recent call last):
  File "C:\Users\user\Desktop\source\lpq_parametric.py", line 58, in lpq
    filterResp=convolve2d(convolve2d(img,w0,convmode),w1,convmode);
  File "C:\Python27\lib\site-packages\scipy\signal\signaltools.py", line 548, in convolve2d
    out = sigtools._convolve2d(in1, in2, 1, val, bval, fillvalue)
ValueError: object of too small depth for desired array

由于我是python中的菜鸟和scipy/numpy,我需要帮助将这个python脚本置于功能状态.你能帮我修一下脚本中的bug吗？

Answer 1

问题似乎是w0,w1和w2.Matlab没有1D数组,所以如果你创建了一个1D数组,它会自动转换为2D数组.然而,numpy确实有1D数组,而arange会生成1D数组.

在这种情况下,从x导出的x和w0,w1和w2都是1D阵列.顾名思义,convolve2d需要2D数组.另一方面,img可能是2D阵列.这就是你得到错误的原因.

解决方案是将您的1D x转换为2D值.然后,这将通过涉及x的其他操作继续进行,使得w0,w1和w2也是2D.

这两种方法中的任何一种都应该有效:

x=np.atleast_2d(np.arange(-r,r))

要么

x=np.arange(-r,r)[np.newaxis]

第一种情况总是使数组2D或更高,而第二种情况增加一个新的维度,无论当前维度是什么.既然您已经了解了维度,那就不是问题了.

虽然由于错误而没有达到它,但是转置操作也不会按预期工作,因为1D数组的转置什么都不做.这也将解决这个问题.

另外,还有三件事:

首先,您不需要在行尾添加分号(;).当它放在Python的一行末尾时它什么都不做.

其次,你正在努力转移.简单的方法是:

filterResp=convolve2d(convolve2d(img,w1.T,convmode),w0,convmode)

第三,可能更容易做类似的事情:

filterResp1=convolve2d(convolve2d(img,w1.T,convmode),w0,convmode)
filterResp2=convolve2d(convolve2d(img,w1.T,convmode),w1,convmode)
filterResp3=convolve2d(convolve2d(img,w1.T,convmode),w2,convmode)
freqResp=np.dstack([np.real(filterResp1),
                    np.imag(filterResp1),
                    np.real(filterResp2),
                    np.imag(filterResp2),
                    np.real(filterResp3),
                    np.imag(filterResp3),
                    np.real(filterResp4),
                    np.imag(filterResp4)])

更新

我注意到阅读代码时还有许多其他问题:

首先,如果您使用的是python 2.x,那么默认情况下,整数除法会返回一个整数.所以,例如,1/2 == 0.要更改这个1/2 ==.5,然后将其添加到现有导入之上:

from __future__ import division

这在python 3.x中不是问题.

接下来,在python np.arange中不包括上端,所以:

x=np.arange(-r,r)

应该:

x=np.arange(-r,r+1)[np.newaxis]

接下来,该行的Matlab版本使用复数:

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

python版本没有.这意味着python和Matlab版本之间的w1和w2都不同.所以:

w1=np.exp(-2*np.pi*x*STFTalpha)

应该:

w1=np.exp(-2*np.pi*x*STFTalpha*1j)

接下来,在matlab版本中你转换第一个w0:

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

在python版本中你没有,所以:

filterResp=convolve2d(convolve2d(img,w0,convmode),w1,convmode);

应该:

filterResp=convolve2d(convolve2d(img,w0.T,convmode),w1,convmode)

接下来,在python版本中,我已经从0开始,因此从中减去1会改变与Matlab版本相比的值.另外,python使用**代表幂,而不是^(在python中是^二进制异或).所以:

for i in range(0, freqNum):
    LPQdesc=LPQdesc+((freqResp[:,:,i])>0)*(2^(i-1))

应该:

for i in range(0, freqNum):
    LPQdesc=LPQdesc+((freqResp[:,:,i])>0)*(2**i)

接下来,在Matlab中,如你所知,arr(:)会使数组变平.在python中不是这种情况.Python使用arr.flatten().并且直方图返回您不想要的bin边缘.所以你需要改变:

LPQdesc=np.histogram(LPQdesc[:],range(256));

至:

LPQdesc=np.histogram(LPQdesc.flatten(),range(256))[0]

一旦你这样做,为了保持与Matlab版本相同的东西,你需要改变:

if mode=='nh':
    LPQdesc[0]=LPQdesc[0]/sum(LPQdesc[0]);

print LPQdesc[0]
return LPQdesc[0]

至:

if mode=='nh':
    LPQdesc=LPQdesc/sum(LPQdesc)

print LPQdesc
return LPQdesc

接下来,虽然这不是错误,但您可以简化:

w0=(x*0+1);

至:

w0=np.ones_like(x);

最后,我看不出这条线可能如何工作:

LPQdesc=uint8(LPQdesc);

您可能没有尝试过mode =='im'.该行应该是:

LPQdesc=np.uint8(LPQdesc)

也不是错误,但你可以使用arr.real和arr.imag来获取python中的实部和虚部,并使用arr.sum()来获得总和.

也不是错误,但由于在python中使用子数组不需要任何额外的内存,这个:

# Read the size of frequency matrix
freqRow,freqCol,freqNum=freqResp.shape

## Perform quantization and compute LPQ codewords
LPQdesc=np.zeros((freqRow,freqCol))

可以变成:

LPQdesc=np.zeros_like(freqResp[:,:,0])

接下来,你永远不会使用你,所以你可以完全摆脱这条线(如果你有你,你将再次需要r + 1):

u=np.arange(1,r)

而你也从未使用过装饰.

接下来,如果freqestim不是1,你永远不会计算w0,w1或w2,所以如果你试图运行除1之外的任何freqstim值你的程序将会崩溃.我无法解决这个问题,因为我不知道该怎么做对于freqstim的其他值.

最后,在Matlab和Python中,使用循环来计算总和是很慢的.Python允许您将低维数组广播到更高的维度,因此您可以执行以下操作:

    inds = np.arange(freqResp.shape[2])[np.newaxis,np.newaxis,:]
    LPQdesc=((freqResp>0)*(2**inds)).sum(2)

所以你的最终python函数看起来应该是这样的:

# -*- coding: cp1253 -*-

from __future__ import division

import numpy as np
from scipy.signal import convolve2d

def lpq(img,winSize=3,freqestim=1,mode='nh'):
    rho=0.90

    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)

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

    img=np.float64(img) # Convert np.image to double
    r=(winSize-1)/2 # Get radius from window size
    x=np.arange(-r,r+1)[np.newaxis] # Form spatial coordinates in window

    if freqestim==1:  #  STFT uniform window
        #  Basic STFT filters
        w0=np.ones_like(x)
        w1=np.exp(-2*np.pi*x*STFTalpha*1j)
        w2=np.conj(w1)

    ## Run filters to compute the frequency response in the four points. Store np.real and np.imaginary parts separately
    # Run first filter
    filterResp1=convolve2d(convolve2d(img,w0.T,convmode),w1,convmode)
    filterResp2=convolve2d(convolve2d(img,w1.T,convmode),w0,convmode)
    filterResp3=convolve2d(convolve2d(img,w1.T,convmode),w1,convmode)
    filterResp4=convolve2d(convolve2d(img,w1.T,convmode),w2,convmode)

    # Initilize frequency domain matrix for four frequency coordinates (np.real and np.imaginary parts for each frequency).
    freqResp=np.dstack([filterResp1.real, filterResp1.imag,
                        filterResp2.real, filterResp2.imag,
                        filterResp3.real, filterResp3.imag,
                        filterResp4.real, filterResp4.imag])

    ## Perform quantization and compute LPQ codewords
    inds = np.arange(freqResp.shape[2])[np.newaxis,np.newaxis,:]
    LPQdesc=((freqResp>0)*(2**inds)).sum(2)

    ## Switch format to uint8 if LPQ code np.image is required as output
    if mode=='im':
        LPQdesc=np.uint8(LPQdesc)

    ## Histogram if needed
    if mode=='nh' or mode=='h':
        LPQdesc=np.histogram(LPQdesc.flatten(),range(256))[0]

    ## Normalize histogram if needed
    if mode=='nh':
        LPQdesc=LPQdesc/LPQdesc.sum()

    print LPQdesc
    return LPQdesc