Image Processing Tutorials with Python &#174

Yann Gavet, Johan Debayle

► To cite this version:

Yann Gavet, Johan Debayle. Image Processing Tutorials with Python &#174. Spartacus-Idh, 316 p., 2019, 978-2-36693-068-9. emse-02469242

HAL Id: emse-02469242

https://hal-emse.ccsd.cnrs.fr/emse-02469242

Submitted on 13 May 2020

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

IMAGE PROCESSING
TUTORIALS
with PYTHON®

<figure><img src="image_3.png" /></figure>

Yann GAVET
Johan DEBAYLE
This is a human-readable summary of (and not a substitute for) the license, available at:
https://creativecommons.org/licenses/by-sa/4.0/legalcode

Your are free to:

**Share** - copy and redistribute the material in any medium or format

**Adapt** - remix, transform or build upon the material for any purpose, even commercially.

The licensor cannot revoke these freedoms as long as you follow the license
terms.

Under the following terms:

<figure><img src="image_4.png" /></figure>

**Attribution** — You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.

**Appropriate credit** — You must mention the authors and their host institution (MINES SAINT-ETIENNE).

No additional restrictions — You may not apply legal terms or technological measures that legally restrict others from doing anything the license permits.
## Credits

Cover image: Pepper & Carrot, David Revoy, www.davidrevoy.com.

## Contributors:

* Victor Rabiet

* Séverine Rivollier

* Place your name here and get special credit if you want to participate to this book.
<figure><img src="image_5.png" /></figure>
# TABLE OF CONTENTS

15

PART I

## Enhancement and Restoration

<table><tr><td>1</td><td>Introduction to image processing</td><td>17</td></tr><tr><td>2</td><td>Image Enhancement</td><td>35</td></tr><tr><td>3</td><td>2D Fourier Transform</td><td>47</td></tr><tr><td>4</td><td>Introduction to wavelets</td><td>57</td></tr><tr><td>5</td><td>Image restoration: denoising</td><td>73</td></tr><tr><td>6</td><td>Image Restoration: deconvolution</td><td>87</td></tr><tr><td>7</td><td>Shape From Focus</td><td>103</td></tr><tr><td>8</td><td>Logarithmic Image Processing (LIP)</td><td>113</td></tr><tr><td>9</td><td>Color LIP</td><td>123</td></tr><tr><td>10</td><td>GANIP</td><td>133</td></tr><tr><td>11</td><td>Image Filtering using PDEs</td><td>145</td></tr><tr><td>12</td><td>Multiscale Analysis</td><td>155</td></tr><tr><td>13</td><td>Introduction to tomographic reconstruction</td><td>165</td></tr></table>

173

PART II

## Mathematical Morphology

<table><tr><td>14</td><td>Binary Mathematical Morphology</td><td>175</td></tr><tr><td>15</td><td>Morphological Geodesic Filtering</td><td>187</td></tr><tr><td>16</td><td>Morphological Attribute Filtering</td><td>195</td></tr><tr><td>17</td><td>Morphological skeletonization</td><td>201</td></tr><tr><td>18</td><td>Granulometry</td><td>209</td></tr></table>
PART III

Registration and Segmentation

<table><tr><td>19</td><td>Histogram-based image segmentation</td><td>217</td></tr><tr><td>20</td><td>Segmentation by region growing</td><td>229</td></tr><tr><td>21</td><td>Hough transform and line detection</td><td>235</td></tr><tr><td>22</td><td>Active contours</td><td>243</td></tr><tr><td>23</td><td>Watershed</td><td>251</td></tr><tr><td>24</td><td>Segmentation of follicles</td><td>259</td></tr><tr><td>25</td><td>Image Registration</td><td>267</td></tr></table>

281

PART IV

Stochastic Analysis

<table><tr><td>26</td><td>Stochastic Geometry / Spatial Processes</td><td>283</td></tr><tr><td>27</td><td>Boolean Models</td><td>299</td></tr><tr><td>28</td><td>Geometry of Gaussian Random Fields</td><td>307</td></tr><tr><td>29</td><td>Stereology and Bertrand's paradox</td><td>317</td></tr><tr><td>30</td><td>Convex Hull</td><td>335</td></tr><tr><td>31</td><td>Voronoï Diagrams and Delaunay Triangulation</td><td>343</td></tr></table>

351

PART V

Image Characterization and Pattern Analysis

<table><tr><td>32</td><td>Integral Geometry</td><td>353</td></tr><tr><td>33</td><td>Topological Description</td><td>361</td></tr><tr><td>34</td><td>Image Characterization</td><td>367</td></tr><tr><td>35</td><td>Shape Diagrams</td><td>375</td></tr><tr><td>36</td><td>Freeman Chain Code</td><td>385</td></tr><tr><td>37</td><td>Machine Learning</td><td>395</td></tr><tr><td>38</td><td>Harris corner detector</td><td>403</td></tr><tr><td>39</td><td>Local Binary Patterns</td><td>409</td></tr></table>

417

PART VI

Exams
<table><tr><td>40 Practical exam 2016</td><td>419</td></tr><tr><td>41 Theoretical exam 2016</td><td>421</td></tr><tr><td>42 Theoretical exam 2017</td><td>425</td></tr></table>
<figure><img src="image_6.png" /></figure>
# About the authors

## Yann GAVET

received his "Ingénieur Civil des Mines de Saint-Etienne" diploma in 2001. He then obtained a Master of Science and a PhD thesis on the segmentation of human corneal endothelial cells (in 2004 and 2008). He is now an assistant professor at the Saint-Etienne School of Mines, where he teaches computer science and image processing to engineering and master students. He is a member of the PMDM Department of the LGF Laboratory, UMR CNRS 5307, dedicated to granular media analysis and modelisation.

He is particularly interested in the world of free (LIBRE) software in computer science. His research interests include image processing and analysis, stochastic geometry and numerical simulations. He published more than 70 papers in international journals and conference proceedings. He is a member of the Institute of Electrical and Electronics Engineers (IEEE), the International Association for Pattern Recognition (IAPR), International Society for Stereology and Image Analysis (ISSIA). He has worked for CS-SI (Toulouse, France) as an IT engineer, and for Thalès-Angénieux (Saint-Héand, France) as an image processing expert.

## Johan DEBAYLE

received his M.Sc., Ph.D. and Habilitation degrees in the field of image processing and analysis, in 2002, 2005 and 2012 respectively. Currently, he is a Full Professor at the Ecole Nationale Supérieure des Mines de Saint-Etienne (ENSM-SE) in France, within the SPIN Center and the LGF Laboratory, UMR CNRS 5307, where he leads the PMDM Department interested in image analysis of granular media. In 2015, he was a Visiting Researcher for 3 months at the ITWM Fraunhofer / University of Kaisersleutern in Germany. In 2017 and 2019, he was invited as Guest Lecturer at the University Gadjah Mada, Yogyakarta, Indonesia. He was also Invited Professor at the University of Puebla in Mexico in 2018 and 2019. He is the Head of the Master of Science in Mathematical Imaging and Spatial Pattern Analysis (MISPA) at the ENSM-SE.

His research interests include image processing and analysis, pattern recognition and stochastic geometry. He published more than 120 international papers in international journals and conference proceedings and served as Program committee member in several international conferences (IEEE ICIP, MICCAI, ICIAR...). He has been invited to give a keynote talk in several international conferences (SPIE ICMV, IEEE ISIVC, SPIE-IS&T EI, SPIE DCS...) He is Associate Editor for 3 international journals: Pattern Analysis and Applications (Springer), Journal of Electronic Imaging (SPIE) and Image Analysis and Stereology (ISSIA).

He is a member of the International Society for Optics and Photonics (SPIE), International Association for Pattern Recognition (IAPR), International Society
for Stereology and Image Analysis (ISSIA) and Senior Member of the Institute of
Electrical and Electronics Engineers (IEEE).

## ÉCOLE NATIONALE SUPÉRIEURE DES MINES DE SAINT-ÉTIENNE

One of the missions of École des Mines de Saint-Étienne, France, is scientific research at the highest level and contributions to companies' competitiveness. It aims at conveying the economical politics of the country and speeding up the sustainable industrial development by innovation and efficient contributions.

This high level scientific research leads to publications recognized by the inter-national scientific community. Research and teaching are very closely interwoven and the consequence of this is the attractiveness of our Master's Degree courses and our renowned Doctoral School.

<figure><img src="image_7.png" /></figure>

MINES
Saint-Étienne

Une école de l'IMT
# Liminar considerations

This book is presented as a collection of tutorials. Each chapter presents different objectives, usually with practical applications, in order to develop the image processing and analysis skills by its own. For each tutorial, you will first find in the statement different questions that will guide you to the complete results. Our suggestion is to start by the first tutorials, that acts as a round-up stage. Then, choose the topic that interests and motivates you, and start to work on your code.

The statements does not always contain the necessary theoretical background that you need to answer to the questions. Find resources on the net. Wikipedia is usually a very good start. Go to the university library!

The correction gives you a proposition of solution: it is not unique and it might not be perfect (or it might unfortunately contain errors). This correction is also available in a dedicated website (links will be provided along with the tutorials), we suggest you have a look at it only when you encounter blocking problems or fastidious lines of code to program. Moreover, you can evaluate the processing time of your version of the code and compare it to our proposed version, which usually uses fast and optimal functions, except for readability purposes.

With Python, you can use the following code:

<figure><img src="image_8.png" /></figure>

1. `import time`

3. `time_start = time.time()`

```
# run your code
```

5. `time elapsed = (time.time() - time_start)`

## CC-BY license?

If you manage to code a drastically faster method, if you notice errors or mistakes, if you want to add precisions, remember that this book as well as the code is published under a FREE CC-BY license (as in FREE speech). Contact us and we will be really happy to introduce modifications and insert you in the credits, or more...

<figure><img src="image_9.png" /></figure>

creativecommons

<figure><img src="image_10.png" /></figure>

The book is available in a high quality printed format, obviously not for free, but the pdf format is free (as in FREE beer). This is due to the fact that we (Johan Debayle and Yann Gavet) are employed by the French government, and we are paid to teach and disseminate informations to students. Our work is thus belongs to the society, and it seems a normal process to give the results back to the society. Feel free to use, modify and distribute these tutorials as you wish. Just remember to
keep the appropriate credits (our names and MINES Saint-Etienne). If you want to
express your gratitude, send us a postcard at:

Yann GAVET or Johan Debayle
MINES Saint-Etienne
CS 62362
42023 SAINT-ETIENNE cedex 2 - FRANCE

You can also buy the printed version of the book, which will make our editor
really happy.

## Images

The images belong to their authors, unless otherwise stated. If by mistake we used
images without giving the appropriate credit, please contact us.

## Softwares

This book is available for two programming languages, MATLAB® and Python.
MATLAB® is a proprietary software dedicated to scientific computing. Functions
are grouped into so-called toolboxes. The documentation and the compatibility of
the functions are very good, but the price is a consequence of this. Python is built
upon open-source and Free softwares. Scientific modules are numerous, but other
general purposes functions can also be found. Documentation and code may be of
inequal qualities, but major scientific modules (numpy, scipy, opencv...) are really
well presented and optimized.

Portions of code will be highlighted by the use of special boxes, like:

<figure><img src="image_11.png" /></figure>

% this is matlab code

<figure><img src="image_12.png" /></figure>

1 # this is python code

## Time to spend on a tutorial

The tutorials are almost independent. To evaluate the difficulty of the tutorial, stars are present at the header of each one. Depending on your knowledge, you will have to spend between 2 hours and 10 hours per tutorial.

* One star tutorial should be a relatively easy tutorial,
* two stars tutorials introduce some difficulties, either in programming or in the theoretical concepts,

* three stars tutorials are difficult both in theory and in programming.

They are divided into 6 parts, namely:

* Enhancement and Restoration: dedicated to image processing methods employed to eliminate noise and improve vision of data.

* Mathematical Morphology: introduction to basic operations of mathematical morphology.

* Registration and Segmentation: methods to segment, i.e. detect objects in images.

* Stochastic Analysis: method based on random processes.

* Characterization and Pattern Analysis: measures performed on objects in order to perform analysis and recognition.

* Exams: uncorrected problems and questions.

Enjoy!

You will find online corrections on the editor website and at page: http://iptutorials.science.
You will also find qrcodes for each correction that points to code and images.
<figure><img src="image_13.png" /></figure>
# Part I Enhancement and Restoration
<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /></figure>

# 1 Introduction to image processing

In this tutorial, you will discover the basic functions in order to load, manipulate and display images. The main informations of the images will be retrieved, like size, number of channels, storage class, etc. Afterwards, you will be able to perform your first classic filters.

The different processes will be realized on the following images:

<figure><img src="image_16.png" /><figcaption>Figure 1.1: Image examples.</figcaption></figure>

## 1.1 First manipulations

<figure><img src="image_17.png" /></figure>

Image loading can be made by the use of the python function `imageio.read`. The visualization of the image in the screen is realized by using the module ` matplotlib .pyplot`.

?

* Load and visualize the first image as below. Notice the differences.

* Look at the data structure of the image *I* such as its size, type....

* Visualize the green component of the image. Is it different from the red one? What is the most contrasted color component? Why?

* Enumerate some digital image file formats. What are their main differences? Try to write images with the JPEG file format with
different compression ratios (0, 50 and 100), as while as the lossless
compression, and compare.

1.2 Color quantization

Color quantization is a process that reduces the number of distinct colors used in an image, usually with the intention that the resulting image should be as visually similar as possible to the original image. In principle, a color image is usually quantized with 8 bits (i.e. 256 gray levels) for each color component.

<figure><img src="image_18.png" /></figure>

* By using the gray level image 'muscle', reduce the number of gray levels to 128, 64, 32, and visualize the different resulting images.

• Compute the different image histograms and compare.

1.3 Image histogram

An image histogram represents the gray level distribution in a digital image. The histogram corresponds to the number of pixels for each gray level.

The python function has the following prototype:

<figure><img src="image_19.png" /></figure>

def histogram(I):

<figure><img src="image_20.png" /></figure>

Compute and visualize the histogram of the image 'muscle.jpg'.

1.4 Linear mapping of the image intensi-
ties

The gray level range of the image 'cellules_cornee.jpg' can be enhanced by a linear
mapping such that the minimum (resp. maximum) gray level value of the resulting
image is 0 (resp. 255). Mathematically, it consists in finding a function $f(x) = ax+b$
such that $f(min) = 0$ and $f(max) = 255$.

<figure><img src="image_21.png" /></figure>

* Load the image and find its extremal gray level values.

* Adjust the intensities by a linear mapping into [0, 255].

* Visualize the resulting image and its histogram.

## 1.5 Aliasing effect

<figure><img src="image_22.png" /></figure>

* Create an image (as right) that contains rings as sinusoids. The function takes two input parameters: the sampling frequency and the signal frequency.

* Look at the influence of the two varying frequencies. What do you observe? Explain the phenomenon from a theoretical point of view.

<figure><img src="image_23.png" /></figure>

## 1.6 Low-pass filtering

The different processes will be realized on the following images:

<figure><img src="image_24.png" /><figcaption>(a) osteoblast cells</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(b) blood cells</figcaption></figure>
Low-pass filtering aims to smooth the fast intensity variations of the image to be processed.

<figure><img src="image_26.png" /></figure>

Test the low-pass filters 'mean', 'median', 'min', 'max' and 'gaussian' on the noisy image 'blood cells'.

<figure><img src="image_27.png" /></figure>

The module `scipy.ndimage filters` contains a lot of classical filters.

<figure><img src="image_28.png" /></figure>

Which filter is suitable for the restoration of this image?

## 1.7 High-pass filtering

High-pass filtering aims to smooth the low intensity variations of the image to be processed.

<figure><img src="image_29.png" /></figure>

* Test the high-pass filters `HP` on the two initial images in the following way: $HP(f) = f - LP(f)$ where `LP` is a low-pass filtering (see the previous exercise).

* Test the Laplacian (high-pass) filter on the two initial images with the following convolution mask:

$$ \begin{bmatrix} -1 & -1 & -1 \\ -1 & +8 & -1 \\ -1 & -1 & -1 \end{bmatrix} $$

## 1.8 Derivative filters

Derivative filtering aims to detect the edges (contours) of the image to be processed.
<figure><img src="image_30.png" /></figure>

* Test the Prewitt and Sobel derivative filters (corresponding to first order derivatives) on the image 'blood cells' with the use of the following convolution masks:

$$ \begin{bmatrix} -1 & 0 & +1 \\ -1 & 0 & +1 \\ -1 & 0 & +1 \end{bmatrix} \begin{bmatrix} -1 & -1 & -1 \\ 0 & 0 & 0 \\ +1 & +1 & +1 \end{bmatrix} $$

$$ \begin{bmatrix} -1 & 0 & +1 \\ -2 & 0 & +2 \\ -1 & 0 & +1 \end{bmatrix} \begin{bmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ +1 & +2 & +1 \end{bmatrix} $$

* Look at the results for the different gradient directions.

* Define an operator taking into account the horizontal and vertical directions.

Remark : the edges could be also detected with the zero-cressings of the Laplacian filtering (corresponding to second order derivatives).

## 1.9 Enhancement filtering

Enhancement filtering aims to enhance the contrast or accentuate some specific image characteristics.

<figure><img src="image_31.png" /></figure>

* Test the enhancement filter E on the image 'osteoblast cells' defined as: $E(f) = f + HP(f)$ where HP is a Laplacian filter (see tutorial 2.).

* Parameterize the previous filter as: $E(f) = \alpha f + HP(f)$, where $\alpha \in \mathbb{R}$.

## 1.10 Open question

Find an image filter for enhancing the gray level range of the image 'osteoblast cells'.
<figure><img src="image_32.png" /></figure>

## 1.11. Python correction

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

```
# display images
```

2 import matplotlib.pyplot as plt

```
# ndimage defines a few filters
```

4 from scopy import misc, ndimage

```
# numeric calculation
```

6 import numpy as np

```
# measure time
```

8 import time

```
# read and save images
```

10 import imageio

### 1.11.1 First manipulations

Open, write images

The following file loads the ascent image and display it. The `print` function is optional.

<figure><img src="image_35.png" /></figure>

```
# load file ascent
```

2 ascent = misc.ascent()

```
# load file cerveau.jpg
```

4 brain = imageioandravel('cerveau.jpg')

```
print(type(brain))
```

6 print(brain.shape, brain da type)

```
# save file
```

8 imageio.imwrite('test.png', brain)

Display images

You can modify to previous example to add the following lines:

<figure><img src="image_36.png" /></figure>

```
plt.imshow(ascent)
```

2 plt.show()

Notice that you have to close the image window to write commands again. Also, the colormap is not the good one by default (see Fig. 1.2).
<figure><img src="image_37.png" /></figure>

```python
plt.figure(figsize=(10, 3.6))

2 # first subplot
pltsubplot(131)

4 plt.imshow(ascent)
# second subplot
6 pltsubplot(132)
plt.imshow(ascent, cmap= plt.cm .gray)

8 plt .axis('off')
# third subplot (zoom)
10 pltsubplot(133)
plt ..imshow(ascent [200:220, 200:220],  cmap= plt.cm .gray, interpolation='nearest')

12 plt .subplots Adjust (wspace=0, hspace=0.,
                         top=0.99, bottom=0.01,
                         left=0.05, right=0.99)

14 plt.show()
```

<figure><img src="image_38.png" /><figcaption>Figure 1.2: Displaying images with an adapted colormap.</figcaption></figure>

Color channels

A color image is constituted of (generally) three channels. This representation
follows the human visual perception principles: in the human retina, the sensitive
cells (the cones) react to specific wavelength that correspond to red, green and blue.
The sensors technology adopted the *same* characteristics and a so-called Bayer filter
has 2 green filters for 1 red and 1 blue. Consequently, the green channel presents a
better resolution than the other channels.
Figure 1.3: The green channel presents the best contrasts in the case of retina images.

<figure><img src="image_39.png" /></figure>

## Image Resizing

The number of pixels is reduced by subsampling the image. Notice that there is no anti-aliasing filter applied to the image before reducing its size. The result is presented in Fig.1.4 with images of the same size, and in Fig.1.5 with images at the same resolution.

Figure 1.4: Reduction of the scale of the image on each axis, showing a so-called *pixellisation* effect. Represented at the same size, the density of pixels is thus reduced.

<figure><img src="image_40.png" /></figure>

Figure 1.5: At a constant resolution, images with different definitions are represented with different sizes.

<figure><img src="image_41.png" /></figure>
Color quantization

The following code uses the properties of integer operations to round values to the nearest integer. Illustration is presented Fig. 1.6.

<figure><img src="image_42.png" /></figure>

1 q4 = ascent // 4*4;

q16= ascent //16*16;

3 q32= ascent //32*32;

<figure><img src="image_43.png" /><figcaption>Figure 1.6: Reduction of the number of gray levels (quantization).</figcaption></figure>

JPEG file format

In order to test the effect of jpeg compression, one can use the parameter quality.
For the lowest quality, the loss of informations is really important (see Fig.1.7).

<figure><img src="image_44.png" /></figure>

# test jpeg quality

2 imageio.imwrite("a_25.python.jpeg", ascent, quality=25);
imageio.imwrite("a_100.python.jpeg", ascent, quality=100);

4 imageio.imwrite("a_50.python.jpeg", ascent, quality=50);
imageio.imwrite("a_75.python.jpeg", ascent, quality=75);

6 imageio.imwrite("a_1.python.jpeg", ascent, quality=1);
<figure><img src="image_45.png" /><figcaption>Figure 1.7: Different quality used to compress jpeg files.</figcaption></figure>

### 1.11.2 Histogram

To plot the histogram of an image, you can use the hist function (see result in Fig. 1.8).

<figure><img src="image_46.png" /></figure>

1. `hist(ascent flatten() , 256)`

2. `plt.show()`

<figure><img src="image_47.png" /><figcaption>Figure 1.8: Histogram of ascent image.</figcaption></figure>

You can also write your own code. The execution time is higher for these two functions than for the previous one.
<figure><img src="image_48.png" /></figure>

# Histogram function with 2D image

2 def compute_histogram(image):
    tab = np.zeros((256, ), dtype='I')

4 X, Y = image.shape
for i in range(X):

6 for j in range(Y):
    tab(image[i, j]]+=1

8

return tab

<figure><img src="image_49.png" /></figure>

1 # Histogram function with flatten image (vector)
def compute_histogram2(image):

3 im = image. flatten ()
tab = np zeros((256, ), ndtype='I')

5 for i in im:
    tab[i] += 1

7

return tab

The following code presents a comparison of the different methods for histogram computation.
<figure><img src="image_50.png" /></figure>

```python
1 # load ascent image and compute histograms
    ascent = misc.ascent()
3 t0 = time.clock()
    h = compute_histogram(ascent)
5 t1 = time.clock()
    h2=compute_histogram2(ascent)
7 t2 = time.clock()
    # .... plots
9 print "execution time 2D:%g" %(t1-t0)
    plt Subplot (131)
11 plt .plot (h)
    plt .title ( '2D function ')
13 print "execution time 1D:%g" %(t2-t1)
15 plt .subplot (132)
    plt .plot (h2)
17 plt .title ( '1D function ')
    # last plot : with matplotlib function
19 plt .subplot (133)
    t3 = time .clock()
21 plt .hist (ascent .flatten() , 256)
    t4 = time .clock()
23 print "execution time matplotlib:%g" %(t4-t3)
    # display
25 plt.show()
```

The console outputs the following computation durations:

<figure><img src="image_51.png" /></figure>

```
1 execution time 2D: 1.40284 s
    execution time 1D: 1.27649 s
3 execution time matplotlib : 0.261095 s
```

### 1.11.3 Linear mapping of the image intensities

The linear mapping is a simple method that stretches linearly the histogram. If displayed with `matplotlib`, the images is linearly stretched, thus the modification cannot be observed.
<figure><img src="image_52.png" /></figure>

def image_stretch(image):

2 # return image with new maximum and minimum at 255 and 0
minimum = image.min();

4 maximum = image.max();

a = 255/(maximum-minimum);

6 b = -255*minimum/(maximum-minimum);
return a*image+b;

1.11.4 Aliasing effect

<figure><img src="image_53.png" /></figure>

1 # aliasing effect (Moire)
def circle(fs, f):

3 # Generates an image with aliasing effect
# fs: sample frequency

5 # f : signal frequency

t = np.arange(0,1,1./fs);

7 ti, tj = np.meshgrid(t, t);

C = np.sin(2* np.pi * f * np.hypot(ti, tj))^2;

9 return C

The image of Fig. 1.9 is generated with the following code.

<figure><img src="image_54.png" /></figure>

1 C = circle(300,50);
 plt.imshow(C, cmap= plt.cm .gray);

3 plt.show()
imageio .imwrite(' moire.png ', C);
<figure><img src="image_55.png" /><figcaption>Figure 1.10: Mean filters.</figcaption></figure>

## Gaussian filter

The gaussian filter is presented in Fig. 1.11.

<figure><img src="image_56.png" /></figure>

1. `# ascent image`
   `ascent = misc.ascent()`

3. `# Gaussian filter`
   `gaussian = ndimage.filters.gaussian_filter(ascent, 5)`

5.

<figure><img src="image_57.png" /><figcaption>Figure 1.11: Gaussian filter of size 5.</figcaption></figure>

## 1.11.6 Derivative filters

Derivative filters (Prewitt, Sobel...) use a finite derivation approximation. They are very sensitive to noise (as every system using a derivation). The gradient is defined
as a vector, and a norm should be used to display a resulting image. Notice that with these filters, the connectivity of the contours is not preserved. Illustration is proposed in Fig. 1.12.

<figure><img src="image_58.png" /></figure>

1 # ascent image
ascent = misc.ascent()
3 ascent .astype('int32');
5 # Prewitt filter
prewitt0 = ndimage .filters .prewitt(ascent, axis=0)
7 prewitt1 = ndimage .filters .prewitt(ascent, axis=1)
9 # Sobel filter
dy = ndimage .filters .sobel(ascent, axis=0) # vertical
11 dx = ndimage .filters .sobel(ascent, axis=1) # horizontal
mag = np.hypot(dx, dy) # magnitude
13 sobel = mag * 255.0 / mag.max() # normalize (Q&D)
15 # display results
plt .subplot(131)
17 plt .imshow(prewitt0, cmap= plt .cm .gray)
plt .axis('off')
19 plt .title('Prewitt filter axis 0')
21 plt .subplot(132)
23 plt .imshow(prewitt1, cmap= plt .cm .gray)
25 plt .axis('off')
27 plt .title('Prewitt filter axis 1')
29 plt .subplot(133)
31 plt .imshow(sobel, cmap= plt .cm .gray)
33 plt .axis('off')
35 plt .title('Sobel filter')
37 plt .show()
<figure><img src="image_59.png" /><figcaption>Figure 1.12: Prewitt filter.</figcaption></figure>
<figure><img src="image_60.png" /></figure>
★ 2 Image Enhancement

The objective of this tutorial is to implement some image enhancement methods, based on intensity transformations or histogram modifications. It will make use of statistical notions like probability density functions or cumulative distribution functions.

<figure><img src="image_61.png" /><figcaption>Figure 2.1: The different processes of this tutorial will be applied on these images.</figcaption></figure>

## 2.1 Intensity transformations (LUT)

Two transformations will be studied: the $\gamma$ correction and the contrast stretching. They enable the intensity dynamics of the gray tone image to be changed. These two operators are based on the following Look Up Tables (LUT):

<figure><img src="image_62.png" /></figure>

<figure><img src="image_63.png" /></figure>
<figure><img src="image_64.png" /></figure>

1. Test the transformation 'γ correction' on the image 'osteoblast'.

2. Implement the operator 'contrast stretching' with the following LUT (also called cumulative distribution function cdf), with m being the mean gray value of the image, and r being a given gray value:

$$s = T(r) = \frac{1}{1 + (m/r)^E}$$

3. Test this transformation with different values of E on the image 'osteoblast'.

<figure><img src="image_65.png" /></figure>

The module skimage.exposure contains different methods for contrast enhancement, among them adjustgamma.

## 2.2 Histogram equalization

The objective is to transform the image so that its histogram would be constant (and its cumulative distribution function would be linear). The notations are:

* I is the image of n pixels, with intensities between 0 and L.

* h is the histogram, defined by:

$$h_I(k) = p(x = k) = \frac{n_k}{n}, \quad 0 \le k < L$$

The following transformation T(I) is called histogram equalization.

$$T(x_k) = (L - 1) \sum_{j=0}^{k} p(x_j)$$

where

$$cdf_I(k) = \sum_{j=0}^{k} p(x_j)$$

is the cumulative distribution function (cumulative histogram).

<figure><img src="image_66.png" /></figure>

1. Compute and visualize the histogram of the image 'osteoblast'.
2. Test this histogram equalization transformation on the image 'os-teoblast' (with builtin functions) and visualize the resulting histogram.

3. Code your own function.

4. The corresponding LUT to this transformation is the cumulative sum of the normalized histogram. Evaluate and visualize this intensity transformation.

<figure><img src="image_67.png" /></figure>

See the numpy.histogram and skimage.exposure functions.

### 2.3 Histogram matching

The objective is to enhance the original image by matching its histogram with a modeled one. The principle is to transform the gray value $x_1$ of first image into $x_2$: $T(x_1) = x_2$ (see Fig.2.2). Based on the fact that $cdf_1(x_1) = cdf_2(x_2)$, the formula to find $x_2$ is:

$$x_2 = cdf_2^{-1}(cdf_1(x_1)).$$

As $x_1$ and $x_2$ are discrete values, this requires an interpolation.

<figure><img src="image_68.png" /><figcaption>Figure 2.2: Histogram matching principle.</figcaption></figure>

<figure><img src="image_69.png" /></figure>

1. Visualize the histogram of the image 'phobos'.

2. Make the histogram equalization and visualize the resulting image.

3. Construct a bi-modal histogram (for example) and code your own function for histogram matching.
<figure><img src="image_70.png" /></figure>

See numpy.interp for interpolation and LUT application.
<figure><img src="image_71.png" /></figure>

## 2.4. Python correction

<figure><img src="image_72.png" /></figure>

<figure><img src="image_73.png" /></figure>

fromigliomisc

2import matplotlib .pyplot as plt
fromskimageimport exposure

4import numpy as np
import sys

### 2.4.1 Intensity transformations

γ correction

The function `skimage.exposure-adjustgamma` is used in order to adjust the gamma of the image. It is before normalized (values are float between 0 and 1). The results are illustrated in Fig.2.3.

<figure><img src="image_74.png" /></figure>

I=imageio.read("osteoblaste.png")

2I = I / np.max(I);

4gamma=2;

I2=exposure Adjustgamma(I, gamma);

6 plt .imshow(I2);
plt .show();
<figure><img src="image_75.png" /><figcaption>Figure 2.3: Gamma transforms.</figcaption></figure>

Contrast stretching

The LUT look-up-table is simply a function applied to the original image. In order to avoid division by zero, the smallest float value is introduced. Fig.2.4 illustrates the results.

<figure><img src="image_76.png" /></figure>

1. def contrast_stretching (I, E):
   epsilon = sys. float_info .epsilon ;

3 m = np.mean(I);

I = I.astype("float");

5 Ar = 1. / (1.+(m/(I+epsilon ))**E);

2.4.2 Histogram equalization

The histogram is displayed with the following function:
<figure><img src="image_77.png" /><figcaption>Figure 2.4: Contrast stretching.</figcaption></figure>

<figure><img src="image_78.png" /></figure>

def displaySaveHisto(I, filename=None):

Display and save pdf (if filename provided) of histogram of image I

If values are between 0 and 1, they are multiplied by 255

if np.max(I) <=1:

I = 255 * I;

hist, bins = np histogram(I. flatten(), 256, range=(0,255))

fig = plt. figure ()

plt. bar(bins[:-1], hist, width=1);

plt. show();

if filename!=None:

fig. savefig(filename, bbox_inches='tight')

The histogram equalization is done via the skimage.exposure. equalize_histogram function.
The results are illustrated in Fig.2.5.
<figure><img src="image_79.png" /></figure>

# histogram equalization

2 I2 = exposure. equalize_histogram(I);

plt.imshow(I2);

4 plt.show()

The look-up-table (LUT) or cumulative distribution function (cdf) is computed
and used in the next function:

<figure><img src="image_80.png" /></figure>

def histeq(I):

###

2 histogram equalization, version with look-up-table
I: original image, with values in 8 bits integer

4

###

6 hist, bins= np histogram(I, flatten(), 256, range=(0,255))

cdf = hist.cumsum();

8 cdf = (cdf / cdf[-1]);

return cdf[I];
<figure><img src="image_81.png" /><figcaption>Figure 2.5: Histogram equalization.</figcaption></figure>

### 2.4.3 Histogram matching

This method is only an approximate method. It requires an interpolation of the cumulative distribution functions in order to find the transformation. Results are shown in Fig.2.6
<figure><img src="image_82.png" /></figure>

1. `def hist_match(I, cdf_dest):`
2. """
3. Histogram matching of image I, with cumulative histogram cdf_dest
4. This should be normalized, between 0 and 1.
5. This version uses interpolation
6. """
7. `imhist`, `bins` = np.histogram(I, flatten(), len(cdf_dest), density=True)
8. `cdf` = `imhist`.cumsum() # cumulative distribution function
9. `cdf` = (cdf / cdf[-1]) # normalize between 0 and 1
10. # first : histogram equalization
11. `im2` = np.interp(I, flatten(), bins[:-1], cdf)
12. # 2nd: reverse function
13. `im3` = np.interp(im2, cdf_dest, bins[:-1])
14. # reshape into image
15. `imres` = im3.reshape(I.shape)
16. return imres;
<figure><img src="image_83.png" /><figcaption>Figure 2.6: Histogram matching.</figcaption></figure>
<figure><img src="image_84.png" /></figure>
<figure><img src="image_85.png" /></figure>

3 2D Fourier Transform

The main objective of this tutorial is to study image filters applied in the
frequencial or spatial domain with the Fourier transform.

<figure><img src="image_86.png" /></figure>

Use module np. fft for Fourier Transform functions (function fft2 ), phase
and abs for phase and amplitude.

The image to be used comes from a human cornea endothelium observed ex vivo
by optical microscopy.

Figure 3.1: Image of human cornea endothelium observed by optical microscopy.
The ophthalmologists would like to know the cell density, without manually count-
ing all the cells.

(a) Human corneal endothelium, extracted
from a donor and observed here before graft-
ing.

<figure><img src="image_87.png" /></figure>

(b) Human eye (from Wikipedia, authors: Rhcastilhos and Jmarchn, CC-By-SA).

<figure><img src="image_88.png" /></figure>
## 3.1 Fourier transform

?

1. Load an image and visualise it.

2. Compute the Fourier Transform by the fft algorithm.

3. Visualise the images of the phase and amplitude of the Fourier Transform.

## 3.2 Inverse Fourier transform

In this exercise, it can be interesting to consider different images, like a Lena picture for example.

?

1. Apply the inverse Fourier transform on the Fourier transform to find the original image.

2. Now, apply the inverse Fourier transform on the phase information only (without using the frequency informations).

3. In a same spirit, apply the inverse Fourier transform on the frequency informations only (without the phase).

## 3.3 Low-pass and high-pass filtering

?

1. Modify the Fourier transform of the image to

* keep only low frequencies,

* keep only high frequencies.

2. Apply the inverse Fourier transform on both and comment.
3.4 Application: evaluation of cellular den-
sity

The ophthalmologists would like to evaluate the cell density of the cornea endothelium observed in Fig. 3.1.

<figure><img src="image_89.png" /></figure>

1. Compute the Fourier transform of the image. If one considers that the cells constitute a repeated pattern on the whole image, locate the repetition frequency on the amplitude image.

2. Can this frequency be linked to the cell density ?

The answer to this last question is obviously yes. Here follows a simple method to evaluate the cell density (or the mean cell radius), if these cells a considered as circular [37].

<figure><img src="image_90.png" /></figure>

1. The amplitude information is very noisy. First of all, a gaussian filter should be applied (see functions fspecial and imfilter ).

2. Find a simple way to evaluate the mean radius of the cells.
<figure><img src="image_91.png" /></figure>

## 3.5. Python correction

<figure><img src="image_92.png" /></figure>

<figure><img src="image_93.png" /></figure>

import numpy as np

2 fromigliimport misc, ndimage #read/write images

import matplotlib .pyplot as plt # plots

### 3.5.1 Introduction and utilities

To display a spectrum, the following function can be useful:

<figure><img src="image_94.png" /></figure>

1 # Displayspectrum and phase in an image (grayscale)

def viewSpectrumPhase(amplitude, phase):

3 plt. figure ()

5 plt.subplot(1,2,1)

5 plt.imshow(np.log(1+amplitude), plt.cm.grayscale);

7 plt.subplot(1,2,2)

mmax = np.max(phase);

9 mmin = np.min(phase);

if (mmax == mmin):

11 B=0;

else :

13 B = 255*(phase-mmin)/(mmax-mmin);

15 plt.imshow(B, cmap= plt.cm .gray);

### 3.5.2 Fourier transform

Results (amplitude and phase) are represented in Fig. 3.4.

<figure><img src="image_95.png" /></figure>

cornea = imageio .imread('cornee.png')

2 print type(cornea)

print cornea.shape, cornea .dtype

4

 plt.subplot(131)

6 plt.imshow(cornea, cmap= plt.cm .gray);
<figure><img src="image_96.png" /></figure>

### Fourier transform

2. `# result is complex`

3. `# fftshift` is used by convention the get frequency 0 at center of image
4. `spectre = npsons.fshift(npsons.fshift2(cornea))`

6. `A = abs(spectre);`

7. `G = npangle(spectre);`

8. `viewSpectrumPhase(A, G);`

### 3.5.3 Inverse Fourier Transform

<figure><img src="image_97.png" /></figure>

1. `## inverse Fourier transform`
   `cornee2 = np.real(npsons.fshift(npsons.fshift2(npsons.fshift(spectre))));`
2. `plt.figure()`
3. `plt.imshow(cornee2, cmap= plt.cm .gray);`
4. `plt .title('inverse FT');`

The Fourier transform, although very powerful, is really difficult to interpret in 2D. The main information is not contained in the amplitude, but in the phase. The following reconstructions will illustrate it (see Fig.3.2).

<figure><img src="image_98.png" /></figure>

1. `## inverse Fourier transform, without the phase`
   `cornee_amplitude = np.real(npsons.fshift(npsons.fshift2(npsons.fshift(A))));`
2. `plt .figure()`
3. `plt .show(cornee_amplitude, cmap= plt.cm .gray);`
4. `plt .title('inverse FT on amplitude');`

<figure><img src="image_99.png" /></figure>

1. `## inverse Fourier transform on phase only`
   `complex_phase = np.exp(1j*G);`
2. `cornee_phase = np.real(npsons.fshift(npsons.fshift2(complex_phase)));
   plt .figure()`
3. `plt .show(cornee_phase, cmap= plt.cm .gray);`
4. `plt .title('inverse FT on phase');`
<figure><img src="image_100.png" /><figcaption>Figure 3.2: Reconstruction of partial informations (phase or amplitude only).</figcaption></figure>

### 3.5.4 Low-pass and high-pass filtering

The results are illustrated in Fig. 3.3. Two functions are defined.

A low pass-filter consists in the suppression of values for frequencies lower than a cut-off frequency.

<figure><img src="image_101.png" /><figcaption>Figure 3.3: Fourier basic filtering of cornea image.</figcaption></figure>
<figure><img src="image_102.png" /></figure>

```python
def LowPassFilter(spectrum, cut):
    """Low pass filter of the FFT (spectrum)
    """
    The shape of this filter is a square. fftshift has been applied so that
    frequency 0 lays at center of spectrum image
    @param spectrum: FFT2 transform
    @param cut : cut value of filter (no physical unit, only number of pixels)
    """
    X,Y = spectrum.shape;
    mask = np.zeros((X,Y), "int");
    mx = X/2; my = Y/2;
    mask[mx-cut:mx+cut, my-cut:my+cut] = 1;
    f = spectrum * mask;
    plt.figure
    plt.imshow(abs(f)); plt.title('Low pass filter')
    return f;
```

A high pass filter is exactly the opposite: suppression of the values under the cut-off frequency.

<figure><img src="image_103.png" /></figure>

```python
1 def HighPassFilter(spectrum, cut):
    """High pass filter of the FFT (spectrum)
    """
    The shape of this filter is a square. fftshift has been applied so that
    frequency 0 lays at center of spectrum image
    @param spectrum: FFT2 transform
    @param cut : cut value of filter (no physical unit, only number of pixels)
    """
    X,Y = spectrum.shape;
    mask = np.ones((X,Y), "int");
    mx = X/2; my = Y/2;
    mask[mx-cut:mx+cut, my-cut:my+cut] = 0;
    f = spectrum * mask;
    plt.figure
    plt.imshow(abs(f)); plt.title('High pass filter')
    return f;
```

In the following application, the image is loaded and the effects of a low-pass and high-pass filters are illustrated.
<figure><img src="image_104.png" /><figcaption>Figure 3.4: The two peaks represent the frequency of repetition of the cellular pattern, i.e. it can be used to compute the cellular density.</figcaption></figure>
<figure><img src="image_105.png" /></figure>
<figure><img src="image_106.png" /></figure>

# 4 Introduction to wavelets

This tutorial introduces practically the basic wavelet decomposition and reconstruction algorithms. The objectives are to code some basic programs that will decompose and reconstruct 1D and 2D signals.

## 4.1 Introduction

Wavelets are based on a mother wavelet $\Psi$ ($s$ is a scale parameter, $\tau$ is the time translation factor ($s, \tau) \in \mathbb{R}_+^* \times \mathbb{R}$):

$$ \forall t \in \mathbb{R}, \psi_{s,\tau}(t) = \frac{1}{\sqrt{s}} \Psi\left(\frac{t-\tau}{s}\right) $$

The continuous wavelet transform is written as follows, where $\psi^*$ means the complex conjugate of $\psi$ and $< \cdot, \cdot >$ is the complex scalar product (Hermitian form):

$$ g(s, \tau) = \int_{-\infty}^{\infty} f(t) \psi_{s,\tau}^{*}(t) ds d\tau = \langle f, \psi_{s,\tau} \rangle $$

The reconstruction is defined by:

$$ f(t) = \frac{1}{C} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{1}{|s|^2} g(s, \tau) \psi_{s,\tau}(t) ds d\tau $$

with

$$ C = \int_{-\infty}^{\infty} \frac{|\hat{\Psi}(\omega)|^2}{|\omega|} d\omega $$

and $\hat{\Psi}$ is the Fourier Transform of $\Psi$.

The discrete wavelet transform corresponds to a sampling of the scales. To compute the different scales, one has to introduce a “father” wavelet, that similarly defines a family of functions orthogonal to the family $\psi_{s,\tau}$.

## 4.2 Fast discrete wavelet decomposition / reconstruction

A simple algorithm (cascade algorithm, from Mallat) is defined as two convolutions (by a lowpass *ld* filter for the projection on the $\psi$ family, and a high pass *hd* filter for the orthogonal projection) followed by a subsampling (see Fig. 4.1).
<figure><img src="image_107.png" /><figcaption>Figure 4.1: Algorithm for wavelet decomposition. First, a convolution is performed (square node), then, a subsampling. $a_i$ stands for decomposition at scale i, $d_i$ stands for detail at scale i.</figcaption></figure>

### 4.2.1 Simple 1D example

Let's consider the signal [4; 8; 2; 3; 5; 18; 19; 20]. We will use the Haar wavelets defined by $ld = [1; 1]$ and $hd = [-1; 1]$. Basically, these filters perform a mean and a difference (see Table 4.1). The result of the decomposition in 3 scales with these wavelets is the concatenation of the details and the final approximation

$$C = \{[-4; -1; -13; -1]; [7; -16]; [-45]; [79]\}$$

<table><caption>Table 4.1: Illustration of Haar decomposition in a simple signal.</caption><thead><tr><th>Scale <i>i</i></th><th>Approximation <i>a<sub>i</sub></i></th><th>Details <i>d<sub>i</sub></i></th></tr></thead><tbody><tr><td>0 (original signal)</td><td>[4; 8; 2; 3; 5; 18; 19; 20]</td><td>[-4; -1; -13; -1]</td></tr><tr><td>1</td><td>[12; 5; 23; 39]</td><td>[7; -16]</td></tr><tr><td>2</td><td>[17; 62]</td><td>[-45]</td></tr><tr><td>3</td><td>[79]</td><td></td></tr></tbody></table>

<figure><img src="image_108.png" /></figure>

In the case of these Haar wavelets, code a function that performs the decomposition for a given number of scales. The prototype of this function will be:

1. def simpleWaveDec(signal, nb_scales):

3. wavelet decomposition of <signal> into <nb_scales> scales
This function uses Haar wavelets for demonstration purposes.

5.
### 4.2.2 Reconstruction

To reconstruct the original signal (see Fig. 4.2), we need the definition of two reconstruction filters, *hr* and *lr*. For the sake of simplicity, we use $lr = ld/2$ and $hr = -hd/2$. These filter will perform an exact reconstruction of our original signal.

Figure 4.2: Reconstruction algorithm. The oversampling is done by inserting zeros.

<figure><img src="image_109.png" /></figure>

Using this algorithm, you can now go on to the next exercise.

Code a function that performs the reconstruction of the signal. The prototype of this function will be, with the previous definition of C:

<figure><img src="image_110.png" /></figure>

1. def simpleWaveRec(C):

2.     wavelet simple reconstruction function of a 1D signal
3.     C: Wavelet coefficients

5.     """

## 4.3 2D wavelet decomposition

Let A be the matrix of an image. We consider that A is of size $2^n \times 2^n$, $n \in \mathbb{N}$. We consider, as for the 1D transform, the filters *ld* and *hd*.

The wavelet decomposition is as follows:

* Apply *ld* and *hd* on rows of A. Results are denoted *ld*<sub>*r*</sub>A and *hd*<sub>*r*</sub>A, of size $2^n \times 2^{n-1}$, with *r* standing for row.

* Then, apply *ld* and *hd* again, to get the four new matrices: *ld*<sub>*c*</sub>*ld*<sub>*r*</sub>A, *ld*<sub>*c*</sub>*hd*<sub>*r*</sub>A, *hd*<sub>*c*</sub>*ld*<sub>*r*</sub>A and *hd*<sub>*c*</sub>*hd*<sub>*r*</sub>A, of sizes $2^{n-1} \times 2^{n-1}$, with *c* standing for column. The matrix *ld*<sub>*c*</sub>*ld*<sub>*r*</sub>A is the approximation, and the other matrices are the details.
?

* Code a function to perform the wavelet decomposition and for a given number of scales (lower than *n*). Test it on images of size $2^n \times 2^n$.

* Code the reconstruction function.

## 4.4 Built-in functions

Let us explore some built-in functions.

### 4.4.1 Continuous wavelet decomposition

One has to make the difference between the continuous wavelet transform and the discrete transform.

<figure><img src="image_111.png" /></figure>

#### Informations

In `scipy` signal you can find function to perform wavelets decomposition, and among them `cwt` for the continuous wavelet transform. There is also the module `pywt` that may be more developed (see also `cwt`).

<figure><img src="image_112.png" /><figcaption>Figure 4.3: Continuous wavelet transform.</figcaption></figure>
4.4.2 Discrete decomposition

?

1. Generate the signal $sin(2π * f₁ * t) + sin(2π * f₂ * t)$, with $f₁ = 3Hz$ and $f₂ = 50Hz$.

2. Apply the wavelet transformation for a given wavelet (like 'db4').

3. Display the 5th approximation of the wavelet decomposition.

4. Display the 4th detail coefficients of the wavelet decomposition.

<figure><img src="image_113.png" /></figure>

Look at `pywt.downcoef` for the decomposition at a given level and `pywt.wavedec` for the entire decomposition.
<figure><img src="image_114.png" /></figure>

## 4.5. Python correction

<figure><img src="image_115.png" /></figure>

<figure><img src="image_116.png" /></figure>

import numpy as np

2 from scipy import misc

import matplotlib .pyplot as plt

### 4.5.1 1D signals

Two functions are required: a function (simpleWaveDec) that loops over the different scales and calls the second function (waveSingleDec) that performs the single step wavelet decomposition. Notice that the Haar wavelet is defined here with integer values (see Fig.4.4), so that the mental computation can be done easily.

<figure><img src="image_117.png" /><figcaption>Figure 4.4: Haar wavelets. From wikipedia, author Omegatron.</figcaption></figure>
## Simple 1D decomposition

<figure><img src="image_118.png" /></figure>

```python
def simpleWaveDec(signal, nb_scales):
    """
wavelet decomposition of <signal> into <nb_scales> scales
This function uses Haar wavelets for demonstration purposes.
    """
    # Haar Wavelets filters for decomposition and reconstruction
    ld = [1, 1];
    hd = [-1, 1];

    # transformation
    C=[];
    A = signal; # approximation
    for i in range nb_scales:
        A, D = waveSingleDec(A, ld, hd);
        # get the coefficients
        C.append(D);

    C.append(A);
    return C;
```

2

4

6

8

10

12

14

16

18

<figure><img src="image_119.png" /></figure>

1

3

5

7

9

11

13

15

## Simple 1D reconstruction

The reconstruction starts from the highest scale and computes the approximation signal with the given details.
<figure><img src="image_120.png" /></figure>

1 def simpleWaveRec(C):

###

3 wavelet simple reconstruction function of a 1D signal

### C: Wavelet coefficients

5 The Haar wavelet is used

###

7 ld = np.array([1, 1]);

9 hd = np.array([-1, 1]);

9 lr = ld/2;

hr = -hd/2;

11

A = C[-1];

13

for scale in reversed(C[:-2]):

### A = waveSingleRec(A, scale, lr, hr);

15

return A;

<figure><img src="image_121.png" /></figure>

1 def waveSingleRec(a, d, lr, hr):

###

3 1D wavelet reconstruction at one scale

### a: vector of approximation

5 d: vector of details

### lr: low pass filter defined by wavelet

7 hr: high pass filter defined by wavelet

This is Mallat algorithm.

9 NB: to avoid side effects, the convolution function does not use the
'same' option

###

11

approx = np zeros((len(a)*2,));

13

approx[:2] = a;

approx = np convolve(approx, lr);

15

detail = np zeros((len(a)*2,));

17

detail[:2] = d;

detail = np convolve( detail , hr);

19

# sum up approximation and details to reconstruct signal at lower scale

approx = approx + detail ;

21

#get rid of last value

approx = np.delete(approx, -1)

23

return approx

## Results

This is the result for the decomposition of the vector with 3 scales.
<figure><img src="image_122.png" /></figure>

1. `s = [4, 8, 2, 3, 5, 18, 19, 20];`
   `print(s)`

2. `C = simpleWaveDec(s, 3);`
   `print(C)`

3. `srec = simpleWaveRec(C);`
   `print(srec);`

<figure><img src="image_123.png" /></figure>

4. `[4, 8, 2, 3, 5, 18, 19, 20]`
2. `[array([-4, -1, -13, -1]), array([ 7, -16]), array([-45]), array([79])]`

### 4.5.2 2D signals

#### Decomposition

The `simpleImageDec` function is the main interface. It takes the image as first parameter, and the number of scales of decomposition. It makes a call to `decWave2D` for the decomposition at one given scale. The latter uses the previous 1D decomposition method.
<figure><img src="image_124.png" /></figure>

def decWave2D(image, ld, hd):

###

2

% wavelet decomposition of a 2D image into four new images.

4

% The image is supposed to be square, the size of it is a power of 2 in the % x and y dimensions.

###

6

# Decomposition on rows

8

sx, sy = image.shape;

10

LrA = np zeros((sx, int(sy/2)));

12

HrA = np zeros((sx, int(sy/2)));

14

for i in range(sx):

A, D = waveSingleDec(image[i:i], ld, hd);

16

LcLrA[i] = A;

HcLrA[i:;] = D;

18

# Decomposition on cols

20

LcLrA = np zeros((int(sx/2), int(sy/2)));

HcLrA = np zeros((int(sx/2), int(sy/2)));

22

LcHrA = np zeros((int(sx/2), int(sy/2)));

HcHrA = np zeros((int(sx/2), int(sy/2)));

24

for j in range(int(sy/2)):

A, D = waveSingleDec(LrA[j:j], ld, hd);

26

LcLrA[j:j] = A;

HcLrA[j:j] = D;

28

A, D = waveSingleDec(HrA[j:j], ld, hd);

30

LcHrA[j:j] = A;

HcHrA[j:j] = D;

32

return LcLrA, HcLrA, LcHrA, HcHrA
<figure><img src="image_125.png" /></figure>

```python
def simpleImageDec(image, nb_scales):
    """
wavelet decomposition of <image> into <nb_scales> scales
This function uses Haar wavelets for demonstration purposes.
"""

#Haar Wavelets filters for decomposition and reconstruction
ld = [1,1];
hd = [-1, 1];

#transformation
C=[];
A = image; # first approximation

coeffs = [];
for i in range nb_scales:
    [A, HcLrA, LcHrA, HcHrA] = decWave2D(A, ld, hd);
    coefs.append(HcLrA);
    coefs.append(LcHrA);
    coefs.append(HcHrA);
    #set the coefficients
    C.append(coefs.copy());
    coefs.clear();
    C.append(A);
    return C;
```

## 2D reconstruction

The `simpleImageRec` function performs the reconstruction of a multiscale wavelet decomposition. The `recWave2D` performs the reconstruction of one scale.
<figure><img src="image_126.png" /></figure>

```python
def recWave2D(LcLrA, HcLrA, LcHrA, HcHrA, lr, hr):
    """
```

Reconstruction of an image from `lr` and `hr` filters and from the wavelet
$\leftrightarrow$ decomposition.

A: resulting (reconstructed) image

NB: This algorithm supposes that the number of pixels in x and y dimensions is
$\leftrightarrow$ a power of 2.

```
sx, sy = LcLrA.shape;
```

# Allocate temporary matrices

```
LrA = np.zeros((sx*2, sy));
```

```
HrA = np zeros((sx*2, sy));
```

```
A = npENG zeros((sx*2, sy*2));
```

#Reconstruct from cols

for j in range(sy):

```
LrA[:,j] = waveSingleRec(LcLrA[:,j], HcLrA[:,j], lr, hr);
```

```
HrA[:,j] = waveSingleRec(LcHrA[:,j], HcHrA[:,j], lr, hr);
```

# Reconstruct from rows

for i in range(sx*2):

```
A[i::] = waveSingleRec(LrA[i::], HrA[i::], lr, hr);
```

return A;

<figure><img src="image_127.png" /></figure>

```python
def simpleImageRec(C):
    """
```

wavelet reconstruction of an image described by the wavelet coefficients C

"""

#The Haar wavelet is used

```
ld = np.array([1, 1]);
```

```
hd = np.array([-1, 1]);
```

```
lr = ld/2;
```

```
hr = -hd/2;
```

```
A = C[-1];
```

for scale in reversed(C[:-1]):

```
A = recWave2D(A, scale[0], scale[1], scale[2], lr, hr);
```

return A;
### Results

The illustration Fig. 4.5 is obtained by the following code. The useful functions are presented below.

<figure><img src="image_128.png" /></figure>

```python
I = imageio.imread('lena256.png');
2 C = simpleImageDec(I, 3);
Irec = simpleImageRec(C);
4 plt.imshow(Irec);
plt.show()
```

The image is recursively split into 4 areas, where the left upper corner is the approximation, and the three others are the details. As the details can have negative values, the intensities are adjusted in order to display the image correctly.

<figure><img src="image_129.png" /></figure>

```python
1 def displayImageDec(C):
    """
3 Construct a single image from a wavelet decomposition
    C: the decomposition
        """
5 n, m = C[0][0].shape;
7 A = np.zeros((2*n, 2*m));
    prev = C[-1];
9 for s, scales in reversedEnumerate(C[:-1]):
    11 ns = n / 2**(s-1);
        ms = m / 2**(s-1);
    13 T = imdec2im(prev, scales);
        A[0:int(ns), 0:int(ms)] = T;
    15 prev = A[0:int(ns), 0:int(ms)];
    return prev
```

These two functions `adjust` and `reversedEnumerate` are used to simplify the notations. The first one performs a linear stretching of the image intensities, and the second one allows an enumeration of a list in a reverse order.
<figure><img src="image_130.png" /></figure>

def adjust(I):

2

### Adjust

simple image intensity stretching

4

return I

### Return

6

I = I - np.min(I);

I = I / np.max(I);

8

return I;

10

def reversedEnumerate(l):

### Reversed

Utility function to perform reverse enumerate of a list

12

returns zip

### Return

14

return zip (range(len(1)-1, -1, -1), 1 [:: -1]);

<figure><img src="image_131.png" /></figure>

1

def imdec2im(LcLrA, lvLC):

### Imdec

3

constructs a single image from:

LcLrA: the approximation image

5

lvLC: the wavelet decomposition at one level

7

for display purposes

### Display

9

HcLrA=lvLC[0];

11

LcHrA=lvLC[1];

HcHrA=lvLC[2];

13

n, m = HcLrA.shape;

15

A = np(zeros((2*n, 2*m)));

17

# Approximation image can be with high values when using Haar coefficients

A[0:n, 0:m] =adjust(LcLrA);

19

# details are low, and can be negative

21

A[0:n, m:2*m] = adjust(HcLrA);

A[n:2*n, 0:m] = adjust(LcHrA);

23

A[n:2*n, m:2*m] = adjust(HcHrA);

25

return A;
<figure><img src="image_132.png" /><figcaption>Figure 4.5: (Haar) Wavelet decomposition of the Lena image.</figcaption></figure>

4.5.3 Built-in functions

An interesting module is `pywt`. This is illustrated by the following example. This gives the same results as previously, with the multiplication by $1/\sqrt{2}$.

<figure><img src="image_133.png" /></figure>

import pywt

2 cA, cD = pywt.dwt(s, 'haar');
print(cA, cD);

<figure><img src="image_134.png" /></figure>

1 [ 8.48528137 3.53553391 16.26345597 27.57716447]
[-2.82842712 -0.70710678 -9.19238816 -0.70710678]

The continuous wavelet transform is applied in the following code and the
result is displayed in Fig.4.6.
<figure><img src="image_135.png" /></figure>

```python
t = np.linspace(-1, 1, 1000, endpoint=False)

f1 = 3;
f2 = 50;

sig = np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t);
widths = np.arange(1, 129)

coef, freqs = pywt.cwt(sign, widths, 'morl')
fig = plt.figure();

plt.imshow(coef, extent=[-1, 1, 1, 31], cmap='PRGn', aspect='auto',
           vmax=abs(coef).max(), vmin=-abs(coef).max())

plt.show()
fig savefig('cwt.python.pdf', bbox_inches='tight')
```

<figure><img src="image_136.png" /><figcaption>Figure 4.6: Continuous wavelet decomposition of the sum of sinusoids.</figcaption></figure>
<figure><img src="image_137.png" /></figure>

## 5 Image restoration: denoising

This tutorial aims to study some random noises and to test different image restoration methods (image denoising).

The different processes will be applied on the following MR image.

<figure><img src="image_138.png" /><figcaption>Figure 5.1: Leg.</figcaption></figure>

### 5.1 Generation of random noises

Some random noises are defined with the given functions (corresponding to specific distributions):

* uniform noise: $R = a + (b - a) * U(0, 1)$.

* Gaussian noise: $R = a + b * N(0, 1)$.

* Salt and pepper noise: $R : \begin{cases} 0 \le U(0,1) \le a & \mapsto 0 \\ a < U(0,1) \le b & \mapsto 0.5 \\ b < U(0,1) \le 1 & \mapsto 1 \end{cases}$

* Exponential noise: $R = -\frac{1}{a} * \ln(1 - U(0, 1))$

<figure><img src="image_139.png" /></figure>

Generate sample images with the four kinds of random noise and visualize their histograms with the built-in functions. Pay attention to the intensity range of the resulting images when calculating the histograms
## 5.2 Noise estimation

The objective is to evaluate the characteristics of the noise in a damaged/noisy image (in a synthetic manner in this tutorial).

<figure><img src="image_140.png" /></figure>

1. Visualize the histogram of a Region Of Interest (ROI) of the image of Fig.5.1. The ROI should be extracted from a uniform (intensity) region.

2. Add an exponential noise to the original image and visualize the histogram of the selected ROI.

3. Add a Gaussian noise to the original image and visualize the histogram of the selected ROI.

## 5.3 Image restoration by spatial filtering

<figure><img src="image_141.png" /></figure>

1. Add a salt-and-pepper noise to the image of Fig.5.1.

2. Test the 'min', 'max', 'mean' and 'median' image filters.

<figure><img src="image_142.png" /></figure>

See the module `scipy.ndimage`

The median filter is efficient in the case of salt-and-pepper noise. However, it replaces every pixel value (first problem) by the median value determined at a given
scale (second problem). In order to avoid these two problems, Gonzalez and Woods proposed an algorithm [11].

**Data:** Input (noisy) image *I*

**Data:** Maximal scale $S_{max}$

**Result:** Filtered image *F*

**Main Function amf (I, S<sub>max</sub>):**

forall pixels (i, j) do

$S \leftarrow 1$

while isMedImpulseNoise(I, i, j, S) AND $S \le S_{max}$ do

$S \leftarrow S + 1$

$med \leftarrow Med(I, i, j, S)$

end

if $I(i, j) = \text{Min}(I, i, j, S)$ OR $I(i, j) = \text{Max}(I, i, j, S)$ OR

$S = S_{max}$ then

$F(i, j) \leftarrow Med(I, i, j, S)$

else

$F(i, j) \leftarrow I(i, j)$

end

end

**Algorithm 1:** Adaptive median filter [11]. The main function takes two arguments: the image *I* to be filtered and the maximal size of neighborhood $S_{max}$. For each pixel (*i*, *j*), the good scale is the scale where the median value is different from the min and the max (median value is thus not an impulse noise). Then, if the pixel value is an impulse noise or the maximal scale has been reached, it must be filtered. Otherwise, it is kept untouched. Min, Max and Med functions compute the minimum, maximum and the median value in a neighborhood of size *S* centered at a pixel (*i*, *j*) in image *I*.

**Function isMedImpulseNoise(I, i, j, S):**

$med \leftarrow Med(I, i, j, S)$

if $med = \text{Max}(I, i, j, S)$ OR $med = \text{Min}(I, i, j, S)$ then

return True

else

return False

end

**Algorithm 1:** End.
?

Implement the adaptive median filter from the following algorithm (Alg.1)
based on two steps (the operator acts on an operational window of size $k \times k$
where different statistics are calculated: median, min or max).
<figure><img src="image_143.png" /></figure>

## 5.4. Python correction

<figure><img src="image_144.png" /></figure>

<figure><img src="image_145.png" /></figure>

import matplotlib, plot as plt

2 import numpy as np
from scipy import ndimage

4 from scipy import misc

### 5.4.1 Generation of random noises

In order to convert stretch the images into the range [0;1], the following function can be used:

<figure><img src="image_146.png" /></figure>

def hist_stretch(I):

2 # histogram stretching
# returns values of I linearly stretched to range [0;1]

4 I = I - np.min(I);
I = I / np.max(I);

6 return I;

Random noise illustrations are proposed in Figs.5.2 and 5.3, and a size S = 32 is used to generate the noisy images.

#### Uniform noise

The python `rand` function generates values between 0 and 1 with uniform distribution.

<figure><img src="image_147.png" /></figure>

S=32

2 a=0; b=255;
R1 = a + (b-a) * np.randomrand(S, S);

#### Gaussian noise

The python `randn` function generates values with normal centered distribution.
<figure><img src="image_148.png" /><figcaption>Figure 5.2: Resulting noise images.</figcaption></figure>

<figure><img src="image_149.png" /></figure>

a=0; b=1;

2 R2 = a + (b-a)*np.random.random(S);
## Salt and pepper noise

<figure><img src="image_150.png" /></figure>

a = 0.05; b = 0.1;

2 R3 = 0.5 * np.ones((S,S));
X = np.randomrand(S,S);

4 R3[X<=a] = 0;
R3[(X>a) & (X <= b)] = 1;

6 R3 = hist_stretch(R3);

## Exponential noise

<figure><img src="image_151.png" /></figure>

a=1;

2 R4 = -1/a * np.log(1-np.randomrand(32, 32));
R4 = hist_stretch(R4);

<figure><img src="image_152.png" /><figcaption>Figure 5.3: Histograms of the noise images generated with 32 × 32 pixels.</figcaption></figure>
### 5.4.2 Noise estimation

In order to estimate the noise, a ROI of visually constant gray level is chosen, and its histogram is displayed. This is simulated by the following code, the result is displayed in Fig.5.4:

<figure><img src="image_153.png" /></figure>

1. `# noise estimation`

fig = plt.figure ()

2. `A = imageio.read('jambe.png');`

roi = A[160:200, 200:240];

3. `plt hist(roi. flatten(), 255);`

fig. savefig("histo_roi_leg.pdf", bbox_inches='tight');

Then, some noise is added to the image, and the histogram of the ROI is displayed. The results are observed in Figs.5.5 and 5.6.

<figure><img src="image_154.png" /></figure>

1. `# add exponential noise to image`

2. `nx, ny = A.shape;`

expnoise = -1./5 * np.log(1-np.randomrand(nx, ny));

3. `expnoise = expnoise / np.max(expnoise);`

B = A + 255*expnoise;

4. `B = hist_stretch(B);`

fig = plt. figure ()

5. `plt.imshow(B, cmap='gray');`

imageio imwrite("leg_exponential.png", B);

10. `# plt. show()`

<figure><img src="image_155.png" /></figure>

1. `# add gaussian noise to image`

2. `nx, ny = A.shape;`

gaussnoise = 50* np.randomrand(nx, ny);

3. `B = A + gaussnoise;`

B = hist_stretch(B);

4. `fig = plt. figure ()`

plt..imshow(B, cmap='gray');

5. `imageio. imwrite("leg_gaussian.png", B);`
<figure><img src="image_156.png" /><figcaption>Figure 5.4: Histogram of the Region of Interest.</figcaption></figure>

In the case of exponential and Gaussian noise added to the image, The histograms are displayed in Figs.5.5 and 5.6.

<figure><img src="image_157.png" /><figcaption>Figure 5.6: Gaussian noise.</figcaption></figure>

<figure><img src="image_158.png" /></figure>
### 5.4.3 Image restoration by spatial filtering

The following code is used to filter the images. The results are displayed in Fig.5.7. Salt and Pepper noise is added to the image.

Figure 5.7: Different filters applied to the noisy image. The median filter is particularly adapted in the case of salt and pepper noise (impulse noise), but still destroys the structures observed in the images.

<figure><img src="image_159.png" /><figcaption>(a) Original image.</figcaption></figure>

<figure><img src="image_160.png" /><figcaption>(b) Noisy image (salt and pepper).</figcaption></figure>

<figure><img src="image_161.png" /><figcaption>(c) Maximum filter.</figcaption></figure>

<figure><img src="image_162.png" /><figcaption>(d) Minimum filter.</figcaption></figure>

<figure><img src="image_163.png" /><figcaption>(e) Mean filter.</figcaption></figure>

<figure><img src="image_164.png" /><figcaption>(f) Median filter.</figcaption></figure>
<figure><img src="image_165.png" /></figure>

```python
a = 0.05; b = 0.05;
2 spnoise = 0.5 * np.ones((nx,ny));
X = np.randomrand(nx,ny);
4 B = A.copy();
B[X<=a] = 0;
6 B[(X>a) & (X <= (a+b))] = 255;
fig = plt.figure()
8 plt.imshow(B, cmap='gray');
imageio imwrite("leg_sp.png", B);
10 plt.show()
```

Then, different filters are illustrated. Notice that the contours informations are somehow damaged. The median filter corresponds to the best filter in this case.

<figure><img src="image_166.png" /></figure>

### filtering by convolution

2 **average filtering**

B1 = ndimage. uniform_filter (B,5);
4 imageio.savwrite("leg_uniform.png", B1);
B2 = ndimage. minimum_filter(B,3);
6 imageio.savwrite("leg_minimumpng", B2);
B3 = ndimage. maximum_filter(B,3);
8 imageio.savwrite("leg_maximumpng", B3);
B4 = ndimage. median_filter (B, 7);
10 imageio.savwrite("legarfng", B4);
B5 = amf(B, 7);
12 plt. والسB5, cmap='gray');
imageio.savwrite("legamf.png", B5);

What can be noticed is that min and max filters are unable to restore the image (opening and closing filters, from the mathematical morphology, could be a solution to explore). The mean filter is a better solution, but an average value is highly modified by an impulse noise. The median filter is the optimal solution is order to suppress the noise, but fine details are lost. An interesting solution is to apply an adaptive median filter.

## Adaptive median filter

The following code is a simple implementation of the algorithm previously presented. It has not been optimized in order to have a simple presentation. A more sophisticated version can be found in [11]. The results are illustrated in Fig.5.8 for $S_{max} = 7$.
<figure><img src="image_167.png" /></figure>

```python
def amf(I, Smax):
    """
    Adaptive median filter
    I: grayscale image
    Smax: maximal size of neighborhood. Limits the effect of median filter
    """
    f = np.copy(I);
    nx, ny = I.shape;
    sizes = np.arange(1, Smax, 2);
    zmin = np(zeros((nx, ny, len(sizes)));
    zmax = np(zeros((nx, ny, len(sizes)));
    zmed = np(zeros((nx, ny, len(sizes)));
    for k, s in enumerate(sizes):
        zmin[k] = ndimage minimum filter(I, s);
        zmax [k] = ndimage maximum_filter(I, s);
        zmed [k] = ndimage median_filter (I, s);
    isMedImpulse = np. logical_or (zmin==zmed, zmax==zmed);
    for i in range(nx):
        for j in range(ny):
            k = 0;
            while k < len(sizes) - 1 and isMedImpulse[i, j, k]:
                k += 1;
                if I[i, j] == zmin[i, j, k] or I[i, j] == zmax[i, j, k] or k == len(sizes):
                    f[i, j] = zmed[i, j, k];
    return f;
```
<figure><img src="image_168.png" /><figcaption>Figure 5.8: Illustration of the adaptive median filter.</figcaption></figure>
<figure><img src="image_169.png" /></figure>
<figure><img src="image_170.png" /></figure>

## 6 Image Restoration: deconvolution

This tutorial aims to study and test different image restoration methods for blurred and noisy images. Although the different software tools contain functions for each filter the objective is to exhaustively code these methods to fully understand their principles. The reader can refer to [36] for more informations on the algorithms.

Some of the images are based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA). You may access to these resources at https://hla.stsci.edu/

The different processes will be applied on the following images.

### 6.1 Damage modeling

A damage process can be modeled by both a damage function $D$ and an additive noise $n$, acting on an input image $f$ for producing the damaged image $g$.

$$g = D(f) + n \qquad (6.1)$$

Knowing $g$, the objective of restoration is to get an estimate $\hat{f}$ (the restored image) of the original image $f$. If $D$ is a linear process, spatially invariant, then the equation can be simplified as:

$$g = h * f + n \qquad (6.2)$$

where $h$ is the spatial representation of the damage function (called the Point Spread Function - PSF), and $*$ denotes the convolution operation. In general, the more knowledge you have about the function $H$ and the noise $n$, the closer $\hat{f}$ is to $f$.

This equation can be written in the frequency domain:

$$G = H.F + N \qquad (6.3)$$

where the letters in uppercase are the Fourier transforms of the corresponding terms in the eq. 6.2.
Figure 6.1: Images that can be used for testing the different restoration algorithms. Intensity levels are represented in false colors.

<figure><img src="image_171.png" /><figcaption>(b) PSF for Jupiter image</figcaption></figure>

<figure><img src="image_172.png" /><figcaption>(c) Saturn, HST.</figcaption></figure>

<figure><img src="image_173.png" /><figcaption>(d) PSF for Saturn image</figcaption></figure>

<figure><img src="image_174.png" /></figure>

<figure><img src="image_175.png" /></figure>

• Generate the 'chessboard' image.

• Generate a PSF corresponding to a blur (motion or isotropic).
<figure><img src="image_176.png" /><figcaption>Figure 6.2: Chessboard image.</figcaption></figure>

* Create a damaged image (a circular option is used, because the FFT implies that functions are periodical) and add a Gaussian noise.

* Visualize the different images.

## 6.2 Simple case: no noise

The objective is now to restore the damaged image.

For the first steps, no noise will be added. The degradation functions is $g = h \ast f$. Knowing $H$ (Fourier Transform of the psf $h$) and ignoring the noise in the damage model, a simple estimate of the initial image can be done by the inverse filtering:

$$f = FT^{-1}\left(\frac{G}{H}\right) \qquad (6.4)$$

where $FT^{-1}$ denotes the inverse Fourier transform.

<figure><img src="image_177.png" /></figure>

* Try the inverse filter in the Fourier space. Why is this filter not working?

* Try to prevent division by 0, with $f = FT^{-1}\left(\frac{G}{H+\alpha}\right)$, $\alpha$ being a small constant value (e.g. $\alpha = 0.1$).

## 6.3 Wiener filter

The Wiener filter is an optimal filter that minimizes the following error:

$$e^2 = E\{(f - \hat{f})^2\}$$
where *E* denotes the expected (mean) value, *f* is the original image and $\hat{f}$ is the restored image.

The solution to this expression, within the frequency domain, is given by:

$$ \hat{F} = \underbrace{\left( \frac{1}{H} \frac{|H|^2}{|H|^2 + R} \right)}_{H_w} G $$

The value of *R* can be arbitrarily fixed. Another way is to define $R = S_n/S_f$, where $S_n$ and $S_f$ denote the power spectrum of the noise *n* and the original image *f*, i.e. $|N|^2$ and $|F|^2$ (matrices), respectively. This ratio can be defined globally (as the constant) or locally (i.e. it is computed for each pixel). If the ratio $S_n/S_f$ (function) is non null, it can be passed to the Wiener filtering process.

### 6.3.1 Simple case: no noise

In the noiseless case, the Wiener filter reduces to the inverse filter:

$$ H_w = \begin{cases} \frac{1}{H(u, v)} & \text{if } |H(u, v)| \neq 0 \\ 0 & \text{otherwise} \end{cases} $$

<figure><img src="image_178.png" /></figure>

Try a Wiener filter for the noiseless case.

### 6.3.2 Noisy images

Generally, the ratio $S_n/S_f$ used in the Wiener filter is replaced by a constant value equal to the ratio of the mean power spectrum:

$$ R = \frac{\frac{1}{PQ} \cdot \sum_{u} \sum_{v} S_n(u, v)}{\frac{1}{PQ} \cdot \sum_{u} \sum_{v} S_f(u, v)} \qquad (6.5) $$

where *PQ* denotes the matrix size (the number of elements).

<figure><img src="image_179.png" /></figure>

1. Calculate the constant ratio *R* and code the Wiener filter.
2. Test the different algorithms on the image 'brain' damaged by a blur and an additive noise.

# 6.4 Iterative filters for noisy images

## 6.4.1 Van Cittert iterative filter

The VanCittert iterative filter is based on the following iterative relation :

$$f_{k+1} = f_k + \beta(g - h * f_k)$$

with $f_0 = g$, the damaged image.

<figure><img src="image_180.png" /></figure>

Code and test on different cases a function that performs the VanCittert restoration filter.

It should have the following prototype:

<figure><img src="image_181.png" /></figure>

function [ fr ] = vca( g, psf, n_iter )

2 % Van Cittert iteration algorithm for deconvolution %

4 % g: degraded image % psf: point spread function

6 % n_iter : number of iterations

<figure><img src="image_182.png" /></figure>

def vca(g, psf, n_iter, beta=.01):

2 Van Cittert iterative filter

4 g: noisy image psf: point spread function

6 n_iter : number of iterations beta: Jansson parameter

8 """
### 6.4.2 Lucy-Richardson filtering

We are now going to use the nonlinear restoration filtering of Lucy-Richardson that is based on a maximum likelihood formulation in which the image is modeled by Poisson statistics. This optimization leads to the solution if the following iterative equation converges:

$$f_{k+1}(x, y) = f_k(x, y) \cdot \left( h(-x, -y) * \frac{g(x, y)}{h(x, y) * f_k(x, y)} \right) \quad (6.6)$$

<figure><img src="image_183.png" /></figure>

Code and test your own function for the Lucy-Richardson filtering process. It should have the same prototype as the VanCittert function. Remember that `*` is the classical multiplication (point to point) and that `*` is the convolution operation. You then have the choice to make a direct call to the convolution function, or to perform this operation in the Fourier space.

## 6.5 Blind deconvolution: restoration with no a priori

If the PSF is unknown, one of the main difficulties in image restoration is to get an accurate estimate of this PSF in order to use it in the deconvolution algorithms. These methods are called blind deconvolution methods. Generally, the PSF is estimated in an iterative way from an initial PSF.

<figure><img src="image_184.png" /></figure>

1. Generate a damaged image of 'chessboard' with a Gaussian blur and an additive Gaussian noise.

2. Restore from a blind deconvolution the damaged image and compare the results with the previous ones.

<figure><img src="image_185.png" /></figure>

Informations

Use the function `skimage restoration unsupervised_wiener` for example.
6.6 Astronomy images

A classical file format used in astronomy is the FITS format.

<figure><img src="image_186.png" /></figure>

With python, fits images can be manipulated with the astropy module.
You can use the following code to display an image.

<figure><img src="image_187.png" /></figure>

hdu_list = fits.open('ic3g01qlq_drz.fits');

2 image = hdu_list [1].data;

image = np.nan_to_num(image);

4 plt.imshow(image, cmap='gray')

plt.show()
<figure><img src="image_188.png" /></figure>

## 6.7. Python correction

<figure><img src="image_189.png" /></figure>

<figure><img src="image_190.png" /></figure>

from astropy.io import fits

2. import matplotlib.pyplot as plt
import numpy as np

4. from Kapteyn.misc import imresize
from Kapteyn import signal

6. import progressbar

### 6.7.1 Damage modeling

There is no built-in function to generate the checkerboard image. A simple code can be used.

<figure><img src="image_191.png" /></figure>

def checkerboard():

2. arr = np.zeros((8,8), type=int)
   arr[::2,::2] = 1

4. arr[1::2,1::2] = 1

6. out = imresize(arr, 8*arr.shape, interp='nearest') / 255;

return out;

The psf can be generated as a Gaussian psf, or a motion psf. In this tutorial, the motion psf is loaded from a file that comes from the MATLAB® function fspecial.

<figure><img src="image_192.png" /></figure>

C = checkerboard();

2. psf = np.load('psfirc6.mpy');

In order to add Gaussian noise, the following function is used:
<figure><img src="image_193.png" /></figure>

```python
def addGaussianNoise(I, sigma=1000):
    """ Adds Gaussian noise to an image
    I: original image
    sigma: signal/noise ratio
    """
    I2 = I.copy();
    m= np.min(I);
    M= np.max(I);
    N = (M-m)/sigma * np.randomrand(I.shape[0], I.shape[1]);
    I2 = I2 + N;
    I2[ I2>M ] = M;
    I2[ I2<m ]=m;
    return I2, I-I2;
```

The image is blurred and noise is then added.

<figure><img src="image_194.png" /></figure>

1. `Cb = signal(convolve2d(C, psf, boundary='wrap', mode='same');`
   `Cbn, noise = addGaussianNoise(Cb);`
2. `plt.imshow(Cbn, cmap='gray')`
3. `plt.title('Motion blur + Gaussian noise on checkerboard')`
4. `plt.show()`

Optical Transfer Function: OTF

The OTF is the centered Fourier Transform of the PSF (Point Spread Function). The following function is used to get the OTF from the PSF:
<figure><img src="image_195.png" /></figure>

1 def psf2otf(psf, s):
    """

3 Get OTF (Optical Transfer Function) from PSF (Point Spread Function)
    OTF is basically the Fourier Transform of the PSF, centered

5 psf: PSF
s: shape of the result, zero-padding is used to center is Fourier Transform
"""

7 sh = psf.shape;
9 sh = np.array(sh);
s = np.array(s);

11 pad = s - sh;
h = np.pad(psf, ((0, pad[0]), (0, pad[1], mode='constant'));

13 shift = (int(pad[0]/2+1), int(pad[1]/2+1));
h_centered = npRoll(h, tuple(shift), axis=(0,1));

15 h_centered = npfeftshift(h_centered);

17 H = np.real(H);
return H;

6.7.2 Simple case: no noise

For the direct reconstruction, just ensure that there is no division by 0. The different
results are illustrated in Fig.6.3.

<figure><img src="image_196.png" /></figure>

H = psf2otf(psf, C.shape);

2 G = np.fft.fft2(Cb);
alpha = 0.001;

4 F = G / (H+alpha);
fr = np.real(np.fft. ifft2(F));

6 plt.imshow(fr);
    plt.title('Noiseless (only motion blur) direct reconstruction')

8 plt.imsave("cb_reconstruction.png", fr, cmap='gray');
    plt.show()

In the presence of noise, the reconstructed image is not perfect.

<figure><img src="image_197.png" /></figure>

1 H = psf2otf(psf, C.shape);
G = np.fft.fft2(Cbn);

3 alpha = 0.001;
F = G / (H+alpha);

5 fr = np.real(np.fft. ifft2(F));

The automatic evaluation of the noise parameter is done in the Wiener filter by:
<figure><img src="image_198.png" /></figure>

1. `SpecPuissNoise=mp.abs(np.fft.fft2(noise))**2;`
   PuissMoyNoise=mp.mean(SpecPuissNoise);

3. `SpecPuissImageOrig=mp.abs(np.fft.fft2(C))**2;`
   PuissMoyImageOrig=mp.mean(SpecPuissImageOrig);

5. `ratio =PuissMoyNoise/PuissMoyImageOrig;`
   `Hw = np.conjugate(H)/(np.abs(H)**2+ratio);`

7. `fr = np.fft. ifft2(Hw* np.fft.fft2(Cbn));`

### 6.7.3 Iterative filters

#### Van-Cittert iterative filter

The parameter (Jansson parameter) controls the precision of convergence, but a small value requires a high number of iterations (and thus a high computation time).

<figure><img src="image_199.png" /></figure>

1. `def vca(g, psf, n_iter, beta=.01):`
   """

3. `Van Cittert iterative filter`
   `g: noisy image`

5. `psf: point spread function`
   `n_iter: number of iterations`

7. `beta: Jansson parameter`
   """

9. `H = psf2otf(psf, g.shape);`

11. `# init iterations`

`fr = g.copy();`

13. `bar = progressbar PROGRESSBAR();`

`for iter in bar(range(n_iter)):`

15. `estimation = np.real(np. ifft2 (H* np.fft.fft2 (fr)));
   fr = fr + beta* (g - estimation);`

17. `return fr;`

#### Richardson-Lucy iterative filter

This algorithm is probably the most famous one in the domain. Remember that the symmetry in the space domain becomes the complex conjugate in the spectral domain (after Fourier transform).
<figure><img src="image_200.png" /></figure>

```python
def rla(g, psf, n_iter):
    """
    Richardson Lucy algorithm
g: noisy image
psf: point spread function
n_iter: number of iterations
    """
    H = psf2otf(psf, g.shape);
    # init iterations
    fr = g.copy();

    bar = progressbar PROGRESSBAR ()
    for iter in bar(range(n_iter)):
        estimation = np(np,fft, ifft2 (H * np,fft,fft2 (fr)));
        M = np.real(np,fft, ifft2 (np.conjugate(H) * np,fft,fft2 (g / estimation)));
        fr = fr * M;
        fr[fr<0] = 0;

    return fr;
```
Figure 6.3: Different deconvolution methods applied on the synthetic image 'checkerboard'.

<figure><img src="image_201.png" /><figcaption>(b) Reconstruction of the checkerboard without noise.</figcaption></figure>

<figure><img src="image_202.png" /><figcaption>(c) Reconstruction of the checkerboard with noise.</figcaption></figure>

<figure><img src="image_203.png" /><figcaption>(d) Wiener deconvolution.</figcaption></figure>

<figure><img src="image_204.png" /><figcaption>(e) Van-Cittert iterative deconvolution with 100 iterations, and Jansson parameter at 0.01.</figcaption></figure>

<figure><img src="image_205.png" /><figcaption>(f) Richardson-Lucy iterative deconvolution with 1000 iterations.</figcaption></figure>

<figure><img src="image_206.png" /><figcaption>(g) Landweber iterative deconvolution with 1000 iterations and Jansson parameter at 1.</figcaption></figure>

<figure><img src="image_207.png" /><figcaption>(h) Poisson Maximum A Posteriori iterative deconvolution with 1000 iterations.</figcaption></figure>

<figure><img src="image_208.png" /></figure>
### 6.7.4 Astronomy images

These real images come from the Hubble Space Telescope (HST, see credits), shown in Fig.6.4. The results for the Richardson-Lucy, Van-Cittert and Landweber algorithms are shown in Fig.6.5. The direct deconvolution does not give a correct result.

<figure><img src="image_209.png" /><figcaption>Figure 6.4: Original images and PSF, from the HST.</figcaption></figure>
<figure><img src="image_210.png" /><figcaption>Figure 6.5: Different deconvolution methods applied for the astronomy images.</figcaption></figure>
<figure><img src="image_211.png" /></figure>
<figure><img src="image_212.png" /></figure>

7 Shape From Focus

This tutorial introduces the basic shape from focus concepts. The objective is to reconstruct a focused image from a serie (generally a stack) of images with inhomogene focus (see Fig. 7.1).

<figure><img src="image_213.png" /><figcaption>Figure 7.1: Different images of the stack, from a corneal endothelium in optical microscopy (with 3D microscope).</figcaption></figure>

<figure><img src="image_214.png" /></figure>

Introduction to classical methods

An optic system has a limited depth of field. When observing non plane surfaces with a microscope, some parts of the observation may be blurred as well as some others may be correctly focused.
To overcome this problem and reconstruct an all-focused image as well as a surface (see Fig. 7.2), an algorithm will look at every pixel of the images in the stack and select the most focused one, by the way of a focus measure. This tutorial proposes to test some classical focus measures. You can have a look at [7] and [8] to see real applications.

<figure><img src="image_215.png" /><figcaption>Figure 7.2: Reconstruction of the surface and the texture.</figcaption></figure>

Practically, TIF files can handle stacks of images. Here is a way to open such a file:

<figure><img src="image_216.png" /></figure>

1. from skimage import io

2. I = io .imread('volume.tif');
3. I = I .astype('float');

<figure><img src="image_217.png" /></figure>

4. % store the images into stack
stack=zeros(info(1).Height, info(1).Width, num_images);
6. for k = 1:redimages
    stack(:,k)=imread stackfile, k);
8. % you can now do something with each image
end

In the following, each of the proposed methods will compute a focus measure
layer by layer, and will maximize this measure over the stack to find the most
focused layer.
### 7.1.1 Sum of Modified Laplacian

One of the first methods had been proposed by [25]. It is based on the second derivatives, specifically the Laplacian operator:

$$ \Delta^2 I = \frac{\partial^2 I}{\partial x^2} + \frac{\partial^2 I}{\partial y^2} $$

The problem with this operator is that the second derivatives in the x and y dimensions can have opposite signs. One way to overcome this problem is to introduce the modified Laplacian as follows:

$$ \Delta_M^2 I = \left| \frac{\partial^2 I}{\partial x^2} \right| + \left| \frac{\partial^2 I}{\partial y^2} \right| $$

The discrete approximation of the modified Laplacian is computed as:

$$ ML(x, y) = \left| 2I(x, y) - I(x - 1, y) - I(x + 1, y) \right| + \left| 2I(x, y) - I(x, y - 1) - I(x, y + 1) \right| \quad (7.1) $$

Then, the focus measure based on the modified Laplacian is:

$$ F_{ML} = \sum_{(x,y) \in \omega} ML(x, y) $$

### 7.1.2 Variance

The measure of focus based on the variance is today the mainly used method [15, 41]. It is based on the computation of the variance in a window $\omega$, with $N = \#\omega$ being the size of the window:

$$ F_v = \sum_{\omega} \left( I - \frac{1}{\#\omega} \underbrace{\sum_{\omega} I}_{\text{Mean of I}} \right)^2 $$

### 7.1.3 Tenengrad

In [21], we can find the definition of the tenengrad [42]. Let $S(x, y)$ be the norm of the Sobel gradient of image $I$.

$$ F_t(I) = \sum_{\omega} S^2 $$

Notice that the original definition requires a threshold value that requires heuristic choices, which is out of the topic of this tutorial.
### 7.1.4 Variance of Tenengrad

We define the variance of Tenengrad measure of focus by:

$$F_{vt}(I) = F_v(S)$$

## 7.2 Texture and surface reconstruction

<figure><img src="image_218.png" /></figure>

A set of images (Vickers indentation test [7], and a human corneal endothelium [8]) are proposed. For all the detailed methods:

* reconstruct the surface and the texture of the image (an image focused on its all field of view).

### 7.2.1 Open question

<figure><img src="image_219.png" /></figure>

Several methods have been implemented in order to perform the 3D surface/-texture reconstruction. Propose a numerical measure that could compare the results and measure the efficiency of these methods? What is the best method?
<figure><img src="image_220.png" /></figure>

## 7.3. Python correction

<figure><img src="image_221.png" /></figure>

<figure><img src="image_222.png" /></figure>

1. import numpy as np
import matplotlib .pyplot as plt

3. from mpl_toolkits.mplot3d import Axes3D
from skimage import io

5. from scipy import ndimage
from matplotlib import cm

### 7.3.1 Main function

This code is used as the main function to perform shape-from-focus reconstruction.

<figure><img src="image_223.png" /></figure>

volumes = ['cornee', 'vickers'];

2. for v in volumes:
   I = io .imread('volume_'+v+' .tif');

4. I = I .astype('float');
F = np .zeros(I.shape)

6. N = 11;

8. myfunctions = [variance, tenengrad, sml];

for f in myfunctions:
    for i, im in enumerate(I):

10. F[i] = f(im, N);

12. # Evaluates altitudes and textures

14. Z = np .argmax(F, axis =0);

16. Z = ndimage .minimum_filter(Z, size =5);

T = extractTexture(I, Z);

The extraction of the texture from the altitudes (indexes) in the stack of images is performed by the following function:
<figure><img src="image_224.png" /></figure>

```python
def extractTexture (I, Z):
    """
    Extract texture from stack of images I, where Z is the index (altitude)
    I: stack of images, of shape (n, X, Y)
    Z: index of SFF maximum (of shape (X,Y)), values are between 0 and n-1
    returns basically I[Z(i,j), i, j]] for all (i,j)
    """
m,n = I.shape[1:]
ii, jj = np.o grid [:m,:n]
T = I[Z, ii, jj];
return T;
```

### 7.3.2 Sum of Modified Laplacian

Results are illustrated in Fig.7.3.

<figure><img src="image_225.png" /></figure>

```python
def sml(I, N):
    """
    SFF measure, SUM of modified Laplacian
    I: image
    N: neighborhood size
    returns: SFF measure for each pixel
    """
h = np.array([[ -1,  2,  -1]]);
ML = np.abs( ndimage convolve(I,h) ) + np abs( ndimage convolve(I, np transpose(h)) ) ;
S = ndimage uniform_filter (ML, N);
return S;
```
<figure><img src="image_226.png" /><figcaption>Figure 7.3: Texture and altitude reconstruction with the SML method.</figcaption></figure>

### 7.3.3 Variance

The focus measure based on the variance i a really simple method that works in most cases, see Fig.7.4.

<figure><img src="image_227.png" /></figure>

1. def variance(I, N):
    """

2. I: image
3. N: neighborhood size

5. returns: SFF measure for each pixel, results is the same shape as I
    """

7. M = ndimage uniform_filter(I, N);

9. D2 = (I - M) ** 2;

10. V = ndimage uniform_filter(D2, N);

11. return V;

### 7.3.4 Tenengrad

The tenengrad method is based on a Sobel filter, see Fig.7.5.
<figure><img src="image_228.png" /><figcaption>Figure 7.4: Texture and altitude reconstruction with the variance method.</figcaption></figure>

<figure><img src="image_229.png" /></figure>

```python
def tenengrad(I, N):
    """

    SFF measure, Tenengrad method
    I: image
    N: neighborhood size
    returns: SFF measure for each pixel, results is the same shape as I
    """
    Sx = ndimage.sobel(I, axis=0);
    Sy = ndimage.sobel(I, axis=1);
    S = np.hypot(Sx, Sy);
    T = ndimage uniform_filter (S, N);
    return T;
```

### 7.3.5 Variance of Tenengrad

The variance of Tenengrad is an improvement of the Tenengrad method, see Fig.7.6.
<figure><img src="image_230.png" /></figure>

def varianceTenengrad(I, N):

2

"
""

4

I: image
N: neighborhood size

6

returns: SFF measure for each pixel
"
""

8

Sx = ndimage.sobel(I, axis=0);

Sy = ndimage.sobel(I, axis=1);

10

S = np.hypot(Sx, Sy);

vt = variance(S, N);

12

return vt;

<figure><img src="image_231.png" /><figcaption>Figure 7.5: Texture and altitude reconstruction with the SML method.</figcaption></figure>
<figure><img src="image_232.png" /><figcaption>Figure 7.6: Texture and altitude reconstruction with the variance of Tenengrad method.</figcaption></figure>
<figure><img src="image_233.png" /></figure>

8

Logarithmic Image Processing
(LIP)

This tutorial aims to study the LIP model, which is a vector space of gray-level images, consistent with the physical laws of Weber and Fechner as well with the visual perception laws. More particularly, the LIP dynamic expansion, used for image enhancement, will be implemented as well as the Sobel LIP filter for edge detection.

The different processes will be applied on a mammographic image.

<figure><img src="image_234.png" /><figcaption>Figure 8.1: Breast</figcaption></figure>

8.1

**Introduction**

The LIP model (Logarithmic Image processing) has been introduced in the mid 1980s. It defines a mathematical framework for image processing in a bounded interval. It is mathematically rigorous and physically justified. For more informations, refer to [28, 16, 29].

An image is represented by its gray-tone function *f*, defined on a spatial domain *D* ⊂ ℝ², and with values into [0, *M*[, where *M* > 0. For all gray-tone functions *S* defined on *D*, a vector space is defined by the operations of addition △ and
multiplication $\triangle$, and by extension with the operations of subtraction $\triangle$ and negation:

$$ \forall f, g \in S, f \triangle g = f + g - \frac{fg}{M} \qquad (8.1) $$

$$ \forall f \in S, \forall \lambda \in \mathbb{R}, \lambda \triangle f = M - M\left(1 - \frac{f}{M}\right)^{\lambda} \qquad (8.2) $$

$$ \forall f \in S, \triangle f = \frac{-Mf}{M - f} \qquad (8.3) $$

$$ \forall f, g \in S, f \triangle g = M \frac{f - g}{M - g} \qquad (8.4) $$

The vector space $S$ of gray-tone functions is algebraically and topologically isomorph to the classical vector space, defined by the isomorphism $\varphi$:

$$ \forall f \in S, \varphi(f) = -M\ln\left(1 - \frac{f}{M}\right) \qquad (8.5) $$

The inverse isomorphism $\varphi^{-1}$ is then defined by:

$$ f = \varphi^{-1}(\varphi(f)) = M\left(1 - \exp\left(-\frac{\varphi(f)}{M}\right)\right) \qquad (8.6) $$

This fundamental isomorphism allows the definition of the following operations on gray-tone functions:

* Dot product: $\forall f, g \in S, \langle f, g \rangle_{\triangle} = \int \varphi(f)\varphi(g)$

* Euclidean norm: $\forall f \in S, \|f\|_{\triangle} = |\varphi(f)|_{\mathbb{R}}$, with $|\cdot|_{\mathbb{R}}$ is the classical absolute value.

## 8.2 Elementary LIP operations

Classically, gray-tones are coded with 8 bits, and thus one can choose $M = 256$. Be careful that there might exist sampling issues, and it might be required to considere images in the LIP space with a double precision.

<figure><img src="image_235.png" /></figure>

1. Implement the elementary LIP operations $\triangle$, $\triangle$ and $\triangle$.

2. Test these operators on the image 'breast'.
## 8.3 LIP dynamic expansion

The LIP model enables to define an image transformation that expands, in an optimal way, the overall dynamic range (gray-tones) while preserving a physical or visual sense. Let a graytone function denoted $f$ be defined on the spatial support $D \subset \mathbb{R}^2$. Its upper and lower bounds are denoted $f_u = \sup_{x \in D} f(x)$ and $f_l = \inf_{x \in D} f(x)$, with $0 < f_l < f_u < M_0$.

The dynamic range of $f$ on $D$ is defined by $R(f) = f_u - f_l$. The LIP scalar multiplication of $f$ by a real number $\lambda > 0$ yields to the following dynamic :

$$R(\lambda \triangleq f) = \lambda \triangleq f_u - \lambda \triangleq f_l \quad (8.7)$$

It has been shown [17] that there exists an optimal value $\lambda_0(f)$ that maximizes the dynamic range, i.e.:

$$R(\lambda_0(f) \triangleq f) = \max_{\lambda > 0} R(\lambda \triangleq f). \quad (8.8)$$

<figure><img src="image_236.png" /></figure>

1. Determine the explicit (analytic) expression of the parameter $\lambda_0(f)$.

2. Calculate the value of the parameter $\lambda_0(f)$ for the image 'breast'.

3. Enhance the image and compare the result with the histogram equalization.

<figure><img src="image_237.png" /></figure>

In python, one can use the scikit-image module, and the skimage.exposure equalize_histogram function.

## 8.4 LIP edge detection

<figure><img src="image_238.png" /></figure>

1. Implement the Sobel filter within the LIP framework.

2. Compare the results with the classical Sobel filter.
<figure><img src="image_239.png" /></figure>

## 8.5. Python correction

<figure><img src="image_240.png" /></figure>

### 8.5.1 Elementary operations

The most important function to code is the graytone transformation function. It considers M as the maximum value (the absolute white). This code means that the absolute white cannot be reached, and thus $f \in [0; M[$. After discretizing the gray values, $F \in [0; M - 1]$.

<figure><img src="image_241.png" /></figure>

### file LIP.py

2 LIP simple module

4 USAGE: import LIP

6 """

import numpy as np

8

def graytone(F, M):

10 \# graytone function transform

\# M: maximal value

12 \# F: image function

f = M-np.finfo(float).eps-F;

14 return f;

<figure><img src="image_242.png" /></figure>

def phi(f, M):

2 \# LIP isomorphism
\# f: graytone function

4 \# M: maximal value
l = -M * np.log(-f/M + 1);

6 return l

<figure><img src="image_243.png" /></figure>

def invphi(l, M):

2 \# inverse isomorphism
f = M*(1-np.exp(-l/M));

4 return f
<figure><img src="image_244.png" /></figure>

def plusLIP(a, b, M):

2 # LIP addition
z = a+b-(a*b)/M;

4 return z;

<figure><img src="image_245.png" /></figure>

def timesLIP(alpha, x, M):

2 # LIP multiplication by a real
z = M-M*(np.ones(x.shape)-x/M)**alpha;

4 return z;

8.5.2 LIP dynamic expansion

The optimal value for dynamic expansion is given by $\lambda_0$:

$$\lambda_0(f) = \arg\max_{\lambda} \{\max(\lambda \triangleq f) - \min(\lambda \triangleq f)\}$$

Let $A(\lambda) = \max(\lambda \triangleq f) - \min(\lambda \triangleq f) = \lambda \triangleq \max(f) - \lambda \triangleq \min(f)$ and $B = \ln(1 - \min(f)/M)$ et $C = \ln(1 - \max(f)/M)$.

$$\begin{align*}
A'(\lambda) &= 0 \leftarrow [(M - M \exp(\lambda C)) - (M - M \exp(\lambda B))]' = 0 \\
&\leftarrow [\exp(\lambda B) - \exp(\lambda C)]' = 0 \\
&\leftarrow B \exp(\lambda B) - C \exp(\lambda C) = 0 \\
&\leftarrow \ln(B) + \lambda B = \ln(C) - \lambda C \\
&\leftarrow \lambda = \frac{\ln(C) - \ln(B)}{B - C} \\
&\leftarrow \lambda = \frac{\ln(C/B)}{B - C}
\end{align*}$$

Thus, yielding to:

$$\lambda_0(f) = \frac{\ln\left(\frac{\ln(1 - \max(f)/M)}{\ln(1 - \min(f)/M)}\right)}{\ln\left(\frac{M - \min(f)}{M - \max(f)}\right)}$$

This is coded in python by:
<figure><img src="image_246.png" /></figure>

def computeLambda(f, M):

2

# compute optimal value for dynamic expansion

B = np.log(1 - f.min() / M);

4

C = np.log(1 - f.max() / M);

l = np.log(C/B)/(B-C);

6

return l;

### 8.5.3 Complete code

This code uses a transform called histogram equalization. This version proposes our own code for the histogram equalization (see also 2.

<figure><img src="image_247.png" /></figure>

import LIP

2

import numpy as np

from scipy import misc

4

import matplotlib .pyplot as plt

import skimage

6

### Histogram equalization

###

8

def histeq(im,nbrbins=256):

### get image histogram

10

imhist, bins = np.histogram(im. flatten() ,nbr bins, normed = True)

cdf = imhist .cumsum() #cumulative distribution function

12

cdf = 255 * cdf / cdf[ -1] #normalize

### use linear interpolation of cdf to find new pixel values

14

im2 = np . interp(im. flatten() ,bins [: -1], cdf)

16

return im2 . reshape(im.shape), cdf

It is compared to the LIP dynamic enhancement. The results are illustrated in Fig. 8.2.
<figure><img src="image_248.png" /><figcaption>Figure 8.2: Comparison of the LIP enhancement with the classical histogram equalization method.</figcaption></figure>
<figure><img src="image_249.png" /></figure>

1 # contours detection

# Sobel detection

3 def Sobel(I):
    sobx = cv2.Sobel(I, -1, 1, 0);
    soby = cv2.Sobel(I, -1, 0, 1);
    sob = np.sqrt(sobx**2 + soby**2);
    return sob;

9 # go into LIP space
tonelip = LIP.phi(tone, M);

# apply Sobel filter

13 sobellip = LIP.grytone(LIP.invphi(Sobel(tonelip), M), M);
imageio imwrite("sobellip.png", sobellip);

15 plt.figure();
    pltsubplot(1,2,1); plt.imshow(sobellip, cmap= plt.cm .gray);

17 plt. title ( 'LIP Sobel edge detection ')

19 # apply Sobel filter in the classic space
sobel = 255-Sobel(B);

21 plt.subplot(1,2,2); plt. imshow(sobel, cmap= plt.cm .gray);
    plt. title ( 'Sobel edge detection ')

23 imageio.savwrite("sobel.png", sobel);

<figure><img src="image_250.png" /><figcaption>Figure 8.3: Sobel edge detection</figcaption></figure>
<figure><img src="image_251.png" /></figure>
<figure><img src="image_252.png" /></figure>

9

Color Logarithmic Image Process-
ing (CoLIP)

The objective of this tutorial is to have an overview of the Color Logarithmic Image Processing (CoLIP) framework [11, 13, 8, 10]. The CoLIP is a mathematical framework for the representation and processing of color images. It is psychophysically well justified since it is consistent with several human visual perception laws and characteristics. It allows to consider color images as vectors in an abstract linear space, contrary to the classical color spaces (e.g., RGB and $L^*a^*b^*$). The purpose of this tutorial is to present the mathematical fundamentals of the CoLIP by manipulating the color matching functions and the chromaticity diagram.

9.1 Definitions

This diagram summarizes the transitions between the classical RGB space and the CoLIP space. The transfer matrices are presented in the following.

<figure><img src="image_253.png" /></figure>

9.1.1 RGB to XYZ

Here, we choose to keep the conventions of the CIE 1931, i.e. a 2° 1931 observer, with
the equal energy illuminant E and the sRGB space. Depending on the illuminant,
there exist many conversion matrices.

<figure><img src="image_254.png" /></figure>

Use the module skimage.color.
### 9.1.2 XYZ to LMS

The matrix $M_{HPE}$ (Hunt, Pointer, Estevez) is defined by:

$$ \begin{pmatrix} L \\ M \\ S \end{pmatrix} = M_{HPE} \times \begin{pmatrix} X \\ Y \\ Z \end{pmatrix}, \qquad (9.1) $$

with

$$ M_{HPE} = \begin{pmatrix} 0.38971 & 0.68898 & -0.07868 \\ -0.22981 & 1.18340 & 0.04641 \\ 0.00000 & 0.00000 & 1.00000 \end{pmatrix}. \qquad (9.2) $$

### 9.1.3 LMS to Lms

$$ \forall c \in \{l, m, s\}, C \in \{L, M, S\}, c = M_0 \left(1 - \frac{C}{C_0}\right), \qquad (9.3) $$

with $C_0$ is the maximal transmitted intensity value. $M_0$ is arbitrarily chosen at normalized value 100. Notice that $C \in [0; C_0]$ and $c \in [0; M_0[$.

### 9.1.4 CoLIP homeomorphism

$$ \forall f \in S, \varphi(f) = -M_0 \ln \left( 1 - \frac{f}{M_0} \right). $$

The logarithmic response of the cones, as in the LIP theory, is modeled through the isomorphism $\varphi$:

$$ \text{for } c \in \{l, m, s\}, \tilde{c} = \varphi(c) = -M_0 \ln \left(1 - \frac{c}{M_0}\right), \qquad (9.4) $$

where $(\tilde{l}, \tilde{m}, \tilde{s})$ are called the logarithmic chromatic tones.

### 9.1.5 Opponent process

$$ \begin{pmatrix} \tilde{a} \\ \tilde{r}g \\ \tilde{y}b \end{pmatrix} = P \times \begin{pmatrix} \tilde{l} \\ \tilde{m} \\ \tilde{s} \end{pmatrix}, \qquad (9.5) $$

with

$$ P = \begin{pmatrix} 40/61 & 20/61 & 1/61 \\ 1 & -12/11 & 1/11 \\ 1/9 & 1/9 & -2/9 \end{pmatrix}. \qquad (9.6) $$
9.1.6 Bounded vector space

Three channels $\hat{f} = (\hat{a}, \hat{r}g, \hat{y}b)$ are defined by:

$$ \hat{a} = a \qquad (9.7) $$

$$ \hat{r}g = \begin{cases} rg & \text{if } rg \ge 0 \\ - \Delta rg & \text{if } rg < 0 \end{cases} \qquad (9.8) $$

$$ \hat{y}b = \begin{cases} yb & \text{if } yb \ge 0 \\ - \Delta yb & \text{if } yb < 0 \end{cases} \qquad (9.9) $$

9.2 Applications

9.2.1 CMF

<figure><img src="image_255.png" /></figure>

Represent the classical chromaticity by using the color matching functions.

<figure><img src="image_256.png" /></figure>

The color matching functions in the different spaces, as well as the wavelength and the purple line, are provided in the matrix 'cmf.npz'.
Load it with numpy.load.
9.2.2 CoLIP chromaticity diagram

Represent the chromaticity diagram in the plane $(\hat{r}g, \hat{y}b)$ as well as the
Maxwell triangle of the RGB colors, as in Fig.9.1.

<figure><img src="image_257.png" /><figcaption>Figure 9.1: Chromaticity diagram in the plane $(\hat{r}g, \hat{y}b)$.</figcaption></figure>
<figure><img src="image_258.png" /></figure>

## 9.3. Python correction

<figure><img src="image_259.png" /></figure>

<figure><img src="image_260.png" /></figure>

import numpy as np

2 import matplotlib .pyplot as plt

import skimage.color as color

First of all, the maximal value $M_0$ is arbitrarily fixed at 100.

<figure><img src="image_261.png" /></figure>

1 def getColipM0():

# return M0 value

3 return 100;

### 9.3.1 LMS tones

This is the difficult part of LIP and CoLIP. Be careful with the use of the function `eps` that returns the precision at a given double value.

<figure><img src="image_262.png" /></figure>

def lmstone(LMS):

2 M0= getColipM0();
return (M0-np.finfo( float ) .eps)*(1-LMS/M0);

### 9.3.2 Isomorphism

The isomorphism is the conversion into/back from the logarithmic space.

<figure><img src="image_263.png" /></figure>

1 def phi(f, M):
    # LIP isomorphism

3 # f: graytone function

# M: maximal value

5 l = -M * np.log(-f/M+1);
return l
<figure><img src="image_264.png" /></figure>

```python
def invphi(l, M):
    # inverse isomorphism
    f = M*(1-np.exp(-1/M));
    return f
```

### 9.3.3 XYZ to LMS

This conversion uses the HPE matrix (Hunt, Pointer, Estevez).

<figure><img src="image_265.png" /></figure>

```python
def XYZ2LMS(XYZ):
    """
    conversion function
    XYZ: data in XYZ space
    """
    U(np.array([[0.38971, 0.68898, -0.07869],
                 [-0.22981,1.18340,  0.04641],
                 [0,         0,         1    ]]));
    m, n = XYZ [:, :0]. shape;
    print (m,n)
    XYZ = XYZ.reshape((m*n, 3)). transpose ();
    print (XYZ.shape)
    LMS = np.matmul(U, XYZ);

    return LMS transpose() . reshape((m, n, 3));
```

### 9.3.4 Conversions into Colip space

The conversion into $(\hat{a}, \hat{r}g, \hat{y}b)$ goes first into $(\tilde{a}, \tilde{r}g, \tilde{y}b)$. It consists on a conversion matrix followed by the absolute value and the operation to get symmetric channels in the colip space.

<figure><img src="image_266.png" /></figure>

```python
def LMS2ARGYBhat(LMS):
    """
    ARGYBtilde = LMS2ARGYBtilde(LMS);
    """
    return ARGYBtilde2ARGYBhat(ARGYBtilde);
```
<figure><img src="image_267.png" /></figure>

```python
def LMS2ARGYBtilde(LMS):
    """
    conversion function
    LMS: data in LMS space
    """
    M0 = getColipM0();
    m, n, p = LMS.shape;
    LMStone = lmstone(LMS);
    LMStilde = phi(LMStone, M0);
    P = [[40/61,20/61,1/61], [1, -12/11,1/11], [1/9,1/9, -2/9]];

    LMStilde = LMStilde. reshape((m*n, 3)). transpose ();
    ARGYBtilde = np.matmul(P, LMStilde). transpose ();
    return ARGYBtilde. reshape((m, n, 3));
```

<figure><img src="image_268.png" /></figure>

```python
def ARGYBtilde2ARGYBhat(ARGYBtilde):
    """
    conversion into ARGYB tilde space
    ARGYBtilde: data in ARGYB tilde space
    """
    M0 = getColipM0();
    ARGYBhat = np zeros(ARGYBtilde.shape);
    ARGYBhat = invphi(ARGYBtilde, M0);
    for c in (1,2):
        tmp = np abs(ARGYBtilde[::,c]);
        ARGYBhat[::,c] = np sign(ARGYBtilde[::,c]) * invphi(tmp, M0);
    return ARGYBhat;
```

### 9.3.5 CMF

The color matching functions are provided for convenience. There exist many resources on the internet where they can be found. The classical diagram in the xy space (the horseshoe) is shown in Fig.9.2.
<figure><img src="image_269.png" /></figure>

with np.load('cmf.npz') as data:

2
cmap = data['cmap'];
pourpresLMS = data['pourpresLMS'];

4
SpecXYZ = data['SpecXYZ'];
SpecLMS = data['SpecLMS'];

6
SpecRGB = data['SpecRGB'];
l = data['l'];

8
# display classical CMF in xy

10
xn = SpecXYZ [::,0]/ np.sum(SpecXYZ, axis=2);
yn = SpecXYZ [::,1]/ np.sum(SpecXYZ, axis=2);

12
zn = 1-xn-yn;

plt . scatter (xn, yn, c=cmap);

<figure><img src="image_270.png" /><figcaption>Figure 9.2: Color matching functions in the xy space.</figcaption></figure>

To display the CMF and the cube of all RGB colors in the (rhat, yhat) space, the
following code is used: The result is shown in Fig.9.3.

<figure><img src="image_271.png" /></figure>

1 #Color matching functions in colip space

ARGYB_hat = LMS2ARGYBhat(SpecLMS);

3 plt . figure ();

plt . scatter (ARGYB_hat [::,0,1], ARGYB_hat [::,0,2], c =cmap);

5 # purple line
purpleARGYB_hat = LMS2ARGYBhat(pourpresLMS);

7

    plt . scatter (purpleARGYB_hat [::,0,1], purpleARGYB_hat [::,0,2], c='black');
The RGB cube is first generated and then converted into the colip space.

<figure><img src="image_272.png" /></figure>

1. `# number of color in each direction`
   `N=10;`

3. `cols = np.linspace(0, 255, num=N);`
   `R, G, B = np.meshgrid colsin, colsin, colsin);`

5. `# reshape the cube for manipulation`

7. `R = R.reshape((R.size, 1));`

G = G.reshape((G.size, 1));

9. `B = B.reshape((B.size, 1));`

NN = R.size;

11. `colRGB = np.concatenates((R, G, B), axis=1);`

cubeRGB = colRGB.reshape((NN, 1, 3));

13. `colRGB = colRGB / 255;`

15. `# conversions`

cubeXYZ = colorIFT2XYZ cubeRGB;

17. `cubeLMS = XYZ2LMS cubeXYZ);`

cubeARGYBhat = LMS2ARGYBhat (cubeLMS);

19. `plot scatter (cubeARGYBhat [1:0,1], cubeARGYBhat [1:0,2], c=colRGB);`

The color matching functions in the colip space are represented in Fig.9.3. A lot could be said about this diagram. For examples, one can notice the approximations in the high (red) wavelengths: this comes from the numerical approximations and the lack of data of the original CMFs. Moreover, the purple line (represented here in black) is not a line in this diagram. This raises the question of the choice of this line in (xy), which was made only by simplicity.

<figure><img src="image_273.png" /><figcaption>Figure 9.3: Color matching functions in the $(\hat{r}g, \hat{y}b)$ space.</figcaption></figure>
<figure><img src="image_274.png" /></figure>
<figure><img src="image_275.png" /></figure>

## 10 GANIP

This tutorial aims to test the elementary operators of the GANIP framework. The General Adaptive Neighborhood Image Processing (GANIP) is a mathematical framework for adaptive processing and analysis of gray-tone and color images. An intensity image is represented with a set of local neighborhoods defined for each pixel of the image to be studied. Those so-called General Adaptive Neighborhoods (GANs) are simultaneously adaptive with the spatial structures, the analyzing scales and the physical settings of the image to be addressed and/or the human visual system.

The following figure illustrates the GANs of two points on an image of retinal vessels.

<figure><img src="image_276.png" /></figure>

The GANs are then used as adaptive operational windows for local image transformations (morphological filters, rank/order filters...) and for local image analysis (local descriptors, distance maps...).

The different processes will be applied on the following gray-tone images:

### 10.1 GAN

Let $f$ be a gray-tone image. The GAN of a point $x$ using the luminance criterion and the homogeneity tolerance $m$ within the CLIP framework is defined as:

$$V_m^f(x) = C_{\{y; |f(y)-f(x)| \le m\}}(x) \qquad (10.1)$$

where $C_X(x)$ denotes the connected component of $X$ holding $x$.
<figure><img src="image_277.png" /><figcaption>Figure 10.1: Two example images for testing the GAN framework.</figcaption></figure>

<figure><img src="image_278.png" /></figure>

1. Load the image 'lena' and compute the GAN of a selected point in the image.

2. Look at the influence of the homogeneity tolerance.

3. Compute the GAN of different points and comment.

## 10.2 GAN Choquet filtering

The Choquet filters generalize the rank-order filters. In this exercise, we are going to implement some GAN Choquet filters such as the GAN mean operator. For each point $x$ of the image, the mean value of all the intensities of the points inside the GAN of $x$ is computed. So, a first algorithm consists in making a loop on the image points for computing the different GAN and the mean intensity values, but it is time consuming. Nevertheless, by using some properties of the GAN (particularly the one giving that iso-valued points can have exactly the same GAN), it is possible to create a second algorithm by making a loop on the gray-tone range:
**Data:** original (8-bit) image *f*, homogeneity tolerance *m*

**Result:** GAN mean filtered image *g*

for *s* = 0 to 255 do

*seed* = points *x* with intensity *f*(*x*) = *s*;

*thresh* = points *y* with intensity satisfying *s* − *m* ≤ *f*(*y*) ≤ *s* + *m*;

*threshGAN* = connected components of *thresh* holding *seed*;

foreach *label* of *threshGAN* do

*currentLabel* = current connected component of *threshGAN*;

*meanValue* = intensity mean value of the image points inside
*currentLabel*;

set *g*(*x*) = *meanValue* to the points *x* of *seed* belonging to
*currentLabel*;

end

end

<figure><img src="image_279.png" /></figure>

1. Implement the proposed algorithm.

2. Test this operator on the 'lena' image with different homogeneity tolerances.

3. Compare the result with a classical mean filtering.

<figure><img src="image_280.png" /></figure>

See spipy.ndimage.filters.
## 10.3 GAN morphological filtering

As previously explained, it is also possible to compute the GAN dilation and GAN erosion by making a loop on the gray-tone range. The algorithm for the GAN dilation is:

**Data:** original (8-bit) image $f$, homogeneity tolerance $m$

**Result:** GAN dilated image $g$

$set g(x) = 0$ for all points $x$;

for $s = 0$ to 255 do

$seed$ = points $x$ with intensity $f(x) = s$;

$thresh$ = points $y$ with intensity satisfying $s - m \le f(y) \le s + m$;

$threshGAN$ = connected components of $thresh$ holding $seed$;

foreach label of $threshGAN$ do

$currentLabel$ = current connected component of $threshGAN$;

$maxValue$ = intensity maximum value of the image points inside

$currentLabel$;

$set g(x) = max(g(x),MaxValue)$ to the points $x$ belonging to

$currentLabel$;

end

end

**Algorithm 2:** GAN dilation

The algorithm for the GAN erosion is similar.

<figure><img src="image_281.png" /></figure>

1. Implement the GAN dilation and GAN erosion.

2. Test this operator on the 'lena' image with different homogeneity tolerances.

3. Compare the results with the classical morphological dilation and erosion.

4. Compute, test, compare and check the properties (idempotence, extensivity/anti-extensivity) of the GAN opening and GAN closing.

<figure><img src="image_282.png" /></figure>

Compare with the functions dilation and erosion from
skimage.morphology.
<figure><img src="image_283.png" /></figure>

## 10.4. Python correction

<figure><img src="image_284.png" /></figure>

<figure><img src="image_285.png" /></figure>

import numpy as np

2 import skimage.morphology

import skimage.io

4 import skimage.measure

import matplotlib .pyplot as plt

6 import scpy.ndimage. filters

### 10.4.1 GAN

In order to show the General Adaptive Neighborhood (GAN) of one image point, the following function is implemented. Note that the function is given within the Classical Linear Image Processing (CLIP) framework, i.e. with the usual addition, subtraction and scalar multiplication. But the function could be generalized to other models such as the Logarithmic Image Processing (LIP).

<figure><img src="image_286.png" /></figure>

def GAN(A, p, m):

2

### GAN construction

4

A: grayscale image

p: point

6

m: tolerance

### GAN construction

8

RES = np. zeros(A.shape);

RES[p[0], p[1]] = 1;

10

s = A[p[0], p[1]];

thresh = np. logical_and(A>= s-m, A<=s+m);

12

SE = np. ones((3,3));

RES = skimage.morphology.reconstruction(RES, thresh, selem=SE);

14

return RES;

In this way, we can visualize the GAN of any point of the 'Lena' image:

<figure><img src="image_287.png" /></figure>

I = skimage.io .imread(' lena.png ');

2

p=[100,200];

for m in [5,50,75]:

4

RES = GAN(I, p, m);

plt .imshow(RES, cmap='gray');

6

plt.show()
### 10.4.2 GAN Choquet filtering

In order to compute the GAN mean filtering, the basic idea is to make a loop on the image points, compute the GAN of the current point and then calculate the mean intensity of the points within the GAN. But this is time computing in the sense that two points with the same intensity can have exactly the same GAN (when one point is included in the GAN of the second point with same intensity). Therefore, a more effective way is to make a loop on the gray levels of the image. The GAN mean filtering is then implemented as:

<figure><img src="image_288.png" /><figcaption>Figure 10.2: GAN of a specific point of the 'Lena' image using different homogeneity tolerances *m* within the CLIP framework.</figcaption></figure>
<figure><img src="image_289.png" /></figure>

```python
def GANmean(A, m):
    """
    GAN mean filter
    Apply mean value of the adaptive neighborhood to result
    A: original image
    m: homogeneity tolerance
    return: GAN mean filter
    """
    RES = np.zeros(A.shape).astype('float');
    SE = np.ones((3,3));
    for s in np.arange(256):
        thresh = np logical and (A>=s-m, A<=s+m);
        seed = A==s;
        thresh = skimage.morphology.reconstruction(seed, thresh, selem=SE);
        L,n = skimage measure.label(thresh, connectivity =2, return_num=True);
        for l in np.arange(1, n+1):
            currentLabel = L==l;
            values = A(currentLabel);
            meanValues = np.mean(values);
            result = meanValues * np logical and (seed, currentLabel).astype('float
                                               c) );
            RES = RES + result;
    return RES;
```

Looking at this function, for the first loop on the gray levels `s`, we detect the GANs of the pixels with gray level `s`. So 'thresh' contains all these connected components (GANs), where some points with identical gray levels `s` can have the same GAN. 'label' enables the different GANs to be labeled. For the second loop on these GANs, we extract for each GAN its gray levels which are stored in 'values'. Afterwards, we take the mean value 'meanValue' that is stored as the resulting gray level for pixels inside the GAN with the same gray level `s`.

We can compare this GAN filtering with the classical one using a fixed operational window.

<figure><img src="image_290.png" /></figure>

```python
B scipy.ndimage. filters uniform_filter(I, size =5);
```

```python
2 m=30;
C = GANmean(I, m);
```

You can see the blurring effect caused by the classical filtering, contrary to the adaptive GAN filtering where the transitions are much more preserved while smoothing the image.
<figure><img src="image_291.png" /><figcaption>Figure 10.3: Classical ($r = 5$) vs. GAN ($m = 30$) mean filtering of the 'Lena' image.</figcaption></figure>

GAN morphological filtering

The following codes enable the GAN dilation and erosion to be computed. As previously mentioned, it is based on a loop on the gray levels and not on the image points. Note that the GANs are computed on the criterion image (it is important for satisfying the good properties of opening and closing, such as idempotence). Please notice that this function as well as the following will have the same parameters:

<figure><img src="image_292.png" /></figure>

1

### A: original image

3

### criterion : criterion image

5

### m: tolerance

<figure><img src="image_293.png" /></figure>

1

### def GANdilation(A, criterion, m):

3

### RES = np.zeros(A.shape).astype('float');

5

### SE = np.ones((3,3));

7

### for s in np.arange(256):
    thresh = nplogical and (criterion >=s-m, criterion <=s+m);
    seed = criterion ==s;
    thresh = skimage.morphology.reconstruction(seed, thresh, selem=SE);
    L,n = skimage.measure.label(thresh, connectivity =2, return_num=True);

9

### for l in np浓ate(1, n+1):
    currentLabel = L==l;

11

### values = A(currentLabel);
    values. sort();

13

### result = currentLabel.astype('float') * values[-1];
    RES = npmaximum(RES, result);

15

### return RES;
<figure><img src="image_294.png" /></figure>

```python
def GANerosion(A, criterion, m):
    RES = 255 * np.ones(A.shape).astype('float');
    SE = np.ones((3,3));

    for s in np.arange(256):
        thresh = np logical_and (criterion >=s-m, criterion <=s+m);
        seed = criterion ==s;
        thresh = skimage.morphology.reconstruction(seed, thresh, selem=SE);
        L,n = skimage.measure.label(thresh, connectivity =2, return_num=True);
        for l in np浓ate(1, n+1):
            currentLabel = L==l;
            values = A(currentLabel);
            values.sort();
            result = currentLabel.astype('float') * values[0] + 255* np.logical_not
                  (currentLabel.astype('float'));
            RES = np minimum(RES, result);
        return RES;
```

The following script shows the difference between classical and adaptive morphology.

<figure><img src="image_295.png" /></figure>

```python
Cdil = GANdilation(I, m);

2 Cero = GANerosion(I, m);
    se =skimage.morphology.disk(2);

4 Bdil = skimage.morphology.dilation(I, selem=se);
    Bero = skimage.morphology.erosion(I, selem=se);
```

As previously mentioned with the Choquet filtering, you can see the blurring effect caused by the classical filtering, contrary to the adaptive GAN filtering where the transitions are much more preserved while smoothing the image.

Now, we can compute the GAN opening and closing as combined operators of dilation and erosion.

<figure><img src="image_296.png" /></figure>

```python
def GANopening(A, criterion, m):
    """
    """
    temp = GANerosion(A, criterion, m);
    RES = GANdilation(temp, criterion, m);
    return RES;
```
<figure><img src="image_297.png" /><figcaption>Figure 10.4: Classical ($r = 2$) vs. GAN ($m = 30$) morphological filtering of the 'Lena' image.</figcaption></figure>
<figure><img src="image_298.png" /></figure>

def GANclosing(A, criterion, m):

2

### GANClosing

### GANclosing

4

temp = GANdilation(A, criterion, m);

RES = GANerosion(temp, criterion, m);

6

return RES;

It is very important to use the same criterion *h* when combining these two operators of GAN dilation and erosion. It means that we compute the GANs at the beginning and thereafter we use these same GANs for computing dilation and erosion. It enables the idempotence, extensivity/anti-extensivity of the GAN opening and GAN closing to be satisfied.

<figure><img src="image_299.png" /></figure>

Copen=GANopening(I, I, m);

2

Cclose=GANclosing(I, I, m);

4

### # classical open close

Bopen = skimage.morphology.opening(I, selem=se);

6

Bclose= skimage.morphology.closing(I, selem=se);

The result of such morphological filters is presented in Fig.10.5.

<figure><img src="image_300.png" /><figcaption>Figure 10.5: GAN ($m = 30$) opening and closing of the 'Lena' image.</figcaption></figure>

We can check some properties of these morphological filters such as extensivity/anti-extensivity and idempotence:
<figure><img src="image_301.png" /></figure>

```python
Copen2 = GANopening(Copen, I, m);

Cclose2= GANclosing(Cclose, I, m);
if np.max(np.abs(Copen2-Copen))<=0:
    print("OK: idempotence in opening");
else:
    print("error in opening");

if np.max(np.abs(Cclose2-Cclose))<=0:
    print("OK: idempotence in closing");
else:
    print("error in opening");
```

<figure><img src="image_302.png" /></figure>

```python
#%% check extensivity and anti-extensivity
if np. all (Copen<=I):
    print("OK: extensivity ");
else:
    print("error in extensivity ");

if np. all (Cclose>=I):
    print("OK: anti-extensivity ");
else:
    print("error in anti-extensivity ");
```

Result verifies these properties.

<figure><img src="image_303.png" /></figure>

```python
OK: idempotence in opening

OK: idempotence in closing

OK: extensivity

OK: anti-extensivity
```
<figure><img src="image_304.png" /></figure>

# 11 Image Filtering using PDEs

This tutorial aims to process images with the help of partial differential equations.

The different processes will be applied on the following MR image Fig. 11.1.

<figure><img src="image_305.png" /><figcaption>Figure 11.1: MRI image of a human brain.</figcaption></figure>

## Notations

The different operators used here are:

$$ \vec{A} = \begin{pmatrix} A_x \\ A_y \end{pmatrix} \qquad (11.1) $$

$$ \mathrm{div} \vec{A} = \frac{\partial A_x}{\partial x} + \frac{\partial A_y}{\partial y} \qquad (11.2) $$

$$ \vec{\mathrm{grad}} u(x, y) = \nabla u = \begin{pmatrix} \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial y} \end{pmatrix} \qquad (11.3) $$

The Laplacian operator is denoted $\Delta u = \nabla^2 u = \mathrm{div} \mathrm{grad} u$. A numerical approximation can be used, but is not recommended in this tutorial.
## 11.1 Linear diffusion

The heat equation is defined as follows:

$$ 
\left\{ 
\begin{aligned}
\frac{\partial u}{\partial t}(x, y, t) &= \operatorname{div}(\nabla u(x, y, t)) \\
&= \frac{\partial^2 u}{\partial x^2}(x, y, t) + \frac{\partial^2 u}{\partial y^2}(x, y, t) \\
u(x, y, 0) &= f(x, y)
\end{aligned}
\right. 
\quad (11.4) $$

If $f(x, y)$ is an image (with $(x, y) \in D \subset \mathbb{R}^2$, $D$ is the spatial support), this defines a filtering method that is equivalent to the convolution of $f$ by a Gaussian function [19]. This equation can be solved by a finite difference numerical method.

### 11.1.1 Numerical scheme

The following notations are employed. $N, S, E, W$ stand for north, south, east and west.

$$ 
\begin{aligned}
+\delta^N u &= \frac{u(x,y+h)-u(x,y)}{h} \\
-\delta^S u &= \frac{u(x,y-h)-u(x,y)}{h} \\
-\delta^E u &= \frac{u(x-h,y)-u(x,y)}{h} \\
+\delta^W u &= \frac{u(x+h,y)-u(x,y)}{h}
\end{aligned} 
\quad (11.5) $$

Thus, the numerical scheme is (h has value 1):

$$ \frac{u^{t+1} - u^t}{\delta t} = \frac{1}{h} \{ \delta^N u - \delta^S u + \delta^E u - \delta^W u \} \quad (11.6) $$

<figure><img src="image_306.png" /></figure>

1. Code the numerical scheme. It should use two parameters: the number of iterations $n$, and the step time $\delta t$.

2. Filter the original image by varying the parameters of the discrete diffusion process.

3. Comment the filtering results. Compare the results with the Gaussian filter.

## 11.2 Nonlinear diffusion

Image filtering by nonlinear diffusion reduces noise in a controled way. The diffusion coefficient is locally adapted, becoming negligible as object boundaries are
approached.

The heat equation is replaced by the following PDE:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial u}{\partial t}(x, y, t) = \operatorname{div} (c(\|\nabla u(x, y, t)\|) \cdot \nabla u(x, y, t)) \\ u(x, y, 0) = f(x, y) \end{array} \right. \qquad (11.7) $$

with c the diffusion function satisfying the following properties:

* $c(0) = 1$,

* $\lim_{s \to +\infty} sc(s) = 0$,

* $\forall s > 0, c'(s) \le 0$.

Perona and Malik have proposed 2 diffusion coefficients [27, 5]:

* $c_1 : \exp\left(-\left(\frac{s}{\alpha}\right)^2\right)$,

* $c_2 : \frac{1}{1 + \left(\frac{s}{\alpha}\right)^2}$.

### 11.2.1 Numerical scheme

The numerical scheme is the following:

$$ \frac{u^{t+1} - u^t}{\delta t} = \frac{1}{h^2} \left\{ \begin{aligned} & c(|\delta^N u|) \cdot \delta^N u + c(|\delta^S u|) \cdot \delta^S u \\ & + c(|\delta^E u|) \cdot \delta^E u + c(|\delta^W u|) \cdot \delta^W u \end{aligned} \right\} \qquad (11.8) $$

<figure><img src="image_307.png" /></figure>

1. Code the numerical scheme

2. Filter the original image by varying the parameters of the discrete nonlinear diffusion process.

3. Comment and compare the filtering results with the previous scheme.
## 11.3 Degenerate diffusion

The following degenerate PDEs have been shown to be equivalent to the morphological operators of dilation and erosion (using a disk as structuring element, see tutorial on mathematical morphology):

$$ \left\{ \begin{array}{l} \frac{\partial u}{\partial t}(x, y, t) = +\|\nabla u(x, y, t)\| \\ u(x, y, 0) = f(x, y) \end{array} \right. \qquad (11.9) $$

$$ \left\{ \begin{array}{l} \frac{\partial u}{\partial t}(x, y, t) = -\|\nabla u(x, y, t)\| \\ u(x, y, 0) = f(x, y) \end{array} \right. \qquad (11.10) $$

### 11.3.1 Numerical schemes

The previous numerical scheme is easy to find, but the results present large differences with the dilation and erosion operators (shocks due to discontinuities in the original image). This is why the following numerical scheme is preferred [23]:

For the erosion:

$$ \frac{u^{t+1} - u^t}{\delta t} = \frac{1}{h^2} \sqrt{ \max^2(0, -\delta^N u) + \min^2(0, -\delta^S u) \dots \\ + \max^2(0, -\delta^W u) + \min^2(0, -\delta^E u) } \qquad (11.11) $$

For the dilation:

$$ \frac{u^{t+1} - u^t}{\delta t} = \frac{1}{h^2} \sqrt{ \min^2(0, -\delta^N u) + \max^2(0, -\delta^S u) \dots \\ + \min^2(0, -\delta^W u) + \max^2(0, -\delta^E u) } \qquad (11.12) $$

<figure><img src="image_308.png" /></figure>

1. Code the numerical scheme.

2. Filter the original image while varying the parameters of the discrete degenerate diffusion process.

3. Comment the morphological filtering results and compare the numerical scheme to the approach with operational windows.
<figure><img src="image_309.png" /></figure>

## 11.4. Python correction

<figure><img src="image_310.png" /></figure>

<figure><img src="image_311.png" /></figure>

import numpy as np

2 from scipy import ndimage, misc

import matplotlib, pyplot as plt

### 11.4.1 Linear diffusion

The linear diffusion filter is equivalent to the classical gaussian filter. The gradients are decomposed into the 4 directions so that the next developments (non linear diffusion filters) are made easy.

<figure><img src="image_312.png" /></figure>

1 def linearDiffusion (I, nbIter, dt): # linear diffusion # I: image # nbIter: number of iterations # dt: step ↪ → time

3 hW = np.array ([[[1, -1, 0]]) ;

hE = np.array ([[[0, -1,1]]) ;

5 hN = npiranpose (hW);

hS = npiranpose (hE);

7 Z = I.copy();

for i in range (nbIter): gW = ndimage .convolve (Z, hW, mode ='constant');

9 gE = ndimage .convolve(Z, hE, mode ='constant');

11 hN = ndimage .convolve(Z, hN, mode ='constant');

gS = ndimage .convolve(Z, hS, mode ='constant');

13 Z = Z + dt*(gW+gE+gN+gS);

return Z

The code to filter images is presented as follows. Results are in Fig. 11.2.

<figure><img src="image_313.png" /></figure>

1 I = imageio .imread("cerveau.png") /255.;

2 F = linearDiffusion (I, 50, .05); F2 = linearDiffusion (I, 200, .05);

4 plt .subplot (1,3,1); plt ..imshow(I, cmap = plt.cm .gray);

6 plt .subplot (1,3,2); plt ..imshow(F, cmap = plt.cm .gray);

8 plt .subplot (1,3,3); plt ..imshow(F2, cmap = plt.cm .gray);
Figure 11.2: Linear diffusion filter. Contours are not preserved: this is equivalent to a Gaussian filter.

<figure><img src="image_314.png" /><figcaption>(b) Linear diffusion filter, with $\alpha = 0.1$, $dt = 0.05$ and 10 iterations.</figcaption></figure>

<figure><img src="image_315.png" /><figcaption>(c) Linear diffusion filter, with $\alpha = 0.1$, $dt = 0.05$ and 50 iterations.</figcaption></figure>

<figure><img src="image_316.png" /></figure>

## 11.4.2 Nonlinear diffusion

The nonlinear diffusion is an adaptation of the diffusion to the informations contained in the images. These informations are high frequency components. The following diffusion coefficient is proposed by Perona and Malik.

<figure><img src="image_317.png" /></figure>

```python
def c(I, alpha):
    # diffusion coefficient
    # I: image
    # alpha: diffusion parameter
    return np.exp(-(I/alpha)**2);
```

The gradients are evaluated in the 4 directions, a specific coefficient is thus applied to each one.
<figure><img src="image_318.png" /></figure>

1 def nonlinearDiffusion(I, nbIter, alpha, dt):
    # linear diffusion

3
# I: image
# nbIter: number of iterations

5
hW = np.array([[1, -1, 0]]);

7
hE = np.array([[0, -1, 1]]);
hN = npiranpose(hW);

9
hS = npiranpose(hE);

11
Z = I.copy();

13
for i in range (nbIter):
    #print "%d" % i

15
gW = ndimage(convolve(Z, hW, mode='constant'));

gE = ndimage(convolve(Z, hE, mode='constant'));

17
gN = ndimage(convolve(Z, hN, mode='constant'));
gS = ndimage(convolve(Z, hS, mode='constant'));

19
Z = Z + dt * (c(np.abs(gW), alpha) * gW + c(np.abs(gE), alpha) * gE
             + c(np abs(gN), alpha) * gN + c(np abs(gS), alpha) * gS);

23
return Z

The script to filter the images, for α = 0.1, dt = 0.05 and for 10 and 50
iterations, is the following. Results are presented in Fig. 11.3.

<figure><img src="image_319.png" /></figure>

1 alpha = 0.1;
dt = .05;

3 I = imageio.read("cerveau.png") /255.;

5 F = nonlinearDiffusion (I, 10, alpha, dt);
F2= nonlinearDiffusion (I, 50, alpha, dt);

11.4.3 Degenerate diffusion

Erosion and dilation can theoretically be coded with this numerical scheme. How-
ever, some numerical problems appear in the first version, which are corrected in
the second version.

First version

This first version is the direct transcription of the numerical scheme. As observed
in Fig.11.4, shocks (peaks) appear after a few iterations.
<figure><img src="image_320.png" /><figcaption>Figure 11.3: Nonlinear diffusion filter. Contours are preserved.</figcaption></figure>

<figure><img src="image_321.png" /></figure>

```python
def degenerateDiffusion1(image, nbIter, dt):
    # degenerate diffusion
    # I: image
    # nbIter: number of iterations
    # dt: step time
    h = np.array([[ -1,0,1] ]) ;
    ht= npzigzag(h);
    Zdilation = image;
    Zerosion = image;
    for i in range nbIter :
        gH = ndimageconvolve(Zdilation, h, mode='constant');
        gV = ndimageconvolve(Zdilation, ht, mode='constant');
        jH = ndimageconvolve(Zerosion, h, mode='constant');
        jV = ndimageconvolve(Zerosion, ht, mode='constant');
        Zdilation = Zdilation + dt*npsqrt(gV**2 + gH**2);
        Zerosion = Zerosion + dt*npsqrt(jV**2 + jH**2);
    return ( Zdilation, Zerosion)
```

<figure><img src="image_322.png" /><figcaption>Figure 11.4: Mathematical morphology operations by diffusion.</figcaption></figure>
More sophisticated version

Another scheme, more numerically stable, can be preferred to the previous one.
Results are presented in Fig. 11.5.

<figure><img src="image_323.png" /></figure>

```python
def degenerateDiffusion2(image, nbIter, dt):
    # degenerate diffusion
    # I: original image
    # nbIter: number of iterations
    # dt: step time
    hW = np.array([[1, -1, 0]]);
    hE = np.array([[0, -1, 1]]);
    hN = npurlapose(hW);
    hS = npurlapose(hE);
    Zdilation = image;
    Zerosion = image;

for i in range (nbIter):
    gW = ndimagee convolve(Zdilation, hW, mode='constant');
    gE = ndimagee convolve(Zdilation, hE, mode='constant');
    gN = ndimagee convolve(Zdilation, hN, mode='constant');
    gS = ndimagee convolve(Zdilation, hS, mode='constant');

    jW = ndimagee convolve(Zerosion, hW, mode='constant');
    jE = ndimagee convolve(Zerosion, hE, mode='constant');
    jN = ndimagee convolve(Zerosion, hN, mode='constant');
    jS = ndimagee convolve(Zerosion, hS, mode='constant');

    g = np. sqrt ( np. minimum(0, -gW)**2 + np. maximum(0, gE)**2 + np. minimum(0,
    -> -gN)**2 + np. maximum(0, gS)**2);
    j = np. sqrt ( np. maximum(0, -jW)**2 + np. minimum(0, jE)**2 + np. maximum(0,
    -> -jN)**2 + np. minimum(0, jS)**2 );

Zdilation = Zdilation + dt * g;
Zerosion = Zerosion - dt * j;

return Zdilation, Zerosion
```
<figure><img src="image_324.png" /><figcaption>Figure 11.5: Mathematical morphology operations by diffusion, sophisticated version.</figcaption></figure>
<figure><img src="image_325.png" /></figure>

## 12 Multiscale Analysis

This tutorial aims to study some multiscale image processing and analysis methods, based on pyramidal and scale-space representations.

The different processes will be applied on the following MR image.

<figure><img src="image_326.png" /><figcaption>Figure 12.1: Brain MR image.</figcaption></figure>

<figure><img src="image_327.png" /></figure>

# 12.1 Pyramidal decomposition and reconstruction

## 12.1.1 Decomposition

**Data:** image $A_0$

**Result:** pyramid of approximations $\{A_i\}_i$, pyramid of details $\{D_i\}_i$

for $i=1$ to 3 do

filtering: $F = \text{filt}(A_{i-1})$;

subsampling: $A_i = \text{ech}(F, 0.5)$;

details: $D_i = A_{i-1} - \text{ech}(A_i, 2)$;

end

**Algorithm 3:** The algorithm of the pyramidal decomposition

The initial image (with the highest resolution), denoted $A_0$, represents the level 0 of the pyramid. To build the pyramid, we successively carry out the following steps:

1. a Gaussian filtering,

2. a subsampling,

3. a calculation of the details (residues).
<figure><img src="image_328.png" /></figure>

Make a pyramidal decomposition with 4 levels of the image 'brain'.

### 12.1.2 Reconstruction

To reconstruct the original image, we carry out the following steps at each level of the pyramid:

1. an oversampling,

2. an addition of the details.
**Data:** image $A_3$, pyramid of details{$D_i$}$_i$

**Result:** reconstructed pyramid {$B_i$}$_i$

initialization : $B_3 = A_3$;

for i=3 to 1 do

1 | oversampling: $R = ech(B_i, 2)$;

2 | adding details: $B_{i-1} = R + D_i$

end

**Algorithm 4:** The algorithm of the pyramidal reconstruction

<figure><img src="image_329.png" /></figure>

1. Reconstruct the original image from the last level of the pyramid.

2. Make the same reconstruction without adding the details.

3. Calculate the resulting error between the reconstructed image and the original image.

<figure><img src="image_330.png" /></figure>

# 12.2 Scale-space decomposition and multiscale filtering

We are going to decompose (without any sampling) the image with the highest resolution with a morphological operator, dilation or erosion. The resulting images will have the same size.

## 12.2.1 Morphological multiscale decomposition

<figure><img src="image_331.png" /></figure>

Build the two scale-space decompositions of the 'brain' image, with a disk of increasing radius as a structuring element

<figure><img src="image_332.png" /></figure>

### Informations

See the functions `morphology.erosion` and `morphology.dilation` from `skimage`.
### 12.2.2 Kramer and Bruckner multiscale decomposition

Now, we are going to use the iterative filter of Kramer and Bruckner [20] defined as:

$$MK_B^n(f) = K_B(MK_B^{n-1}(f)) \qquad (12.1)$$

where $B$ denotes a disk of a fixed radius $r$, and:

$$K_B(f)(x) = \begin{cases} D_B(f)(x) & \text{if } D_B(f)(x) - f \leq E_B(f)(x) \\ E_B(f)(x) & \text{otherwise} \end{cases} \qquad (12.2)$$

where $D_B(f), E_B(f)$ represent respectively the dilation and the erosion of the image $f$ with the structuring element $B$.

<figure><img src="image_333.png" /></figure>

Implement and test this filter on the image 'brain' for different values of $n$.
<figure><img src="image_334.png" /></figure>

## 12.3. Python correction

<figure><img src="image_335.png" /></figure>

### 12.3.1 Pyramidal decomposition and reconstruction

#### Decomposition

The following function makes the decomposition of the Laplacian and Gaussian pyramids at the same time. The Laplacian pyramid can be reconstructed without any additional information. This is illustrated in Fig. 12.2.

<figure><img src="image_336.png" /></figure>

```python
def LaplacianPyramidDecomposition(Image, levels, interp='bilinear'):
    """
    Laplacian / Gaussian Pyramid
    The last image of the laplacian pyramid allows a full reconstruction of the
    ↪→ original image.
    Image: original image, float32
    levels : number of levels of decomposition
    interp : interpolation mode for downsizing the image
    """
    returns : pyrL, pyrG: Laplacian and Gaussian pyramids, respectively, as a list
              ↪→ of arrays
    """
    pyrL = [];
    pyrG = [];
    """
sigma = 3.;
for l in range(levels):
    prevImage = Image.copy();
    g = ndimage.gaussian_filter(Image, sigma);
    print(g.dtype)
    """
    Image = misc.imresize(g, .5, interp=interp, mode='F');
    primeImage = misc.imresize(Image, prevImage.shape, interp=interp, mode='F')
              ↪→ ;
    """
    pyrL.append(prevImage - primeImage);
    pyrG.append(prevImage);
    """
pyrL.append(Image);
pyrG.append(Image);
return pyrL, pyrG;
```
<figure><img src="image_337.png" /></figure>

Informations

Notice that with the `misc.imresize` function, the float mode of the operation has to be specified.

## Reconstruction

The reconstruction is straightforward and exact because of the construction of the residue. The details can be filtered (removed for example), thus giving the following result Fig. 12.3.

<figure><img src="image_338.png" /></figure>

def LaplacianPyramidReconstruction(pyr, interp='bilinear'):

Reconstruction of the Laplacian pyramid, starting from the last image
pyr: pyramid of images (list of arrays)
 interp: interpolation mode, for upsizing the image
returns: Image, reconstructed image

Image = pyr[-1];
for i in range(len(pyr)-2, -1, -1):
    Image = pyr[i] + misc.imresize(Image, pyr[i].shape, interp=interp, mode='→ F');

return Image;

### 12.3.2 Scale-space decomposition and multiscale filtering

The pyramid of erosions and dilations is illustrated in Fig.12.4.
Figure 12.2: Gaussian and Laplacian pyramids, for 3 levels of decomposition. The Laplacian pyramid in addition to the last level of the Gaussian pyramid is required to exactly reconstruct the original image.

<figure><img src="image_339.png" /></figure>
<caption>Figure 12.3: Reconstruction of the Laplacian pyramid.</caption>

<figure><img src="image_340.png" /><figcaption>(b) Reconstruction of the pyramid with all the details.</figcaption></figure>

<figure><img src="image_341.png" /></figure>

<figure><img src="image_342.png" /></figure>

```python
1 def morphoMultiscale(I, levels):
2     """
3     Morphological multiscale decomposition
4     I: original image, float32
5     levels: number of levels, int
6     """
7     returns: pyrD, pyrE: pyramid of Dilations / Erosions, respectively
8     """
9     pyrD=[];
10    pyrE=[];
11   for r in np.arange(1, levels):
15        se = morphology.disk(r);
18        pyrD.append(morphology.dilation(I, selem=se));
19        pyrE.append(morphology.erosion(I, selem=se));
20    return pyrD, pyrE;
```

### 12.3.3 Kramer and Bruckner multiscale decomposition

The results are illustrated in Fig.12.5.
<figure><img src="image_343.png" /><figcaption>Figure 12.5: Kramer and Bruckner multiscale decomposition, with $r = 5$.</figcaption></figure>

<figure><img src="image_344.png" /></figure>
<figure><img src="image_345.png" /></figure>

1 def kb(I, r):
    """
3 Elementary Kramer/Bruckner filter. Also called toggle filter .
4 I: image
5 r: radius of structuring element (disk), for max/min evaluation
    """
7 se = morphology.disk(r);
8 D=morphology.dilation(I, selem=se);
9 E=morphology.erosion(I, selem=se);
10 difbool = D-I < I-E;
11 k = D*difbool + E * (~difbool);
12 return k;

<figure><img src="image_346.png" /></figure>

1 def KBmultiscale(I, levels, r=1):
    """
2 Kramer and Bruckner multiscale decomposition
4 I: original image, float32
5 pyrD: pyramid of Dilations
6 pyrE: pyramid of Erosions
8 returns: MKB: Kramer/Bruckner filters
    """
10 MKB = [];
12 MKB.append(I);
14 for i in range(levels):
    MKB.append(kb(MKB[i-1], r));
    return MKB
<figure><img src="image_347.png" /></figure>

## 13 Introduction to tomographic reconstruction

### 13.1 X ray tomography

<figure><img src="image_348.png" /><figcaption>Figure 13.1: Radon transform notations.</figcaption></figure>

An X-ray beam (considered as monochromatic) going through objects is attenuated, following a logarithmic law ($f$ is the linear attenuation, $v$ an elementary volume, $I$ the intensity of the beam entering the volume $dv$ and $I + dI$ the intensity exiting the volume):

$$ \log \frac{I + dI}{I} = -f dv $$

By integrating this formula on a line $D$ (direction of the beam), yields:

$$ I = I_0 \exp \left( \int_D f(x, y) dv \right) $$

where $f(x, y)$ is the attenuation at point $(x, y)$. This principle is used in scanners, scintigraphy, PET...
## 13.2 Acquisition simulation

This first exercise will allow us to simulate the projections. The mathematical operator behind this is the Radon transform.

<figure><img src="image_349.png" /></figure>

* Generate a “phantom” image of size 256 × 256. Display it. Notice that any image may be good for testing the simulation and reconstruction via backprojection.

* To simulate the attenuation process, consider the sum of all gray levels of all pixels. Angles of projections will be between 0 and 180 degrees with a unit step. The result is presented as a sinogram $p_{\theta}(u)$ (see Fig. 13.2), with $u$ the length (in ordinates), and $\theta$ the projection angle (in abscissa).

<figure><img src="image_350.png" /><figcaption>Figure 13.2: Simulated sinogram of the phantom image.</figcaption></figure>

## 13.3 Backprojection algorithm

A non exact reconstruction is the backprojection operator $B$ (it is not the inverse Radon transform). It attributes the value $p_{\theta}(u)$ to each point of the projection line that gave this attenuation, and sum up all contributions from all projections (at the different angles).

$$B[p](x, y) = \int_{0}^{\pi} p_{\theta}(x \cos \theta + y \sin \theta) d\theta$$
A simple (but slow) method is described to compute the backprojection. $p_{\theta}$ is the attenuation vector for a given angle $\theta$.

* Use a command to replicate a matrix to obtain the contribution of each line $D$, like on Fig. 13.3a.

<figure><img src="image_351.png" /><figcaption>(a) Contribution of a line $D$ before rotation.</figcaption></figure>

<figure><img src="image_352.png" /><figcaption>(b) Reconstruction after projection every 45 degrees (4 projections).</figcaption></figure>

<caption>Figure 13.3: Backprojection</caption>

* Each contribution line $D$ should be turned of a certain angle to be added to the overall contribution. This backprojection is really fuzzy compared to the original phantom image (see Fig. 13.3b).

## 13.4 Filtered backprojection

It can be shown that backprojection is in fact the convolution of $f$ by a filter. A solution would be to make a deconvolution in the Fourier domain. The solution used here is the filtered backprojection.

Numerous filters can be used. Among them, the Ram-Lak filter $RL$ is approximated by:

$$ 
\begin{aligned}
k \in [-B; B] \subset \mathbb{N}, RL(k) &= \frac{\pi}{4} \text{ if } k = 0 \\
&= 0 \text{ if } |k| \text{ is even} \\
&= \frac{-1}{\pi k^2} \text{ if } |k| \text{ is odd}
\end{aligned}
$$
?

* Write a function that generates this filter (see prototype). The parameter *width* is the number of points of the resulting vector.

* Modify the backprojection algorithm to filter (by convolution of the projection vector) each contribution line before the summation. Observe the improved result.
<figure><img src="image_353.png" /></figure>

## 13.5. Python correction

<figure><img src="image_354.png" /></figure>

<figure><img src="image_355.png" /></figure>

1. import numpy as np
import matplotlib.pyplot as plt

3. from skimage.io import readimage

5. from skimage transform import radon, rescale, rotate
from scipy import signal import convolve

7. import skimage.io

### 13.5.1 Acquisition simulation

The phantom image is first loaded and displayed.

<figure><img src="image_356.png" /></figure>

1. # Read image
image = readimage("phantom.png", as_gray=True)

3. image = image [::,0];
plt.imshow(image, cmap='gray');

The simulation of the projection is simply an addition of all gray-levels of the pixels, after rotating the image in order to simulate the rotation of the object (or of the sensor). See Fig. 13.3a

<figure><img src="image_357.png" /></figure>

```python
def simuProjection(I, theta):
    """
    simulation of the generation of a sinogram
    I : original image (phantom for example)
    theta: angles of projection
    """
    N = I.shape[1];
    M = len theta;
    S = np.zeros((N,M));

    for i, ang in enumerate theta:
        image1 = rotate(I, ang);
        S [i, i] = np.sum(image1, axis=1);

    return S;
```
### 13.5.2 Backprojection algorithm

The backprojection algorithm will sum-up all the contributions of each projection.

<figure><img src="image_358.png" /></figure>

```python
def backprojection(P, theta, filtre):
    """
    Backprojection of
    P: image
    theta: list of projection angles
    filtre: bool, True if filtered
    """
    N = P.shape[0];
    R = np(zeros((N,N)));

    # in case of filtered back-projection
    if filtre:
        h = RamLak(31);

    # loops over all angles
    for i, ang in enumerate theta:
        proj = P[:, i];

        # filtered back-projection
        if filtre:
            proj = convolve FouilleKron proj, h, mode='same');
        else:
            proj2 = np.matlab.repmat FouilleKron proj2, h, ang);
            R = R + proj2;
    return R transpose();
```

The results is better in the case of a filtered backprojection. The RamLak function is provided and illustrated in Fig.13.4.
<figure><img src="image_359.png" /></figure>

```python
def RamLak(width):
    """
    Ramlak filter of size width
    width must be odd
    """
    ramlak = np.zeros((2*width+1,));
    for indice, val in enumerate(np.arange(-width, width+1)):
        val = np.abs(val);

    if val==0: # center
        ramlak[index]=np.pi/4;
    else:
        if val%2==1: # even indices
            ramlak[index]=-1/(np.pi*val**2);
        else: # odd indices
            ramlak[index]=0;

    return ramlak;
```

<figure><img src="image_360.png" /><figcaption>Figure 13.4: RamLak function.</figcaption></figure>

The reconstruction of the original image is obtained by the following code:
<figure><img src="image_361.png" /></figure>

2 # Performs unfiltered backprojection

3 ubp = backprojection (P, theta, False);

4 plt.imshow(ubp, cmap='gray');

6

8 # Performs filtered backprojection

9 fbp = backprojection (P, theta, True);

10 plt.imshow(fbp, cmap='gray');

<figure><img src="image_362.png" /><figcaption>(a) Original phantom image.</figcaption></figure>

<figure><img src="image_363.png" /><figcaption>(b) Unfiltered backprojection.</figcaption></figure>

<figure><img src="image_364.png" /><figcaption>(c) Filtered backprojection.</figcaption></figure>

<caption>Figure 13.5: Reconstruction by backprojection.</caption>
# Part II Mathematical Morphology
<figure><img src="image_365.png" /></figure>
<figure><img src="image_366.png" /></figure>

## 14 Binary Mathematical Morphology

The objective of this tutorial is to process binary images with the elementary operators of mathematical morphology. More particularly, different image transformations, based on the morphological reconstruction, will be studied (closing holes, removing small objects...).

The different transformations will be applied on the following images Fig. 14.1:

<figure><img src="image_367.png" /><figcaption>Figure 14.1: Images to use for this tutorial.</figcaption></figure>

### 14.1 Introduction to mathematical morphology

Mathematical morphology started in the 1960s with Serra and Matheron [38]. It is based on Minkowski addition of sets. The main operators are erosion and dilation, and by composition, opening and closing. More informations can be found in [40].

The erosion of a binary set A by the structuring element B is defined by:

$$ \varepsilon_B(A) = A \ominus B = \{z \in A | B_z \subseteq A\} \quad (14.1) $$

where, $B_z = \{b + z | b \in B\}$.
The dilation can be obtained by:

$$ \delta_B(A) = A \oplus B = \{z \in A | (B^s)_z \cap A \neq \emptyset\} \quad (14.2) $$

where $B^s = \{x \in E | -x \in B\}$ is the symmetric of B.

By composition, the opening is defined as: $A \circ B = (A \ominus B) \oplus B$. The closing is defined as: $A \bullet B = (A \oplus B) \ominus B$.

## 14.2 Elementary operators

### Elementary operators

Test the functions of dilation, erosion, opening and closing on the image $I_2$ by varying:

1. the shape of the structuring element,

2. the size of the structuring element.

The python functions come from the python module `scipy.ndimage.morphology`. Useful functions are `binary_dilation`, `binary_erosion`, `binary_opening` and `binary_closing`. The function `ndimageleastelfstructure` creates a structuring element.

<figure><img src="image_368.png" /></figure>

## 14.3 Morphological reconstruction

The operator of morphological reconstruction $\rho$ is very powerful and largely used for practical applications. The principle is very simple. We consider two binary images: $I_1$ (the studied binary image) and $M$ (the marker image). The objective is to reconstruct the elements of $I_1$ marked by $M$ as illustrated in the Fig. 14.2.

To do this, we iteratively dilate the marker $M$ while being included in $I_1$ (see Eq.14.3). In order to guarantee this inclusion in $A$, we keep from each dilated set its intersection with $I_1$. The algorithm is stopped when the process dilation-intersection is equal to the identity transformation (convergence).

Let $\delta_I^c(M) = \delta_{B_1}(M) \cap I$ be the dilation of the marker set $M$ constrained to the set $I$. Then, the morphological reconstruction is defined as:

$$ \rho_I(M) = \lim_{n \to \infty} \underbrace{\delta_I \circ \dots \circ \delta_I(M)}_{n \text{ times}} \quad (14.3) $$
<figure><img src="image_369.png" /><figcaption>Figure 14.2: Illustration of morphological reconstruction of $I_1$ by $M$.</figcaption></figure>

**Data:** image *I* and marker *M*

**Result:** reconstructed image $rec(I, M)$

$r = area(M);$

$s = 0;$

while $r \neq s$ do

$s = r;$

$M = I \cap (M \oplus B_1);$

$r = area(M);$

end

$rec(I, M) = M;$

**Algorithm 5:** The algorithm of this morphological reconstruction

<figure><img src="image_370.png" /></figure>

1. Implement the algorithm.

2. Test this operator with the images $I_1$ and $M$.

<figure><img src="image_371.png" /></figure>

Evaluate the area of a set may be done by counting the number of its pixels.

14.4

Operators by reconstruction

<figure><img src="image_372.png" /></figure>

Using the reconstruction operator, implement the 3 following transformations:

1. removing the border objects,

2. removing the small objects,
3. closing the object holes.

Test these operators on the image $I_2$.

<figure><img src="image_373.png" /></figure>

The morphological reconstruction function to use is binaryPropagation in the module ndimage.morphology.

14.5

Cleaning of the image of cells

<figure><img src="image_374.png" /></figure>

1. Threshold the image of cells (Fig. 14.1d).

2. Process the resulting binary image with the 3 cleaning processes of the previous question.
<figure><img src="image_375.png" /></figure>

## 14.6. Python correction

<figure><img src="image_376.png" /></figure>

<figure><img src="image_377.png" /></figure>

fromigliрован и從scipy中import ndimage, misc

2import numpy as np

4import matplotlib, plt

### 14.6.1 Elementary operators

First of all, one have to create a structuring element. The following function can be
used to create a circular structuring element. A square can also be used.

<figure><img src="image_378.png" /></figure>

def disk radius:

2 # defines a circular structuring element with radius given by 'radius'
x = np浓括(-radius, radius + 1, 1);

4
xx, yy = np.濃括grid(x, x);
d = np.√(xx**2 + yy**2);

6
return d <= radius;

The following code presents basic mathematical morphology operations. First
of all, declare the structuring element.

<figure><img src="image_379.png" /></figure>

# read binary image (and ensure binarization)

2 B = imageio.imead("B.jpg");
B = B>100;

4 # Structuring element
square = np.ones((5,5));

Erosion and dilation are the two elementary function of mathematical morphology.
<figure><img src="image_380.png" /></figure>

1 # Erosion

Bsquare_erode = ndimage.morphology-binary_erosion(B, structure =square);

3 plt .subplot (231);

4 plt .imshow(Bsquare_erode); plt .title("erosion")

5 imageio .imwrite('erosion.png', Bsquare_erode);

# Dilation

7 Bsquare_dilate = ndimage.morphology-binary_dilation(B, structure =square);

8 plt .subplot (232);

9 plt .imshow(Bsquare_dilate); plt .title("dilation")

imageio .imwrite('dilation.png', Bsquare_dilate);

Opening and closing are a combination of the two previous functions.

<figure><img src="image_381.png" /></figure>

# Opening

2 Bsquare_open = ndimage.morphology-binary_opening(B, structure =square);

3 plt .subplot (233);

4 plt .imshow(Bsquare_open); plt .title("opening")

imageio .imwrite('open.png', Bsquare_open);

6 # Closing

Bsquare_close = ndimage.morphology-binary_closing(B, structure =square);

8 plt .subplot (234);

9 plt .imshow(Bsquare_close); plt .title("closing")

10 imageio .imwrite('close.png', Bsquare_close);

The results are presented in Fig. 14.3.
<figure><img src="image_382.png" /><figcaption>Figure 14.3: Basic mathematical morphology operations.</figcaption></figure>

## 14.6.2 Morphological reconstruction

The algorithm of morphological reconstruction is coded like this in python:

<figure><img src="image_383.png" /></figure>

def reconstruct(image, mask):

2 # should be binary images

4 M = np minimum mask, image);

4
area = ndimage.measurements.sum(M);

6 s=0

8 se = np.array ([ [ 0, 1, 0 ], [ 1, 1, 1 ], [ 0, 1, 0 ] ] );

10 while (area != s):
    s = area;
    M = np minimum(image, ndimage.morphology.binary_dilation(M, structure = se))
    ⊲ ;

12 area = ndimage.measurements.sum(M);

14 return M

The Fig. 14.4 illustrates the morphological reconstruction.
<figure><img src="image_384.png" /></figure>

A=imageio.read('A.jpg');
2 A = A > 100;
M=imageio.read('M.jpg');
4 M = M > 100;
# reconstruction de A par M
6 AM=reconstruct(A, M);

8 # display results
    plt Subplot (1, 3, 1);
10 plt .imshow(A);
    plt .subplot (1, 3, 2);
12 plt .imshow(M);
    plt .subplot (1, 3, 3);
14 plt .imshow(AM);
    plt .show();

<figure><img src="image_385.png" /><figcaption>Figure 14.4: Morphological reconstruction.</figcaption></figure>

14.6.3 Operators by reconstruction

Remove objects touching the borders

The Fig. 14.5 illustrates the suppression of objects touching the borders of the
image. If B represents the border of the image (create an array of the same size
as the image, with zeros everywhere and ones at the sides), then this operation is
defined by:

$$
\text{killBorders}(I) = I \setminus \rho_I(\mathcal{B})
$$
<figure><img src="image_386.png" /><figcaption>Figure 14.5: Suppress border objects.</figcaption></figure>

<figure><img src="image_387.png" /></figure>

def killBorders (A):

2 # remove cells touching the borders of the image
m, n = A.shape

4 M = np.zeros((m,n));
M[0] = 1;

6 M[m-1, :] = 1;
M[:, 0] = 1;

8 M[:,n-1] = 1;
M = reconstruct(A, M);

10 return A-M

## Remove small objects

This is illustrated in Fig. 14.6. It consists in an erosion followed by a reconstruction. The structuring element used in the erosion defines the objects considered as “small”.

$$ \text{killSmall}(I) = \rho_I(\varepsilon(I)) $$

<figure><img src="image_388.png" /></figure>

def killSmall (A, n):

2 # destroy small objects ( size smaller than n)
se = np.ones((n, n));

4 M = ndimage.morphology.binary_erosion(A, structure =se);
return reconstruct(A, M);
<figure><img src="image_389.png" /><figcaption>Figure 14.6: Small objects removal.</figcaption></figure>

### Close holes in objects

This is illustrated in Fig. 14.7. The operation is given by the following equation, with $B$ the border of image $I$, and $X^C$ denoting the complementary of set $X$:

$$ removeHoles(I) = \{\rho_{IC}(B)\}^C $$

<figure><img src="image_390.png" /></figure>

1. def closeHoles(A):
   # close holes in objects

3. Ac = ~A;
m,n = A.shape;

5. M = np.zeros((m,n));
M[0] = 1;

7. M[m-1, :] = 1;
M[0, :] = 1;

9. M[1, n-1] = 1;
M = reconstruct(Ac, M);

11. return ~M

### 14.6.4 Application

The application shown in Fig. 14.8 presents the segmentation of blood cells after removing small cells, closing holes and removing the cells touching the borders of the image.
<figure><img src="image_391.png" /><figcaption>Figure 14.7: Hole filling.</figcaption></figure>

<figure><img src="image_392.png" /></figure>

1 cells = imageio.read('cells.jpg')<98;
imageio.write('cellsbw.png', cells);

3 B = closeHoles( cells );

B = killBorders (B);

5 B = killSmall (B, 5);

plt .subplot (1,2,1);

7 plt .imshow (cells);

plt .subplot (1,2,2);

9 plt .imshow(B);

plt .title('clean image')

11 plt .show()

imageio.write('clean.png', B);

<figure><img src="image_393.png" /><figcaption>Figure 14.8: Segmentation of the cells.</figcaption></figure>
<figure><img src="image_394.png" /></figure>
<figure><img src="image_395.png" /></figure>

## 15 Morphological Geodesic Filtering

This tutorial aims to test different morphological filters and particularly the geodesic ones (using reconstruction) for gray-level images.

The different processes will be applied on the following images:

<figure><img src="image_396.png" /><figcaption>(a) Lena</figcaption></figure>

<figure><img src="image_397.png" /><figcaption>(b) 2-D section of a cement paste (X-ray tomography)</figcaption></figure>

<figure><img src="image_398.png" /></figure>

### 15.1 Morphological centre

The morphological centre is an auto-dual filter using a family of operators {$\psi_i$}:

$$C(f) = (f \lor (\land\{\psi_i(f)\}) \land (\lor\{\psi_i(f)\})) \quad (15.1)$$

<figure><img src="image_399.png" /></figure>

* Implement this transformation with the family {$\gamma\phi\gamma, \phi\gamma\phi$} where $\gamma$ denotes the morphological opening and $\phi$ the morphological closing.

* Test this operator by varying the size of the structuring element.

* Add a 'salt and pepper' noise to the 'Lena' image and compare the morphological center with the median filtering.

* Try to reduce the noisy image 'cement paste'.
## 15.2 Alternate sequential filters (ASF)

ASF can be defined from a family of openings and closings:

$$N_i(f) = \gamma_i \phi_i \circ \gamma_{i-1} \phi_{i-1} \dots \gamma_2 \phi_2 \circ \gamma_1 \phi_1(f) \quad (15.2)$$

$$M_i(f) = \phi_i \gamma_i \circ \phi_{i-1} \gamma_{i-1} \dots \phi_2 \gamma_2 \circ \phi_1 \gamma_1(f) \quad (15.3)$$

where $\gamma_k$ (resp. $\phi_k$) denotes the opening (resp. closing) with a structuring element of size $k$.

<figure><img src="image_400.png" /></figure>

* Implement these two operators.

* Test these transformations on the noisy image by varying the parameter $i$ of the filter.

## 15.3 Reconstruction filters

The geodesic dilation of size 1 and $n$ are respectively defined as:

$$\delta_f(g) = \wedge(\delta_{B_1}(g), f) \quad (15.4)$$

$$\delta_f^n(g) = \delta_f(\delta_f \dots (\delta_f(g))) \quad (15.5)$$

where $\delta_{B_1}$ denotes the morphological dilation with a disk of radius 1 as structuring element. The opening $(\gamma_k^{rec}(f))$ and closing $(\phi_k^{rec}(f))$ by reconstruction are then defined as:

$$\gamma_k^{rec}(f) = \bigvee \{\delta_f^n(\epsilon_{B_k}(f)), n > 0\} \quad (15.6)$$

$$\phi_k^{rec}(f) = M - \gamma_k^{rec}(M - f) \quad (15.7)$$

where $\epsilon_{B_k}$ denotes the morphological erosion with a disk $B$ of radius $k$ as structuring element.

<figure><img src="image_401.png" /></figure>

* Implement these two operators $\gamma_k^{rec}(f)$ and $\phi_k^{rec}(f)$.

* Test these transformations by varying the parameter $k$ of the filter.

* Implement and test (on the noisy image) the filters of morphological center and ASF using the geodesic operators. Compare with the classical ones.
<figure><img src="image_402.png" /></figure>

## 15.4. Python correction

<figure><img src="image_403.png" /></figure>

<figure><img src="image_404.png" /></figure>

from skimage.util import random.noise

2 from scipy import misc

import matplotlib.pyplot as plt

4 from skimage import morphology as m

from skimage import filters

6 import numpy as np

### 15.4.1 Morphological center

The noisy image is obtained with the function `random.noise` from `skimage.util`.

<figure><img src="image_405.png" /></figure>

L = imageio.imread('lena512.bmp');

2 A = random.noise(L, mode='s&p', amount=.04);

Following the definition, the morphological center is obtained with the following code, and illustrated in Fig.15.1:

<figure><img src="image_406.png" /></figure>

def morphoCenter(I, c, o, selem=m.disk(1)):

2

**coco**

4 coc = c(o(c(I, selem=selem), selem=selem), selem=selem);

oco = o(c(o(I, selem=selem), selem=selem), selem=selem);

6 cMin = npella minimum(oco, coc);

cMax = npella maximum(oco, coc);

8 F = npella minimum(npella maximum(A, cMin), cMax);

return F;

10

B = morphoCenter(A, m.closing, m.opening);

12 Bmed= filters .median(A, selem= np. ones ((3,3) ))
<figure><img src="image_407.png" /><figcaption>Figure 15.1: Morphological center compared to the classical median filter.</figcaption></figure>

## 15.4.2 Alternate sequential filters

The order of these filters are often chosen empirically. The results are illustrated in Fig.15.2.

<figure><img src="image_408.png" /></figure>

```python
def asf_n(I, order=3):
    F = I.copy();
    for r in np.arange(1, order+1):
        se = m.disk(r);
        F = m.opening(m.closing(F, selem=se), selem=se);
    return F;
```

<figure><img src="image_409.png" /></figure>

```python
def asf_m(I, order=3):
    F = I.copy();
    for r in np.arange(1, order+1):
        se = m.disk(r);
        F = m.closing(m.opening(F, selem=se), selem=se);
    return F;
```
<figure><img src="image_410.png" /><figcaption>Figure 15.2: Alternate Sequential Filters compared to original image.</figcaption></figure>

15.4.3 Geodesic reconstruction filters

These two functions are simply implemented using erosion and reconstruction
operators. Notice the duality property, that is used to code *closerec*. In this example,
8 bits images (unsigned) are considered, and the results are illustrated in Fig.15.3.

<figure><img src="image_411.png" /></figure>

def openrec(I, selem=m.disk(1)):

2 B = m.erosion(I, selem=selem);

F = m.reconstruction (B, I);

4 return F;

<figure><img src="image_412.png" /></figure>

def closerec(I, selem=m.disk(1)):

2 F = 255-openrec(255-I, selem=selem);

return F;
<figure><img src="image_413.png" /><figcaption>Figure 15.3: Opening and closing by reconstruction.</figcaption></figure>

### 15.4.4 ASF by reconstruction

This is an example of a 3rd order alternate sequential filter, illustrated in Fig.15.4.

<figure><img src="image_414.png" /></figure>

```python
1 def asfrec(I, order=3):
2     A = I.copy();
3     for r in np.arange(1, order+1):
4         se = m.disk(r);
5         A = closerec(openrec(A, selem=se), selem=se);
6     return A;
```

<figure><img src="image_415.png" /><figcaption>Figure 15.4: Alternate Sequential Filtering by reconstruction.</figcaption></figure>

### 15.4.5 Morphological center by reconstruction

Morphological center by reconstruction replaces the opening and closing operations by their equivalent by reconstruction (see Fig.15.5).
<figure><img src="image_416.png" /></figure>

```python
def centerrec(I, selem=m.disk(1)):
    2 """
    4 B = morphoCenter(I, closerec, openrec, selem=selem);
    return B;
```

<figure><img src="image_417.png" /><figcaption>Figure 15.5: Morphological center by reconstruction.</figcaption></figure>
<figure><img src="image_418.png" /></figure>
<figure><img src="image_419.png" /></figure>

## 16 Morphological Attribute Filtering

This tutorial aims to test the attribute morphological filters for gray-level images. The objective of such filters is to remove the connected components of the level sets which do not satisfy a specific criterion (area, elongation, convexity...).

Fig.16.1 are presented some results of attribute filtering applied on the grayscale image.

<figure><img src="image_420.png" /><figcaption>Figure 16.1: Attribute filtering examples.</figcaption></figure>

A binary connected opening $\Gamma_x$ transforms a binary image $f$ given a pixel $x$, by keeping the connected component that contains $x$ and removing all the others. Binary trivial opening $\Gamma_T$ operates on a given connected component $X$ according to an increasing criterion $T$ applied to the connected component. If the criterion is satisfied, the connected component is preserved, otherwise it is removed according to:

$$ \Gamma_T(X) = \begin{cases} X & \text{if } T(X) = \text{true} \\ \emptyset & \text{if } T(X) = \text{false} \end{cases} \qquad (16.1) $$
In general, one or more features of the connected component, on which the filter is applied, are compared to a given threshold defined by the rule.

Binary attribute opening $\Gamma^T$, given an increasing criterion $T$, is defined as a binary trivial opening applied on all the connected components of $F$. This can be formally represented as:

$$ \Gamma^T(F) = \bigcup_{x \in F} \Gamma_T(\Gamma_x(F)) \qquad (16.2) $$

If the criterion $T$ evaluated in the transformation is not increasing, e.g., when the computed attribute is not dependent itself on the scale of the regions (e.g. shape factor, orientation, etc.), the transformation also becomes not increasing. Even if the increasingness property is not fulfilled, the filter remains idempotent and anti-extensive. For this reason, the transformation based on a non-increasing criterion is not an opening, but a thinning. In this case, one speak about a binary attribute thinning.

Attribute openings and thinnings introduced for the binary case can be extended to grayscale images employing the threshold decomposition principle. A grayscale image can be represented by thresholding the original image at each of its graylevels. Then, the binary attribute opening can be applied to each binary image and the grayscale attribute opening is given by the maximum of the results of the filtering for each pixel and can be mathematically represented as:

$$ \gamma^T(f)(x) = \max\{k; x \in \Gamma^T(Th_k(f))\} \qquad (16.3) $$

where $Th_k(f)$ is the binary image obtained by thresholding $f$ at graylevel $k$.

The extension to grayscale of the binary attribute thinning is not straightforward and not unique. For example, a possible definition of a *grayscale attribute thinning* can be analogously defined as Equation (3). However, other definitions are possible, according to the filtering rule considered in the analysis (max, min, subtractive, direct...). Again, if the criterion $T$ is increasing (i.e. the transformation is an opening), all the strategies lead to the same result.

Attribute openings are a family of operators that also includes openings by reconstruction. If we consider a binary image with a connected component $X$ and the increasing criterion "the size of the largest square enclosed by $X$ must be greater than $\lambda$", the result of the attribute opening is the same as applying an opening by reconstruction with a squared structuring element of size $\lambda$.

## 16.1 Binary attribute opening/thinning

We are going to implement and test some binary attribute filters according to different criteria (area, elongation, convexity).
?

* Regarding the area criterion, implement the associated binary attribute opening. Test on the binary image.

* Regarding the elongation criterion, implement the associated binary attribute thinning (non-increasing criterion). Test on the binary image.

<figure><img src="image_421.png" /></figure>

Code your own function bwFilter based on regionprops from skimage.measure.

16.2

Grayscale attribute opening/thinning

We are going to extend the binary filters to grayscale images.

?

* Implement the decomposition operator into binary sections.

* Conversely, implement the reconstruction operator from binary sections.

* Implement and test the grayscale attribute filters according to the following criteria: area, elongation, convexity.
<figure><img src="image_422.png" /></figure>

## 16.3. Python correction

<figure><img src="image_423.png" /></figure>

<figure><img src="image_424.png" /></figure>

from skimage import morphology as m

2. from scipy import misc
import matplotlib.pyplot as plt

4. from skimage import regionprops
from skimageapon measure label

6. import numpy as np

### 16.3.1 Binary attribute filtering

If `bw` is a binary image, the different attribute are evaluated with `regionprops` and
filtered according to the given upper and lower thresholds.

<figure><img src="image_425.png" /></figure>

def bwFilter(bw, attribute, thresholds):

2
***

binary filtering according to attribute

4

bw: binary image

attribute : string representing attribute , defined by regionprops

6

thresholds : threshold values, inside which objects are removed

returns binary image

8

###

F = bw.copy();

10

L = label (bw);

for region in regionprops(L):

12

a = getattr(region, attribute);

if a < thresholds [1] and a >= thresholds [0]:

14

F[L==region.label] = False;

return F;

### 16.3.2 Grayscale filtering

The grayscale filtering is the previous binary filtering process applied to all level-sets
of the original image.
<figure><img src="image_426.png" /></figure>

```python
1 def grayFilter(I, attribute, thresholds):
    """
    grayscale image filtering by attribute
    for 8 bits unsigned images
    I: original grayscale image (N, M) ndarray
    attribute : string representing attribute, defined by regionprops
    thresholds : threshold values, inside which objects are removed
    returns grayscale filtered image
    """
    N, M = I.shape
    F = np.zeros((N, M, 256));
    for s in range(256):
        F[:,:, s] = s * bwFilter(I>=s, attribute = attribute, thresholds = thresholds);

11

13

15

17

# reconstruction
R = np.amax(F, axis=2);
return R;
```

Then, for each level, the binary set is filtered by some attribute, and the resulting image is reconstructed by taking the maximum value on all levels.

The shape filtering process do not take the results from the `regionprops` function. The shape is filtered by using the morphological function of reconstruction.

<figure><img src="image_427.png" /></figure>

```python
def shapeFilter(I, selem=m.square(25)):
    """
    image filtering when attribute is a shape of a given size, defined by selem
    I: grayscale image
    selem: structuring element
    returns: grayscale filtered image
    """
    N, M = I.shape
    F = np.zeros((N, M, 256));
    for s in range(256):
        F[:,:, s] = s * m.reconstruction(m.opening(I>=s, selem=selem), I>=s);

10

12

14

# reconstruction
R = np.amax(F, axis=2);
return R;
```

The examples are generated via the following code:
<figure><img src="image_428.png" /></figure>

1. F = bwFilter(I>50, 'area', (0,5000));

2. F = grayFilter(I, 'area', (0,1000));

3. F = grayFilter(I, 'eccentricity', (0, 0.75));

4. F = grayFilter(I, 'solidity', (0, 0.75));

5. F = shapeFilter(I);

Results are illustrated in Fig.16.1.
<figure><img src="image_429.png" /></figure>

## 17 Morphological skeletonization

This tutorial aims to skeletonize objects with specific tools from mathematical morphology (thinning, maximum ball...).

The different processes will be applied on the following image:

<figure><img src="image_430.png" /></figure>

<figure><img src="image_431.png" /></figure>

### 17.1 Hit-or-miss transform

The hit-or-miss transformation enables specific pixel configurations to be detected. Based on a pair of disjoint structuring elements $T = (T^1, T^2)$, this transformation is defined as:

$$ \eta_T(X) = \{x, T_x^1 \subseteq X, T_x^2 \subseteq X^c\} \quad (17.1) $$

$$ = \epsilon_{T^1}(X) \cap \epsilon_{T^2}(X^c) \quad (17.2) $$

where $\epsilon_B(X)$ denotes the erosion of $X$ using the structuring element $B$.

<figure><img src="image_432.png" /></figure>

1. Implement the hit-or-miss transform.

2. Test this operator with the following pair of disjoint structuring elements:

<table><tr><td>+1</td><td>+1</td><td>+1</td></tr><tr><td>0</td><td>+1</td><td>0</td></tr><tr><td>-1</td><td>-1</td><td>-1</td></tr></table>

$$ T^1 = \begin{bmatrix} 1 & 1 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix}, $$

$$ T^2 = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 1 & 1 & 1 \end{bmatrix} $$

where the points with +1 (resp. -1) belong to $T^1$ (resp. $T^2$).
## 17.2 Thinning and thickening

Using the hit-or-miss transform, it is possible to make a thinning or thickening of a binary object in the following way:

$$ \theta_T(X) = X \setminus \eta_T(X) \qquad (17.3) $$

$$ \chi_T(X) = X \cup \eta_T(X^c) \qquad (17.4) $$

These two operators are dual in the sense that $\theta_T(X) = (\chi_T(X^c))^c$.

<figure><img src="image_433.png" /></figure>

1. Implement these two transformations.

2. Test these operators with the previous pair of structuring elements.

## 17.3 Topological skeleton

By using the two following pairs of structuring elements with their rotations (90°) in an iterative way (8 configurations are thus defined), the thinning process converges to a resulting object which is homothetic (topologically equivalent) to the initial object.

<table><tr><td>+1</td><td>+1</td><td>+1</td></tr><tr><td>0</td><td>+1</td><td>0</td></tr><tr><td>-1</td><td>-1</td><td>-1</td></tr></table>

<table><tr><td>0</td><td>+1</td><td>0</td></tr><tr><td>+1</td><td>+1</td><td>-1</td></tr><tr><td>0</td><td>-1</td><td>-1</td></tr></table>

<figure><img src="image_434.png" /></figure>

1. Implement this transformation (the convergence has to be satisfied).

2. Test this operator and comment.

## 17.4 Morphological skeleton

A ball $B_n(x)$ with center $x$ and radius $n$ is maximum with respect to the set $X$ if there exists neither index $k$ nor centre $y$ such that:

$$ B_n(x) \subseteq B_k(y) \subseteq X $$
In this way, the morphological skeleton of a set X is constituted by all the centers of maximum balls. Mathematically, it is defined as:

$$S(X) = \bigcup_{r} \epsilon_{B_r(0)}(X) \setminus \gamma_{B_1(0)}(\epsilon_{B_r(0)}(X)) \quad (17.5)$$

<figure><img src="image_435.png" /></figure>

1. Implement this transformation.

2. Test this operator and compare it with the topological skeleton.
<figure><img src="image_436.png" /></figure>

## 17.5. Python correction

<figure><img src="image_437.png" /></figure>

<figure><img src="image_438.png" /></figure>

1. fromigliimport ndimage
import numpy as np

3. import imageio

5. import matplotlib .pyplot as plt

### 17.5.1 Hit or miss transform

The hit-or-miss transform is illustrated in Fig.17.1.

<figure><img src="image_439.png" /></figure>

1. def hitormiss(X, T):
    """

3. hit or miss transform
X: binary image

5. T: structuring element (values -1, 0 and 1)

7. return : result of hit or miss transform (binary image)
    """

9. T1 = (T==1);
T2 = (T==−1);

11. E1 = ndimage.morphology.binary_erosion(X, T1);
E2 = ndimage.morphology.binary_erosion(np. logical_not (X), T2);

13. B = np .minimum(E1, E2);
return B;

### 17.5.2 Thinning and thickening

The code is split into 2 elementary operations of thinning and thickening, so that the thinning consists in iterating this operation. The thickening operation is obtained by thinning the complementary set.
<figure><img src="image_440.png" /><figcaption>Figure 17.1: Elementary functions</figcaption></figure>

<figure><img src="image_441.png" /></figure>

1. def elementary_thinning(X, T):

2. thinning function

X: binary image

5. T: structuring element (values -1, 0 and 1)

7. return: result of thinning

9. B = np minimum(X, np logical not(hitormiss(X, T))); return B;

<figure><img src="image_442.png" /></figure>

def elementary_thickening(X, T):

2.
thickening function

4. X: binary image

6. T: structuring element (values -1, 0 and 1)

8. return: result of thickening

10. B = not element_thinning(not(X), T)); return B;
<figure><img src="image_443.png" /></figure>

def thinning(X, TT):

2

###

morphological thinning

4

TT is a configuration of 8 pairs of structuring elements

###

6

for T in TT:

X = elementary_thinning(X, T);

8

return X;

## 17.5.3 Skeletons

The topological skeleton is the iteration of the operation of thinning, for structuring elements defined like:

<figure><img src="image_444.png" /></figure>

1 TT = [];

TT.append(np.array([[-1, -1, -1],[0,1,0],[1,1,1]]));

3

TT.append(np.array([0, -1, -1],[1,1, -1],[0,1,0]));

TT.append(np.array([[1,0, -1],[1,1, -1],[1,0, -1]]));

5

TT.append(np.array ([0,1,0],[1,1, -1],[0, -1, -1]));

TT.append(np.array ([1,1,1],[0,1,0],[ -1, -1, -1]));

7

TT.append(np.array ([0,1,0],[ -1,1,1],[ -1, -1,0]));

TT.append(np.array([[-1,0,1],[ -1,1,1],[ -1,0,1]]));

9

TT.append(np.array([[-1, -1, 0],[ -1,1,1],[0,1,0]]));

<figure><img src="image_445.png" /></figure>

1 def topological_skeleton (X, TT):

###

Topological skeleton: preserves topology

3

X: binary image to be transformed

5

TT: set of pairs of structuring elements

return skeleton

###

7

B = np. logical_not (np.copy(X));

9

while not(np. all (X == B)):

B = X;

11

X = thinning(X, TT);

return B;

The topological skeleton is the iteration of the thinning with structuring elements in all 8 directions. It has the property of preserving the topology of the
discrete structures, contrary to the morphological skeleton (see Fig.17.2). The mor-
phological skeleton does not preserve the connexity of the branches, but it can
be used to reconstruct the original image. Pay attention to the construction of
structuring elements, which should be homothetic ($B_r = \underbrace{B_1 \oplus \cdots \oplus B_1}_{r \text{ times}}$).

<figure><img src="image_446.png" /></figure>

def morphological_skeleton(X):

'''

2

morphological skeleton

4

X: binary image

using disk structuring elements

6

in order to perform the reconstruction from the skeleton, one has to store
the value of S for each size of structuring element

8

return: S, morphological skeleton, grayscale image
'''

10

strel_size = -1;

12

pred = True;

14

S = np(zeros(X.shape));

se = ndimageurl generate_binary_structure(2, 1);

16

E = np.copy(X);

while pred:

18

strel_size += 1;

E = ndimage.morphology.binary_erosion(X, se);

20

if np. all(E==0):

pred = False;

22

D = ndimage.morphology.binary_dilation(E, se);

S = np.maxim S, (strel_size +1)*np.min X, np.logical_not(D));

24

X = E;

return S;
<figure><img src="image_447.png" /></figure>

```python
def reconstruction_skeleton (S):
    """
    Reconstruction of the original image from the morphological skeleton
    S: Skeleton, as constructed by morphological_skeleton

    return: original image I
    """

X = np.zeros(S.shape).astype("bool");
n = np.max(S);
se = ndimage(generate_binary_structure(2, 1));

for strel_size in range(int(n)):
    Sn = S == strel_size + 1;

    # this is for preserving homothetic structuring elements
    for k in range(strel_size):
        Sn = ndimage.morphology.binary_dilation(Sn, se);

X = np maximum(X, Sn);
return X;
```

<figure><img src="image_448.png" /><figcaption>Figure 17.2: Skeletons. The topology is not preserved in the morphological skeleton, but it can be used to reconstruct the original image. The intensities are inverted in order to facilitate the display.</figcaption></figure>
<figure><img src="image_449.png" /></figure>

# 18 Granulometry

This tutorial aims to compute specific image measurements on digital images. It consists in determining the size distribution of the 'particles' included in the image to be analyzed.

The image measurements will be realized on a simulated image of disks and an image of silicon carbide powder acquired by Scanning Electron Microscopy (SEM, Fig. 18.1).

Figure 18.1: These images are used for computing a granulometry of objects by mathematical morphology.

<figure><img src="image_450.png" /></figure>

## 18.1 Morphological granulometry

### 18.1.1 Introduction to mathematical morphology

The tutorial 14 introduces the basic operations of the mathematical morphology, as well as the morphological reconstruction. More informations can be found in [40].

With Python, the useful functions are `binary_opening` of the module `scipy.ndimage` for the morphological opening, `binary Propagation` for the geodesic reconstruction. For counting the number of objects, one can use `measurements.label` for a labelling and counting algorithm.

<figure><img src="image_451.png" /></figure>
18.1.2 Granulometry

?

* Load and visualize the image 'simulation' Fig.18.1.

* Generate $n$ images by applying morphological openings (with reconstruction) to the binary image with the use of homogeneous structuring elements of increasing size: $1 < \dots < n$.

* Compute the morphological granulometry expressed in terms of surface area density. It consists in calculating the specific surface area of the grains in relation with the size $n$ of the structuring element.

* Compute the morphological granulometry expressed in terms of number density. It consists in calculating the specific number of grains in relation with the size $n$ of the structuring element.

* Compare and discuss.

<figure><img src="image_452.png" /></figure>

Use `generate_binary_structure` and `iterate_structure` from python module `scipy.ndimage` to define homogeneous structuring elements.

18.2 Real application

?

* Load and visualize the image 'powder' of Fig. 18.1b.

* Realize the image segmentation step. I can simply consist in a binarization by a global threshold, followed by hole filling. Noise can also be removed by my mathematical morphology opening.

* Compute the morphological granulometry (surface area, number).

* What can you conclude?
<figure><img src="image_453.png" /></figure>

1. **## Granulometry of synthetic image**
   `# read binary simulated image, normalize it`

3. `I = imageio.read("simulation.png") /255;`
   `I = I [[:,-,2]>.5;`

5. `plt.figure()`;
   `plt.imshow(I);`

7. `plt.show()`;

9. `granulometry(I, 35);`

<figure><img src="image_454.png" /><figcaption>Figure 18.2: Granulometry on simulated image.</figcaption></figure>

### 18.3.2 Powder image and segmentation

First of all, the image must be binarized, i.e. segmented. The following code is a proposition of segmentation, leading to the result presented in Fig. 18.3c.
<figure><img src="image_455.png" /></figure>

1. **## Granulometry of real image**
I = imageio.read("poudre.bmp");

3. ** segmentation**
BW = I>74;

5. **BW = ndimage.morphology-binary_fill_holes(BW);**

7. **suppress small objects**
se = ndimage(generate_binary_structure(2, 1));

9. **m = ndimage.morphology-binary_opening(BW);**

11. **opening by reconstruction**
BW=ndimage.morphology-binary Propagation(m, mask=BW);

13. **imageio.write("segmentation.png", BW);**

15. ** plt.imshow(BW)**

17. ** plt.show()**

17. **granulometry(BW, 20);**

<figure><img src="image_456.png" /><figcaption>Figure 18.3: Results of granulometry analysis of original image from (a).</figcaption></figure>
<figure><img src="image_457.png" /></figure>
# Part III Registration and Segmentation
<figure><img src="image_458.png" /></figure>
<figure><img src="image_459.png" /></figure>

## 19 Histogram-based image segmentation

This tutorial aims to implement some image segmentation methods based on histograms (thresholding and “k-means” clustering).

The different processes will be applied on the following images:

<figure><img src="image_460.png" /><figcaption>(a) Cells.</figcaption></figure>

<figure><img src="image_461.png" /><figcaption>(b) Muscle cells, author: Damien Freyssenet, University Jean Monnet, Saint-Etienne, France.</figcaption></figure>

### 19.1 Manual thresholding

The most simple segmentation method is thresholding.

<figure><img src="image_462.png" /></figure>

* Visualize the histogram of the grayscale image 'cells'.

* Make the segmentation with a threshold value determined from the image histogram.

### 19.2 k-means clustering

Let $X = \{x_i\}_{i \in [1:n], n \in \mathbb{N}}$ be a set of observations (the points) in $\mathbb{R}^d$. The k-means clustering consists in partitioning $X$ in $k$ ($k < n$) disjoint subsets $\tilde{S}$ such that:

$$ \tilde{S} = \arg \min_{S=\{S_i\}_{i \le k}} \sum_{i=1}^{k} \sum_{x_j \in S_i} \|x_j - \mu_i\|^2 \quad (19.1) $$
where $\mu_i$ is the mean value of the elements in $S_i$. The k-means algorithm is iterative. From a set of k initial elements $\{m_i^{(1)}\}_{i \in [1:k]}$ (randomly selected), the algorithm iterates the following (t) steps:

* Each element of X is associated to an element $m_i$ according to a distance criterion (computation of a Voronoi partition):

$$S_i^{(t)} = \{ \mathbf{x}_j : \|\mathbf{x}_j - \mathbf{m}_i^{(t)}\| \le \|\mathbf{x}_j - \mathbf{m}_{i^*}^{(t)}\|, \forall i^* \in [1:k] \} \quad (19.2)$$

* Computation of the new mean values for each class:

$$\mathbf{m}_i^{(t+1)} = \frac{1}{|S_i^{(t)}|} \sum_{\mathbf{x}_j \in S_i^{(t)}} \mathbf{x}_j \quad (19.3)$$

where $|S_i^{(t)}|$ is the number of elements of $S_i^{(t)}$.

## 19.3 Grayscale image, $k = 2$ in one dimension

The objective is to binarize image 'cells', which is a grayscale image. The set $X$ is defined by $X = \{I(p)\}$, for $p$ being the pixels of the image $I$.

<figure><img src="image_463.png" /></figure>

* Implement the algorithm proposed below (Alg. 6).

* Test this operator on the image 'cells'.

* Test another method of automatic thresholding (defined by Otsu in [26]).

* Compare the values of the thresholds (manual and automatic).
**Data:** Original image A

**Data:** Stop condition $\epsilon$

**Result:** thresholded image

Initialize $T_0$, for example at $\frac{1}{2}(\max(A) + \min(A))$;

done ← False;

while NOT done do

Segment the image A with the threshold value $T$;

it generates two classes $G_1$ (intensities $\ge T$) and $G_2$ (intensities $< T$);

Compute the mean values, denoted $\mu_1, \mu_2$, of the two classes $G_1, G_2$, respectively;

Compute the new threshold value $T_i = \frac{1}{2}(\mu_1 + \mu_2)$;

if $|T_i - T_{i-1}| < \epsilon$ then

done ← True

end

end

Segment the image with the estimated threshold value.

**Algorithm 6:** K-means algorithm for automatic threshold computation of grayscale images.

<figure><img src="image_464.png" /></figure>

For use with python, use the module `skimage`. `filter` and the function `threshold_otsu`.

19.4

Simulation example, $k = 3$ in two dimensions

The objective is to generate a set of 2-D random points (within $k = 3$ distinct classes) and to apply the $k$-means clustering for separating the points and evaluating the method (the classes are known!).

<figure><img src="image_465.png" /></figure>

Write a function for generating a set of $n$ random points around a point with coordinates $(x, y)$.

<figure><img src="image_466.png" /></figure>

We will use the python function `randn` from the module ` numpy.random`.

<figure><img src="image_467.png" /></figure>

* Use this function to generate 3 set of points (in a unique matrix) around the points $(0, 0)$, $(3, 4)$ and $(-5, -3)$.
* Use the python function `sklearn.cluster.KMeans` for separating the points. The result is presented in Fig.19.1.

<figure><img src="image_468.png" /></figure>

Verify the utility of the option `n_init=10`.

<figure><img src="image_469.png" /><figcaption>Figure 19.1: Resulting clustering of random points.</figcaption></figure>

## 19.5 Color image segmentation using K-means: $k = 3$ in 3D

The $k$-means clustering is now used for segmenting the color image representing the muscle cells 'Tv16.png'.

?

Which points have to be separated? Transform the original image into a vector of size $N \times 3$ (where $N$ is the number of pixels) which represent the 3 components R, G and B of each image pixel.

<figure><img src="image_470.png" /></figure>

The python function ` numpy.reshape ` does the same transformation.
?

* Visualize the 3-D map (histogram) of all these color intensities.

* Make the clustering of this 3-D map by using the K-means method.

* Visualize the corresponding segmented image.
<figure><img src="image_471.png" /></figure>

## 19.6. Python correction

<figure><img src="image_472.png" /></figure>

### 19.6.1 Manual thresholding

#### Visual analysis of histogram

When manually choosing a threshold value, one has to analysis the histogram (Fig. 19.2).

<figure><img src="image_473.png" /></figure>

```python
import numpy as np

import imageio
import matplotlib.pyplot as plt  # plots

from skimage import filter  # otsu thresholding

# read image
cells =imageio.openread(' cells .png');

# display histogram

fig = plt .figure ();
plt .hist( cells .flatten (), 256)

fig.show();
fig .savefig("histo .pdf");
```

<figure><img src="image_474.png" /><figcaption>Figure 19.2: Original image and its histogram.</figcaption></figure>
Segmentation

<figure><img src="image_475.png" /></figure>

1 fig = plt .figure (); plt .subplot (1,2,1)

3 plt .imshow cells , plt .cm .gray); plt .title ( 'Original image') ; plt .subplot (1,2,2)

5 plt .imshow( cells >80, plt .cm .gray); plt .title ( 'Manual segmentation'); fig .savefig ( "manual.pdf");

19.6.2 Automatic thresholding

<figure><img src="image_476.png" /></figure>

def autothresh(image): """

@param image: image to segment
@return : threshold value

###
6 s = 0.5*( np .amin(image) + np .amax(image));

done = False ;
while ~done:

B = image >=s; 10 sNext = .5*( np .mean(image[B]) + np .mean(image[~B]));

done = abs(s-sNext) <.5; 12 s = sNext;

return s

The results are displayed using the following code (Fig. 19.3):

<figure><img src="image_477.png" /></figure>

1 # Automatic threshold
s_auto= autothresh( cells );

3
# Otsu thresholding

5 s_otsu = filter .threshold_otsu( cells );

7 plt .figure ();
plt .subplot (1,2,1)

9 plt .imshow( cells >s_auto, plt .cm .gray); plt .title ( 'Automatic thresholding')
plt .subplot (1,2,2);

11 plt .imshow( cells >s_otsu, plt .cm .gray); plt .title ( 'Otsu thresholding');
Figure 19.3: Automatic thresholding and thresholding by Otsu. Results are almost identical because threshold values are 105.3 and 105, respectively.

<figure><img src="image_478.png" /></figure>

### 19.6.3 k-means clustering

Different techniques can be found in the scikit documentation. A point cloud will first be generated, from 3 clustered cloud points. The objective is then to segment all the points into their original cluster.

#### Imports

<figure><img src="image_479.png" /></figure>

1. import numpy as np
2. import matplotlib.pyplot as plt

3. import time

5. from sklearn.cluster import KMeans
k-means clustering

<figure><img src="image_480.png" /></figure>

```python
n=3; # number of clusters

# k-means initialization
kmeans = KMeans(init='k-means++', n clusters=n, n_init=10)
t0 = time.time(); # computation time
tmeans.fitalteps); # kmeans segmentation
t_batch = time.time() - t0;

# retrieve results
kmeans labels = kmeans.labels_;
kmeans_cluster centers = kmeans.cluster centers_;

# plot
fig = plt .figure()
colors = ['#4EACC5', '#FF9C34', '#4E9A06']

# k-means
zip aggregates values two by two
for k, col in zip(range(n), colors):
    my_members = kmeans labels == k
    cluster center = kmeans_cluster centers[k]

# display points
plt .plot(pts[my_members,0], pts[my_members,1], 'w',
    markerfacecolor=col, marker='.')

# display centroid
plt .plot( cluster center [0], cluster center [1], 'o',
    markerfacecolor=col, markeredgecolor='k',
    markersize=6)

# title
fig .title('KMeans')
fig .show()
fig .savefig("kmeans.pdf");
```
<figure><img src="image_481.png" /></figure>

### 19.6.4 Color image segmentation

Three different colors can be observed in the image. The objective is to separate the 3 colors with the help of the K-means algorithm. Thus, the segmentation is performed in the RGB color space, and each pixel is represented by a point in this 3D space. Initialization steps are identical to previous code. The data is converted from a color image (of size (n, m, 3)) to a vector (of size (n × m, 3)), done by the reshape function of numpy.

<figure><img src="image_482.png" /></figure>

```python
# load color image
```

2 cells =imageio.openread('Tv16.png');
[nLines,nColsm channels] = cells.shape

4 # reshape data
    data = np.reshape cells, (nLines*nColsm channels);

6 kmeans.fit(data);
# convert result to an image

8 # as we got labels, we expand the dynamic (multiply by 70)
segmentation = 70*len(np.reshape(kmeans.labels_, (nLines, nColsm channels)));

10
fig = plt .figure();

12 plt ..imshow(segmentation, cmap = plt.cm .gray);
imageio .imwrite("segmentation_kmeans.png", segmentation);

#### 3D scatter plot

This is a method to display colors in the RGB cube. This method is really slow, depending on your GPU.
<figure><img src="image_483.png" /><figcaption>Figure 19.5: Segmentation result.</figcaption></figure>

<figure><img src="image_484.png" /></figure>

1. from mpl_toolkits.mplot3d import Axes3D # 3D scatter plot
   # plot

3. colors = ['#4EACC5', '#FF9C34', '#4E9A06']

5. fig = plt.figure ()
   ax = fig.add_subplot(111, projection='3d')

7.
# Plot scatter points

9. for k, col in zip(range(n), colors):
   my_members = kmeans_labels == k
   cluster_center = kmeans_cluster_centrerk]
11. ax. scatter(data[my_members, 0], data[my_members, 1],
               data[my_members, 2], c=col)
13. ax. scatter.cluster_center[0], cluster_center[1],
15. cluster_center[2], s=30, c=col)
<figure><img src="image_485.png" /></figure>

## 20 Segmentation by region growing

This tutorial proposes to program the region growing algorithm.

### 20.1 Introduction

The region growing segmentation method starts from a seed. The initial region first contains this seed and then grows according to

* a growth mechanism (in this tutorial, the $N_8$ will be considered)

* an homogeneity rule (predicate function)

The algorithm is simple and barely only needs a predicate function:

```
Data: I: image
Data: seed: starting point
Data: queue: queue of points to consider
Result: visited: boolean matrix, same size as I
begin
    queueenqueue(seed);
    while queue is not empty do
        p = queuedequeue();
        foreach neighbor of p do
            if not visited(p) and neighbor verifies predicate then
                queueenqueue(neighbor);
                visited(neighbor) = true;
            end
        end
    end
    return visited
end
```

**Algorithm 7:**

To get the coordinate of the mouse click, use the `matplotlib` connect utilities.
Define a function `def onpick`.