A hitchhiker guide to Matching Pursuit

Laurent Perrinet

The Matching Pursuit algorithm is popular in signal processing and applies well to digital images.

I have contributed a python implementation and we will show here how we may use that for extracting a sparse set of edges from an image.

this will be exposed in the following book chapter (see also http://invibe.net/LaurentPerrinet/Publications/Perrinet15bicv ), to be published by the end of this year:

@inbook{Perrinet15bicv,
    title = {Sparse models},
    author = {Perrinet, Laurent U.},
    booktitle = {Biologically-inspired Computer Vision},
    chapter = {13},
    editor = {Keil, Matthias and Crist\'{o}bal, Gabriel and Perrinet, Laurent U.},
    publisher = {Wiley, New-York},
    year = {2015}
}

setting things up¶

We will take some "baby steps first to show how it works, and then apply that to an image .

At first, we will define the SparseEdges framework which will create the necessary processing steps, from the raw image, to creating the multi-resolution scheme and the Matching Pursuit algorithm:

In [1]:

from SparseEdges import SparseEdges
mp = SparseEdges('https://raw.githubusercontent.com/bicv/SparseEdges/master/default_param.py')
mp.pe.N = 4
mp.pe.do_mask = False
mp.pe.MP_alpha = 1.
mp.pe.do_whitening = False

This defines the following set of parameters:

In [2]:

print('Default parameters: ', mp.pe)

Default parameters:  {'verbose': 0, 'N_image': None, 'seed': 42, 'N_X': 256, 'N_Y': 256, 'noise': 0.33, 'do_mask': False, 'mask_exponent': 3.0, 'do_whitening': False, 'white_name_database': 'kodakdb', 'white_n_learning': 0, 'white_N': 0.07, 'white_N_0': 0.0, 'white_f_0': 0.4, 'white_alpha': 1.4, 'white_steepness': 4.0, 'white_recompute': False, 'base_levels': 1.618, 'n_theta': 24, 'B_sf': 0.4, 'B_theta': 0.17453277777777776, 'N': 4, 'MP_alpha': 1.0, 'd_width': 45.0, 'd_min': 0.5, 'd_max': 2.0, 'N_r': 6, 'N_Dtheta': 24, 'N_phi': 12, 'N_scale': 5, 'loglevel_max': 7, 'edge_mask': True, 'do_rank': False, 'scale_invariant': True, 'multiscale': True, 'kappa_phase': 0.0, 'weight_by_distance': True, 'svm_n_jobs': 1, 'svm_test_size': 0.2, 'N_svm_grid': 32, 'N_svm_cv': 50, 'C_range_begin': -5, 'C_range_end': 10.0, 'gamma_range_begin': -14, 'gamma_range_end': 3, 'svm_KL_m': 0.34, 'svm_tol': 0.001, 'svm_max_iter': -1, 'svm_log': False, 'svm_norm': False, 'figpath': 'results', 'do_edgedir': False, 'edgefigpath': 'results/edges', 'matpath': 'data_cache', 'edgematpath': 'data_cache/edges', 'datapath': 'database', 'figsize': 14.0, 'figsize_hist': 8, 'figsize_cohist': 8, 'formats': ['png', 'pdf', 'jpg'], 'dpi': 450, 'use_cache': True, 'scale': 0.8, 'scale_circle': 0.08, 'scale_chevrons': 2.5, 'line_width': 1.0, 'line_width_chevrons': 0.75, 'edge_scale_chevrons': 180.0}

The useful imports for a nice notebook:

In [3]:

import matplotlib
matplotlib.rcParams.update({'text.usetex': False})
%matplotlib inline
%config InlineBackend.figure_format='retina'
import matplotlib.pyplot as plt

import numpy as np
np.set_printoptions(precision=4)#, suppress=True)

fig_width = 14

In [4]:

print('Range of spatial frequencies: ', mp.sf_0)

Range of spatial frequencies:  [0.4545 0.2935 0.1895 0.1224 0.079  0.051  0.0329 0.0213 0.0137 0.0089
 0.0057]

In [5]:

print('Range of angles (in degrees): ', mp.theta*180./np.pi)

Range of angles (in degrees):  [-82.5 -75.  -67.5 -60.  -52.5 -45.  -37.5 -30.  -22.5 -15.   -7.5   0.
   7.5  15.   22.5  30.   37.5  45.   52.5  60.   67.5  75.   82.5  90. ]

testing one step of (Matching + Pursuit)¶

Here, we will synthesize an image using 2 log-Gabor filters and check that they are correctly retrieved using a few Matching Pursuit steps. An edge will be represented by a nd.array with respectively [x position, y center position, indice of angle in multi resolution, index of scale in multi resolution scheme] and a coefficient (a complex number whose argument determines the phase of the log-Gabor filter):

In [6]:

# one
edge_in, C_in= [3*mp.pe.N_X/4, mp.pe.N_Y/2, 2, 2], 42
# the second
edge_bis, C_bis = [mp.pe.N_X/8, mp.pe.N_Y/4, 8, 4], 4.*np.sqrt(2)*np.exp(1j*np.pi/4.)

Which results in

In [7]:

# filters in Fourier space
FT_lg_in = mp.loggabor(edge_in[0], edge_in[1], sf_0=mp.sf_0[edge_in[3]],
                         B_sf=mp.pe.B_sf, theta= mp.theta[edge_in[2]], B_theta=mp.pe.B_theta)
FT_lg_bis = mp.loggabor(edge_bis[0], edge_bis[1], sf_0=mp.sf_0[edge_bis[3]],
                         B_sf=mp.pe.B_sf, theta= mp.theta[edge_bis[2]], B_theta=mp.pe.B_theta)
# mixing both and shows one
FT_lg_ = C_in *  FT_lg_in + C_bis * FT_lg_bis
image = mp.invert(FT_lg_)

_ = mp.show_FT(FT_lg_, axis=True)

We may start by computing the linear coefficients in the multi resolution scheme

In [8]:

C = mp.linear_pyramid(image)

These coefficients correspond to a measure of how well the image matches with the filters in the multi-resolution scheme. These will be used in the Matching step. We then detect the best match corresponding to the coefficient with maximum absolute activity:

In [9]:

edge_star = mp.argmax(C)
print('Coordinates of the maximum ', edge_star, ', true value: ', edge_in)
print('Value of the maximum ', C[edge_star], ', true value: ', C_in)

Coordinates of the maximum  (192, 128, 2, 2) , true value:  [192.0, 128.0, 2, 2]
Value of the maximum  (41.99999999999932+4.213693403441275e-13j) , true value:  42

At this particular orientation and scale, the absolute activity looks like:

In [10]:

fig, a1, a2 = mp.show_spectrum(np.absolute(C[:, :, edge_star[2], edge_star[3]]), axis=True)
_ = a2.plot([edge_star[1]], [edge_star[0]], 'r*')

On our image, this looks like:

In [11]:

fig, a1, a2 = mp.show_spectrum(image, axis=True)
_ = a2.plot([edge_star[1]], [edge_star[0]], 'r*')

Let's now extract and show the "winner" of the Matching Step:

In [12]:

FT_star = mp.loggabor(edge_star[0], edge_star[1], sf_0=mp.sf_0[edge_star[3]],
                         B_sf=mp.pe.B_sf, theta= mp.theta[edge_star[2]], B_theta=mp.pe.B_theta)
im_star = mp.invert(FT_star)
_ = mp.show_FT(FT_star, axis=True)

Great, it looks like we have extracted the first edge. Let's first check that the energy of

the sum the extracted component energy,
maximum of its power spectrum (equivalent to the above, but we just check),
half the mean spectrum energy (half? rememeber the trick used in Log-Gabor filters) are all equal to one:

In [13]:

print(np.sum(im_star**2), mp.FTfilter(im_star, FT_star).max(), np.mean(np.abs(FT_star)**2)/2)

1.0000000000000002 1.0000000000000004 1.0000000000000004

Let's overlay the extracted edge on the image:

In [14]:

edge_star_in = np.array([edge_star[0], edge_star[1], mp.theta[edge_in[2]], mp.sf_0[edge_star[3]], np.absolute(C[edge_star]), np.angle(C[edge_star])])
mp.pe.figsize_edges = 9
fig, a = mp.show_edges(edge_star_in[:, np.newaxis], image=image)

Now, we may substract the residual from the image, it is the Pursuit step:

In [15]:

image_res = (image - C[edge_star] * im_star).real 
fig, a1, a2 = mp.show_spectrum(image_res, axis=True)
_ = a2.plot([edge_star[1]], [edge_star[0]], 'r*')

Looks pretty clear now that only the second edge remains. We may now repeat another Matching step on the residual image:

In [16]:

C = mp.linear_pyramid(image_res)
edge_star_bis = mp.argmax(C)
print('Coordinates of the maximum ', edge_star_bis, ', true value: ', edge_bis)
print('Value of the maximum ', C[edge_star_bis], ', true value: ', C_bis)

Coordinates of the maximum  (32, 64, 8, 4) , true value:  [32.0, 64.0, 8, 4]
Value of the maximum  (4+4.000000000000005j) , true value:  (4.000000000000001+4j)

Note: to store the complex value of the coefficient, we use the fact that it is easy to transform a complex number to 2 reals and back:

In [17]:

z = np.sqrt(2)*np.exp(1j*np.pi/4.)
z_, z_p = np.absolute(z), np.angle(z)
print(z, z_*np.exp(1j*z_p))

(1.0000000000000002+1j) (1.0000000000000002+1j)

Note that the linear coefficients corresponding to the first winning filter are canceled:

In [18]:

print('Value of the residual ', C[edge_star], ', initial value: ', C_in)

Value of the residual  (-1.1652065246228562e-14+5.182340906836718e-13j) , initial value:  42

Let's overlay the extracted edge on the image:

In [19]:

fig, a1, a2 = mp.show_spectrum(image_res, axis=True)
_ = a2.plot([edge_star_bis[1]], [edge_star_bis[0]], 'r*')

We have well extracted the two edges:

In [20]:

edge_stars = np.vstack((edge_star_in,
                        np.array([edge_star_bis[0], edge_star_bis[1], mp.theta[edge_star_bis[2]], mp.sf_0[edge_star_bis[3]], np.absolute(C[edge_star_bis]), np.angle(C[edge_star_bis])])))
print(edge_stars)
mp.pe.figsize_edges = 9
fig, a = mp.show_edges(edge_stars.T, image=image)

[[ 1.9200e+02  1.2800e+02 -1.1781e+00  1.8950e-01  4.2000e+01  1.0033e-14]
 [ 3.2000e+01  6.4000e+01 -3.9270e-01  7.9005e-02  5.6569e+00  7.8540e-01]]

testing four steps of Matching Pursuit¶

Let's redo these steps using the run_mp function:

In [21]:

# filters in Fourier space
FT_lg_in = mp.loggabor(edge_in[0], edge_in[1], sf_0=mp.sf_0[edge_in[3]],
                         B_sf=mp.pe.B_sf, theta= mp.theta[edge_in[2]], B_theta=mp.pe.B_theta)
FT_lg_bis = mp.loggabor(edge_bis[0], edge_bis[1], sf_0=mp.sf_0[edge_bis[3]],
                         B_sf=mp.pe.B_sf, theta= mp.theta[edge_bis[2]], B_theta=mp.pe.B_theta)
# mixing both and shows one
FT_lg_ = C_in *  FT_lg_in + C_bis * FT_lg_bis

In [22]:

fig = mp.show_FT(FT_lg_)

In [23]:

image = mp.invert(FT_lg_)
edges, C_res = mp.run_mp(image, verbose=True)

Edge  0 / 4  - Max activity  :  0.000  phase=  98.412  deg,  @  (192, 128, 2, 2)
Edge  1 / 4  - Max activity  :  0.000  phase=  -78.690  deg,  @  (32, 64, 8, 4)
Edge  2 / 4  - Max activity  :  0.000  phase=  67.620  deg,  @  (192, 128, 2, 2)
Edge  3 / 4  - Max activity  :  0.000  phase=  -96.698  deg,  @  (190, 127, 2, 2)

In [24]:

fig, a = mp.show_edges(edges, image=image)

Testing edge detection on a natural image¶

trying out on a flikr image from @Doug88888:

In [25]:

!curl https://farm7.staticflickr.com/6058/6370387703_5e718ea681_q_d.jpg -o ../files/2015-05-22-a-hitchhiker-guide-to-matching-pursuit/6370387703_5e718ea681_q_d.jpg

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100 10602  100 10602    0     0  18814      0 --:--:-- --:--:-- --:--:-- 18831

In [26]:

%ls -l ../files/2015-05-22-a-hitchhiker-guide-to-matching-pursuit/MPtutorial*

-rw-r--r--  1 laurentperrinet  staff   98432 Jun 25 14:17 ../files/2015-05-22-a-hitchhiker-guide-to-matching-pursuit/MPtutorial.npy
-rw-r--r--  1 laurentperrinet  staff   16512 Jun 25 14:17 ../files/2015-05-22-a-hitchhiker-guide-to-matching-pursuit/MPtutorial_MSE.npy
-rw-r--r--  1 laurentperrinet  staff   98432 Jun 18 18:26 ../files/2015-05-22-a-hitchhiker-guide-to-matching-pursuit/MPtutorial_lena.npy
-rw-r--r--  1 laurentperrinet  staff  145757 Jun 26 12:05 ../files/2015-05-22-a-hitchhiker-guide-to-matching-pursuit/MPtutorial_rec.gif
-rw-r--r--  1 laurentperrinet  staff   84303 Jun 26 12:05 ../files/2015-05-22-a-hitchhiker-guide-to-matching-pursuit/MPtutorial_rec.mp4

%rm ../files/2015-05-22-a-hitchhiker-guide-to-matching-pursuit/MPtutorial.npy ../files/2015-05-22-a-hitchhiker-guide-to-matching-pursuit/MPtutorial_MSE.npy%rm ../files/2015-05-22-a-hitchhiker-guide-to-matching-pursuit/MPtutorial_MSE.npy

In [27]:

import numpy as np
from SparseEdges import SparseEdges
mp = SparseEdges('https://raw.githubusercontent.com/bicv/SparseEdges/master/default_param.py')
if not mp.pe.do_whitening: print('\!/ Wrong parameters... \!/')
    
# where should we store the data + figures generated by this notebook
import os
name = 'MPtutorial'
mp.pe.matpath, mp.pe.figpath = '../files/2015-05-22-a-hitchhiker-guide-to-matching-pursuit', '../files/2015-05-22-a-hitchhiker-guide-to-matching-pursuit'
mp.pe.do_mask = True

mp.pe.N = 2048
mp.pe.MP_alpha = 1.
mp.pe.MP_alpha = .8

mp.pe.mask_exponent = 6.

mp.init()

# defining input image (get it @ https://www.flickr.com/photos/doug88888/6370387703)
#image = mp.imread('https://farm7.staticflickr.com/6058/6370387703_5e718ea681_q_d.jpg')
image = mp.imread(os.path.join(mp.pe.matpath, '6370387703_5e718ea681_q_d.jpg'))

white = mp.pipeline(image, do_whitening=True)

import os
matname = os.path.join(mp.pe.matpath, name + '.npy')

try:
    edges = np.load(matname)
except:
    edges, C_res = mp.run_mp(white, verbose=True)
    np.save(matname, edges)    

matname_MSE = os.path.join(mp.pe.matpath, name + '_MSE.npy')
try:
    MSE = np.load(matname_MSE)
except:
    MSE = np.ones(mp.pe.N)
    image_rec = np.zeros_like(image)
    for i_N in range(mp.pe.N):
        MSE[i_N] =  ((white-image_rec)**2).sum()
        image_rec += mp.reconstruct(edges[:, i_N][:, np.newaxis])
        
    np.save(matname_MSE, MSE)

Showing the original image with the edges overlaid:

In [28]:

#edges = np.load(matname)

fig_width_pt = 318.670  # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27               # Convert pt to inches
fig_width = fig_width_pt*inches_per_pt  # width in inches
mp.pe.figsize_edges = .382 * fig_width  # useful for papers
mp.pe.figsize_edges = 9 # useful in notebooks
mp.pe.line_width = 1.
mp.pe.scale = 1.
fig, a = mp.show_edges(edges, image=image, show_phase=True, show_mask=True)
#mp.savefig(fig, name + '_rec')