Guglielmo Camporese (gool-yell-moe)

Research Engineer at The Walt Disney Studios

Ph.D. in Brain, Mind and Computer Science

guglielmocamporese [at] gmail [dot] com

Fourier Representation Meets Images

$$N=1$$
$$N=2$$
$$N=5$$
$$N=10$$
$$N=20$$
$$N=50$$
$$N=100$$
$$N=200$$
$$N=500$$

Open In Colab

In general a signal can be approximated by a Fourier Series. In this blog I will represent the signal coming from the edges of the image containing my face using the Fourier representation. Considering $t \in \mathcal{T} = [t_0, t_1] \subset \mathbb{R}$ be an input variable, $T = t_1 - t_0$ the length of the input domain, and $f(\cdot) : \mathcal{T} \rightarrow \mathbb{C}$ being the function that represents the edges signal in the 2D complex plane, the Fourier representation of the signal can be computed as:

$$f_N(t) = \sum_{n=-N}^N c_n e^{i \frac{2 \pi n}{T} t} = \sum_{n=-N}^N c_n e^{i \theta_n t}\ \ \ \ \ (1)$$
where the coefficients are:
$$c_n = \frac{1}{T}\int_\mathcal{T} f(t) e^{-i \theta_n t} dt \ \ \ \ \ (2).$$

Here $N$ is the number of components used for the approximating the input signal, the larger this number, the better the approximation. $T$ is the length of the input domain of $x$, $T = | \mathcal{T}| $.

Let's make an example by considering an image:


# Imports
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
from sklearn.cluster import KMeans
import cv2
import requests

# Load the image
url = 'https://guglielmocamporese.github.io/static/static/img/fourier/image.png'
image = np.array(Image.open(requests.get(url, stream=True).raw))
    

Now let's compute the edges of the image by low-pass filtering the image and using the Canny algorithm, already implemented in OpenCV.


# Compute the edges
def get_edges(x):
  """Extract edges from an image."""
  x = cv2.GaussianBlur(x, (5, 5), 0)
  x = cv2.Canny(x, 10, 200)
  return x

edges = get_edges(image) 
    

The next step is to reduce the number of points contained in the edges. To this end, I clustered these points reducing them from a number $\sim 13000$ to $n\_cluster=2000$ as follows:


# Cluster the points of the edges
edges_pts = np.stack(np.where(edges), 1)
edges_pts[:, 0] = image.shape[0] - edges_pts[:, 0]
edges_pts = np.stack([edges_pts[:, 1], edges_pts[:, 0]], 1)

kmeans = KMeans(n_clusters=2000, random_state=0).fit(edges_pts)
edges_centers = kmeans.cluster_centers_
    
The image above shows the cluster centers that fit the edges.

For now we have a set of points that builds up the edge signal, but this set is un-ordered since each point does not have the meaning of "neighbors". In the image below you can see this fact by connecting adjacent points in the set.

In order to solve the order issue, I designed a sorting algorithm that returns a set of ordered points using the euclidian distance.


# Sort the points to make a smooth path
dist = lambda p0, p1: (((p0 - p1) ** 2).sum() ** 0.5)

def sort_points(points):
  """Sort points by their proximity under the euclidian distance."""
  ordered_pts = [points[0]]
  points = [p for p in points[1:]]
  while len(points) > 0:
    p0 = ordered_pts[-1]
    dist_01 = [dist(p0, p1) for p1 in points]
    i_min = np.argmin(dist_01)
    closest_pt = points[i_min]
    del points[i_min]
    ordered_pts += [closest_pt]
  return np.array(ordered_pts)

edges_centers_sorted = sort_points(edges_centers)
    
The image above shows the same points of the previous image, but with sorted points.

The equation in (1) returns complex values whereas the edge signal that I'm intered in is a real 2D function. In this case the Fourier Series representation can be rewritten as:

$$ \begin{split} f_N(t) &= \sum_{n=-N}^N c_n e^{i \theta_n t} \\ &= \sum_{n=-N}^N \Big[ Re\{c_n\} cos(\theta_n t) - Im\{c_n\} sin(\theta_n t) \Big] + \\ &\ \ \ \ + i \sum_{n=-N}^N \Big[ - Re\{ c_n \} sin(\theta_n t) + Im\{c_n\} cos(\theta_n t) \Big] \end{split}$$
where the coefficients are:
$$ Re\{c_n\} = \frac{1}{T} \int_\mathcal{T} \Big[ Re\{f(t)\} cos(\theta_n t) + Im\{f(t)\} sin(\theta_n t)\Big]dt$$
$$ Im\{c_n\} = \frac{1}{T} \int_\mathcal{T} \Big[ -Re\{f(t)\} sin(\theta_n t) + Im\{f(t)\} cos(\theta_n t)\Big]dt$$
In the following code snippet I approximated the integrals above in order to evaluate them:

def f_transform(f, N=100):
    """Compute the fourier series of the input signal."""
    y = np.zeros_like(f) # output signal of shape [num_points, 2]
    num_points = f.shape[0]
    t = np.linspace(0, 1, num_points)
    dt = np.abs(t[1] - t[0])
    f_real, f_imag = f[:, 0], f[:, 1]
    for n in range(- N, N + 1):
        theta_n = 2 * np.pi * n
        cos_n = np.cos(theta_n * t)
        sin_n = np.sin(theta_n * t)
        c_real_n = (f_real * cos_n + f_imag * sin_n).sum() * dt
        c_imag_n = (- f_real * sin_n + f_imag * cos_n).sum() * dt
        y_real = c_real_n * cos_n - c_imag_n * sin_n
        y_imag = c_real_n * sin_n + c_imag_n * cos_n
        y += np.stack([y_real, y_imag], 1)
    return y

y = f_transform(edges_centers_sorted, N=550)