将特定的 Matlab 脚本转换为 Python

Question

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

此脚本计算一种称为 LPQ（局部相位量化器）的特征，该特征通常用于人脸识别。

可在此处找到记录 LPQ 特征提取的论文： http://www.ee.oulu.fi/mvg/files/pdf/ICISP08.pdf2

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 脚本处于功能状态。你能帮我修复脚本中的错误吗？

Answer 1

问题似乎出在 w0、w1 和 w2 上。 Matlab 没有一维数组，因此如果您创建应该是一维数组的内容，它会自动转换为二维数组。然而，numpy 确实有一维数组，并且 arange 产生一维数组。

在这种情况下，从 x 派生的 x 和 w0、w1 和 w2 都是一维数组。 convolve2d，顾名思义，需要二维数组。另一方面，img 可能是一个二维数组。这就是您收到错误消息的原因。

解决方案是将 1D x 转换为 2D 值。然后这将通过涉及 x 的其他操作进行，使 w0、w1 和 w2 也成为 2D。

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

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

或

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

第一种情况总是使数组 2D 或更高，而第二种情况会添加一个新维度，无论当前维度是什么。既然你已经知道维度，那不是问题。

虽然由于错误您没有达到它，但您的转置操作也不会按预期工作，因为一维数组的转置没有任何作用。这也将解决这个问题。

还有其他三件事：

首先，您不需要在行尾使用分号 (;)。在 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 版本没有。这意味着 w1 和 w2 在 python 和 Matlab 版本之间是不同的。所以：

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 中 ^ 是二进制 XOR）。所以：

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)

也不是报错，但是你可以在python中使用arr.real和arr.imag得到实部和虚部，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])

接下来，你从不使用 u，所以你可以完全摆脱这一行（如果你有 u，你将需要 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