To code image as edges, for instance in the SparseEdges sparse coding scheme, we use a model of edges in images. A good model for these edges are bidimensional Log Gabor filter. This is implemented for instance in the LogGabor library. The library was designed to be precise, but not particularly for efficiency. In order to improve its speed, we demonstrate here the use of a cache to avoid redundant computations.

Let's first initialize the notebook:

In [1]:
from __future__ import division, print_function
import numpy as np
np.set_printoptions(precision=6, suppress=True)
%load_ext autoreload
%autoreload 2

timing of the library without cache

Let's make calls to the library and record the wall clock timing:

In [2]:
from LogGabor import LogGabor
lg = LogGabor('https://raw.githubusercontent.com/bicv/SparseEdges/master/default_param.py')
lg.pe.use_cache = False
lg.pe.verbose = 100
lg.init()
In [3]:
%%timeit
edge = [3*lg.pe.N_X/4, lg.pe.N_Y/2, 2, 2]
FT_lg = lg.loggabor(edge[0], edge[1], sf_0=lg.sf_0[edge[3]], B_sf=lg.pe.B_sf, theta=lg.theta[edge[2]], B_theta=lg.pe.B_theta)
6.46 ms ± 144 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Note that most of the time, we compute the filter at the origin and that whenever it is the case we avoid performing the translation. This makes the call systematically faster:

In [4]:
%%timeit
edge = [0., 0., 2, 2]
FT_lg = lg.loggabor(edge[0], edge[1], sf_0=lg.sf_0[edge[3]], B_sf=lg.pe.B_sf, theta=lg.theta[edge[2]], B_theta=lg.pe.B_theta)
3.54 ms ± 50.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Using a cache

We will use the fact that the many calls to the logGabor library repeat the same operation. We can cache the computed matrices instead of repeating the operation. In particular, we will take advantage of using scales (bands) and orientation separately, the multiplication being rapid in numpy.

In [5]:
lg = LogGabor('https://raw.githubusercontent.com/bicv/SparseEdges/master/default_param.py')
lg.pe.use_cache = True
lg.pe.verbose = 100
lg.init()
print ('Dictionary that will contain the matrices=', lg.cache)
Dictionary that will contain the matrices= {'band': {}, 'orientation': {}}

In the beginning, the cache is empty but every time with compute one matrix, it gets filled up:

In [6]:
edge = [0., 0., 2, 2]
FT_lg = lg.loggabor(edge[0], edge[1], sf_0=lg.sf_0[edge[3]], B_sf=lg.pe.B_sf, theta=lg.theta[edge[2]], B_theta=lg.pe.B_theta)
print ('Dictionary that contains the matrices=', lg.cache)
doing band cache for tag  0.189503126127_0.4
doing orientation cache for tag  -1.1780972451_0.17453277777777776
Dictionary that contains the matrices= {'band': {'0.189503126127_0.4': array([[ 0.006271,  0.006501,  0.006738, ...,  0.006984,  0.006738,
         0.006501],
       [ 0.006501,  0.00674 ,  0.006988, ...,  0.007244,  0.006988,
         0.00674 ],
       [ 0.006738,  0.006988,  0.007246, ...,  0.007513,  0.007246,
         0.006988],
       ..., 
       [ 0.006984,  0.007244,  0.007513, ...,  0.007792,  0.007513,
         0.007244],
       [ 0.006738,  0.006988,  0.007246, ...,  0.007513,  0.007246,
         0.006988],
       [ 0.006501,  0.00674 ,  0.006988, ...,  0.007244,  0.006988,
         0.00674 ]])}, 'orientation': {'-1.1780972451_0.17453277777777776': array([[  2.857228e+05,   2.536583e+05,   2.249385e+05, ...,
          7.830696e-14,   7.439627e-14,   7.074259e-14],
       [  3.217782e+05,   2.857228e+05,   2.534195e+05, ...,
          7.445581e-14,   7.077041e-14,   6.732618e-14],
       [  3.626499e+05,   3.220804e+05,   2.857228e+05, ...,
          7.079868e-14,   6.732618e-14,   6.407993e-14],
       ..., 
       [  1.720206e+13,   1.639461e+13,   1.561163e+13, ...,
          3.499896e-06,   3.101855e-06,   2.752276e-06],
       [  1.638199e+13,   1.560551e+13,   1.485306e+13, ...,
          3.949803e-06,   3.499896e-06,   3.104815e-06],
       [  1.559949e+13,   1.485306e+13,   1.413020e+13, ...,
          4.454120e-06,   3.946026e-06,   3.499896e-06]])}}
In [7]:
%%timeit
edge = [3*lg.pe.N_X/4, lg.pe.N_Y/2, 2, 2]
FT_lg = lg.loggabor(edge[0], edge[1], sf_0=lg.sf_0[edge[3]], B_sf=lg.pe.B_sf, theta=lg.theta[edge[2]], B_theta=lg.pe.B_theta)
3.54 ms ± 23.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [8]:
%%timeit
edge = [0., 0., 2, 2]
FT_lg = lg.loggabor(edge[0], edge[1], sf_0=lg.sf_0[edge[3]], B_sf=lg.pe.B_sf, theta=lg.theta[edge[2]], B_theta=lg.pe.B_theta)
395 µs ± 11 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

That's a great improvement! Let's now apply that to the Matching Pursuit algorithm implemented in the SparseEdges library:

application to SparseEdges

In [9]:
from SparseEdges import SparseEdges
mp = SparseEdges('https://raw.githubusercontent.com/bicv/SparseEdges/master/default_param.py')
mp.pe.N = 32 # number of edges
mp.pe.use_cache = False
mp.init()

# defining a test image
image = np.zeros((mp.pe.N_X, mp.pe.N_Y))
image[mp.pe.N_X//2:mp.pe.N_X//2+mp.pe.N_X//4, mp.pe.N_X//2:mp.pe.N_X//2+mp.pe.N_X//4] = 1
image[mp.pe.N_X//2:mp.pe.N_X//2+mp.pe.N_X//4, mp.pe.N_X//4:mp.pe.N_X//2] = -1
In [10]:
%%timeit -n1 -r1
edges, C_res = mp.run_mp(image, verbose=False)
1min 32s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
In [11]:
from SparseEdges import SparseEdges
mp = SparseEdges('https://raw.githubusercontent.com/bicv/SparseEdges/master/default_param.py')
mp.pe.N = 32 # number of edges
mp.pe.use_cache = True
mp.init()

# defining a test image
image = np.zeros((mp.pe.N_X, mp.pe.N_Y))
image[mp.pe.N_X//2:mp.pe.N_X//2+mp.pe.N_X//4, mp.pe.N_X//2:mp.pe.N_X//2+mp.pe.N_X//4] = 1
image[mp.pe.N_X//2:mp.pe.N_X//2+mp.pe.N_X//4, mp.pe.N_X//4:mp.pe.N_X//2] = -1
In [12]:
%%timeit -n1 -r1
edges, C_res = mp.run_mp(image, verbose=False)
1min 8s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)

Which shows a performance gain of approximately 25%. These changes are now effective in the code (see this commit).

some book keeping for the notebook

In [13]:
%load_ext watermark
%watermark
2017-10-14T23:43:12+02:00

CPython 3.6.3
IPython 6.2.1

compiler   : GCC 4.2.1 Compatible Apple LLVM 9.0.0 (clang-900.0.37)
system     : Darwin
release    : 17.0.0
machine    : x86_64
processor  : i386
CPU cores  : 8
interpreter: 64bit
In [14]:
%load_ext version_information
%version_information numpy, scipy, matplotlib, sympy, pillow, imageio
Out[14]:
Software Version
Python 3.6.3 64bit [GCC 4.2.1 Compatible Apple LLVM 9.0.0 (clang-900.0.37)]
IPython 6.2.1
OS Darwin 17.0.0 x86_64 i386 64bit
numpy 1.13.3
scipy 0.19.1
matplotlib 2.1.0
sympy 1.1.1
pillow 4.3.0
imageio 2.1.2
Sat Oct 14 23:43:13 2017 CEST