In this notebook, we'll explore the spectral representation of signals and how we use filters to alter it. We'll use the PyGSP, a Python package that eases Signal Processing on Graphs.
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pygsp as pg
#plt.rcParams['figure.figsize'] = (17, 5)
The goal of this section is to get familiar with two representations of a signal:
pg.graphs.Bunny()
.my_graph.plot()
.graph = pg.graphs.Bunny()
fig, ax = graph.plot(title='clean bunny')
ax.axis('off');
The position of each point along the x axis can be viewed as a signal on the bunny graph. As such, we can compute its Fourier transform and look at it in the spectral domain.
graph.coords
.
The Fourier transform is computed with my_graph.gft(my_signal)
.my_graph.e
).graph.compute_fourier_basis()
plt.plot(graph.e, abs(graph.gft(graph.coords[:, 0])))
[<matplotlib.lines.Line2D at 0x7fbb106d4eb8>]
Add some noise, e.g. with np.random.normal()
(a reasonable scale is 4e-3
), to the 3D position of the vertices (stored in graph.coords
).
noise = np.random.normal(0, 5e-3, size=(graph.N, 3))
coords_clean = graph.coords
coords_noisy = graph.coords + noise
graph.coords = coords_noisy
fig, ax = graph.plot(title='noisy bunny')
ax.axis('off');
Look again at the frequency content of the position signal. What changed?
plt.semilogy(graph.e, abs(graph.gft(noise[:, 0])), label='noise');
plt.semilogy(graph.e, abs(graph.gft(graph.coords[:, 0])), label='noisy bunny');
plt.legend();
To denoise, we have to remove the noise while keeping the information, which is the shape of the bunny. Visually, you should see that it is much easier to do that in the spectral domain than in the spatial domain!
Design a filter that will remove the high frequencies, where the noise is concentrated, while keeping the low frequencies, where the information is concentrated.
Try the following filters: pg.filters.Heat
, pg.filters.Rectangular
, and pg.filters.Expwin
. Do play with their parameters!
You can plot a filter's response with my_filter.plot()
.
g = pg.filters.Heat(graph, scale=1)
#g = pg.filters.Rectangular(graph, band_max=0.2)
#g = pg.filters.Expwin(graph, band_max=0.6)
g.plot(eigenvalues=False);
my_filter.filter(my_signal)
.
If using a rectangular filter, additionally set the option method='exact'.graph.coords = g.filter(coords_noisy, method='exact')
plt.semilogy(graph.e, abs(graph.gft(coords_clean[:, 0])), label='original bunny');
plt.semilogy(graph.e, abs(graph.gft(coords_noisy[:, 0])), label='noisy bunny');
plt.semilogy(graph.e, abs(graph.gft(graph.coords[:, 0])), label='denoised bunny');
plt.legend();
Finally, look at the denoised bunny in the spatial domain.
fig, ax = graph.plot(title='cleaned bunny')
ax.axis('off');
The goal of this section is to see how we can use filters to solve the heat equation.
Define a domain on which we'll diffuse heat.
You can try the following graphs: Sensor
, Airfoil
, Comet
.
graph = pg.graphs.Airfoil()
graph.estimate_lmax()
fig, ax = graph.plot()
Define an initial heat distribution.
A good one to try first to get some intuition is a delta: i.e., a signal that has value 1 on a vertex (or multiple vertices) and 0 everywhere else.
Plot that signal on the graph with my_graph.plot(my_signal)
.
DELTA = range(0, 800)
initial = np.zeros(graph.n_vertices)
initial[DELTA] = 1
fig, ax = graph.plot(initial)
Instantiate a heat kernel with pg.filters.Heat
. Plot its frequency response.
heat = pg.filters.Heat(graph, scale=[10, 50, 200])
fig, ax = heat.plot(sum=False)
ax.set_xlim(0, 1)
(0, 1)
Filter your initial distribution with the heat kernel and observe how heat diffuses.
Note that the scale
parameter of the kernel controls the amount of diffusion.
What is the solution after infinite time?
diffused = heat.filter(initial)
fig, axes = plt.subplots(1, len(heat), figsize=(17, 5))
for i, ax in enumerate(axes):
fig, ax = graph.plot(diffused[:, i], ax=ax)
graph = pg.graphs.Bunny()
graph.estimate_lmax()
Create a filterbank of 6 mexican hat wavelets. Visualize the filterbank in the spectral domain.
g = pg.filters.MexicanHat(graph, Nf=6)
fig, ax = g.plot()
Visualize the filterbank in the vertex domain.
Look at some of the localized wavelets.
You can localize a filter with my_filter.localize()
.
DELTA = 20
fig = plt.figure(figsize=(12, 3))
for i in range(3):
ax = fig.add_subplot(1, 3, i+1, projection='3d')
_ = graph.plot_signal(g[i].localize(DELTA), ax=ax)
_ = ax.set_title('Wavelet {}'.format(i+1))
ax.set_axis_off()
Let's now try to estimate the curvature of the underlying 3D model by only using spectral filtering on the nearest-neighbor graph formed by its point cloud. A simple, but not theoretically grounded, way to accomplish that is to use the coordinates map $[x, y, z]$ and filter it using the above defined wavelets. Doing so gives us a set of $N_f$ 3-dimensional signals $[g_i(L)x, g_i(L)y, g_i(L)z], \ i \in [0, \ldots, N_f]$ that describe variation along the 3 coordinates.
signal = graph.coords
signal = g.filter(signal)
The curvature is then estimated by taking the $\ell_1$ or $\ell_2$ norm across the 3D position.
signal = np.linalg.norm(signal, ord=2, axis=1)
Plot the result to observe that we indeed have a measure of the curvature at different scales.
fig = plt.figure(figsize=(10, 7))
for i in range(4):
ax = fig.add_subplot(2, 2, i+1, projection='3d')
title = 'Curvature estimation (scale {})'.format(i+1)
_ = graph.plot_signal(signal[:, i], ax=ax, title=title)
ax.set_axis_off()