Numerical Tours
← All Python tours
Download notebook Open in Colab

Introduction to Image Processing

Explore how arrays become images, how convolution removes fine detail, and how derivatives reveal boundaries. Compare spatial operations with their Fourier-domain counterparts to connect the formulas with visible changes in an image.

Run this tour

Run the cells in order with a Python 3 kernel. The first cell locates the companion data and toolbox and installs missing dependencies when needed. All worked examples include their implementation directly in this notebook. Random seeds make comparisons reproducible; you can change them to explore other samples.

# Locate the companion toolbox locally, or fetch it for a standalone/Colab copy.
from pathlib import Path
import importlib.util
import os
import subprocess
import sys

working = Path.cwd()
candidates = [working, working / "python", working.parent / "python"]
python_dir = next((p for p in candidates if (p / "nt_toolbox").is_dir()), None)
if python_dir is None:
    checkout = working / "numerical-tours-support"
    if not checkout.exists():
        subprocess.run(
            [
                "git",
                "clone",
                "--depth",
                "1",
                "--branch",
                "master",
                "https://github.com/gpeyre/numerical-tours.git",
                str(checkout),
            ],
            check=True,
        )
    python_dir = checkout / "python"
os.chdir(python_dir)
if str(python_dir) not in sys.path:
    sys.path.insert(0, str(python_dir))
requirements = python_dir / "requirements.txt"
if any(
    importlib.util.find_spec(name) is None
    for name in [
        "numpy",
        "scipy",
        "matplotlib",
        "skimage",
        "sklearn",
        "pywt",
        "ipywidgets",
        "cvxpy",
        "skfmm",
        "autograd",
        "progressbar",
        "celer",
    ]
):
    subprocess.run(
        [sys.executable, "-m", "pip", "install", "-r", str(requirements)], check=True
    )

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0)
plt.rcParams.update(
    {
        "figure.figsize": (8, 4),
        "figure.dpi": 100,
        "axes.spines.top": False,
        "axes.spines.right": False,
        "font.size": 11,
        "image.cmap": "gray",
    }
)
%matplotlib inline

This numerical tour explores some basic image processing tasks.

$\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

from nt_toolbox.general import np, transform
from nt_toolbox.signal import grad, imageplot, load_image, np, transform
import numpy as np
import matplotlib.pyplot as plt
from numpy import (
    abs,
    arange,
    array,
    concatenate,
    cos,
    exp,
    linspace,
    log,
    matrix,
    number,
    pi,
    real,
    sum,
    transpose,
)
from matplotlib.pyplot import clf, matplotlib, np
from numpy.fft import fft2, fftshift, ifft2

%matplotlib inline

Image Loading and Displaying

Several functions are implemented to load and display images.

First we load an image.

path to the images

name = "nt_toolbox/data/hibiscus.bmp"
n = 256
M = load_image(name, n)

We can display it. It is possible to zoom on it, extract pixels, etc.

m = int(n / 2)
imageplot(M[m - 25 : m + 25, m - 25 : m + 25], "Zoom", [1, 2, 2])
No description has been provided for this image

Image Modification

An image is a 2D array, that can be modified as a matrix.

imageplot(-M, "-M", [1, 2, 1])
imageplot(M[::-1, :], "Flipped", [1, 2, 2])
No description has been provided for this image

Blurring is achieved by computing a convolution with a kernel.

Compute the low pass Gaussian kernel. Warning, the indexes needs to be modulo $n$ in order to use FFTs.

sigma = 5
t = concatenate((arange(0, n / 2 + 1), arange(-n / 2, -1)))
[Y, X] = np.meshgrid(t, t)
h = exp(-(X**2 + Y**2) / (2.0 * float(sigma) ** 2))
h = h / sum(h)
imageplot(fftshift(h))
No description has been provided for this image

Compute the periodic convolution ussing FFTs

Mh = real(ifft2(fft2(M) * fft2(h)))

Display

imageplot(M, "Image", [1, 2, 1])
imageplot(Mh, "Blurred", [1, 2, 2])
No description has been provided for this image

Several differential and convolution operators are implemented.

G = grad(M)
imageplot(G[:, :, 0], "d/ dx", [1, 2, 1])
imageplot(G[:, :, 1], "d/ dy", [1, 2, 2])
No description has been provided for this image

Fourier Transform

The 2D Fourier transform can be used to perform low pass approximation and interpolation (by zero padding).

Compute and display the Fourier transform (display over a log scale). The function fftshift is useful to put the 0 low frequency in the middle. After fftshift, the zero frequency is located at position $(n/2+1,n/2+1)$.

Mf = fft2(M)
Lf = fftshift(log(abs(Mf) + 1e-1))
imageplot(M, "Image", [1, 2, 1])
imageplot(Lf, "Fourier transform", [1, 2, 2])
No description has been provided for this image

Worked example 1: To avoid boundary artifacts and estimate really the frequency content of the image (and not of the artifacts!), one needs to multiply M by a smooth windowing function h and compute fft2(M*h). Use a sine windowing function. Can you interpret the resulting filter ?

# compute kernel h
t = linspace(-pi, pi, n)
h0 = matrix((cos(t) + 1) / 2)
h = array(transpose(h0) * h0)
# compute FFT
Mf = fft2(M * h)
Lf = fftshift(log(abs(Mf) + 1e-1))
# display
clf
imageplot(M * h, "Image", [1, 2, 1])
imageplot(Lf, "Fourier transform", [1, 2, 2])
No description has been provided for this image

Worked example 2: Perform low pass filtering by removing the high frequencies of the spectrum. What do you oberve ?

k = round(0.8 * n)
k = round(k / 2) * 2  # even number
Mf = fft2(M)
r = int(k / 2)
Mf[m - r + 2 : m + r, m - r + 2 : m + r] = 0
Mh = real(ifft2(Mf))
# display
clf
imageplot(M[m - 20 : m + 20, m - 20 : m + 20], "Image", [1, 2, 1])
imageplot(Mh[m - 20 : m + 20, m - 20 : m + 20], "Low pass filtered", [1, 2, 2])
No description has been provided for this image

References and further reading