This is one of the best books in image processing:

You should buy several copies, for yourself, your friends and your family; and ask your lab to buy several copies for the library. Or even better, since the editor (whose name I have trimmed from the cover image) is a notorious bully, instead of buying it, download a copy from libgen and send the money directly to the authors.

The book is a collection of independent self-contained chapters written by different authors, all of them famous people from the french school of image processing.

The only critique of this book that I can conceive is that the ``practice'' part of the title is not really fullfilled. There is not a single line of real computer code displayed in the book. But giving the codes for the hundreds of experiments of a 500 page book is probably too much to ask. The goal of this document is to provide such a code for a small part of the book.

My favourite chapters are 6 and 7:

**6. A short tour of mathematical morphology on edge and vertex weighted
graphs**, *Laurent Najman and Fernand Meyer*

**7. Partial difference quations on graphs for local and nonlocal image
processing**, *Abrerrahim Elmoataz, Olivier Lézoray, Vinh-Thong Ta and
Sébastien Bougleux*.

And these are the chapters whose implementation I detail below.

The book contains this kind of sentences: $\newcommand{\R}{\mathbf{R}}$

Let $G=(V,E)$ be a graph and $\mathcal{H}(V)$ be the Hilbert space of real-valued functions defined on the vertices of $G$. The space $\mathcal{H}(V)$ is endowed with the usual inner product $\left<f,h\right>_{\mathcal{H}(V)}=\sum_{v_i\in V}f(v_i)h(v_i)$, where $f,h:V\to\R$. Similarly, let $\mathcal{H}(E)$ be the Hilbert space of real-valued functions defined on the edges of $G$, $\ldots$. Now, consider a linear operator between Hilbert spaces $A:\mathcal{H}(V)\to\mathcal{E}(V)$...

While these sentences are crystal clear and very appealing to an audience of mathematicians, I have found them to be intimidating when trying to evangelize people to read the book. Thus, I ``translate'' them into the following kind of language, which is $100\%$ equivalent:

Consider a graph with $n$ vertices and $m$ edges. We will use vectors of length $n$ and $m$ to represent functions defined over the vertices or the edges, respectively. We will also use matrices of size $m\times n$ to represent linear maps between them $A:\R^n\to\R^m$. In octave/matlab:n = 100; # number of vertices in the graph m = 200; # number of edges in the graph x = rand(n,1); # define a random function over the vertices A = rand(m,n); # define a random linear map y = A * x; # obtain a function over the edges

This is easier to interpret thanks to the computer code. Of course, linear maps with random coefficients are silly. We will see more interesting examples below.

You would think that to work with graphs on a computer you need some sort of
library for graphs. Nothing farther from the truth. What you really need is
a library for doing *linear algebra*. In all the examples here we use
octave, but you can translate it easily to python+numpy, which is slightly
more verbose.

In what follows we reserve the letters $n$ and $m$ for the following meanings

$n\, =\ \,$ number of vertices in the graph

$m =\ \,$ number of edges in the graph

For the following graph, we have $n=5$ and $m=6$:

$$ Z = \begin{pmatrix} 1 & 2 \\ 1 & 3 \\ 2 & 3 \\ 3 & 4 \\ 3 & 5 \\ 4 & 5 \end{pmatrix}_{\small 6\times 2} A = \begin{pmatrix} 0 & 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 1 & 1 \\ 0 & 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 1 & 0 \\ \end{pmatrix}_{\small 5\times 5} B = \begin{pmatrix} -1& 1 & 0 & 0 & 0 \\ -1& 0 & 1 & 0 & 0 \\ 0 &-1 & 1 & 0 & 0 \\ 0 & 0 &-1 & 1 & 0 \\ 0 & 0 &-1 & 0 & 1 \\ 0 & 0 & 0 &-1 & 1 \\ \end{pmatrix}_{\small 6\times 5} $$

There are **five** basic matrices associated to a graph, for which we will
always use the same letters:

$Z\ \,=\ \,$ adjacency list, of size $m\times 2$,
list of vertex index pairs

$A\ \,=\ \,$ adjacency matrix, of size $n\times n$,
logical matrix of joined vertices

$B\ \,=\ \,$ incidence matrix, of size $m\times n$,
list of input/output vertices of each edge

$C\ \,=\ \,$ centering matrix, of size $m\times n$,
defined as $C=\frac{1}{2}|B|$

$L\ \,=\ \,$ Laplacian matrix, of size $n\times n$,
defined as $L=-B'B$

It is very important to understand now the meaning of the matrices $Z,A,B$. The matrix $Z$ is the easiest to type by hand in a computer, but it is not very useful for doing algebra with it. All the other matrices are fundamental linear operators over functions defined on the graph. Each of these fives matrices, alone, determines completely the graph (modulo the numbering of the edges, in the case of $A$ and $L$).

The matrices $A$, $B$ and $C$ satisfy the identity: $A=2C^TC-B^TB/2$. This allows to compute $A$ from $B$. To recover $B$ from $A$ we must decide on an ordering for the edges.

The following octave functions allow to convert between each representation:

function A = graph_adjacency_from_list(Z) n = max(Z(:)); # number of vertices U = sparse(Z(:,1), Z(:,2), 1, n, n); # directed graph A = U + U'; # symmetrization end

function B = graph_incidence_from_adjacency(A) [i,j] = find(triu(A)); # find the (i,j) positions n = rows(A); # number of vertices m = rows(i); # total number of edges B1 = sparse(1:m, i, 1, m, n); # matrix for destination vertices B2 = sparse(1:m, j, 1, m, n); # matrix for source vertices B = B1 - B2; # signed incidence matrix end

function Z = graph_list_from_adjacency(A) [i,j] = find(triu(A)); Z = [i j]; end

function A = graph_adjacency_from_incidence(B) A = max(-B'*B,0); # equal to (abs(B'*B) - B'*B) / 2 end

Typically, in the applications, you can often build $B$ directly so that you do not really need these functions. From $B$, the other matrices are easily computed if needed by:

L = -B'*B; A = L > 0; C = abs(B)/2;

The following table summarizes the language that we will use everywhere.

$n$ | number of vertices in the graph |

$m$ | number of edges in the graph |

$u\in\R^n$ | scalar field $u$ |

$\mathbf{v}\in\R^m$ | scalar field $\mathbf{v}$ |

$B:\R^n\to\R^m$ | gradient |

$-B^T:\R^m\to\R^n$ | divergence |

$L:\R^n\to\R^n$ | Laplacian |

$C:\R^n\to\R^m$ | centering operator (from vertices to edges) |

$C^T:\R^m\to\R^n$ | centering operator (from edges to vertices) |

$C^TC:\R^n\to\R^n$ | smoothing operator |

The most important notion is that the matrix $B$ is called the **gradient**. It is a linear operator that maps scalar fields (vectors of
lenght $n$) into vector fields (vectors of length $m$).
This definition is used for an arbitrary graph, but it makes a lot of sense
when the graph is the grid of an image, because in that case the gradient
corresponds exactly to the gradient computed using finite differences.

The pixels of an image are arranged naturally in the shape of a grid. Here, for example, you have the grid of a $4\times 3$ image:

\[\footnotesize \begin{array}{c|cccccccccccc} & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12\\ \hline 1 & -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 2 & 0 & -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 3 & 0 & 0 & -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 4 & 0 & 0 & 0 & 0 & -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 5 & 0 & 0 & 0 & 0 & 0 & -1 & 1 & 0 & 0 & 0 & 0 & 0\\ 6 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 1 & 0 & 0 & 0 & 0\\ 7 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 1 & 0 & 0\\ 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 1 & 0\\ 9 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 1\\ 10 & -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 11 & 0 & -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 12 & 0 & 0 & -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 13 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\ 14 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\ 15 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 1 & 0 & 0\\ 16 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 1 & 0\\ 17 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 1 \end{array} \]

Here you see that the graph has $n=12$ vertices and $m=17$ edges. In general, for an image of size $w\times h$, the graph will have $n=wh$ vertices and $m=(w-1)h+(h-1)w$ edges. The matrix $A$ of such a graph is build by the following octave code

function A = grid_adjacency(w, h) # build a grid graph WxH x = sparse(1:w-1, 2:w, 1, w, w); # path graph of length W y = sparse(1:h-1, 2:h, 1, h, h); # path graph of length H U = kron(y,speye(w)) + kron(speye(h),x); # kronecker sum A = U + U'; # symmetrization end

This works because the grid graph is the product graph of two paths, and the adjacency matrix of a product graph is the Kronecker sum of their matrices.

You can also build the incidence matrix directly, by a similar process

function B = grid_incidence(w, h) # grid graph WxH x = sparse(1:w-1, 2:w, 1, w-1, w) - speye(w-1,w); # path of length W y = sparse(1:h-1, 2:h, 1, h-1, h) - speye(h-1,h); # path of length H B = [ kron(speye(h),x) ; kron(y,speye(w)) ]; # kronecker union end

The graph defines just the *domain* of an image. We still need the
*data*. As a sample image, we will use the amazing portrait of Samuel
Beckett by Jane Bown:

`beckett.png`

The easiest operator to understand is the Laplacian. The following octave code thus reads an image, applies the laplacian operator, and saves the result.

x = imread("i/beckett.png"); # load input image [w, h] = size(x); # extract dimensions x = double(x(:)); # flatten image data into a vector A = grid_adjacency(w,h); # build graph adjacency matrix L = A - diag(sum(A)); # Laplacian matrix y = L * x; # Laplacian of the original image z = uint8(reshape(127-2*y,w,h)); # contrast change and reshape imwrite(z, "beckett-lap.png"); # save output image

`beckett-lap.png`

The typical color coding for looking at a laplacian is such that gray=zero, white=negative, black=positive. As expected, the laplacian enhances the edges and textures while setting the constant regions to zero.

By substracting the laplacian to the image, we ``sharpen'' the original image.

y = x - L * x; # image minus laplacian imwrite(reshape(uint8(y),w,h), "beckett-sharp.png"); # save output image

`beckett-sharp.png`

Conversely, we can smooth the image by adding multiples of the laplacian to it, iteratively. This amounts to approximation the solution of the heat equation on the graph:

S = speye(w*h) + L/4; # smoothing operator y = S^8 * x; # run 8 smoothing steps imwrite(reshape(uint8(y),w,h), "beckett-smooth.png"); # save output image

`beckett-smooth.png`

The morphological operations of **dilation** and **erosion** are
defined for functions over graphs, as the maximum and minimum value of
neighboring vertices. The other morphological operations are all defined in
terms of dilation and erosion:

$u$ | function defined on the vertices of the graph |

$d(u)$ | dilation (max value among neighboring vertices) |

$e(u)$ | erosion (min value among neighboring vertices) |

$d(e(u))$ | opening |

$e(d(u))$ | closing |

$d(u) - u$ | inner morphological gradient |

$u - e(u)$ | outer morphological gradient |

$\frac12(d(u) - e(u))$ | centered morphological gradient |

$\frac12(d(u) + e(u)) - u)$ | morphological laplacian |

$3u - d(u) - e(u)$ | morphological sharpening (image minus laplacian) |

$\frac12(d(u) + e(u))$ | morphological blur (image plus laplacian) |

$u - d(e(u))$ | top hat (image minus opening) |

$e(d(u)) -u$ | bottom hat (closing minus image) |

Due to the inequalities $e(u)\le u\le d(u)$, we can see that all these operations (except the laplacian) produce positive images. The morphological gradients are also called upwind and downwind derivatives.

Thanks to sparse matrices, the implementation of these operations is very easy. The crucial matrix here is the structuring element matrix $E$, defined as the adjacency matrix plus the identity $$ E = A + I_n $$

The matrices $A$ and $E$ define linear operators on the scalar fields of a graph. In the particular case when the graph is the grid of an image, they can be interpreted as shift-invariant filters defined by the following stencils: $$ \mathrm{stencil}(A)= \begin{array}{|c|c|c|}\hline0&1&0\\\hline1&0&1\\\hline0&1&0\\\hline\end{array} \qquad \qquad \mathrm{stencil}(E)= \begin{array}{|c|c|c|}\hline0&1&0\\\hline1&1&1\\\hline0&1&0\\\hline\end{array} $$

Of course, this interpretation is only valid far from the boundary of the image domain. However, bear in mind that the matrix product is a well-defined, so that the boundary condition is dealt with ``automatically'' by the form of the graph.

We start with the implementation of mathematical morphology for binary
images, which are also called *masks*.
This implementation is easier to understand than for the
general gray-scale case.
The **dilation** of a mask $m$ can be computed simply by multiplying the
mask by powers of the adjacency matrix. Actually, the **erosion** and the
**median filter** are also computed by the same operation, and they can be
extracted by thresholding the resulting gray-scale image at its minimum,
maximum or central values.

Notice that if you want to apply erosion a certain number of times, you can multiply by a power of the adjacency matrix! Morphology is a nonlinear operation, but the nonlinearity comes from the threshold, not from the matrix multiplication.

The following code shows the computation of 6 dilations and erosions of the same mask (a binarization of the Beckett portrait):

x = imread("i/beckett.png"); # load image [w,h] = size(x); # extract dimensions m = double(x(:) > 66); # flatten and binarize E = grid_adjacency(w,h) + speye(w*h); # structuring element y = E^6 * m; # apply 6 times y1 = y > 0; # extract dilation y2 = !(y < max(y)); # extract erosion y3 = y > max(y)/2; # extract medfilter f = 255 / max(y); # gray-scale factor imwrite(logical(reshape(m,w,h)), "beckett-bin.png"); # save binary mask imwrite(uint8(reshape(y,w,h)*f), "beckett-6gray.png"); # save gray-scale imwrite(logical(reshape(y1,w,h)), "beckett-6dil.png"); # save dilated mask imwrite(logical(reshape(y2,w,h)), "beckett-6ero.png"); # save eroderd mask imwrite(logical(reshape(y3,w,h)), "beckett-6med.png"); # save median filter

The implementation of morphology for gray-scale images is a bit more complicated. The trick of thresholding the matrix multiplication does not work anymore, because we are multiplying by a vector that has many different values (not only zeros and ones). What we need is a matrix$\times$vector multiplication that takes the maximum, instead of the sums along each row. After a few hours of head-scratching, I realized that this is achieved by multiplying the matrix $E$ by the diagonal matrix defined by the vector $x$ (which results in a sparse matrix) and then taking the maximum along rows of the resulting matrix. The code is simple, but maybe not evident:

function y = dilation(E, x) y = full(max(diag(x)*E))'; # maximum value along sparse rows end

Notice that this only works for positive-valued images (otherwise the "max" is perturbed by the zeros in the sparse matrix, which are not ignored). If we want to compute the minimum, this technique does not work because the minimum on a sparse matrix is always zero! Yet, we can ``cheat'' and implement the erosion by computing the dilation of an image in the negative:

function y = erosion(E, x) m = 1 + max(x); t = m - x; y = m - dilation(E, t); end

And finally, the implementation of all the morphological operations consists simply in copying the table above:

function y = opening(E,x) y = dilation(E,erosion(E,x)); end function y = closing(E,x) y = erosion(E,dilation(E,x)); end function y = egradient(E,x) y = x - erosion(E,x); end function y = igradient(E,x) y = dilation(E,x) - x; end function y = cgradient(E,x) y = (igradient(E,x) + egradient(E,x)) / 2; end function y = mlaplacian(E,x) y = (igradient(E,x) - egradient(E,x)) / 2; end function y = msharpen(E,x) y = x - mlaplacian(E,x); end function y = mblur(E,x) y = x + mlaplacian(E,x); end function y = tophat(E,x) y = x - opening(E,x); end function y = bothat(E,x) y = closing(E,x) - x; end

The following code tests all these operations:

x = imread("i/beckett.png"); [w,h] = size(x); x = double(x(:)); A = grid_adjacency(w,h); E = A + speye(w*h); # structuring element imwrite(uint8(reshape( dilation(E,x) ,w,h)), "beckett-dil.png"); imwrite(uint8(reshape( erosion(E,x) ,w,h)), "beckett-ero.png"); imwrite(uint8(reshape( opening(E,x) ,w,h)), "beckett-ope.png"); imwrite(uint8(reshape( closing(E,x) ,w,h)), "beckett-clo.png"); imwrite(uint8(reshape( 2*igradient(E,x) ,w,h)), "beckett-igrad.png"); imwrite(uint8(reshape( 2*egradient(E,x) ,w,h)), "beckett-egrad.png"); imwrite(uint8(reshape( 2*cgradient(E,x) ,w,h)), "beckett-cgrad.png"); imwrite(uint8(reshape( 127-4*mlaplacian(E,x) ,w,h)), "beckett-mlap.png"); imwrite(uint8(reshape( msharpen(E,x) ,w,h)), "beckett-msharp.png"); imwrite(uint8(reshape( mblur(E,x) ,w,h)), "beckett-mblur.png"); imwrite(uint8(reshape( 6*tophat(E,x) ,w,h)), "beckett-top.png"); imwrite(uint8(reshape( 255-6*bothat(E,x) ,w,h)), "beckett-bot.png");

It is an interesting exercise to look at these images and try to describe the
changes verbally. For example: **dilation** enlarges the light objects,
while **erosion** enlarges the dark ones. The **opening** removes the
small bright spots and ridges, while **closing** removes the small dark
spots and valleys. All the three **gradients** look like the euclidean
norm of the linear gradient. The **inner** and **outer** gradients are
very sharp, but are not centered, each one biased towards the dark or the
light side of the discontinuities. The **centered** gradient is the
average of the two, which is centered but less sharp. The **morphological
laplacian** is the difference of the inner and outer gradients; it looks
almost exactly like the linear laplacian despite being a combination of
nonlinear operations. The **morphological sharpening** is similar to the
linear one, based on the reverse heat equation. The **top-hat** keeps only
the britght spots and ridges (those that were removed by the opening), and
set to zero the rest of the image. The **bottom-hat** keeps only the dark
spots and valleys (those that were removed by the closing), and sets to zero
the rest of the image; thus the bottom-hat produces an image in the negative,
that must be negated to visualize it correctly.