Introduction to Image Processing§

author:David Coeurjolly

Value based Approaches§

Principles§

We only focus on values of the functions \, f:
S\subset\mathbb{Z}^n\rightarrow Q using stochastic representation

Important object: histograms

H: Q \rightarrow [0,1]

H(i) = \frac{Card \{ p\in S \,|\, f(p) =i \}}{Card S}

The image can be seen as a random variable with probability density function H

Exemple: S=[0,185]\times[0,85],\, Q=[0,255]

_images/contourS.png _images/histocontourS.png

Histogram based image processing§

Histograms sum up colorimetric information of the image without considering spatial relationships.

Histogram transformations§

Few examples:

Translation

Histogram remapping

Inversion

More generally

\phi: [0,1] \rightarrow [0,1]

H'(i) = H(\phi(i))

Propagation to image f'(p) = \phi(f(p)) <demo>

Histogram Equalization and Histogram Specifications§

In previous examples, \phi was prescribed. Now, we are looking for \phi to map histogram H (empiric or observed) to H' (model).

Example: Equalization

_images/len_dark.png _images/len_egal.png
_images/histolen_dark.png _images/histolen_egal.png

Equalization (bis)§

Be careful not always a good idea

_images/contourS-egal.png _images/histocontourS-egal.png

Histogram based Segmentation§

Image Segmentation classify pixels into classes such that intra-class pixels share the same visual properties (colorimetric information, geometrical properties …)

Let’s go back to initial image:

_images/contourS.png _images/histocontourS.png
_images/contourS-seuil.png _images/histo-seuil.png

Histogram Thresholding

Color transfer using Optimal Transport§

<slides>

Local Approaches§

Principles§

Take into account spatial relationships between pixel values

Key tool: convolution

Algorithmically

Fourier Space§

_images/len_std.png _images/len_spectrum.png

\mathcal{F}(u)=\int_{-\infty}^\infty f(x)\cdot e^{-i u x}dx
\qquad {f}(x)=\int_{-\infty}^\infty \mathcal{F}(u)\cdot e^{i u x}du

F_k=\sum _{{n=0}}^{{{\mathrm  {N}}-1}}x_n\cdot e^{{-2i\pi kn / N}}\qquad
x_n={\frac  {1}{{\mathrm  {N}}}}\sum _{{k=0}}^{{{\mathrm
{N}}-1}}{ {F}}_k\cdot e^{{2i\pi n k /{N}}}

(for 0\leqslant k<{\mathrm  {N}})

Few Fourier Space Properties§

(magnitude of Fourier transform leads to translation and rotational invariant descriptors)

Illustrations§

(http://www.cs.unm.edu/~brayer/vision/fourier.html)

_images/cosines.png _images/phase.png _images/brks_blks.png
Pure cosine signals Magnitude vs. Phase Real image example

Illustrations§

(http://www.cs.unm.edu/~brayer/vision/fourier.html)

_images/slant.png _images/window.png
Continuous vs. discrete FT Edge effect reduction through convolution

Periodic structure:

_images/slant_group.png

Illustrations§

(http://www.cs.unm.edu/~brayer/vision/fourier.html)

_images/lowpass.gif _images/highpass.gif
Lowpass filter Highpass filter

Elementary Image Filtering§

Be careful

\int_{-\infty}^\infty g(s,t)dsdt = 1

Smoothing / Low-pass filters

M_{3\times 3} = \frac{1}{9} \begin{bmatrix} 1 & 1 & 1\\1 &1 & 1\\1 & 1 & 1\\ \end{bmatrix}

G_{3\times 3} = \frac{1}{16} \begin{bmatrix} 1 & 2 & 1\\2 &4 & 2\\1 & 2 & 1\\ \end{bmatrix}

G_{5\times 5} = \frac{1}{256} \begin{bmatrix} 1 & 4 & 6 &4 & 1\\ 4 & 16 & 24 & 16 &4\\6 & 24 & 36 & 24 & 6\\ 4 & 16 & 24 & 16 &4\\1 & 4 & 6 &4 & 1\\ \end{bmatrix}

Mask Filtering Design§

G : Gaussian kernel approximation function

g(x,y) = \frac{1}{\sqrt{2\pi}\sigma}\exp^{-  \frac{x^2+y^2}{2\sigma} }

_images/Gaussian_2d.png

\Rightarrow binomial coefficients, ….

Examples§

M_{3\times 3} _images/len_std.png _images/len_moy.png _images/len_moybis.png
G_{3\times 3} _images/len_std.png _images/len_G.png _images/len_Gbis.png
G_{5\times 5} _images/len_std.png _images/len_G5.png _images/len_G5bis.png

Noise Models§

Match between filtering and noise model

Gaussian noise _images/len_bruitGaussien.png _images/len_bruitGaussienG.png
Specular noise salt/pepper _images/len_bruit.png _images/len_bruitG.png
Contour spread _images/len_bruit2.png _images/len_bruit2G.png

High pass filter§

Complementary filters to low-pass filters

M_{3\times 3} = \frac{1}{9} \begin{bmatrix} -1 & -1 & -1\\-1 &8 & -1\\-1 & -1 & -1\\ \end{bmatrix}

G_{3\times 3} = \frac{1}{16} \begin{bmatrix} -1 & -2 & -1\\-2 &12 & -2\\-1 & -2 & -1\\ \end{bmatrix}

_images/len_std.png _images/len_hautG.png

Non-linear filters§

Median filter

example

\begin{bmatrix} 12 & 13 & 24\\1 &30 & 43\\3  & 15 & 20\\ \end{bmatrix}

\Rightarrow 15

…but…

complexity for a size m domain?

Example§

_images/len_bruit.png _images/len_bruitMedian.png

Contour Detection and Differential Estimators§

Objectives§

Context

\Rightarrow loci if high gradient vector norm of the function f(x,y)

\Rightarrow Zero-crossing of the Laplacian \Delta f

_images/objects.png _images/NaturalTexture.png

We will focus later to segmentation process.

Image Gradient§

Definition

\nabla f(x,y) = \left (\frac{\partial f}{\partial x}(x,y),
\frac{\partial f}{\partial y}(x,y)\right)

Finite difference approximation

Prewitt_x = \frac{1}{3} \begin{bmatrix} -1 & 0 & 1\\-1 &0 & 1\\-1 & 0 & 1\\ \end{bmatrix}

Sobel_x = \frac{1}{4} \begin{bmatrix} -1 & 0 & 1\\-2 &0 & 2\\-1 & 0 & 1\\ \end{bmatrix}

Kirsch_x = \frac{1}{15} \begin{bmatrix} -3 & -3 & 5\\-3 &0 & 5\\-3 & -3 & 5\\ \end{bmatrix}

(other masks exist along specific directions \theta)

Image Gradient (bis)§

Amplitude = gradient vector norm

\|\nabla f\|_2 = \sqrt{\left(\frac{\partial f}{\partial x}\right)^2 +   \left(\frac{\partial f}{\partial y}\right)^2}

\|\nabla f\|_1 = \left |\frac{\partial f}{\partial x}\right| +   \left |\frac{\partial f}{\partial y} \right|

\|\nabla f\|_\infty = \max\left(\left |\frac{\partial f}{\partial x}\right|, \left |\frac{\partial f}{\partial y}\right|\right)

Orientation

\theta = atan\left( \frac{\frac{\partial f}{\partial y}}{\frac{\partial f}{\partial x}}\right)

\Rightarrow a contour can be characterized has pixels with high gradient vector norm

Example§

_images/objects.png _images/objects-prewittx.png _images/objects-prewitty.png

(Prewitt)

Mask Design§

Gradient of ta filterer image and Gaussian smoothing filters

\nabla (f*g) = \nabla f * g = f*\nabla g

\Rightarrow let us consider Gaussian kernel filter g(x,y) =
\frac{1}{\sqrt{2\pi}\sigma}\exp^{- \frac{x^2+y^2}{2\sigma} }, we have

\frac{\partial g(x,y)}{\partial x} = -\frac{x}{\sigma^2\sqrt{2\pi}\sigma}\exp^{-  \frac{x^2+y^2}{2\sigma} }

_images/gaussien.png _images/gaussien_x.png _images/gaussien_xx.png
g g’ g’‘

Separable Filters§

Let g(x,y) be a convolution filter (\equiv Impulse response of a filter). The filter is separable if

g(x,y) = g_x(x,y)g_y(x,y)

Thus:

f*g = g_x * (g_y*f)

and for partial derivatives:

\frac{\partial (f*g)(x,y)}{\partial x} = f(x,y)* \left( g_x(x)\frac{dg_y}{dy}(y)\right)

\Rightarrow direct consequences when implementing filters (only 1D convolutions)

i.e. Prewitt’s filter corresponds to finite difference approximation of the result of a constant smoothing filtering of f

Theoretical Design of a Contour Detector§

Canny’s criteria (1983)

Noise model + signal model + Criteria modeling + initial conditions \Rightarrow Optimal filter as solution of a PDE

2g(x) - 2\lambda_1g'(x) + 2\lambda_2g''(x) + \lambda_3 = 0

_images/shencastan.png _images/deriche.png
Shen-Castan : c\exp^{-\alpha|x|} Deriche : c(\alpha|x|+1)\exp^{-\alpha|x|}

The Gaussian filter is a good approximation of such optimal filters

Contour Extraction§

Contour extraction requires gradient vector norm (or zero-crossing laplacian) thresholding

or

\|\nabla f(x,y)\| > \sigma_2 and \exists
(x',y') neighbor of \exists (x,y) such that \|\nabla f(x',y')\| > \sigma_1

More Advanced Image Processing§

Bilateral Filter§

[Tomasi & Manduchi, 98]

Non-linear filter with a pair of spatial intensity kernels

I^*(x) = \frac{1}{W} \sum_{y \in \Omega_x}
I(y)f(|I(y)-I(x)|)g(|y-x|)

_images/Bilateral_Filter.jpg

Bilateral Filter (bis)§

_images/Bilateral2.png

(From [Bilateral Filtering: Theory and Applications Sylvain Paris, Pierre Kornprobst, Jack Tumblin, and Frédo Durand Foundations and Trends in Computer Graphics and Vision, 2009])

Question: How to efficiently implement bilateral filtering ?

Variational Approaches: Anisotropic diffusion§

[Perona, Malik]

with for example:

c\left(\|\nabla I\|\right)=e^{{-\left(\|\nabla I\|/K\right)^{2}}}

Anisotropic diffusion process that removes noise preserving edges

Rationale: for some g\mathbb{R}\rightarrow \mathbb{R} we want to minimize

E[I]={\frac  {1}{2}}\int _{{\Omega }}g\left(\|\nabla I(x)\|^{2}\right)\,dx

and thus

{\frac  {\partial I}{\partial t}}=-\nabla E_{I}={\mathrm
{div}}(g'\left(\|\nabla I(x)\|^{2}\right)\nabla I)

Examples§

_images/len_bruitG.png _images/peronna.jpg _images/peronna2.jpg

N.B Bilateral filtering is asymptotically close to Perona-Malik’s diffusion

Mumford-Shah functional§

Idea Variational formulation edge aware smoothing

\mathcal{MS} (K,u)= \underbrace{ \alpha \int_{{M\backslash K}} | u - g |^2 dx}_{\text{attachment term}}  + \underbrace{\int_{{M\backslash K}} | \nabla u |^2 dx}_{\text{smoothness term}} + {\underbrace{\lambda \mathcal{H}^{1}(K \cap M)}_{\text{discontinuities length}}}

But:

Ambrosio-Tortorelli functional§

\Gamma- convergence mathemtatical framework

\mathcal{AT}_\epsilon(u,v) = {\alpha} \underbrace{\int_{{M}} |u - g|^{{2}} dx}_{\text{attachment term}} + \underbrace{\int_{{M}} |{v} \nabla u|^{{2}} dx}_{\text{smoothness term}} + \underbrace{{\lambda}  \int_{{M}} {\epsilon} |\nabla {v}|^{{2}} + \frac{1}{4 {\epsilon}} |1 - {v}|^{{2}} dx}_{\text{discontinuities length}}

:-)

Examples§

_images/AT1.png

Examples§

_images/AT2.png

Examples§

_images/AT3.png

3D Examples§

_images/ATMesh1.png

3D Examples§

_images/ATMesh2.png

3D Examples§

_images/ATMesh3.png

Mathematical Morphology§

Principles§

[Matheron, Serra, …]

Idea

Operators§

Translation

X_p = \{ x + p\,|\, x\in X\}

Dilation by a structuring element B

\delta_B(X) =  X \oplus B = \bigcup_{x\in X} B_x = \bigcup_{b\in B} X_b

Erosion by a structuring element B

\epsilon_B(X) = X \ominus B = \bigcap_{b\in B} X_{-b}

Properties

with \check{B} = \{ -p \,|\, p \in B\} and X^c = E \setminus X

(X \oplus B)^c = X^c \ominus \check{B}

(X \ominus B)^c = X^c \oplus \check{B}

Illustrations§

Structuring element: Euclidean disc

_images/Dilation.png _images/Erosion.png

(blue: set X, gray circle: structuring element and cyan: result of the operators)

More operators§

Opening B

A \circ B = (A \ominus B)\oplus B \quad (\gamma)

Closing

A \bullet B = (A \oplus B)\ominus B \quad (\phi)

_images/Opening.png _images/Closing.png

Examples§

3x3 structuring element

_images/Illustration_morpho.png _images/Illustration_ouverture.png _images/Illustration_fermeture.png

Properties§

Properties

Useful tool for granulometric analysis

Sequence of increasing structuring elements B_k= B\oplus\ldots\oplus B k times

\gamma_k(X) = X \circ B_k

G_k = |\gamma_k(X)|

PS_k = G_k(X) - G_{k+1}(X)

G_k is called the granulometry function of X and PS_k the spectrum

Intuitive explanation

X is defined as the union of grains and G_k is the size of the set \gamma_k(X) defined by grains larger than k

Example§

_images/snapAl.png _images/snapSam.png
_images/al-result.png _images/M2DISCO81.png

Example§

_images/distrib2.png

Generalizations§

Operators on gray-level images

F,G: E\rightarrow E

(F\oplus G)(x) = \sup_{y\in E} \{ F(y) + G(x-y)\}

(F\ominus G)(x) = \inf_{y\in E} \{ F(y) - G(x-y)\}

Example

Generalization (bis)§

\begin{bmatrix}13 &  16 & 17\\
15 & 10 & 13\\
16 & 9 & 15\\
 \end{bmatrix} \oplus B =   \begin{bmatrix}16 &  17 & 17\\
15 & 15 & 13\\
16 & 16 & 15\\
 \end{bmatrix}

Example (input)§

_images/objects-deg.png

Example (output)§

_images/ob-dilation.png _images/erode.png
_images/ob-closing.png _images/ob-opening.png

Morphogical Gradient/Laplacian§

Basic “à-la” finite difference definition

_images/gradient_m.png _images/laplacian_m.png

Mathemtatical model

Operators acting on complete lattices (L, \leq)

Good morpholigical filters§

Principle

Given an specific image

_images/m_input.png _images/m_vert.png _images/m_horiz.png _images/m_union.png

[Soille]

Example§

_images/snap2.png