Converting a specific Matlab Script to Python

Question

The problem seems to be with w0, w1, and w2. Matlab does not have 1D arrays, so if you create what should be a 1D array it is automatically converted to a 2D array. numpy, however, does have 1D arrays, and arange produces 1D array.

In this case, x, and w0, w1, and w2 that are derived from x, are all 1D arrays. convolve2d, as its name implies, requires 2D arrays. img, on the other hand, is likely a 2D array. That is why you are getting the error.

The solution is to convert your 1D x to a 2D value. This will then be carried forward through the other operations involving x, making w0, w1, and w2 also 2D.

Either of these two approaches should work:

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

or

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

The first case always make the array 2D or higher, while the second adds one new dimension no matter what the current dimension is. Since you know the dimensionality already, that is not a problem.

Although you didn't reach it because of the error, your transpose operation also will not work as expected, since the transpose of a 1D array does nothing. This will also fix that problem.

Also, three other things:

First, you don't need the semicolon (;) at the end of the lines. It does nothing when placed at the end of a line in Python.

Second, you are doing transpose the hard way. The easy way is this:

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

Third, it is probably easier to do something like:

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)])

Update

There are a number of other problems that I noticed reading the code:

First, if you are using python 2.x, then by default division of integers returns an integer. So, for example, 1/2==0. To change this so 1/2==.5, then add this above your existing imports:

from __future__ import division

This is not an issue in python 3.x.

Next, in python np.arange does not include the upper end, so:

x=np.arange(-r,r)

Should be:

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

Next, the Matlab version of this line uses complex numbers:

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

The python version doesn't. This means both w1 and w2 will be different between the python and Matlab versions. So:

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

Should be:

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

Next, in the matlab version you transpose the first w0:

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

In the python version you don't, so:

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

Should be:

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

Next, in the python version, i already starts at 0, so subtracting 1 from it changes the values compared to the Matlab version. Also, python uses ** for powers, not ^ (in python ^ is binary XOR). So:

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

Should be:

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

Next, in Matlab, as you know, arr(:) flattens an array. This is not the case in python. Python uses arr.flatten(). And the histogram returns the bin edges, which you don't want. So you need to change:

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

to:

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

Once you do that, to keep things the same as the Matlab version, you would then need to change:

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

print LPQdesc[0]
return LPQdesc[0]

to:

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

print LPQdesc
return LPQdesc

Next, while this is not an error, you can simplify:

w0=(x*0+1);

To:

w0=np.ones_like(x);

Finally, I don't see how this line could possibly work:

LPQdesc=uint8(LPQdesc);

You have probably not tried mode=='im'. The line should be:

LPQdesc=np.uint8(LPQdesc)

Also not an error, but you can use arr.real and arr.imag to get the real and imaginary parts in python, and arr.sum() to get the sum.

Also not an error, but since taking a sub-array in python doesn't take any additional memory, this:

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

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

Can become:

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

Next, you never use u, so you can get rid of this line entirely (if you had u you would need r+1 again):

u=np.arange(1,r)

And you never use decorr either.

Next, you never compute w0, w1, or w2 if freqestim is not 1, so your program will crash if you try to run for any freqstim value other than 1. I can't fix this since I don't know what to do for other values of freqstim.

Finally, in both Matlab and Python, using the loop for the sums is slow. Python lets you broadcast low-dimensional arrays to higher dimensions, so you can do something like:

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

So your final python function should look something like this:

# -*- 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