CODEWITHSUNDEEP
pythonperformancematlabpython-2.7numpy
Victor Pira
Victor Pira

Reputation: 1172

Matlab corr2 equivalent in python - performance issue

I ran into the same problem as presented in this question, i.e. I wanted scalar output of 2D correlation as matlab corr2 but in python. I've solved that by writing my own algorithm, but it has performance issues.

It's based on the formula used in matlab docs:

import numpy as np
from scipy.misc import imread

#Reading the input
a = imread(src[0],flatten=True)
b = imread(src[1],flatten=True)

#Getting shapes and prealocating the auxillairy variables
k = np.shape(a)
H=k[0]
W=k[1]
c = np.zeros((H,W))
d = np.zeros((H,W))
e = np.zeros((H,W))

#Calculating mean values
AM=np.mean(a)
BM=np.mean(b)  

#Calculating terms of the formula
for ii in range(0,H):
  for jj in range(0,W):
    c[ii,jj]=(a[ii,jj]-AM)*(b[ii,jj]-BM)
    d[ii,jj]=(a[ii,jj]-AM)**2
    e[ii,jj]=(b[ii,jj]-BM)**2

#Formula itself
r = np.sum(c)/float(np.sqrt(np.sum(d)*np.sum(e)))

It works, but it's terribly slow. Any ideas how to make it more elegant and faster?

Thanks in advance!

Upvotes: 0

Views: 783

Answers (1)

Divakar
Divakar

Reputation: 221634

This vectorized approach should speed it up for you -

# Vectorized versions of c,d,e
c_vect = (a-AM)*(b-BM)
d_vect = (a-AM)**2
e_vect = (b-BM)**2

# Finally get r using those vectorized versions
r_out = np.sum(c_vect)/float(np.sqrt(np.sum(d_vect)*np.sum(e_vect)))

Runtime tests and verify results

In [70]: M = 100
    ...: N = 100
    ...: a = np.random.rand(M,N)
    ...: b = np.random.rand(M,N)
    ...: 
    ...: k = np.shape(a)
    ...: H=k[0]
    ...: W=k[1]
    ...: AM=np.mean(a)
    ...: BM=np.mean(b)
    ...: 

In [71]: %%timeit
    ...: c = np.zeros((H,W))
    ...: d = np.zeros((H,W))
    ...: e = np.zeros((H,W))
    ...: 
    ...: for ii in range(0,H):
    ...:   for jj in range(0,W):
    ...:     c[ii,jj]=(a[ii,jj]-AM)*(b[ii,jj]-BM)
    ...:     d[ii,jj]=(a[ii,jj]-AM)**2
    ...:     e[ii,jj]=(b[ii,jj]-BM)**2
    ...: 
    ...: r = np.sum(c)/float(np.sqrt(np.sum(d)*np.sum(e)))
    ...: 
10 loops, best of 3: 24.7 ms per loop

In [72]: %%timeit
    ...: c_vect = (a-AM)*(b-BM)
    ...: d_vect = (a-AM)**2
    ...: e_vect = (b-BM)**2
    ...: r_out = np.sum(c_vect)/float(np.sqrt(np.sum(d_vect)*np.sum(e_vect)))
    ...: 
10000 loops, best of 3: 64.7 µs per loop

In [73]: np.allclose(r,r_out)
Out[73]: True

Upvotes: 1

Related Questions