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

```python
# start by displaying a figure,
2 # ask for mouse input ( click )
    fig = plt.figure ();
4 ax = fig.add_subplot (211);
    ax. set_title ( 'Click on a point')
6 # load image
    img = misc.ascent ();
8 ax.imshow (img, picker = True, cmap = plt. gray());
# connect click on image to onpick function
10 fig canvas.mpl_connect ( 'button press event', onpick);
    plt. show ();
12 def onpick(event):
    """ connector """
```

## 20.2 Region growing implementation

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

The seed pixel is denoted s.

* Code the predicate function: for an image f and a pixel p, p is in the same segment as s implies $|f(s) - f(p)| \le T$.

* Code a function that performs region growing, from a starting pixel (seed).

* Try others predicate functions like:

- pixel p intensity is close to the region mean value, i.e.:

$$|I(p) - m_s| \le T$$

- Threshold value $T_i$ varies depending on the region $R_i$ and the intensity of the pixel $I(p)$. It can be chosen this way:

$$T_i = \left(1 - \frac{\sigma_i}{m_i}\right) \cdot T$$
<figure><img src="image_488.png" /></figure>

## 20.3. Python correction

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

### 20.3.1 imports

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

1. `import queue`
   `from scipy import misc`

3. `import matplotlib.plt`
   `import numpy as np`

### 20.3.2 Predicate

This function defines the aggregation condition.

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

```python
def predicate(image, i, j, seed):
    f=image[i, j];
    g=image[seed[0], seed[1]];
    return abs(f-g)<20
```

The following code is used to start the region growing from a pixel manually clicked on an image.

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

```python
# start of code
2 fig = plt.figure();
   ax = fig.add_subplot(211);
4 ax. set_title('Click on a point')

# load image
6 img = misc.ascent();
8 ax.imshow(img, picker=True, cmap= plt .gray());
10 fig .canvas.mpl_connect('button press event', onpick)
   plt.show();
```

And here comes the main function for region growing. The result is illustrated in Fig.20.1.
<figure><img src="image_493.png" /></figure>

```python
def onpick(event):
    """
    this functions gets the event 'click', and starts the region growing algorithm
    Notice that img is a global variable
    """
    print("x, y: ", event.xdata, event.ydata);
    # original pixel
    seed = np.array ([int (event.ydata), int (event.xdata)]);
    q = queue.Queue ();

myfunctions = [ predicate , predicate2 , predicate3 ];

for index, f in enumerate(myfunctions):
    """
    # initializes the queue
    q.put (seed);
    """
    #Visited matrix : result of segmentation
    #this matrix will contain 1 if in the region, -1 if visited but not in the
    #  →  region, 0 if not visited
    visited = np.zeros (img.shape)
    """
    #----------------------
    #Start of algorithm
    visited [seed [0], seed [1]] = 1;

while not q.empty ():
    """
    for i in range (max (0, p[0] - 1), min (img.shape [0], p[0] + 2)):
        for j in range (max (0, p[1] - 1), min (img.shape [1], p[1] + 2)):
            if not visited [i, j]:
                if f (img, i, j, seed, visited ):
                    visited [i, j] = 1;
                    q.put (np.array ([i, j]));
                else :
                    visited [i, j] = -1;
                """
    #visited matrix contains the segmentation result
    #display results
    ax = fig.add_subplot (12+index);
    ax.imshow (visited == 1);

    imageio imsave (f.__name__ +'_seg.python.png ', (visited ==1 ).astype ( 'int'))
    """
    fig .draw ();
    plt.show();
    return ;
```

Notice that values -1 of the visited matrix avoid testing multiple times the same pixel. In the predicate function, the visited matrix is used in case of adapting the
predicate to the current region. In the next case, the candidate pixel's graylevel is
compared to the mean gray value of the region. The results are illustrated Fig.20.1.

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

def predicate2 (image, i, j, seed, visited ):
2     f = int(image[i, j]);
     m = np.mean(image[visited==1]);
4     return abs (f - m) < 20

Another predicate function would be:

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

def predicate3 (image, i, j, seed, visited ):
2     f = int(image[i, j]);
     m = np.mean(image[visited==1]);
4     s = np.std(image[visited==1]);
     return abs (f - m) < 20 * (1 - s/m)
<figure><img src="image_496.png" /><figcaption>Figure 20.1: Region growing illustration for pixel (189,136). The segmentation result highly depend on the order used to populate the queue, on the predicate function and on the seed pixel.</figcaption></figure>
<figure><img src="image_497.png" /></figure>

## 21 Hough transform and line detection

This tutorial introduces the Hough transform. Line detection operators are implemented.

### 21.1 Introduction

This tutorial deals with line detection in an image. For a given point of coordinates $(x, y)$ in $\mathbb{R}^2$, there exists an infinite number of lines going by this point, with different angles $\theta$. These lines are represented by the following equation:

$$ \rho = x \cdot \cos(\theta) + y \cdot \sin(\theta). $$

Thus, for each point $(x, y)$ (Fig. 21.1a) corresponds a curve parametered by $[\theta, \rho]$, where $\theta \in [0; 2\pi]$ (Fig. 21.1b). The intersection of these curves represents a line (in this case, $y = x$).

<figure><img src="image_498.png" /><figcaption>Figure 21.1: Representation of the Hough transform.</figcaption></figure>

### 21.2 Algorithm

The (general and simple) method for line detection is then:
1. Compute contours detections (get a binary image BW).

2. Apply the Hough transform on the contours BW.

3. Detect the maxima of the Hough transform.

4. Get back in the Euclidean space and draw the lines on the image.

Results should look like in Fig. 21.2.

<figure><img src="image_499.png" /><figcaption>Figure 21.2: Lines detection via Hough transform.</figcaption></figure>

## 21.3 Hough transform

?

Code a function that will transform each point of a binary image into a
curve in the Hough space. For each curve, increment each pixel by one in
the Hough space.
**21.4** Maxima detection

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

Use or code a function to detect maxima (regional maxima). For each maximum, keep only one point.

**21.5** Display lines

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

For each maximum, display the corresponding line above the original image.
<figure><img src="image_502.png" /></figure>

## 21.6. Python correction

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

This correction makes use of the python modules numpy, opencv, skimage and matplotlib.

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

1. import numpy as np
import cv2

3. import matplotlib.pyplot as plt
from skimage feature import peak_local_max

### 21.6.1 Contours detection

The contours are detected using the Canny edge detection method. In this code, the method from OpenCV is employed, see Fig.21.3a.

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

```python
img = cv2.imread('TestPR46.png');
2. plt.figure()
    plt.imshow(img)
4
# perform contours detection
6. edges = cv2.Canny(img,100,200);
    plt.figure()
8. plt.imshow eds
```

### 21.6.2 Hough transform

Notice that OpenCV contains Hough fonctions: HoughLines and HoughLinesP. The result (sinogram) of the image is presented Fig.21.3b.
<figure><img src="image_506.png" /><figcaption>(a) Canny edge detection.</figcaption></figure>

<figure><img src="image_507.png" /><figcaption>(b) Representation of the sinogram and detection of the maxima in the Hough space.</figcaption></figure>

Figure 21.3: The algorithm of the Hough line detection consists in detecting the edges, then representing each pixel in the Hough space and finally detecting the maximal value in the sinogram.

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

## Hough transform

2 # size of image

X = img.shape[0];

4 Y = img.shape[1];

6 angular_sampled = 0.01; # angles in radians

8 # initialization of matrix H

rho_max = np.hypot(X, Y);

10 rho = np.arange(-rho_max, rho_max, 1);
theta = np.arange(0, np.pi, angular_sampled);

12 cosTheta = np.cos theta;

14 sinTheta = np.sin theta;

16 H = np(zeros([rho.size, theta.size]));

18 # Hough transform

19 for i in range(X):
    for j in range(Y):

22 R = np.random(R + rho.size /2) .astype(int);
H[R,rangeUCEsize)] += 1;

24 plt.imshow(H);
### 21.6.3 Maxima detection

The function `peak_local_max` from `skimage` is used to detect local maxima in the Hough transform. Matrix H is first smoothed with a Gaussian filter and represented in Fig.21.3b.

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

# Maxima detection

2. `G = cv2.GaussianBlur(H, (5,5), 5);`
   `maxima = peak_local_max(H, 5, threshold_abs=150, numPeaks=5);`

4. `plt.figure()`;
   `plt.imshow(G);`

6. # display maxima on Hough transform image G

8. `plt.scatter(maxima[:,1], maxima[:,0], c='r');`
   `plt.show();`

### 21.6.4 Resulting lines

The result is shown in Fig.21.4.

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

1. # display the results as lines in image
   for i_rho, i_theta in maxima:

3. `print rho[i_rho], theta[i_theta]`

a = np.cos theta[i_theta]

5. `b = np.sin theta[i_theta]`

y0 = a*rho[i_rho]

7. `x0 = b*rho[i_rho]`

y1 = int(y0 + 1000*(-b))

9. `x1 = int(x0 + 1000*(a))`

y2 = int(y0 - 1000*(-b))

11. `x2 = int(x0 - 1000*(a))`

13. `cv2 line(img,(x1,y1),(x2,y2),(0,0,255),2)`

15. # display in window
   `cv2.imshow('hough transform', img);`

17. # write resulting image
   `cv2-imwrite('cv_hough.png', img);`
<figure><img src="image_511.png" /><figcaption>Figure 21.4: Lines detected with the Hough transform.</figcaption></figure>

### 21.6.5 OpenCV builtin function

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

```python
import cv2

2 import numpy as np

4 # read image and convert it to gray
img = cv2.imread('TestPR46.png')

6 gray = cv2.CvtColor(img, cv2.COLOR_BGR2GRAY)
edges = cv2.Canny(grey, 100, 200, apertureSize = 3)

8 # threshold value for lines selection :
10 # lower value means more lines
threshold = 150;

12 # perform lines detection
14 lines = cv2.HoughLines eds, 1, np.pi/180, threshold

16 # display lines
for rho, theta in lines[0]:
    print rho, theta
    a = np.cos theta
    b = np.sin theta
    x0 = a*rho
    y0 = b*rho
    x1 = int(x0 + 1000*(-b))
    y1 = int(y0 + 1000*(a))
    x2 = int(x0 - 1000*(-b))
    y2 = int(y0 - 1000*(a))
    print x1, y1, x2, y2
28 cv2 line (img, (x1, y1), (x2, y2), (0,0,255), 2)
```
<figure><img src="image_513.png" /></figure>
<figure><img src="image_514.png" /></figure>

## 22 Active contours

This tutorials aims at introducing the active contours (a.k.a. snakes) method as originally presented in [18]. It is a segmentation method based on the optimization of a contour that will converge to a specific object.

### 22.1 Definition

A snake is a parametric curve $v(s)$ with $s \in [0; 1[$. The energy functional is represented by:

$$E_{snake} = \int_{0}^{1} E_{int}(v(s))ds + \int_{0}^{1} E_{ext}(v(s))ds.$$

The internal energy is detailed in Eq.22.1. The first order derivation constrains the length of the curve, the second order derivation constrains the curvature. The parameters $\alpha$ and $\beta$ a priori depend on $s$, but for simplicity, constant values will be taken.

$$E_{int}(v(s)) = \frac{1}{2} \left( \underbrace{\alpha(s)|v'(s)|^2}_{\text{length}} + \underbrace{\beta(s)|v''(s)|^2}_{\text{curvature}} \right) \quad (22.1)$$

The snake will be attracted by some edges points from the external energy (the image to be segmented). The, the energy functional will be in a local minimum: it is shown that the Euler-Lagrange equation is satisfied:

$$-\alpha v^{(2)} + \beta v^{(4)} + \nabla E_{ext} = 0$$

with $v^{(2)}$ and $v^{(4)}$ denoting the 2nd and 4th order derivatives. To solve this equation, the gradient descent method is employed: the snake is now transformed into a function of the position $s$ and the time $t$, with $F_{ext} = -\nabla E_{ext}$

$$\frac{\partial s}{\partial t} = \alpha v^{(2)} - \beta v^{(4)} + F_{ext}$$
Figure 22.1: Illustration of the active contours segmentation method. Two energies are at stake: internal energies depend only on the snake shape and control points, external energies are related to the image properties.

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

## 22.2 Numerical resolution

The spatial derivatives are approximated with the finite difference method, and the snake is now composed of $n$ points, $i \in [1; n]$:

$$x^{(2)} = x_{i-1} - 2x_i + x_{i+1}$$

$$x^{(4)} = x_{i-2} - 4x_{i-1} + 6x_i - 4x_{i+1} + x_{i+2}$$

The gradient descent method can be written as, with $\gamma$ being the time step and controls the convergence speed:

$$\frac{x_t - x_{t-1}}{\gamma} = A \cdot x_t + f_x(x_{t-1}, y_{t-1})$$

$$\frac{y_t - y_{t-1}}{\gamma} = A \cdot y_t + f_y(x_{t-1}, y_{t-1})$$

with

$$A = \begin{bmatrix} -2\alpha - 6\beta & \alpha + 4\beta & -\beta & 0 & \dots & 0 & -\beta & \alpha + 4\beta \\ \alpha + 4\beta & -2\alpha - 6\beta & \alpha + 4\beta & -\beta & \dots & 0 & -\beta & 0 \\ -\beta & \alpha + 4\beta & -2\alpha - 6\beta & \alpha + 4\beta & \dots & 0 & \dots & 0 \\ 0 & -\beta & \alpha + 4\beta & -2\alpha - 6\beta & \alpha + 4\beta & -\beta & 0 & \dots \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & \dots & -\beta & \alpha + 4\beta & -2\alpha - 6\beta & \alpha + 4\beta & -\beta & 0 \\ 0 & \dots & 0 & -\beta & \alpha + 4\beta & -2\alpha - 6\beta & \alpha + 4\beta & -\beta \\ -\beta & 0 & \dots & 0 & -\beta & \alpha + 4\beta & -2\alpha - 6\beta & \alpha + 4\beta \\ \alpha + 4\beta & -\beta & 0 & \dots & 0 & -\beta & \alpha + 4\beta & -2\alpha - 6\beta \end{bmatrix}$$

The forces $f_x$ and $f_y$ are the components of $F_{ext}$. For example for an image $I$, if $*$ denotes the convolution and $G_\sigma$ a gaussian kernel of standard deviation $\sigma$:

$$F_{\text{ext}} = -\nabla (\|\nabla (G_{\sigma} * I)\|)$$
?

* Generate a binary image (of size $n \times n$ pixels, with values 0 and 1) containing a disk (of a radius $R$).

* Generate the initial contour as an ellipse, with the same center as the disk.

* Generate the pentadiagonal matrix $A$. The different parameters can be $\alpha = 10^{-5}$ and $\beta = 0.05$.

* Program the iterations and visualize the results. For example, $\gamma = 200$ and 1000 iterations may give an idea of the parameters to use.
<figure><img src="image_516.png" /></figure>

step =.01;

2 x = n/2 + 400 * np.cos(np.arange(0,2* np.pi+step, step));

y = n/2 + 200 * np.sin(np.arange(0,2* np.pi+step, step));

The different parameters are defined by:

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

1 k=.1;
alpha = .0001;

3 beta = 10;
gamma= 100;

5 iterations = 1000;

<figure><img src="image_518.png" /><figcaption>Figure 22.3: Circle.</figcaption></figure>

### 22.3.3 Matrix construction

This is maybe the hardest part of this code, with the use of the `scipy sparse`. `diags` function.

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

N = x.size;

2 X = np.array([-beta, alpha+4*beta, -2*alpha-6*beta, alpha+4*beta, -beta, -beta,
  -> alpha+4*beta, -beta, alpha+4*beta])

A = scipy sparse diags(X, np.array([-2, -1, 0, 1, 2, N-2, N-1, -N+2, -N+1]),
  -> =(N,N)).toarray();

4 AA = np identity (N)-gamma*A;
invAA = np linalg inv(AA);
### 22.3.4 External forces

The external forces are computed with the following code. Notice that axis 0 corresponds to the vertical axis (y) and the axis 1 corresponds to the horizontal axis (x).

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

1. `G` = `scipy.ndimage filters.gaussianGradient_magnitude(I, 30);`

2. `Fy` = `scipy.ndimage.prewitt(G, axis=0);`

3. `Fx` = `scipy.ndimage.prewitt(G, axis=1);`

### 22.3.5 Display results

To enhance the role of the external forces, the arrows showing the force are displayed (quiver function, see Fig.22.4).

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

4. `step` = 20;
5. `hold` on
6. `plot([x;x(1)], [y; y(1)], 'g', 'linewidth', 3);`

7. `# display arrows for external forces`
   `step = 20;`

8. `subx = 1:step: size (I,1);`

9. `suby = 1:step: size (I,2);`

10. `[Xa, Ya] = meshgrid(subx, suby);`
11. `quiver(Xa, Ya, Fx(subx, suby), Fy(subx, suby));`

<figure><img src="image_522.png" /><figcaption>Figure 22.4: External forces that will be applied to the snake.</figcaption></figure>
### 22.3.6 Convergence algorithm

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

1. # interpolation methods to get values of the external forces at the coordinates of the snake

2. sx, sy = I.shape;
   ix = interpolate.interp2d(np.arange(n), np.arange(n), Fx);

5. iy = interpolate.interp2d(np.arange(n), np.arange(n), Fy);
   # loop for convergence of the snake

7. bar = progressbar PROGRESSBAR();
   for index in bar(range(iterations)):

9. fex = np.array([float(ix(XX,YY)) for XX,YY in zip(x,y)]);

fey = np.array([float(iy(XX,YY)) for XX,YY in zip(x,y)]);

11. #print(np.max(fex), np.min(fex))

x = np.matmul invAA, x+gamma*fex;

13. y = np.matmul invAA, y+gamma*fey;

The results are displayed in Fig.22.5.

Figure 22.5: Result of the snake converging toward the disk, after 1000 iterations with the proposed parameters.

<figure><img src="image_524.png" /></figure>
<figure><img src="image_525.png" /></figure>
<figure><img src="image_526.png" /></figure>

## 23 Watershed

This tutorial aims to study the watershed transform for image segmentation. In image processing, an image can be considered as a topographic surface. If we flood this surface from its minima and if we prevent the merging of the water coming from different sources, we partition the image into two different sets: the catchment basins separated by the watershed lines.

The different processes will be applied on the following images:

<figure><img src="image_527.png" /><figcaption>(a) circles</figcaption></figure>

<figure><img src="image_528.png" /><figcaption>(b) gel</figcaption></figure>

### 23.1

#### Watershed and distance maps

The objective is to individualize the disks by disconnecting them with the distance map.

?

* Calculate the distance map of the image 'circles'.

* Take the complementary of this distance map and visualize its minima.

* Calculate the watershed transform of the inverted distance map.

* Subtract the watershed lines to the original image.

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

Look at the documentation of `scipy.ndimage.morphology` for Euclidean or Chamfer distance transform. Reconstruction operator can be found in `skimage.morphology`.
## 23.2 Watershed and image gradients

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

* Calculate the Sobel gradient of the image 'gel'.

* Visualize the minima of the image gradient.

* Apply the watershed transform on the gradient image.

The watershed transform, applied in a direct way, leads to an over-segmentation. To overcome this limitation, the watershed operator can be applied on a filtered image.

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

* Smooth the original image with a low pass filter. To stay in the mathematical morphology field, you can use an alternate morphological filter (opening followed by closing). A Gaussian filter is also a good solution.

* Calculate the gradient operator on the filtered image.

* Calculate the corresponding watershed.

## 23.3 Constrained watershed by markers

In order to individualize the image spots, we have to determine the internal and external markers for the constrained watershed.

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

* Calculate the gradient (Sobel) of the filtered image.

* Calculate the internal markers (minima of the filtered image) and external (watershed of the filtered image) as minima of the gradient image.

* Calculate the corresponding watershed.

* Superimpose the watershed lines of the resulting segmentation to the original image.
<figure><img src="image_533.png" /></figure>

## 23.4. Python correction

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

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

import numpy as np

2 from scipy import misc

import matplotlib .pyplot as plt

4 from scipy import ndimage as ndi

from scapy .ndimage.morphology import distance_transform

6 from scapy .ndimage.morphology import distance_transform_cdt

from skimage import morphology

### 23.4.1 Watershed and distance map

#### Regional maxima

The following code is an implementation of the regional maxima (from mathematical morphology definition). It deals with the case of float images as input.

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

1 def rmax(I):

###

3 Own version of regional maximum

This avoids plateaus problems of peak_local_max

5 I: original image, int values

returns: binary array, with True for the maxima

###

7 I = I .astype( 'float ');

9 I = I / np .max(I) * 2**31;

I = I .astype( 'int32 ');

11 h = 1;

rec = morphology.reconstruction(I, I+h);

13 maxima = I + h - rec;

return maxima>0

#### Distance map followed by a watershed

This method is a classical way of performing the separation of some objects by proximity or influence zones. It is illustrated in Fig.23.1.
<figure><img src="image_537.png" /></figure>

A = imageio.imread(' circles . tif ');

2 # chamfer distance gives integer distances

dm = distance_transform_edt(A);

4 # regional maxima

local_maxi = rmax(dm);

The local maxima will be used as a marker for the watershed operation, in order
to perform the separation of the grains, see Fig.23.1.

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

1 # watershed segmentation for separating the circles

markers = ndi. label (local_maxi, np. ones ((3, 3))) [0]

3 W = morphology.watershed(-dm, markers, watershed_line=True)

# separation of the grains

5 B = A & W==0;

separation = ndi. label (B, np. ones ((3,3) )) [0];

<figure><img src="image_539.png" /><figcaption>Figure 23.1: Steps of the separation of the grains.</figcaption></figure>

23.4.2 Watershed and image gradients

The gradient image amplifies the noise. Thus, the watershed operator directly ap-
plied to the gradient of the image produces an over-segmented image (see Fig.23.2).
<figure><img src="image_540.png" /></figure>

```python
def sobel_mag(im):
    """
    Returns Sobel gradient magnitude
    im: image array of type float
    returns: magnitude of gradient (L2 norm)
    """
    dx = ndi.sobel(im, axis=1)  # horizontal derivative
    dy = ndi.sobel(im, axis=0)  # vertical derivative
    mag = np.hypot(dx, dy)      # magnitude
    return mag;
```

Figure 23.2: Performing the watershed on the gradient image is usually not a good idea.

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

In fact, this method produces as many segments as there are minima in the gradient image Fig.23.3.

<figure><img src="image_542.png" /><figcaption>Figure 23.3: Local minima of the gradient image.</figcaption></figure>

Solution: filtering the image

Before evaluating the gradient, the image is filtered. The number of minima is lower and this leads to a less over-segmented image (Fig.23.4).
<figure><img src="image_543.png" /></figure>

```python
SE = morphology.disk(2);
2 O = morphology.opening(gel, selem=SE);
F = morphology.closing(O, selem=SE).astype('float');
4 g = sobel_mag(F).astype('float');
```

Figure 23.4: Even if the gradient is performed on the filtered image, there is still a high over-segmentation.

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

### 23.4.3 Watershed constrained by markers

The watershed can be constrained by markers: the markers can provide the correct number of regions. This method imposes both the background (external markers) and the objects (internal markers). The results are illustrated in Fig.23.5. The ultimate erosion of the internal markers is used to deconnect these markers from the external markers.

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

```python
local_maxi = rmax(255-F);
2 markers = ndi. label(local_maxi, np.ones((3, 3))) [0]
W = morphology.watershed(F, markers, watershed_line=True)
4
markers2 = local_maxi | (W==0);
6 M = ndi. label/markers2, np.ones((3, 3))) [0]
segmentation = morphology.watershed(g, M, watershed_line=True);
8
gel [segmentation==0] = 255;
10 plt.imshow(gel);
12 plt.show();
```
<figure><img src="image_546.png" /><figcaption>Figure 23.5: Watershed segmentation by markers.</figcaption></figure>
<figure><img src="image_547.png" /></figure>
<figure><img src="image_548.png" /></figure>

## 24 Segmentation of follicles

This practical work aims to investigate image segmentation with a direct application to ovarian follicles. The overall objective is to extract and quantify the granulosa cells and the vascularization of each follicle included in an ewe's ovary.

The image to be processed is a 2D histological image of an ewe's ovary acquired by optical microscopy in Fig.24.1. The presented image contains one entire follicle (the white region and its neighborhood) and a part of a second one (right-upper corner). The follicle is composed of different parts shown in: antrum, granulosa cells and vascularization. The theca is the ring region around the antrum where the follicle is vascularized.

<figure><img src="image_549.png" /><figcaption>Figure 24.1: Different parts of the follicles, to be segmented.</figcaption></figure>

### 24.1 Vascularization

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

* Load and visualize the image.

* Extract the antrum of the follicle.

* Extract the vascularization (inside a ring around the antrum).
## 24.2 Granulosa cells

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

Which kind of processing could be suitable for extracting the granulosa cells?

## 24.3 Quantification

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

Provide some geometrical measurements of the different entities of the follicle (antrum, vascularization, granulosa cells).
<figure><img src="image_553.png" /></figure>

## 24.4. Python correction

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

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

import numpy as np

2 import matplotlib.pylplot as plt
from scipy import misc, ndimage

4 from skimage import morphology

### 24.4.1 Vascularization

#### Antrum segmentation

The first step consists in the segmentation of the antrum by thresholding the blue component. Some post-processes are used to remove artifacts such as holes. This is illustrated in Fig. 24.2.

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

A = imageio.read('follicle.png');

2 plt.imshow(A);
plt.show();

4 ## Antrum
# segmentation by mathematical morphology

6 # manual selection of the antrum
B = A [:, : ,2];

8 antrum = B > 220;
L = morphology.label(antrum, connectivity =2);

10 antrum = L == L [300,300];
antrum = ndimage.morphology-binary_fill_holes(antrum);

12 plt.imshow(antrum);
plt.title('Antrum')

14 plt.show();

#### Theca segmentation

The second step provides the segmentation of the theca which is extracted as a spatial region (corona) adjacent to and outside the antrum. The width of the corona is selected by the user (expert).
<figure><img src="image_557.png" /></figure>

```matlab
se40 = morphology.disk(40);
2 theca = morphology-binary_dilation(antrum, selem = se40);
theca = theca - antrum;
4 plt.imshow(theca); plt.title('Theca')
plt.show()
```

## Vascularization segmentation

The final step extracts the vascularization of the considered follicle. Pixels belonging to the vascularization are considered to have a low blue component and they should also be included in the antrum of the follicle.

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

```matlab
1 vascularization = B < 140;
    vascularization = vascularization * theca;
3 plt.imshow(vascularization);
    plt.title('vascularization')
5 plt.show();
```
## Results

<figure><img src="image_559.png" /><figcaption>Figure 24.2: Extraction of the different components of the follicle.</figcaption></figure>

### 24.4.2 Granulosa cells

#### First solution

The granulosa cells have a low contrast, so it is difficult to use thresholding techniques. But we know they are localized between the antrum and the vascularization. Nevertheless the vascularization is not a closed region outside the antrum. Therefore, the proposed solution consists in first trying to close the vascularization region and then taking the corona between this region and the antrum. The result is shown in Fig. 24.3.
<figure><img src="image_560.png" /></figure>

1. `se10` = morphology.disk(10);
   dil = 1-morphology-binary_closing( vascularization , se10);

3. L = morphology.label(dil, connectivity =1);
   dil = L == L [300,300];

5. `granulosa` = dil - antrum;
   `plt.imshow(granulosa);`

7. `plt.title('granulosa')`
   `plt.show()`

<figure><img src="image_561.png" /><figcaption>Figure 24.3: Extraction of the granulosa cells of the follicle.</figcaption></figure>

Second solution

This first solution is not really robust. The closing of the vascularization region is
not really accurate. A more robust solution consists in using deformable models
to get the corona between the vascularization and the antrum. But this kind of
method is out of the scope of this tutorial.

**24.4.3 Quantification**

The quantification is easy to process. In addition, we can represent the different
extracted components of the follicle in false colors.
<figure><img src="image_562.png" /></figure>

```python
result = antrum + 2*granulosa + 3* vascularization;

2 plt.imshow(result, cmap='jet');
plt.show();

4 plt.imshow("result.png", result, cmap='jet');

6 follicle = antrum + theca;
q_vascularization = np.sum(vascularization) / np.sum(follicle);

8 print('vascularization: ', q_vascularization)
q_granulosa = np.sum(granulosa)/np.sum(follicle);

10 print('Granulosa: ', q_granulosa);
```

This quantification gives:

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

```
vascularization : 0.0428495233516
2 Granulosa: 0.0913707490403
```
<figure><img src="image_564.png" /></figure>
<figure><img src="image_565.png" /></figure>

## 25 Image Registration

This tutorial aims to implement the Iterative Closest Point (ICP) method for image registration. More specifically, we are going to estimate a rigid transformation (translation + rotation without scaling) between two images.

The different processes will be applied on T1-MR images of the brain in Fig.25.1.

<figure><img src="image_566.png" /><figcaption>Figure 25.1: Original images.</figcaption></figure>

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

### 25.1 Transformation estimation

A classical method in image registration first consists in identifying and matching some characteristic points by pairs. Thereafter, the transformation is estimated from this list of pairs (displacement vectors).

#### 25.1.1 Preliminaries

Pairs of points are first manually selected.

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

* Read and visualize the two MR images 'brain1' and 'brain2' (moving and source images).

* Manually select a list of corresponding points.
There is no built-in function for manual selection of pairs of points. You
can use the following points or code your own function, see tutorial 20

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

# define control points

2 A_points = np.array ([[136, 100], [127, 153], [96, 156], [87, 
  ↪ → 99]]);

B_points = np.array ([[144, 99] , [109, 140], [79, 128], [100, 
  ↪ → 74]]);

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

25.1.2 Rigid transformation

With this list of pairs $(p_i, q_i)_i$, this tutorial proposes to estimate a rigid transformation between these points. It is composed of a rotation and a translation (and not scaling). For doing that, we make a Least Squares (LS) optimization, which is defined as follows. The parameters of the rotation $R$ and the translation $t$ minimize the following criterion:

$$
C(R, t) = \sum_{i} \|q_{i} - R.p_{i} - t\|^{2}
$$

**Calculation of the translation :**

The optimal translation is characterized by a null derivative of the criterion:

$$
\frac{\partial C}{\partial t} = -2 \sum_{i} (q_{i} - R.p_{i} - t)^{\mathrm{T}} = 0 \Leftrightarrow \sum_{i} q_{i} - R. \left( \sum_{i} p_{i} \right) = N.t
$$

where $N$ denotes the number of matching pairs. By denoting $\bar{p} = \frac{1}{N} \sum_i p_i$ and $\bar{q} = \frac{1}{N} \sum_i q_i$ the barycenters of the point sets and by changing the geometrical referential: $p'_i = p_i - \bar{p}$ et $q'_i = q_i - \bar{q}$, the criterion can be written as:

$$
C'(R) = \sum_{i} \|q'_{i} - R.p'_{i}\|^{2}
$$

The estimated rotation $\hat{R}$ will provide the expected translation:

$$
\hat{t} = \bar{q} - \hat{R} \cdot \bar{p} \tag{25.1}
$$

**calculation of the rotation by the SVD method:**

We will use the following theorem: Let $U.D.V^T = K$ a singular decomposition of the correlation matrix $K = q'^T.p'$, for which the singular values are sorted in the increasing order. The minimum of the criterion: $C(R) = \sum_i ||q'_i - R.p'_i||^2$ is reached by the matrix :

$$
\hat{R} = U.S.V^T \tag{25.2}
$$

with $S = \text{DIAG}(1, \dots, 1, \det(U), \det(V))$.
?

* Code a function that takes as parameters the pairs of points, and returns the rigid transformation elements of rotation $\hat{R}$ and translation $\hat{t}$ as defined in Eqs.25.1 and 25.2.

* Apply the transformation to the moving image.

* Visualize the resulting registered image.

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

* See svd function from scipy.linalg or numpy.linalg.

* See cv2.warpAffine and cv2 transform for transformation application

## 25.2

# ICP-based registration

When the points are not correctly paired, it is first necessary to reorder them before estimating the transformation. In this way, the ICP (iterative Closest Points) consists in an iterative process of three steps: finding the correspondence between points, estimating the transformation and applying it. The process should converge to the well registered image.

?

To simulate the mixing of the points, randomly shuffle them selected on the first image and perform the registration with the previous method.

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

See np.random.permutation.

?

1. From the list of the characteristic points $p$ of the image 'brain1', find the nearest neighbors $q$ in the image 'brain2'. Be careful to the order of the input arguments.

2. Estimate the transformation $T$ by using the LS minimization coded previously.

3. Apply this transformation to the points $p$, find again the correspondence between these resulting points $T(p)$ and $q$ and estimate a new transformation. Repeat the process until convergence.
4. Visualize the resulting registered image.

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

Look at `scipy.spatial.cKDTree` for nearest neighbor search.

## 25.3 Automatic control points detection

The manual selection of points may be fastidious. Two simple automatic methods for detecting control points are the so-called Harris corners detection or Shi-Tomasi corner detector.

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

The function `goodFeaturesToTrack` returns the corners from the Shi-Tomasi method.

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

Replace the manual selection in the two previous parts of this tutorial.

Notice that the automatic points detection does not ensure to give the same order in the points of the two images. More generally, other salient points detectors do not give the same number of points, and thus the algorithms have to remove outliers (non matching points). This case is not taken into account in this tutorial.
<figure><img src="image_576.png" /></figure>

## 25.4. Python correction

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

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

fromigliomisc

2import matplotlib .pyplot as plt

import numpy as np

4import cv2

fromigliomiscimport cKDTree

### 25.4.1 Transformation estimation based on corresponding points

Image visualization

The following code is used to display the images (see Fig.25.2).

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

1 # Read images and display
A=imageio.openread("brain1.bmp")

3B=imageio.openread("brain2.bmp")
 plt .imshow(A,cmap='gray');

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

7 plt.show();

9 # define control points
A_points = np.array([[136, 100], [127, 153], [96, 156], [87, 99]]);

11 B_points = np.array([[144, 99], [109, 140], [79, 128], [100, 74]]);

<figure><img src="image_580.png" /><figcaption>Figure 25.2: Initial images.</figcaption></figure>
If you want to display and save the fusion of both images, you can use this function:

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

```python
def superimpose(G1, G2, filename=None):
    """
    superimpose 2 images, supposing they are grayscale images and of same shape
    """
    r, c=G1.shape;
    S = np.zeros((r, c, 3));
    S [:, :0] = np maximum(G1-G2, 0)+G1;
    S [:, :1] = np.max imum(G2-G1, 0)+G2;
    S [:, :2] = (G1+G2) / 2;
    S = 255 * S / np.max(S);
    S = S .astype('uint8');
    plt .show();
    plt .show()
    if filename!=None:
        cv2 .imwrite(filename, S);
    return S
```

## Manual selection of corresponding points

With openCV, there is not built-in function to make a manual selection of pairs of control points. The following code uses a global variable I in order to manage the display of the points, which are finally stored into the pts variable. First, the callback function onMouse is defined to handle mouse event.

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

```python
pts = [];

2 def onMouse(event, x, y, flags, param):
    """
callback method for detecting click on image
It draws a circle on the global variable image I
    """
    global pts, I;
    if event == cv2.EVENT_LUP:
        pts.append((x, y));
    cv2 circle (I, (x, y), 2, (0,0,255), -1)
```

The function cpselect allows the selection of multiple points.
<figure><img src="image_583.png" /></figure>

```python
def cpselect():
    """
    method for manually selecting the control points
    It waits until 'q' key is pressed.
    """
    cv2.namedWindow("image")
    cv2.setMouseCallback("image", onmouse)
    print("press 'q' when finished")
    # keep looping until the 'q' key is pressed
    while True:
        """
        # display the image and wait for a keypress
        cv2.imshow("image", I)
        key = cv2.waitKey(1) & 0xFF
        """
        # if the 'c' key is pressed, break from the loop
        if key == ord("q"):
            break
        """
        # close all open windows
        cv2 destoyAllWindows()
        return pts;
```

## Transformation estimation

The rigid transformation is estimated from the corresponding points by the following function:
<figure><img src="image_584.png" /></figure>

```python
def rigid Registration (data1, data2):
    """
    Rigid transformation estimation between n pairs of points
    This function returns a rotation R and a translation t
    data1 : array of size nx2
    data2 : array of size nx2
    returns transformation matrix T of size 2x3
    """
    data1 = np.array(data1);
    data2 = np.array(data2);

    # computes barycenters, and recenters the points
    m1 = np.mean(data1,0);
    m2 = np.mean(data2,0);
    data1_inv_shifted = data1-m1;
    data2_inv_shifted = data2-m2;

    # Evaluates SVD
    K = np.matmul(npurlap( data2_inv_shifted ), data1_inv_shifted );
    U,S,V = np. linalg.svd(K);

    # computes Rotation
    S = np. eye(2);
    S [1,1] = np. linalg. det(U)*np. linalg. det(V);
    R = np. matmul(U,S);
    R = np. matmul(R, np. transpose(V));

    # computes Translation
    t = m2-np. matmul(R, m1);

    T = np. zeros((2,3));
    T [0:2,0:2] = R;
    T [0:2,2] = t;
    return T;
```

Then, you can apply this function to the manually selected points.
<figure><img src="image_585.png" /></figure>

```matlab
% 1st case, rigid registration, with pairs of points in the correct order
2 T = rigid Registration (A_points, B_points);
4 % Apply transformation on control points and display the results
    data_dest = applyTransform(A_points, T);
6 I = B.copy();
    for pb, pa in zip(data_dest, B_points):
        cv2 circle (I, totuple (pa), 1, (255,0,0), -1)
        cv2 circle (I, totuple (pb), 1, (0,0,255), -1)
10 plt ..imshow(I);
12 plt .show();
12 % Apply transformation on image
    rows, cols = B.shape;
14 dst = cv2.warpAffine(A,T, (cols, rows));
    plt ..imshow(dst);
16 superimpose(dst, B, "rigid_manual.png");
```

<figure><img src="image_586.png" /><figcaption>Figure 25.3: Result of the registration for the manually selected control points.</figcaption></figure>

The result is good, because the manual selection of the points is good (the points are given in the correct order for both images).

## 25.4.2 ICP registration

### Random permutation of points

The following code randomly shuffles the points of the first vector. The result is of course a wrongly registered image, see Fig.25.4.
<figure><img src="image_587.png" /></figure>

# random permutation of the points

2 p = np.random.permutation(np.arange(4));

A_points = A_points[p];

4 T = rigidRegistration (A_points, B_points);

Figure 25.4: Result of the registration for when the control points are not found in
the same order.

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

ICP

The previous operations simulates a general non-manual selection of the control
points: there is not reason to finding the points by matching pairs. Thus, a reordering
is necessary. This propositions implies that the number of points is exactly the same
in order to perform the registration process, and that these are matching points
(every point has a matching point in the other image). The ICP method (see code
for icp_transform) reorders the points via a nearest neighbor rule.
<figure><img src="image_589.png" /></figure>

```python
def icp_transform(dataA, dataB):
    """
    Find a transform between A and B points
    with an ICP (Iterative Closest Point) method.
    dataA and dataB are of size nx2, with the same number of points n
    returns the transformation matrix of shape 2x3
    """
    data2A = dataA.copy();
    data2B = np.zeros(data2A.shape);
    T = np zeros ((2,3) );
    T [0:2,0:2] = np .eye(2) ;
    T [0:2,2] = 0;

    nb_loops=5;
    tree = cKDTree( dataB );

    for loop in range nb_loops:
        """
        # search for closest points and reorganise array of points accordingly
        d, inds = tree .query( data2A);

        data2B = dataB[inds ,:] ;
        """
        # find rigid registration with reordered points
        t_loop = rigid _registration (data2A, data2B);

        T = composeTransform(t_loop, T);
        """
        # evaluates transform on control points, to make the iteration
        data2A = applyTransform(dataA, T);

    return T;
```

Automatic extraction of corner points

Generally, the points are automatically detected, and thus, there is no warranty
that they are found in the same order, nor that each pair of point correspond to
matching points (some points -called outliers- need to be eliminated to compute
the correct transformation). In this tutorial, we do not address the problem of
outliers. Please notice that with these parameters and images, by chance, the same
points are detected in the correct order.
Figure 25.5: Result of the registration for when the control points are not found in the same order. The ICP algorithm reorders the points and gives a good result.

(a) Random shuffle of points and direct rigid transformation estimation.

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

(b) ICP registration on the same points.

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

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

1. Automatic extraction of corner points

# Harris corners detection or Shi and Tomasi

3. This method does not ensure the correct order in the points (A_points and B_points)

5. `nb_points = 4; # number of points to detect`
   corners = cv2.goodFeaturesToTrack(A,nb_points,0.01,10,blockSize=3,
                                     → useHarrisDetector=False)

7. `corners = np.int0(corners)`

a, b, c = corners.shape;

9. `A_points = np.reshape(corners, (a,c));`

11. `corners = cv2.goodFeaturesToTrack(B,nb_points,0.01,10,blockSize=3,
                                     → useHarrisDetector=False)`

corners = np.int0(corners)

a, b, c = corners.shape;

B_points = np.reshape(corners, (a,c));

The Fig.25.6 displays the results with the automatic detection of corners (in this case, the method from Shi and Tomasi). Notice that by chance, the points in both images match. In other cases, one would have to remove outliers.
<figure><img src="image_593.png" /><figcaption>Figure 25.6: Shi and Tomasi corners detection. By chance, the points correspond and the ICP method gives a correct result.</figcaption></figure>
<figure><img src="image_594.png" /></figure>
# Part IV Stochastic Analysis
<figure><img src="image_595.png" /></figure>
<figure><img src="image_596.png" /></figure>

# 26 Stochastic Geometry / Spatial Processes

This tutorial aims to simulate different spatial point processes.

## 26.1 Poisson processes

This process simulates a conditional set of point$_{nb}$ points in a spatial domain $D$ defined by the values $x_{min}, x_{max}, y_{min}, y_{max}$. In order to simulate a non conditional Poisson point process of intensity $\lambda$ within a domain $S$, it is necessary to generate a random number of points according to a Poisson law with the parameter $\lambda S$ : point$_{nb}$ = poisson($\lambda S$) (i.e. the number of points is a random variable following a Poisson Law).

The coordinates of each point follow a uniform distribution.

<figure><img src="image_597.png" /><figcaption>Figure 26.1: Conditional Poisson point process, with 200 points.</figcaption></figure>

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

1. Simulate a conditional Poisson point process on a spatial domain $D$ with a fixed number of points (see Fig. 26.1).

2. Simulate a Gaussian distribution of points around a given center.
<figure><img src="image_599.png" /></figure>

Import spipy.stats to use the function poisson and np.random for more classical stochastic functions.

## 26.2 Neyman-Scott processes

This process simulates aggregated sets of points within a spatial domain $D$ defined by the values $x_{min}, x_{max}, y_{min}, y_{max}$. For each aggregate ($n_{par}$), we first generate the random position of the 'parent' point. Then, the 'children' points are simulated in a neighborhood (within a square box of size $r_{child}^2$) of the 'parent' point. The number of points for each aggregate is either fixed or randomized according to a Poisson law of parameter $n_{child}$ (see Fig. 26.2).

<figure><img src="image_600.png" /><figcaption>Figure 26.2: Neyman-Scott point process with $h_{child} = 10$ and $n_{par} = 10$</figcaption></figure>

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

1. Simulate a process of 3 aggregates with 10 points.

2. Simulate a process of 10 aggregates with 5 points.

## 26.3 Gibbs processes

The idea of Gibbs processes is to spatially distribute the points according to some laws of interactions (attraction or repulsion) within a variable range.
Such a process can be defined by a cost function $f(d)$ that represents the cost associated to the presence of 2 separated points by a distance $d$ (see Fig. 26.3). For a fixed value $r$, if $f(d)$ is negative, there is a high probability to find 2 points at a distance $d$ (attraction). Conversely, if $f(d)$ is positive, there is a weak probability to find 2 points at a distance $d$ (repulsion).

* Code such a function, with prototype function e=f(d) or def energy(d):, where

$$f(d) = \begin{cases} 50 & \text{if } 0 < d \le 15 \\ -10 & \text{if } 15 < d \le 30 \\ 0 & \text{otherwise} \end{cases}$$

<figure><img src="image_602.png" /><figcaption>Figure 26.3: Energy function $f$.</figcaption></figure>

The generation process reorganizes an initial Poisson point process within the spatial domain according to a specific cost piecewise constant function. The reorganization consists in a loop of $nb_{iter}$ iterations. For each iteration, we calculate the total energy:

$$e = \sum_{(i,j), i \neq j} f(dist(x_i, x_j)) \qquad (26.1)$$

The objective is to reduce this energy iteratively. At each iteration step, we try to replace a point by four (for example) other random points and we calculate the energy for each new configuration. The initial point is either preserved or replaced (by one of the four points) according to the configuration of minimal energy (see Fig.26.4).

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

1. Code a function with the following prototype (the function energy-Function is previously noticed $f$). It computes the energy between all the points present in an array of coordinates $[x, y]$ (the points that
<figure><img src="image_604.png" /><figcaption>Figure 26.4: Illustration of iterative construction of the Gibbs point process.</figcaption></figure>

do not move) and point $[x_k, y_k]$ (the new point). The energyFunction converts a distance to an 'energy', reflecting attractivity or repulsivity.

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

def energy(P, eFunction=exampleEnergyFunction):

2. Simulate a regular point process by choosing a specific energy function.

3. Change the cost function and simulate a few aggregated point processes (see Fig. 26.5).

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

Informations

Use the function `pdist` from `scipy.spatial.distance`, which computes pairwise distances of all points.
<figure><img src="image_607.png" /><figcaption>Figure 26.5: Gibbs aggregated point process and its energy function.</figcaption></figure>

## 26.4 Ripley function

The Ripley function $K(r)$ characterizes the spatial distribution of the points. For a Poisson process of density $\lambda$, $\lambda K(r)$ is equal to the mean value of the number of neighbors at a distance lower than $r$ to any point. In the case where the process is not known ($\lambda?$), the Ripley function has to be estimated (approximated) with the unique known realization. It is the first estimator of $K(r)$:

$$K(r) = \frac{1}{\lambda N} \sum_{i=1}^{N} \sum_{i \neq j} k_{ij} \qquad (26.2)$$

where $N$ is the number of points in the studied domain $D$, $\lambda = N/D$ is the estimator of the density of the process and $k_{ij}$ takes the value 1 if the distance between the point $i$ and the point $j$ is lower than $r$, and 0 in the other case.

We denote:

$$L(r) = \sqrt{K(r)/\pi} \qquad (26.3)$$

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

1. Code a function to compute $K$ (function K=ripley(points, box, r)), with points being the considered points, box the simulation domain, and $r$ the distance (or an array of all distances).

2. Calculate the Ripley function for an aggregated point process, a Poisson point process and a regular point process.

3. Display the value $L(r) - r$ as in Fig 26.6.
<figure><img src="image_609.png" /><figcaption>Figure 26.6: Ripley’s function $L(r) - r$.</figcaption></figure>

## 26.5 Marked point processes

To simulate a complex point process, it can be useful to associate several random variables (marks) for each point.

?

1. Simulate a disk process, where the points (disk centers) are defined according to a Poisson process and the radii are selected with a Gaussian law.

2. Add a second mark for allocating a color to each disk (uniform law).

An example result is shown in Fig. 26.7

<figure><img src="image_610.png" /><figcaption>Figure 26.7: A Poisson point process is used to generate the center of the circles. The radii are chosen according to a Gaussian law, and the color according to a uniform law.</figcaption></figure>
<figure><img src="image_611.png" /></figure>

## 26.6. Python correction

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

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

import numpy as np;

2 import matplotlib .pyplot as plt;

from scipy .spatial .distance import pdist;

4 from scipy .stats import poisson

### 26.6.1 Conditional Poisson Process

The conditional Poisson Point Process uses a given number of point, contrary to the Poisson Process where the number of points follows a Poisson law. The result is illustrated in Fig.26.8.

<figure><img src="image_614.png" /><figcaption>Figure 26.8: Conditional Poisson point process, with N=100 points.</figcaption></figure>

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

def cond_Poisson (nb_points, xmin, xmax, ymin, ymax):

2 # Conditional Poisson Point Process

# uniform distribution

4 # nb_points: number of points

# xmin, xmax, ymin, ymax: defined the domain (window)

6 x =xmin+(xmax-xmin)*np.random .rand (nb_points)

y =ymin+(ymax-ymin)*np.random .rand (nb_points);

8 return x,y

10 def test_ppp () :

# testing function

12 x,y =cond_Poisson (100, 0, 100, 0, 100);

plt .plot (x,y,'+');
## 26.6.2 Normal distribution

The normal distribution is illustrated in Fig.26.9. Each marginal distribution (distribution on each axis) follows the normal distribution.

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

1. def normal_distribution (nb_points, mu, sigma): # Normal distribution centered around the point mu with stdev sigma

3. x = mu[0] + sigma[0]*np.random.random nb_points;

y = mu[1] + sigma[1]*np.random.random nb_points;

5. return x, y;

<figure><img src="image_617.png" /><figcaption>Figure 26.9: Normal distribution of points around (0, 0), with $\sigma = (10, 20)$.</figcaption></figure>

## 26.6.3 Neyman-Scott Process

Neymann-Scott point process is an aggregated Poisson point process. It is illustrated in Fig.26.10. It consists on a “sub” point processes generated at locations corresponding to a point process.
<figure><img src="image_618.png" /></figure>

```R
def neyman_scott(nRoot, xmin, xmax, ymin, ymax, lambdaS, rSon):
    # Neyman-Scott process simulation
    # nRoot: number of aggregates
    # xmin, xmax, ymin, ymax: domain
    # lambdaS: number of points. lambda is a density, S is the spatial domain
    # rSon: radius around aggregate (points are distributed in a square)

    # number of sons
    nSons = poisson.rvs(lambdaS, size=nRoot);
    # results
    x = [];
    y = [];

    # father points coordinates
    xf, yf = cond_Poisson(nRoot, xmin, xmax, ymin, ymax);
    for i in range(nRoot):
        # loop over all aggregates
        xs, ys = cond_Poisson(nSons[i], xf[i]-rSon, xf[i]+rSon, yf[i]-rSon, yf[i]+
                                i % rSon);
        x = np.concatenate((x,xs), axis=0);
        y = np.concatenate((y,ys), axis=0);

    return x,y
```

<figure><img src="image_619.png" /><figcaption>Figure 26.10: Neyman-Scott point process, with lambdaS=10 and rSon=10.</figcaption></figure>

### 26.6.4 Gibbs Point Process

Gibbs point process allows attraction and repulsion at different distances. It is illustrated in Fig.26.11. The attraction/repulsion law is given by the following code for regular or aggregated point process.
<figure><img src="image_620.png" /></figure>

```python
def gibbs nb_points, xmin, xmax, ymin, ymax, nbiter, eFunction=
    ↪→ exampleEnergyFunction):
        """
        Gibbs point process
        xmin, xmax, ymin, ymax represents the spatial window
        nb_points: number of generated points
        nbiter: number of iterations
        returns (x,y) coordinates of the points
        """
        # start with a Poisson Point Process
        x,y = spat_pp.cond_Poisson nb_points, xmin, xmax, ymin, ymax;
        nb_move = 0;
        e_prev = energy(np.vstack((x,y)), eFunction);
        print(" initial energy: {0: f}".format(e_prev));

        for i in range (nbiter):
            # choose a random point
            j = np.randomondeint(0, nb_points);
            x2 = np.delete(x, j);
            y2 = np.delete(y, j);

            P = np.vstack((x2, y2));
            e1 = energyFromPoint([x[j], y[j]], P, eFunction);

            for m in range(10):
                zm, ym = spat_pp.cond_Poisson(1, xmin, xmax, ymin, ymax);
                e2 = energyFromPoint([zm, ym], P, eFunction);
                if e2<e1:
                    nb_move+=1;
                    x[j]=zm;
                    y[j]=ym;
                    e1=e2;

        print("Number of moves: " + str (nb_move));
        print("Final energy: " + str (e1));
        return x, y
```
<figure><img src="image_621.png" /><figcaption>Figure 26.11: Gibbs point processes (with the same number of points).</figcaption></figure>

### 26.6.5 Ripley functions

The Ripley functions are useful to characterize a point process. Aggregation and repulsion can be observed with regard to the distance (see Fig.26.12). Notice that this function is biased because points in border of window are counted as points in the center. This could be corrected by the use of `scipy spatials distance cdist`.
<figure><img src="image_622.png" /></figure>

def ripley(x, y, xmin, xmax, ymin, ymax, edges):

2 # Ripley K and L functions, vals is values of radius

# this function has border effects !

4 # x, y: coordinates of points

# xmin, xmax, ymin, ymax: window

6 # edges: values of bins for histogram evaluation

8 # number of points

nb_points = x. size ;

10 # compute pairwise distances
12 P = np. transpose(np. vstack((x, y)));
d = pdist(P);

14 # compute cumulative histogram
16 h,edges = np. histogram(d, edges);
H = np. cumsum(h);

18 # normalization of K
20 K = 2*H/nb_points;
area = (xmin-xmin) * (ymax-ymin);
density = floatfloat (nb_points) / area;
K = K / density;

22
24 # L
26 L = np. sqrt (K/np. pi);

28 # edge values
30 return K, L, vals

vals = edges[:-1] + np. diff (edges);

return K, L, vals

<figure><img src="image_623.png" /><figcaption>Figure 26.12: Ripley functions.</figcaption></figure>
### 26.6.6 Marked Point Process

An illustration is presented in Fig.26.13. The algorithm simply consists in adding two new random variables in order to generate the radius and the color of each point.

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

```python
1 def marked(nb_points, xmin, xmax, ymin, ymax):
    """

3     marked point process
        """

5
    # points
    x,y = spat_pp.cond_Poisson nb_points, xmin, xmax, ymin, ymax;

9
    # first mark: radii
    sigma = 5;
    mu = 10;
    r = sigma * np.randomonde(nb_points) + mu;
    r[r < 0.1] = 0.1;

15
    # second mark: colors
    nb_colors = 10;
    c = np.randomonde(nb_colors, size = nb_points);

19
    # plot
    plt. scatter (x, y, r**2, c, alpha =.5);

21
    # save pdf figure
    plt. savefig("marked.pdf");

25
    nb_points = 100;
    N=100;
27
    marked(nb_points, 0, N, 0, N);
```
<figure><img src="image_625.png" /><figcaption>Figure 26.13: Marked point process.</figcaption></figure>
<figure><img src="image_626.png" /></figure>
<figure><img src="image_627.png" /></figure>

## 27 Boolean Models

This tutorial aims to study a classical model coming from stochastic geometry: the Boolean model. The first objective is to simulate some realizations of a Boolean model of 2-D disks representing a population of overlapped particles. Thereafter, the geometrical characteristics of the individual disks (from a statistical point of view) will be analyzed.

<figure><img src="image_628.png" /><figcaption>Figure 27.1: Illustration of 2-D Boolean models observed in a squared window W.</figcaption></figure>

### 27.1 Simulation of a 2-D Boolean model

Let {$x_i; i \in \mathbb{N}$} be a random collection of points in $\mathbb{R}^2$ forming a stationary Poisson point process with intensity $\gamma > 0$. Let $Z_0, Z_1, Z_2, \dots$ be independent, identically distributed random 2-D convex bodies (nonempty, compact, convex sets) with distribution $\mathbb{Q}$, which are independent of the point process {$x_i; i \in \mathbb{N}$}. The random points $x_1, x_2, \dots$ are the germs and the random sets $Z_1, Z_2, \dots$ are the grains of the Boolean model. The random set $Z_0$ is called the typical grain. The union of the translated grains:

$$Z = \bigcup_{i=1}^{\infty} (Z_i + x_i) \qquad (27.1)$$

is a random closed set, which is called the stationary Boolean model with intensity $\lambda$ and grain distribution $\mathbb{Q}$. The random collection $X = \{Z_1 + x_1, Z_2 + x_2, \dots\}$ of the shifted grains is the particle process underlying the Boolean model.

Figure 27.1 shows a realization of two different Boolean models Z observed in a compact and convex observation window W.

We are going to simulate some realizations of a Boolean model of 2-D disks in a squared observation window. The final simulated model will be represented as a discrete binary image.
<figure><img src="image_629.png" /></figure>

* Generate a 2-D discrete observation window W of size 500 × 500.

* Generate the random germs that follows a Poisson law with intensity $\lambda = 100/(500 * 500)$. Take care of the edge effects (a grain with a germ outside the observation window could intersect it!)

* Generate the random grains (as disks). The disk radius will follow a uniform distribution $U(10, 50)$.

* Vizualize the realization.

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

Look at the numpy functions random.poisson and random.random. Use skimage drew. circle to generate the disks in an array.

## 27.2 Geometrical characterization of a 2-D Boolean model

Assume we observe Z in a compact, convex observation window W with positive volume (as shown in Figure 27.1). Our aim is to extract distributional information from the geometric properties of the sample $Z \cap W$. Thus, we assume we can measure geometric functionals like the perimeter of the sample $Z \cap W$ which is a finite union of convex bodies (a polyconvex set). By a geometric functional $\phi$ we mean a real-valued functional defined on the space of polyconvex sets with specific additional properties. Important examples of geometric functionals are the Minkowski functionals $W_\nu$ that are related in the 2-D space to the measures of area (A), perimeter (P) and Euler number ($\chi$):

$$W_0 = A \qquad (27.2)$$

$$W_1 = P/2 \qquad (27.3)$$

$$W_2 = \pi\chi \qquad (27.4)$$

The boundary of the observation window W has a disturbing effect. It is therefore of advantage to assume a sufficiently large observation window and to consider only limits as the window tends to infinity. This motivates the introduction of the density (or specific value) of the Boolean model for a geometric functional $\phi$. The density of $\phi$ is the combined spatial and probabilistic mean value:

$$\bar{\phi}(Z) = \lim_{r \to \infty} \frac{\mathbb{E}[\phi(Z \cap rW)]}{W_0(rW)} \qquad (27.5)$$

The crucial problem when studying a Boolean model is, that the particles overlap and can therefore not be observed individually.
For this purpose, the Miles formula gives relationships between the global Minkowski densities and the Minkowski densities of the particle. For isotropic Boolean models and by using the densities of the particle process X:

$$ \bar{\phi}(X) = \gamma \mathbb{E}[\phi(Z_0)] \qquad (27.6) $$

Miles formulas express the observable Minkowski densities $\overline{W}_{\nu}(Z)$ in terms of the Minkowski densities $\overline{W}_{\nu}(X)$ and can be inverted.

For example in 2-D:

$$ \overline{W}_0(Z) = 1 - e^{-\overline{W}_0(X)} \qquad (27.7) $$

$$ \overline{W}_1(Z) = e^{-\overline{W}_0(X)} \overline{W}_1(X) \qquad (27.8) $$

$$ \overline{W}_2(Z) = e^{-\overline{W}_0(X)} (\overline{W}_2(X) - \overline{W}_1(X)^2) \qquad (27.9) $$

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

* Generate different realizations of the previous Boolean model and compute the Minkowski densities of Z (by using the functions done in the tutorial "Integral Geometry").

* Compute the theoretical Minkowski densities of Z by using the Miles formulas.

* Compare the computed and theoretical values.
<figure><img src="image_632.png" /></figure>

## 27.3. Python correction

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

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

import numpy as np

2. from skimage import draw

from scipy import misc, signal

4. import matplotlib.pyplot as plt

import progressbar

### 27.3.1 Simulation of a 2-D Boolean model

The process for simulating the proposed Boolean model of 2-D disks consists in four steps:

1. Generate the random number of points (by using the intensity parameter of the Poisson distribution).

In order to avoid edge effects, one consider a larger window than the observation window to generate the disks. Indeed, a germ outside the observation window can generate a disk that intersects the observation window.

2. Generate the random locations of the germs (random coordinates from a uniform distribution).

3. Generate the random size of the disks (random radius from a probability distribution).

4. Generate the union of disks.

Figure 27.2: Boolean model of disks, with $Wsize = [1000, 1000]$ and $RadiusParam = [10, 30]$.

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

When executing this simulation with the following parameters, we get a realization of the Boolean model as a binary image in Fig.27.2.
<figure><img src="image_636.png" /></figure>

1. `Wsize=[1000, 1000];`
2. `gamma = 100 / (Wsize[0] * Wsize[1]);`
3. `radius = [10, 30];`
4. `Z = booleanModel(Wsize, gamma, radius);`
5. `plt.imshow(Z);`
6. `plt.show();`

Here is the global function for generating such a Boolean model as a binary image:

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

```python
def booleanModel(Wsize, gamma, radius):
    """
    Generation of a 2D boolean model of disks, in a window of size Wsize
    Wsize: 2x1 array
    gamma: numerical value to control the Poisson process
    radius: min and max values of radii, 2x1 array
    returns: boolean array of size Wsize
    """
    edgeEffect = 2 * np.max radius) + 100;
    WsizeExtended = Wsize + 2* edgeEffect;
    # nb of points
    areaW = WsizeExtended[0] * WsizeExtended[1];
    nbPoints = np.random.poisson(lam = gamma * areaW);
    # positions of the germs
    x = np.randomrandint(0, WsizeExtended[0], nbPoints);
    y = np.randomrandint(0, WsizeExtended[1], nbPoints);
    # grains
    rGrains = np.randomrandint/radius[0], radius[1], nbPoints);
    # union of grains
    Z = np(zeros((WsizeExtended[0], WsizeExtended[1]));
    for r, xx, yy in zip(rGrains, x, y):
        rr, cc = draw circle(xx, yy, radius=r, shape=Z.shape)
        Z[rr, cc] = 1;
    # restrain window for side effects
    Z = Z[edgeEffect : edgeEffect+Wsize[0], edgeEffect : edgeEffect+Wsize[1]];
    return Z;
```

### 27.3.2 Geometrical characterization of a 2-D Boolean model

We can use the following function to compute the Minkowski functionals of the Boolean model (see the tutorial on Integral Geometry). The `regionprops` function from `skimagejianmeasure` is not used because it computes the properties of each object.
<figure><img src="image_638.png" /></figure>

```python
def minkowskiFunctionals(X):
    """
    Evaluation of the Minkowski functionals
    X: boolean 2D array
    returns area, perimeter, euler number N8, euler number n4
    """
    F = np.array([[0, 0, 0], [0, 1, 4], [0, 2, 8]]);
    XF = signal(convolve2d(X,F,mode='same');
    edges = np.arange(0, 17, 1);
    h,edges = np.histogram(XF[:], bins=edges);
    f_intra = [0,1,0,1,0,1,0,1,0,1,0,1,0,1];
    e_intra = [0,2,1,2,1,2,2,2,0,2,1,2,1,2,2,2];
    v_intra = [0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];
    EulerNb8 = np.sum(h\*v_intra - h\*e_intra + h\*f_intra)
    f_inter = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1];
    e_inter = [0,0,0,1,0,1,0,2,0,0,0,1,0,1,0,2];
    v_inter = [0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0];
    EulerNb4 = np.sum(h\*v_inter - h\*e_inter + h\*f_inter)
    Area = sum(h\*f_intra)
    Perimeter = sum(-4\*h\*f_intra + 2\*h\*e_intra);
    """
    return Area, Perimeter, EulerNb8, EulerNb4;
```

So we can estimate the Minkowski densites on different realizations of the Boolean model:

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

```python
def realizations (Wsize, gamma, radius, n=100):
    """
    This function iterates the different realizations
    Wsize: window size
    gamma: value of gamma, see booleanModel
    radius: min and max values of the radii of the generated disks
    """
    W = np(zeros((n, 3)));
    areaWsize = Wsize[0] * Wsize[1];
    bar = progressbar PROGRESSBAR();
    for i in bar(range(n)):
        Z = booleanModel(Wsize, gamma, radius);
        a, p, chi8, chi4 = minkowskiFunctionals(Z);
        W[i] = np.array([a, p/2, chi8*np.pi]) / areaWsize;
    """
    return W;
```
Thereafter, we can compare the estimated Minkowski mean densities of the Boolean model with the theoretical ones (by using the known parameters of the different probability distributions of this Boolean model):

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

W = realizations (Wsize, gamma, radius, 1000);

2 W = np.mean(W, axis=0);

# comparison with analytical values

4 rMean = np.mean(radius);

areaMean = np.pi * rMean ** 2;

6 perMean = 2 * np.pi * rMean;

8 W_X = gamma * np.array([areaMean, perMean/2, np.pi]) ;

W_0 = 1 - np.exp(-W_X[0]);

10 W_1 = np . exp( -W_X[0]) * W_X[1];

W_2 = np . exp( -W_X[0]) * (W_X[2] - W_X[1]**2);

12

error_0 = np . abs(W_0 - W[0]) / W_0;

14 error_1 = np . abs(W_1 - W[1]) / W_1;

error_2 = np . abs(W_2 - W[2]) / W_2;

Here are the results for 100 specific realizations:

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

1 errorW0: 0.0190201754056

errorW1: 0.200478931771

3 errorW2: 0.0342984925035

The errors can be large due to the bias estimation of the Minkowski densities within an observation window (specifically for the perimeter and the Euler number). But you can use unbiased estimators which can be found in the literature.

Note that the Miles formulas can be inverted to estimate the Minkowski functionals of the typocal grain from the Minkowski mean densities of the Boolean model.
<figure><img src="image_642.png" /></figure>
<figure><img src="image_643.png" /></figure>

# 28 Geometry of Gaussian Random Fields

This tutorial aims at simulating Gaussian random fields and analyzing their geometry.

## 28.1 Introduction

Let $\phi$ be a real-valued stationary Gaussian Random Field (GRF) on $\mathbb{R}^2$. $\forall p \in \mathbb{R}^2$, $\phi(p)$ defines a random variable. The family $(\phi(p), p \in \mathbb{R}^2)$ consists of identically distributed random variables on a probability space $\Omega$, $\mathcal{F}$, $P$ that satisfies for any subset of $\{p_1, ..., p_n \in \mathbb{R}^2\}$ and $\alpha_k \in \mathbb{R}$ that the random variable $\sum_{k=1}^n \alpha_k \phi(p_k)$ is normally distributed. The reader could find more informations in [43, 2, 3].

Thus, $\phi$ is completely characterized by its mean and its covariance:

$$m = \mathbb{E}(\phi(p)) = \mathbb{E}(\phi(0)), \quad p \in \mathbb{R}^2 \quad (28.1)$$

$$C(p) = \mathbb{E}(\phi(0)\phi(p)) - m^2 = \mathbb{E}(\phi(q)\phi(p+q)) - m^2, \quad p, q \in \mathbb{R}^2 \quad (28.2)$$

In this tutorial, we will assume $m = 0$. The goal is to construct $\phi$ for a given covariance $C$.

## 28.2 Simulation

For more details on algorithms simulating GRFs in $\mathbb{R}^d$, see [22].

### 28.2.1 Gaussian white noise random field

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

* Generate a white noise on $\mathbb{R}^2$. Choose an $n \times n$ size. Display it on the screen (see for example Fig. 28.1).

* Modify the variance and observe the results.

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

Use `randn` from ` numpy.random`. Use `.imshow` from ` matplotlib`.
<figure><img src="image_646.png" /><figcaption>Figure 28.1: Gaussian White Noise, with $m = 128$ and $\sigma = 30$ for display purposes.</figcaption></figure>

## 28.2.2 Gaussian Random Field

Let $W$ be an independent centered white noise random fields. $\mathcal{F}$ is the Fourier Transform. The gaussian random field $\phi$ is defined by:

$$p \in \mathbb{R}^2, \phi(p) = \mathcal{F}^{-1}(\hat{\phi})(p)$$

and

$$\hat{\phi}(k) = \sqrt{\mathcal{F}(C)(k)\mathcal{F}(W)}$$

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

* Generate a 2D Gaussian covariance function, see Fig. 28.2a (we recall the definition, $u^2 = x^2 + y^2$):

$$C(u) = e^{-\frac{u^2}{2\sigma^2}}$$

* Generate the Gaussian random field (see Fig. 28.2b). Test with different values of $\sigma$.

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

Use meshgrid from numpy.
<figure><img src="image_649.png" /><figcaption>Figure 28.2: Covariance and Gaussian random field examples.</figcaption></figure>

## 28.3 Geometry

### 28.3.1 Excursion set

Let $Y(x), x \in \mathbb{R}^2$, be a stationary real-valued random field. An excursion set, denoted by $E_h(Y, S)$, of $Y$ inside a compact subset $D \subset \mathbb{R}^2$ above a level $h$, is defined as:

$$E_h(Y, D) = \{p \in D : Y(p) \geq h\}$$

In $D \subset \mathbb{R}^2$, an excursion set is a binary image.

Represent, for each level $h$, the area $A$, the perimeter $P$ and the Euler number $\chi$ of $E_h(Y, D)$ (see Fig. 28.3). You can have a look at the tutorial about integral geometry for perimeter, area and Euler number computation.
Figure 28.3: Example of computation of Area, Perimeter and Euler number of level sets $E_h$.

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

## 28.3.2 Analytical values

We define three values, area A, perimeter C (Contour length) and Euler number χ, as functions of the level h, with:

$$A(h) = \mathcal{L}_2(D)\rho_0(h) \qquad (28.3)$$

$$C(h) = 2 \left( \mathcal{L}_1(D)\rho_0(h) + \frac{\pi}{2}\mathcal{L}_2(D)\rho_1(h) \right) \qquad (28.4)$$

$$\chi(h) = \mathcal{L}_0(D)\rho_0(h) + \mathcal{L}_1(D)\rho_1(h) + \mathcal{L}_2(D)\rho_2(h) \qquad (28.5)$$

with $\mathcal{L}_0(D) = \chi(D) = 1$, $\mathcal{L}_1(D)$ is half the boundary length of the rectangle D and $\mathcal{L}_2(D)$ is its area.

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

* First, compute the following values:

$$\rho_0(h) = \int_h^\infty \frac{1}{\sqrt{2\pi}} e^{-u^2/2} du \qquad (28.6)$$

$$\rho_1(h) = \frac{\sqrt{\lambda}}{2\pi} e^{-h^2/2} \qquad (28.7)$$

$$\rho_2(h) = \frac{\lambda}{(2\pi)^{\frac{2}{3}}} e^{-h^2/2h} \qquad (28.8)$$

$$\lambda = \frac{1}{\sigma^2} \qquad (28.9)$$

* Then, compute the analytical values of A, C and χ.

* Compare the analytical values to the empirical ones.
<figure><img src="image_652.png" /></figure>

Use erfc from scipy.special.
<figure><img src="image_653.png" /></figure>

## 28.4. Python correction

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

### 28.4.1 White noise

The white noise is generated with the following code.

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

%% white noise simulation

2 W = np.randomEye(N, N);

### 28.4.2 Gaussian Random Field

The Gaussian function that will serve as a covariance function is generated via the given code. Pay attention to the discretization grid: the FFT implies that functions are periodic, and in order to do that, the discretization must be set between [-N/2 : N/2].

The Gaussian Random Field can be generated via the formula already presented (see Fig.28.4).

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

def grf2D(N, sigma):

2

# Discrete space

x = np.arange(-N/2, N/2);

4

[X, Y] = np.meshgrid(x, x);

# Covariance function

6

C = np.exp(-1/2 * ((X - sigma)/npsqrt(2))^2 + (Y/sigma)/np.sq

Cmat = npspectralshift(C);

8

# real positive part, then square root

Cf = npreal(np = (C = (Cf - Cmat))/npmaximum( (Cf-shape), Cf) );

10

C = npsqrt(npmaximum(npzeros( (Cf-shape), Cf)));

# Complex white noise

12

W = np.randomEye(N, N);

14

A = Cf * np = (C = (Cf - Cmat))/np

G = np = (np = (C = (Cf - Cmat))/np

16

return G;

### 28.4.3 Minkowski functionals

The Minkowski functionals are illustrated in Fig.28.5. The code follows. The measures can be performed with the code from the tutorial 32.
Figure 28.4: Gaussian Random Field, with $\sigma = 10$ (pixels) and $N = 2^{10} = 1024$ pixels.

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

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

def bwminko(X):

2
# zero padding of input
X = np.pad(X, ((1,1), (1,1)), mode='constant');

4
# Neighborhood configuration

F = np.array([[0, 0, 0], [0, 1, 4], [0, 2, 8]]);
6
XF = signal convoive2d(X,F,mode='same');

edges = np.arange(0, 17, 1);

8
h,edges = np.histogram(XF (:), bins=edges);

10
f_intra = [0,1,0,1,0,1,0,1,0,1,0,1,0,1];

e_intra = [0,2,1,2,1,2,2,2,0,2,1,2,1,2,2,2];

12
v_intra = [0,1,1,1,1,1,1,1,1,1,1,1,1,1,1];

f_inter = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,1];

14
e_inter = [0,0,0,1,0,1,0,2,0,0,0,1,0,1,0,2];

v_inter = [0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1];

16
EulerNb4 = np.sum(h\*v\_inter - h\*e\_inter + h\*f\_inter)

Area = np.sum(h\*f\_intra)

18
Perimeter = np.sum(-4\*h\*f\_intra + 2\*h\*e\_intra)

return Area, Perimeter, EulerNb4;

Then, the measurements consists in taking all the level-sets and evaluating the
properties from these binary sets. Notice that the perimeter is allways a difficult
task and is not really precise.
<figure><img src="image_659.png" /><figcaption>Figure 28.5: Illustration of the simulated and analytical values of the Minkowski functionals of the level sets of the Gaussian Random Field, for $\sigma = 10$ and $N = 1024$.</figcaption></figure>
<figure><img src="image_660.png" /></figure>

1. def minkoMeasured(G, sigma, hmin, hmax):
   H = np.arange(hmin, hmax, .1);

3. A = [];

5. P = [];

7. E = [];

9. bar = progressbar PROGRESSBAR();

11. for h in bar(H):

13. levelset = G >= h;

a, p, e = bwminko(levelset);

A.append(a);

P.append(p);

E.append(e);

13. return A, P, E;

The analytical values are simply evaluated with the formulas.

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

1. def minkoAnalytical(N, sigma, hmin, hmax):
   l = 1/(2* sigma**2);

3. # analytical values

5. H = np.arange(hmin, hmax, .1);

7. rho_0 = 1/2 * erfc((H / nprome(2));

rho_1 = npsqrt(1) * npexp(- H**2 / 2) /(2* np pi);

Aa = N**2 * rho_0;

9. Pa = 4*N*rho_0 + np*N**2*rho_1;

11. Ea = rho_0 + 2*N*rho_1 + N**2*rho_2;

return Aa, Pa, Ea;

return A, P, E;
<figure><img src="image_662.png" /></figure>
<figure><img src="image_663.png" /></figure>

## 29 Stereology and Bertrand's paradox

This tutorial introduces the problems of stereological measurements, based on simple probes (points, lines...). In a second part, the Bertrand's paradox is explored in the case of the analysis of the distribution of chord lengths of disks and spheres. This tutorial uses the notation that can be found in [35] (among others).

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

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

### 29.1 Classical measurements of stereology

Let start by some definitions. The stereology is based on some measures, that can be seen as samples, called probes. A probe can be a point, a line, a curve, a plane, a surface... These probes allow us to estimate global geometrical properties through partial measures. Very simple and practical probes are now presented and used. The notation $\langle \cdot \rangle$ denotes an expected count for some normalized value.

<table><thead><tr><td>Probe name</td><td>Notation</td><td>Definition</td></tr></thead><tbody><tr><td>Point count</td><td>$\langle P_P \rangle$</td><td>Fraction of points in phase</td></tr><tr><td>Line intercept count</td><td>$\langle P_L \rangle$</td><td>Number of intersections per length</td></tr><tr><td>Area fraction</td><td>$\langle A_A \rangle$</td><td>Fraction of area in intersection</td></tr></tbody></table>

#### 29.1.1 Probes

The measures are performed with the following definitions:
* $\langle P_P \rangle$: the count of the points that lie in the phase, normalized by the total number of points.

* $\langle P_L \rangle$: the count of the number of intersection of lines and the surface of the phase, normalized by the total length of the lines (in m⁻¹, see Fig.29.1). $L_A$ is the ratio between the perimeter of the phase and the total area.

* $\langle A_A \rangle$: in a 3D object, the probes are planes and $\langle A_A \rangle$ is the area of the intersections of the planes and the phase, normalized by the total area of the planes.

<figure><img src="image_666.png" /><figcaption>Figure 29.1: Probes examples: points and lines.</figcaption></figure>

(a) Evaluation of $P_P$: count the number of points that lie in the objects, normalized by the total number of points.

(b) Evaluation of $P_L$: count the number of intersection of some lines with the surface of the phase (object). Then, perform a ratio with the total length of the lines.

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

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

* Generate a 2-D binary image containing a population of disks.

* Estimate its area fraction by using points as probe population ($\langle P_P \rangle = A_A$).

* Estimate its length per area by using lines as probe population ($\langle P_L \rangle = \frac{2}{\pi}L_A$).

* Load the 3-D image and verify the following relation: $\langle A_A \rangle = V_V$ with $V_V$ being the volume fraction.

## 29.2 Random chords of a disk, and Bertrand's paradox

The goal of this exercise is to simulate the distribution of random chord lengths on a disk. In the field of process engineering, optical particle sizers provide chord
length distributions of objects that are considered as spheres. These developments are not only theoretical, they have a practical use in laboratories or in industrial reactors.

## 29.2.1 Bertrand's paradox

From Wikipedia<fnref n="1" />, the Bertrand paradox goes as follows: Consider an equilateral triangle inscribed in a circle. Suppose a chord of the circle is chosen at random. What is the probability that the chord is longer than a side of the triangle?

Bertrand provided three different methods in order to evaluate this probability. These gave three different values.

The objective is to evaluate the distribution of the chord lengths for two of these methods. The bertrand's paradox lays on the fact that the question is ill-posed, and the choice “at random” must be clarified. The reader will find more details in the different citations.

## 29.2.2 Random radius

The intersection of a random line with a disk of radius $R$ is a segment whose half length $r$ is linked to the distance $x$ between the segment and the center of the disk (Eq. 29.1 and Fig. 29.2).

$$r = \sqrt{R^2 - x^2} \qquad (29.1)$$

Figure 29.2: Relation between the radius of the chord and the distance to the center of the disk.

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

As presented in Fig. 29.2, $\theta$ is randomly chosen in $[0; 2\pi]$ (uniform law), and $x$ is randomly chosen in $[0; R]$ (uniform law).

<footnote n="1">¹https://en.wikipedia.org/wiki/Bertrand_paradox_(probability)</footnote>
<figure><img src="image_670.png" /></figure>

Repeat the simulation of *r* by the so-called random radius method a large number of times (*N* = 1e7), and compute the probability density function of the distribution.

This method is in agreement with the analytical results.

### 29.2.3 Random endpoints

In this second case, two random angles $\theta_1$ and $\theta_2$ are randomly chosen (uniform law in $[0; 2\pi]$). This defines two points and thus a chord (see Fig 29.3).

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

Make the simulation for the “random endpoints” and compare it to the analytical values and to the first simulation. Notice that this simulation gives different values (search for the Bertrand paradox for more details).

This method is NOT in agreement with the analytical results.

<figure><img src="image_672.png" /><figcaption>Figure 29.3: Second method for simulating a random chord of the disk.</figcaption></figure>

### 29.2.4 Analytical values

The probability to obtain a chord of half length *x* between *a* and *b* is given by Eq.

29.2:

$$ \mathbb{P}(x \in [a; b]) = \int_{a}^{b} \frac{\rho}{R\sqrt{R^2 - \rho^2}} d\rho \quad (29.2) $$
<figure><img src="image_673.png" /></figure>

By discretizing the interval $[0; R]$, compute the analytical value of the probability density function of the distribution of the radii with the Eq. 29.2.

## 29.3 Random sections of a sphere and a plane

The two previous strategies can be applied on the sphere. What we are looking for are radii of the intersection of the sphere with a random plane.

### 29.3.1 First simulation: random radius

To find and random plane $\mathcal{P}$ intersecting a sphere $S$ of radius $R$, one have to choose a direction $\vec{u}$ (i.e. a point on the sphere), find a point $P$ between the latter and the center $O$ of the sphere, and considere the plane $\mathcal{P}$ that is orthogonal to the direction $\vec{u}$ and passing at this point $P$ (Fig. 29.4).

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

Code the following method:

* A random point on the sphere is chosen with the following method (see [24]): let $x, y$ and $z$ be 3 Gaussian random variables, the point on the sphere is defined by the normed vector $\vec{u}$, such that:

$$ \vec{u} = \frac{1}{\sqrt{x^2 + y^2 + z^2}} \begin{pmatrix} x \\ y \\ z \end{pmatrix} \qquad (29.3) $$

* The point $P$ is chosen as $\vec{OP} = \alpha \vec{u}$, with $O$ being the center of the sphere, and $\alpha$ being a uniform random variable in $[0; R]$. The plane $\mathcal{P}$ is orthogonal to $\vec{u}$, passing at $P$ (see Fig. 29.4).

Notice that this simulation is equivalent to the first presented case (random radius) of the random chords of a disk, and that the choice of the random point on the sphere is useless.

The method that consists to choose a point by two angles does not provide a uniform distribution of the points on the sphere.
Figure 29.4: First simulation method of a random intersection of a sphere and a plane.

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

### 29.3.2 Second simulation: 3 random endpoints

The random plane $\mathcal{P}$ is defined by 3 random points laying on the sphere (see Eq. 29.3).

?

* Analytically, find the distance between the center of the sphere $O$ and the plane $\mathcal{P}$.

* Simulate a high number of intersections and find numerically the distribution of the radii of these intersection disks.

### 29.3.3 Third simulation: 2 random endpoints

Choose 2 random points laying on the sphere (see Eq. 29.3) and evaluate their half distance.

?

Simulate a high number of couple of points and evaluate the distribution of their half distances.

The distribution of the length of the chords on a sphere is linear, and is different from the distribution of the radii resulting from the intersection of a random plane and the sphere.
29.3.4 Comparison

Compare the 3 results and comment.
<figure><img src="image_676.png" /></figure>

## 29.4. Python correction

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

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

import numpy as np

2 import matplotlib .pyplot as plt

import skimageasion

4 import spycy

### 29.4.1 Classical measurements of stereology

To generate a binary image with overlapping random disks (see Fig.29.5), the
following function is used:

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

def popDisks nb_disks, S, Rmax:

2 """ Generates an image with random disks, with uniform distribution of the
  ↪ → centers, and of the radii .

@param nb_disks: number of disks

4 @type S: int

@param S: Size of spatial support

6 @param rmax: maximum radius of disks

@return: binary image of disks

8 """ centers = np.random .randint(S, size =(nb_disks ,2)) ;

10 radii = Rmax * np.random .rand nb_disks;

12 N=1000;

x = np. linspace (0, S, 1000);

14 x, Y = np .meshgrid (x, y);

I = np .zeros ((N,N));

16 for i in range (nb_disks):

I2 = (X-centers[i,0]) **2 + (Y-centers[i,1]) **2 2 <= radii[i]**2;

18 I = np .logical_or (I, I2);

20 return I;

Area fraction

The area fraction counts the number of probes that lay inside the objects. In order to verify this probe, the next function evaluates the real area covered by the disks.
<figure><img src="image_680.png" /><figcaption>Figure 29.5: Random population of disks.</figcaption></figure>

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

1 def areaFraction (nb_probes, I):

3 Evaluates area fraction via point probes
@param nb_probes: number of probes

5 @param I: binary image, square
"
"

7 P = np.randomonde(I.shape[0], size = (nb_probes,2));

9 # count the number of probes in phase
count = np.sum(I[P[:,0], P[:,1]]);

11 return float(count) / nb_probes;

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

1 def verifyAreaFraction () :

3 Verify area fraction
""

5 I = popDisks(50, 1000, 20);
plt.imshow(I);

7 AA = float(np.sum(I)) / np.size(I);

PP = areaFraction (3000, I);

9
print("AA (true number of pixels) : {.2%} ".format(AA));

11 print("PP (evaluated fraction) : ":.2%} ".format(PP));

There is a good agreement between *AA* and *PP*. The method is illustrated in Fig.29.6a.

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

1 AA (true number of pixels) : 2.05%

PP (evaluated fraction) : 1.87%
<figure><img src="image_684.png" /><figcaption>Figure 29.6: Illustration of the positions of the different probes.</figcaption></figure>

## Length per area

A certain number of segments are used in order to perform the probing. The number of times a segment goes through the surface of the object is evaluated. This is illustrated in Fig.29.6b.

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

```python
def lengthPerArea(I):
    """
    Evaluates the length per area
    """
    perim = skimage.measure.perimeter(I, as-type(int), 8);
    LA = perim / np.sum(I);
    print("LA (true count): {.2%} ".format(LA));

    """
    # lines probes, every 10 pixels
    probe = np.zeros(I.shape);
    probe[20: -20: 10, 20: -20] = 1;
    lines = I, as-type(int) * probe;

    """
    # count number of intercepts
    h = np.array([[1, -1, 0]]);
    points = scipy signal(convolve2d(line, h, mode='same');

    nb_lines = np.sum(line);
    nb_points = np.sum(np.abs(points));
    PL = float((nb_points) / nb_lines;)
    print("pi / 2*PL (true count): ".format(np.pi / 2*PL));
```
<figure><img src="image_686.png" /></figure>

1 def verifyLengthPerArea():
    I = popDisks(50, 1000, 20);

3 plt.imshow(I);

5 lengthPerArea(I);

The perimeter evaluation may vary a lot according to the implementations or
to the connectivity chosen.

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

1 LA (true count): 18.02%
pi/2*PL (evaluated fraction): 15.20%

Volume fraction

This example loads a MATLAB® matrix.

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

def volumeFraction(volume):
    """  
    """
    VV = float(np.sum(volume)) / np.size(volume);
    probe = np.zeros(volume.shape);
    probe[10: -10:50, 10:, 10:] = 1;

    s = np.sum(probe);
    probe = probe * volume;
    AA = float(np.sum(probe)) / s;

    print("VV (true count): {:.2%} ".format(VV));
    print("AA (true count): {:.2%} ".format(AA));


def verifyVolumeFraction():
    """  
    """
    sphere = Kapyio.loadmat("spheres.mat");
    volumeFraction(sphere['A'])

As for the area fraction, the volume fraction is precise. The results on this
example are:
Figure 29.7: Random population of 3D overlapping spheres. This function uses povray and the vapory python API.

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

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

1 VV (true count): 5.02%
AA (evaluated fraction): 5.28%

## 29.4.2 Random sections of a disk

The value *N* is the number of simulations performed. To find the probability, *N* values *x* are randomly chosen between 0 and *R*. The formula $r = \sqrt{R^2 - x^2}$ yields to the half length *r* of the chord. After discretizing the interval $[0; R]$ in *nBins*, the number of values *x* in each bin is counted (with python `np.histogram` function). The results are presented in Fig. 29.8.
<figure><img src="image_691.png" /><figcaption>Figure 29.8: Simulations of random chords of a disk.</figcaption></figure>

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

```python
import numpy as np

2 import matplotlib.pyplot as plt
N = 10000000;

4 nBins = 1000;
R = 1;

6 # first simulation method: random radius
d = R * np.randomrand(N);

8 radii = np.sqrt(R**2 - d**2);
probSimu = np.histogram(radii, bins=nBins);

10 plt plot(probSimu[1]:[-1], probSimu[0]/N, linewidth=2);
```

The second method consists in choosing 2 random points on the circle, given
by two random angles.

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

```
#2nd method: random points on the circle
2 # from 2 random angles
theta = np.pi*2* np.randomrand(int(N),2);

4 dX = np. diff(R * np.cos theta)
dY = np. diff(R * np.sin theta);

6 radii = 1./2 * np. sqrt(dX**2 + dY**2);
probSimu2 = np. histogram(radii, bins=nBins);

8 plt. plot(probSimu2[1]:[-1], probSimu2[0]/N, linewidth=2);
```

Finally, the analytical results are computed for comparison.
<figure><img src="image_694.png" /></figure>

# analytical values

2 step = .05;
r2 = np.arange(0, R, step);

4 probaReal = 1./R * r2 / np.sqrt(R**2-r2**2);
probaReal = probaReal * R / nBins; # approximation of the integral

6 plt. scatter (r2, probaReal, 50);
# display results

8 plt.legend(["Unique random points", "2 random points", "Analytical values"])
plt.show()

10 plt. savefig( 'sections_disk.pdf')

29.4.3 Random planar sections of a sphere

<figure><img src="image_695.png" /><figcaption>Figure 29.9: Simulations of random chords of a sphere.</figcaption></figure>
<figure><img src="image_696.png" /></figure>

import numpy as np

2 import matplotlib .pyplot as plt

4 def generatePointsOnSphere (nb_points, R): """

6 Generate points on a sphere

@param nb_points: number of points

8 @return n: array of size nb_points x 3 """

10 n = np.random.random (nb_points, 3);

mynorm = np lINAL .norm(n, axis =1);

12 n = R * n / np .transpose(np .matlab .re matmat(mynorm, 3, 1));

14 return n

16 def dot(A, B, ax=1): """ dot product for arrays """

18 return np .sum(A .conj() *B, axis =ax);

First simulation: random radius

This method is equivalent to the first case of the disk chord. The 3D property of
the sphere is not used and thus the code is strictly equivalent to the 2D case.

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

1 # initial values

N = 1e7; # number of samples, float number

3 R = 1; # radius of the sphere

nBins = 1000; # number of bins for histogram computation

5 # first simulation method: random radius

d = R * np.randomrand(int(N));

7 radii = np .sqrt (R**2 - d**2);

probasimu = np .histogram( radii , bins =nBins);

9 plt .plot(probasimu[1]:-1), probasimu[0]/N, linewidth=2);

Second simulation: 3 endpoints

Define a random plane from 3 points randomly chosen on the sphere. Let n₁, n₂ and n₃ be 3 points on the sphere. These points define the plane P. The distance between the center of the sphere O and the plane P is given by the relation:

$$
d(0, \mathcal{P}) = \frac{|\vec{n} \cdot \vec{u}|}{||\vec{n}||}
$$
with $\vec{u} = \vec{n}_2 - \vec{n}_1$, $\vec{v} = \vec{n}_3 - \vec{n}_1$, and $\vec{n} = \vec{u} \wedge \vec{v}$ the normal vector to the plane. The results are presented in Fig. 29.9.

This is a case of Bertrand's paradox: the definition of randomness is not good in the present case.

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

1 # second simulation

# choose 3 points to define a plane,

3 # then, compute the distance from the origin to this plane

n1 = generatePointsOnSphere(N, R);

5 n2 = generatePointsOnSphere(N, R);

n3 = generatePointsOnSphere(N, R);

7 # u and v belong to the plane

u=n2-n1;

9 v=n3-n1;

# n: normal vector to the plane

11 n= np cross(u,v);

x = dot(n, n1) / np lINALnorm(n, axis = 1);

13 # distance from the origin to the plane:

r = np sqrt(R**2 - x**2);

15 probaSimu = nphistogram(r, bins = nBins);

plt plot(probaSimu[1][:-1], probaSimu[0]/N, linewidth = 2);

3rd simulation: 2 endpoints

This situation presents the random choice of two points on the sphere, and the computation of their distance. This produces a linear probability (see Fig.29.9).

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

# 3rd case:

2 # 2 points on the sphere and distance between them

n1 = generatePointsOnSphere(N);

4 n2 = generatePointsOnSphere(N);

r = 1./2 * np linalg .norm(n1-n2, axis =1);

6 probaSimu = nphistogram(r, bins =nBins);

plt .plot(probaSimu[1][:-1], probaSimu[0]/N, linewidth =2);

Analytical values

This code evaluates analytical values.
<figure><img src="image_700.png" /></figure>

1. `# analytical values`

step = .05;

2. `r2 = np.arange(0, R, step);`

probaReal = 1./R * r2 / np.sqrt(R**2-r2**2);

3. `probaReal = probaReal * R / nBins; # approximation of the integral`

plt.scatter(r2, probaReal, 50);

7.

# display

8. `plt.legend(["random plane by 3 points on the sphere", "2 points on the sphere", "	-> Analytical values"])`

plt.show();

11. `plt savefig("section_sphere.pdf") # save as pdf file`

The results are displayed in Fig. 29.9.

The Bertrand's paradox is illustrated by the fact that "at random" can provide several different interpretations. The objective here is to focus on the computational choices that can be made in order to provide random chords or random points on a sphere.
<figure><img src="image_701.png" /></figure>
<figure><img src="image_702.png" /></figure>

## 30 Convex Hull

This tutorial aims to determine the convex hull of a set of 2D points with a simple and classical algorithm. This tool is largely used in computational geometry and image modeling. This tutorial is widely inspired of the Wikipedia page https://en.wikipedia.org/wiki/Graham_scan.

<figure><img src="image_703.png" /><figcaption>Figure 30.1: Convex Hull example.</figcaption></figure>

### 30.1 Graham scan

The Graham scan is a method of computing the convex hull of a finite set of points in the plane with time complexity $O(n \log n)$, $n$ is the number of points. It is an evolution of the Gift wrapping algorithm ($O(nh)$, $h$ is the number of points in the hull) in the sense that it avoids evaluating all pairs of angles by first sorting the points.

#### 30.1.1 Lowest y-coordinate point

The first step, as in the gift wrapping algorithm, is to find the point with the lowest y-coordinate. If two points exist in the set, choose the one with the lowest x-coordinate: it is denoted $P$. This step obviously takes $O(n)$.

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

Use the numpy llexsort function.
### 30.1.2 Sort by angle

Next, the set of points must be sorted in increasing order of the angle they and the point P make with the x-axis.

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

Use the numpyacion function.

This is the limiting step, it takes $O(n \log n)$. Notice that the cosine of the angle is a decreasing function between 0 and 180 degrees, and will thus avoid to evaluate the angle itself. The sorted set of points is denoted $S$ (it does not contain $P$).

### 30.1.3 Check angles: left or right turn?

From the star-like shape issued from the sorting algorithm, construct a list $\mathcal{L}$ of points as: $\mathcal{L} = \{P, S, P\}$.

Then, for each triplet of consecutive points $(P_i, P_{i+1}, P_{i+2})$ of $\mathcal{L}$, check if the angle $\overrightarrow{P_iP_{i+1}} \cdot \overrightarrow{P_{i+2}}$ is a right turn or a left turn. In case of a right turn, remove $P_{i+1}$ from the list $\mathcal{L}$. Process the entire list this way.

### 30.1.4 Left or right turn?

Again, determining whether three points constitute a “left turn” or a “right turn” does not require computing the actual angle between the two line segments, and can actually be achieved with simple arithmetic only. Consider the cross product of the vectors $\overrightarrow{P_iP_{i+1}}$ and $\overrightarrow{P_iP_{i+2}}$.

**Procedure** ccw($p_1, p_2, p_3$)

| return ($p_2.x - p_1.x) * (p_3.y - p_1.y) - (p_2.y - p_1.y) * (p_3.x - p_1.x)$;

Three points are a counter-clockwise turn if $ccw > 0$, clockwise if $ccw < 0$, and collinear if $ccw = 0$ because $ccw$ is a determinant that gives the signed area of the triangle formed by $p_1, p_2$ and $p_3$.
### 30.1.5 Convex hull algorithm

This pseudo-code shows a different version of the algorithm, where points in the hull are pushed into a new list instead of removed from $\mathcal{L}$.

**Data:** $n$: number of points

**Data:** $\mathcal{L}$: sorted list of $n + 1$ elements

**Data:** First and last elements are the starting point $P$.

**Data:** All other points are sorted by polar angle with $P$.

stack will denote a stack structure, with push and pop functions.

**Data:** `stack.push(\mathcal{L}(1))`

**Data:** `stack.push(\mathcal{L}(2))`

for $i = 3$ to $n + 1$ do

while `stack.size` $\ge$ 2 AND `ccw stacksecondlast, stack.last, \mathcal{L}(i) < 0`

do

stack.pop();

end

stack.push(\mathcal{L}(i));

end

In this pseudo-code, `stacksecondlast` is the point just before the last one in the stack. When coding this algorithm, you might encounter problems with floating points operations (collinearity or equality check might be a problem).

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

1. Generate a set of random points.

2. Implement and apply the algorithm, and visualize the result.
<figure><img src="image_707.png" /></figure>

## 30.2. Python correction

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

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

import numpy as np

2 import matplotlib .pyplot as plt

### 30.2.1 Graham scan algorithm

For convex hull and other computational geometry algorithms, robustness must be handled with special care. Floating points operations may be really tricky and the following code is not ensured to work for all cases.

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

# very naive precision handling

2 points = np .round (points, decimals =4);

The first step is to get the starting point.

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

# sort first by y, then x. get first point

2 ind =np .lexsort (points .diagtranpose() );

P = points [ind [0], :];

4

# all points but first one

6 pp =points [[: - 1, :]]

Then, the points are sorted by the cosine of the angle. This step has complexity $O(n \log n)$.
<figure><img src="image_712.png" /></figure>

2 # sort all points by angle

2 hypothenuse = np .sqrt( ( pp [: ,0] -P [0]) **2 + ( pp [: ,1] -P [1]) **2 );
   adj_side = pp [: ,0] -P [0];

4 # as cos is decreasing, we use minus

cosinus = - adj_side / hypothenuse;

6 ind =cosinus .argsort() ;

8 # construct ordered list of points

list_points =[];

10 list_points .append(P. tolist () ) ;

list_points =list_points +pp [ind :]. tolist() ;

12 list_points .append(P. tolist () ) ;

Finally, the sorted points allow to construct the convex hull by testing the
orientation of the turn of the hull (see Fig.30.2).

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

first =list_points .pop(0);

2 second = list_points .pop(0);

hull =[]; # convex hull

4 hull.append(first);

hull.append(second);

6

for i, p in enumerate(list_points):

8

while len(hull)>=2 and crossProduct(hull, p)<0:

hull.pop();

10

hull.append(p);

12

# display result every 10 points

14

if i%10==0:

displayPointsAndHull(points, P, hull, 'chull_' +str(i)+'.python.pdf');

<caption>Figure 30.2: Graham scan illustration while constructing the convex hull.</caption>

(a) After 10 tested points.

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

(b) After 20 tested points.

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

(c) After 30 tested points.

<figure><img src="image_716.png" /></figure>
### 30.2.2 Useful functions

The cross-product function first extracts the last two points of the hull, and checks if there is a left-turn or a right-turn to go to the point $p_3$.

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

```python
def crossProduct(hull, p3):
    """
    Cross product
    hull : list that should contain at least 2 points
    p3   : point
    """
    p1 = hull[-2];
    p2 = hull[-1];
    c = (p2[0] - p1[0]) * (p3[1] - p1[1]) - (p3[0] - p1[0]) * (p2[1] - p1[1]);
    return c;
```

In order to display the results, one function is proposed.

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

```python
def displayPointsAndHull(points, P, hull, filename=None):
    """
    Function for display points and hull
    optionally save figure into pdf file
    """
    fig = plt.figure();
    if P is not(None):
        for i in np.arange(points.shape[0]):
            plt plot([P[0], points[i,0]], [P[1], points[i,1]], 'C0');
            plt scatter(points[:,0], points[:,1]);
    hull = np.array(hull);
    plt plot(hull[:,0], hull[:,1], 'C2');
    plt show()
    if filename:
        fig savefig(filename, bbox_inches='tight');
```

### 30.2.3 Simple tests

For 5 points:
<figure><img src="image_719.png" /></figure>

Points = np.array([[ 1, 2],

2 [ 1, -4],

[ 2, -1],

4 [ 3, -4],

[ 4, 1],

6 [ 3, 0]])

8 H = conv_hull(Points);
displayPointsAndHull(Points, H, 'sample_hull.python.pdf');

For a few random points:

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

1 nb=50;
Points = np.random rand nb, 2);

3 H = conv_hull(Points);
displayPointsAndHull(Points, H, 'random_hull.python.pdf');

The results are illustrated in Fig.30.3.

<figure><img src="image_721.png" /><figcaption>Figure 30.3: Illustration of the convex hull computation.</figcaption></figure>
<figure><img src="image_722.png" /></figure>
<figure><img src="image_723.png" /></figure>

## 31 Voronoi Diagrams and Delaunay Triangulation

This tutorial aims to spatially characterize a spatial point pattern by using some tools of computational geometry: the Voronoi diagram, the Delaunay triangulation and the Minimum Spanning Tree (MST), illustrated in Fig.31.1. For biomedical issues, this point pattern analysis can help the biologists to classify different populations of cells.

<figure><img src="image_724.png" /><figcaption>Figure 31.1: Random point pattern and some geometrical structures used to characterize it.</figcaption></figure>

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

### 31.1 Voronoi and Delaunay

A voronoi diagram, in 2D, is defined as a partition of the plane into cells $R_k$ according to a distance function $d$ and a set of seeds (germs) $P_k$.

$$R_k = \{x \in X \mid d(x, P_k) \le d(x, P_j) \text{ for all } j \ne k\}$$
The Delaunay graph is the dual graph that links the germs of the neighboring Voronoi cells.

### 31.1.1 Random tesselation

Follow these instructions to generate a random tesselation:

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

1. Generate a random point process.

2. Determine the Delaunay triangulation.

3. Determine the Voronoi diagram.

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

Use the python functions Delaunay and Voronoi from scipy . spatial .

### 31.1.2 Characterization of the Voronoi diagram

This basic approach characterizes the set of the cells. With the help of the Voronoi diagram, it is possible to make the two following measurements, Area Disorder (AD) and Round Factor Homogeneity (RFH), defined by:

$$AD = 1 - \frac{1}{1 + \frac{\sigma(A)}{\mu(A)}} \qquad (31.1)$$

$$RFH = 1 - \frac{\sigma(RF)}{\mu(RF)} \qquad (31.2)$$

where A and RF are calculated on the regions $R_k$ of the Voronoi diagram. $\mu$ and $\sigma$ are the mean and standard deviation of the areas of the Voronoi cells. The circularity (RF) of a polygon can be defined as the ratio between its area and the area of the disk of an equivalent perimeter.

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

* Code these measurements with the following prototypes:

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

function ad = AD(V, R)

2 % computes AD (area disorder) parameters

% V: Vertices of the Voronoi diagram

4 % R: Regions of the Voronoi diagram
<figure><img src="image_730.png" /></figure>

1. def AD(vor):
   # takes a voronoi diagram to compute area disorder

In order to evaluate the area of each Voronoi cell, transform each cell
to a polygon.

* Represent the couple `(ad, rfh)` in a graph, which gives a characterization of the Voronoi diagram.

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

See shapely Official Documentation for evaluating the area of a polygon.

31.1.3 Characterization of the Delaunay graph

If *L* denotes the set of the edge lengths of the Delaunay triangulation, the mean and the standard deviation of *L* can also give informations on the graph.

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

* Compute and display in a graph the point of coordinates $(\mu(L), \sigma(L))$, with $\mu$ representing the mean and $\sigma$ the standard deviation.

31.2

**Minimum spanning tree**

Definition from Wikipedia: a minimum spanning tree (MST) or minimum weight
spanning tree is a subset of the edges of a connected, edge-weighted undirected
graph that connects all the vertices together, without any cycles and with the
minimum possible total edge weight. That is, it is a spanning tree whose sum of
edge weights is as small as possible.

One of the methods to compute the MST is the Kruskal algorithm.

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

* Compute the MST.

* Compute $(\mu(L^*), \sigma(L^*))$ where $L^*$ denotes the set of the edge lengths of the MST.
31.3 Characterization of various point patterns

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

1. Generate $n$ condition Poisson point processes of 100 points each. For each realization, calculate the parameters $(AD, RFH)$, $(\mu(L), \sigma(L))$ and $(\mu(L^*), \sigma(L^*))$. Display these $n$ points in a 2D diagram in order to analyze the robustness of the quantification.

2. Generate 3 different point processes with regular, uniform and Gaussian dispersion. Display the different diagrams. Which one is the most discriminant?
<figure><img src="image_735.png" /></figure>

## 31.4. Python correction

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

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

from scpy. spatial import Voronoi, voronoi_plot_2d, Delaunay, distance

2 import numpy as np

4 import matplotlib. pyplot as plt
from shapely. import geometry # for polygons, area and perimeter

6 from scpy. sparse. import csgraph # minimum spanning tree

### 31.4.1 Random tessellations

Random tessellations are generated following normal standard and uniform distribution. A regular pattern is also employed. They are illustrated in Fig.31.2.

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

def dist_poisson (N=100):
    points = np. random. rand(N, 2)
    return points

4
def dist_gaussienne (N=100):
    points = np. random. randn(N, 2)
    return points

8
def dist fu regular (N=100):
10
c = np. floor (np. sqrt (N));
x2, y2 = np. meshgrid (range (int (c)), range (int (c)));
12 points = np. vstack ([x2. ravel () , y2. ravel ()])
return points. transpose ();

<figure><img src="image_739.png" /><figcaption>Figure 31.2: Different point patterns.</figcaption></figure>
### 31.4.2 Voronoi diagram and analysis

The Voronoi diagram is simply generated via the following command, with `points` being generated with the previous functions:

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

```python
vor = Voronoi(points);
```

The two characterization functions RFH and AD are defined on the Voronoi cells.

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

```python
1 def RFH(vor):
    """
3     Evaluates Round Factor Homogeneity from voronoi diagram
    """
5     rfs = [];
    for cell in vor.regions:
        if cell and -1 not in cell:
            poly = geometry.Polygon([((voranny vertices [p-1]) for p in cell)]);
        rfs.append(4*np.pi*poly.area/(poly.length**2));
    res = 1 - np.std(rfs) / np.mean(rfs);
    return res;
```

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

```python
1 def AD(vor):
    """
3     Evaluates Area Disorder from voronoi diagram
    """
5     areas = [];
    for cell in vor.regions:
        if cell and -1 not in cell:
            poly = geometry.Polygon([((vor.anny vertices [p-1]) for p in cell)]);
        areas.append(poly.area);
    res = 1 - 1/(1+np.std(areas) / np.mean(areas));
    return res;
```

### 31.4.3 Delaunay triangulation and minimum spanning tree

The Delaunay triangulation is computed with
<figure><img src="image_743.png" /></figure>

1 tri = delaunay(points);

Then, the characterization of the triangulation is done by measuring the dis-
tances of the edges.

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

1 def characterization ( tri ) : """ Characterization of the Delaunay triangulation (mean and std dev of edges) """ 5 M = triToMat( tri ); m = np.mean(M[M>0]); 7 s = np.std(M[M>0]); return m, s; 9 def triToMat( tri , value =0.) : """ Transforms the triangulation into a matrix representation , for simplicity """ 13 M = np. full (( tri .npoints , tri .npoints ), value ); 15 d = distance .pdist( tri .points ); distances = distance .squareform(d); 17 for s in tri .simplices : M[s[0], s [1]] = distances [s [0], s [1]]; 19 M[s[1], s [2]] = distances [s [1], s [2]]; M[s[2], s [0]] = distances [s [2], s [0]]; 21 return M;

31.4.4 Characterization of different realizations

n realizations of the two distributions (uniform and Gaussian) are simulated. Then,
the Voronoi diagram and the Delaunay triangulation are computed and character-
ized. The results are presented in Fig.31.3.
Figure 31.3: Characterization of several point processes. Each color represent a different process, and each point represent one realization. These characterizations can enhance a difference between the spatial distributions of the points.

<figure><img src="image_745.png" /><figcaption>(b) Characterization of the minimum spanning tree of the Delaunay triangulation by the mean and standard deviation of the lengths of the MST.</figcaption></figure>

<figure><img src="image_746.png" /><figcaption>(c) AD and RFH on the Voronoi diagram.</figcaption></figure>

<figure><img src="image_747.png" /></figure>
# Part V Image Characterization and Pattern Analysis
<figure><img src="image_748.png" /></figure>
<figure><img src="image_749.png" /></figure>

32 Integral Geometry

This tutorial aims to characterize objects by measurements from integral
geometry.

The different processes will be applied on the following synthetic image:

<figure><img src="image_750.png" /><figcaption>Figure 32.1: X</figcaption></figure>

32.1 Cell configuration

The spatial support of an image can be covered by cells associated with pixels. A cell (square of interpixel distance size) is composed of 1 face, 4 edges and 4 vertices. Either a cell is centered in a pixel (intrapixel cell) or a cell is constructed by connecting pixels (interpixel cell). The following figure shows the two possible representations of an image with 16 pixels:

<figure><img src="image_751.png" /><figcaption>Figure 32.2: Pixels representation.</figcaption></figure>

We respectively denote *f*<sup>intra</sup> (resp. *f*<sup>inter</sup>), *e*<sup>intra</sup> (resp. *e*<sup>inter</sup>) and *v*<sup>intra</sup> (resp. *v*<sup>inter</sup>) the number of faces, edges and vertices for the intrapixel (resp. interpixel) cell configuration Using these two configurations, the measurements from integral
geometry (area A, perimeter P, Euler number $\chi_8$ or $\chi_4$) can be computed as:

$$A = f^{\text{intra}} = v^{\text{inter}} \qquad (32.1)$$

$$P = -4f^{\text{intra}} + 2e^{\text{intra}} \qquad (32.2)$$

$$\chi_8 = v^{\text{intra}} - e^{\text{intra}} + f^{\text{intra}} \qquad (32.3)$$

$$\chi_4 = v^{\text{inter}} - e^{\text{inter}} + f^{\text{inter}} \qquad (32.4)$$

?

1. Count manually the number of faces, edges and vertices of the image X for the two configurations (Intra- and Inter-pixel) of Fig. 32.1.

2. Deduce the measurements from integral geometry (Eq. 32.1-32.4).

## 32.2 Neighborhood configuration

In order to efficiently calculate the number of vertices, edges and faces of the object, the various neighborhood configurations (of size 2x2 pixels) of the original binary image X are firstly determined. Each pixel corresponds to a neighborhood configuration $\alpha$. Thus, sixteen configurations are possible, presented in Tab. 32.1.

<table><caption>Table 32.1: Neighborhood configurations.</caption><tr><td></td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>1</td></tr><tr><td></td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td><td>1</td></tr><tr><td>α</td><td>0</td><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td></tr><tr><td></td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td></td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td></tr><tr><td>α</td><td>8</td><td>9</td><td>10</td><td>11</td><td>12</td><td>13</td><td>14</td><td>15</td></tr></table>

Thereafter, each configuration contributes to a known number of vertices, edges and faces (Tab. 32.2). To determine the neighborhood configurations of all the pixels, an efficient algorithm effective consists in convolving the image X by a mask F, whose values are powers of two, and whose origin is the top-left pixel:

$$F = \begin{pmatrix} 1 & 4 \\ 2 & 8 \end{pmatrix}$$

The resulting image is $X * F$. Notice that this image X has no pixel touching the borders, your code should ensure this (for example by padding the array with zeros). In this way, the histogram h of $X * F$ gives the distribution of the neighborhood configurations from image X. And each configuration contributes to a known
number of vertices, edges and faces:

$$v = \sum_{\alpha=0}^{15} v_{\alpha}h(\alpha) \qquad (32.5)$$

$$e = \sum_{\alpha=0}^{15} e_{\alpha}h(\alpha) \qquad (32.6)$$

$$f = \sum_{\alpha=0}^{15} f_{\alpha}h(\alpha) \qquad (32.7)$$

The following table gives the values of $v_\alpha$, $e_\alpha$ and $f_\alpha$ for each cell configuration.

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

1. Compute the distribution of the neighborhood configurations from image X.

2. Deduce the number of vertices, edges and faces for each cell representation and compare these values with the previous (manually computed) results.

<table><caption>Table 32.2: Contributions of the neighborhood configurations to the computation of $v$, $e$ and $f$</caption><thead><tr><th colspan="15">intrapixel cells</th></tr><tr><th>α</th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th><th>14</th><th>15</th></tr></thead><tbody><tr><th>f<sub>α</sub></th><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td></tr><tr><th>e<sub>α</sub></th><td>0</td><td>2</td><td>1</td><td>2</td><td>1</td><td>2</td><td>2</td><td>2</td><td>0</td><td>2</td><td>1</td><td>2</td><td>1</td><td>2</td><td>2</td><td>2</td></tr><tr><th>v<sub>α</sub></th><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr></tbody></table>

<table><thead><tr><th colspan="15">interpixel cells</th></tr><tr><th>α</th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th><th>14</th><th>15</th></tr></thead><tbody><tr><th>f<sub>α</sub></th><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td></tr><tr><th>e<sub>α</sub></th><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>2</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>2</td></tr><tr><th>v<sub>α</sub></th><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td></tr></tbody></table>

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

## Crofton perimeter

The Crofton perimeter could be computed from the number of intercepts in different random directions. In discrete case, only the directions $0, \pi/4, \pi/2$ and $3\pi/4$ are considered; they are selected according to the desired connectivity.

The number of intercepts are denoted $i_0, i_{\pi/4}, i_{\pi/2}$ and $i_{3\pi/4}$ for the orientation angles $0, \pi/4, \pi/2$ and $3\pi/4$ respectively.
In this way, the Crofton perimeter (in 4 and 8 connexity) is defined in discrete case as:

$$P_4 = \frac{\pi}{2} (i_0 + i_{\pi/2}) \qquad (32.8)$$

$$P_8 = \frac{\pi}{4} \left( i_0 + \frac{i_{\pi/4}}{\sqrt{2}} + i_{\pi/2} + \frac{i_{3\pi/4}}{\sqrt{2}} \right) \qquad (32.9)$$

These perimeter measurements can be computed from the neighborhood configurations of the original image:

$$P_4 = \sum_{\alpha=0}^{15} P_{\alpha}^4 h(\alpha) \qquad (32.10)$$

$$P_8 = \sum_{\alpha=0}^{15} P_{\alpha}^8 h(\alpha) \qquad (32.11)$$

with the following weights $P_{\alpha}^4$ and $P_{\alpha}^8$ of these linear combinations (Tab. 32.3):

<table><caption>Table 32.3: Weights for the computation of the Crofton perimeter</caption><tr><td>α</td><td>0</td><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td></tr><tr><td>P<sup>4</sup><sub>α</sub></td><td>0</td><td>π/2</td><td>0</td><td>0</td><td>0</td><td>π/2</td><td>0</td><td>0</td></tr><tr><td>P<sup>8</sup><sub>α</sub></td><td>0</td><td>π/4 (1 + 1/√2)</td><td>π/4√2</td><td>π/2√2</td><td>0</td><td>π/4 (1 + 1/√2)</td><td>0</td><td>π/4√2</td></tr><tr><td>α</td><td>8</td><td>9</td><td>10</td><td>11</td><td>12</td><td>13</td><td>14</td><td>15</td></tr><tr><td>P<sup>4</sup><sub>α</sub></td><td>π/2</td><td>π</td><td>0</td><td>0</td><td>π/2</td><td>π</td><td>0</td><td>0</td></tr><tr><td>P<sup>8</sup><sub>α</sub></td><td>π/4</td><td>π/2</td><td>π/4√2</td><td>π/4√2</td><td>π/4</td><td>π/2</td><td>0</td><td>0</td></tr></table>

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

Evaluate the perimeters $P_4$ and $P_8$ with the previous formula on Fig.32.1
<figure><img src="image_755.png" /></figure>

## 32.4. Python correction

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

### 32.4.1 Cell configurations

The following values are reported:

$$ 
\begin{aligned}
f^{\text{intra}} &= 50 \\
e^{\text{intra}} &= 158 \\
v^{\text{intra}} &= 107
\end{aligned}
 $$

$$ 
\begin{aligned}
f^{\text{inter}} &= 4 \\
e^{\text{inter}} &= 42 \\
v^{\text{inter}} &= 50
\end{aligned}
 $$

Then, it is easy to compute the following values:

$$ 
\begin{aligned}
A &= f^{\text{intra}} = 50 \\
P &= -4f^{\text{intra}} + 2e^{\text{intra}} = 116 \\
\chi_8 &= v^{\text{intra}} - e^{\text{intra}} + f^{\text{intra}} = -1 \\
\chi_4 &= v^{\text{inter}} - e^{\text{inter}} + f^{\text{inter}} = 12
\end{aligned}
 $$

### 32.4.2 Neighborhood configuration

The configuration is computed using the convolution function `scipy.signal convolve2d`. The histogram of the different configurations is presented in Fig.32.3.

Be aware that this algorithm works if there is no object pixel touching the borders of the image. The example image is in this case, but you can ensure this property by padding 0 values around the image:

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

X = np.pad(I, ((1,1), (0,0)), mode='constant');
<figure><img src="image_758.png" /></figure>

1. # Neighborhood configuration

F = np.array([[0, 0, 0], [0, 1, 4], [0, 2, 8]]);

3. XF = signal(convolve2d(X,F,mode='same'));
edges = np.arange(0, 17, 1);

5. h,edges = np.histogram(XF[:], bins=edges);
plt.figure()

7. plt.bar eds[0:-1], h);

7. plt. title ("Histogram of the different configurations ")

9.

 plt.show()

<figure><img src="image_759.png" /><figcaption>Figure 32.3: Distribution of the different neighborhood configurations.</figcaption></figure>

The Minkowski functionals are computed using the cells contributions:

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

1. # Computation of the functionals

f_intra = [0,1,0,1,0,1,0,1,0,1,0,1,0,1];

3 e_intra = [0,2,1,2,1,2,2,0,2,1,2,1,2,2,2];

v_intra = [0,1,1,1,1,1,1,1,1,1,1,1,1,1,1];

5 EulerNb8 = np.sum(h*v_intra - h*e_intra + h*f_intra)

f_inter = [0,0,0,0,0,0,0,0,0,0,0,0,0,1];

7 e_inter = [0,0,0,1,0,1,0,2,0,0,0,1,0,1,0,2];

v_inter = [0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1];

Then, the values are easily verified.

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

EulerNb4 = np.sum(h*f_inter - h*e_inter + h*v_inter)

2 Area = sum(h*f_intra)

Perimeter = sum(-4*h*f_intra + 2*h*e_intra)

4 print("E_4 :{0}, A :{1}, P :{2} ".format(EulerNb4, Area, Perimeter));
<figure><img src="image_762.png" /></figure>

E_4:12, A:50, P:116

32.4.3 Crofton perimeter

The Crofton perimeter is a good approximation of a perimeter. One should notice
that there is no definitive solution to perimeter evaluation. The Crofton perimeter
is approximated in 2 or 4 directions, denoted $P_4$ and $P_8$ with a reference to the
connectivity.

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

1 # Crofton perimeter

P4 = [0, np.pi /2,0,0,0, np.pi /2,0,0, np.pi /2, np.pi ,0,0, np.pi /2, np.pi ,0,0];

3 Perimeter4 = sum(h\*P4)

P8 = [0, np.pi /4*(1+1/(np Pag(2))), np.pi /(4* np. sqrt (2) ), np pi /(2* np. sqrt (2) ), 0, np.
    ↪→ pi /4*(1+1/( np. sqrt (2) )) ,0, np. pi /(4* np. sqrt (2) ), np. pi /4, np. pi /2, np. pi /(4*
    ↪→ np. sqrt (2) ), np. pi /(4* np. sqrt (2) ), np. pi /4, np. pi /2,0,0];

5 Perimeter8 = sum(h\*P8)

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

1 Perimeter4: 91.10618695410399, Perimeter8: 77.76399477870015
<figure><img src="image_765.png" /></figure>
<figure><img src="image_766.png" /></figure>

33 Topological Description

The objective of this tutorial is to classify all foreground points of a binary
image according to their topological signification: interior, isolated, border....
The reader can refer to [6]for more details.

The different processes will be realized on the following binary image.

<figure><img src="image_767.png" /><figcaption>Figure 33.1: Test</figcaption></figure>

**Preliminary definitions:**

• y is 4-adjacent to x if |y1 - x1| + |y2 - x2| ≤ 1.

• y is 8-adjacent to x if $\max(|y_1 - x_1|, |y_2 - x_2|) \le 1$.

• $V_4(x) = \{y : y$ is 4-adjacent to $x\}; V_4^*(x) = V_4(x)\setminus\{x\}$.

• $V_8(x) = \{y : y$ is 8-adjacent to $x\}; V_8^*(x) = V_8(x)\setminus\{x\}$.

Figure 33.2: Different neighborhoods. By convention, pixels in white are of value 1, in black of value 0.

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

• a n-path is a point sequence ($x_0, \dots, x_k$) with $x_j$ n-adjacent to $x_{j-1}$ for $j = 1, \dots, k$.

• two points $x, y \in X$ are n-connected in $X$ if there exists a n-path ($x = x_0, \dots, x_k = y$) such that $x_j \in X$. It defines an equivalence relation.

• the equivalence classes of the previous binary relation are the n-components of X.
## 33.1 Connectivity numbers

Let $Comp_n(X)$ be the number of n-components ($n = 4$ or $n = 8$ within the selected topology $V_4$ or $V_8$) of the set $X$ of foreground points (object). We define the following set:

$$CAdj_n(x, X) = \{C \in Comp_n(X) : C \text{ is n-adjacent to } x\}$$

We select the 8-connectivity for the set $X$ of foreground points (object) and the 4-connectivity for the complementary $\bar{X}$.

Warning: the definition of $CAdj_n$ introduces the n-adjacency to the central pixel $x$. In the case of the following configuration (Fig.33.3), $T_8 = 2$, $\bar{T}_8 = 2$ and $TT_8 = 3$.

<figure><img src="image_769.png" /><figcaption>Figure 33.3: The pixel in the bottom right corner is not C-adjacent-4 to $x$.</figcaption></figure>

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

1. Create a function for determining the connectivity number: $T_8(x, X) = #CAdj_8(x, V_8^*(x) \cap X)$,

2. Create a function for determining the connectivity number: $\bar{T}_8(x, X) = #CAdj_4(x, V_8^*(x) \cap \bar{X})$,

3. Create a function for determining the number: $TT_8(x, X) = #(V_8^*(x) \cap X)$.

Test these functions on some foreground points of the image 'test'.

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

Use `scipy.ndimage.measurements`.

## 33.2 Topological classification of binary points

From the connectivity numbers $T_8(x, X)$, $\bar{T}_8(x, X)$ and $TT_8(x, X)$, it is possible to classify a foreground point $x$ within the binary image $X$ according to its topological
signification:

<table><thead><tr><th>T<sub>8</sub>(x, X)</th><th>T&#773;<sub>8</sub>(x, X)</th><th>TT<sub>8</sub>(x, X)</th><th>Type</th></tr></thead><tbody><tr><td>0</td><td>1</td><td></td><td>isolated point</td></tr><tr><td>1</td><td>0</td><td></td><td>interior point</td></tr><tr><td>1</td><td>1</td><td>&gt; 1</td><td>border point</td></tr><tr><td>1</td><td>1</td><td>1</td><td>end point</td></tr><tr><td>2</td><td>2</td><td></td><td>2-junction point</td></tr><tr><td>3</td><td>3</td><td></td><td>3-junction point</td></tr><tr><td>4</td><td>4</td><td></td><td>4-junction point</td></tr></tbody></table>

The following Fig. 33.4 shows the classification of 4 points.

<figure><img src="image_772.png" /><figcaption>Figure 33.4: Points configurations.</figcaption></figure>

With the help of this table, classify the points of the image 'test'.
<figure><img src="image_773.png" /></figure>

## 33.3. Python correction

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

### 33.3.1 Topological description

The operations are not difficult, except that the *CAdj*<sub>4</sub> should be coded carefully.

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

def nc(A):

# A : block 3x3, binary

# complementary set of A

invA=1-A;

# neighborhoods

V8= np. ones((3,3)).astype(int);

V8_star = np. copy(V8);

V8_star [1,1] = 0;

V4 = np. array([ [0, 1, 0], [1, 1, 1], [0, 1, 0] ] ). astype( int );

# intersection is done by the min operation

X1 = np. minimum(V8_star, A);

TT8 = np. sum(X1);

L, T8 = mes. label (X1, structure =V8);

# The C-ajd-4 might introduce some problems if a pixel is not 4-connected

# to the central pixel

X2 = np. minimum(V8, invA);

Y = np. minimum(X2, V4);

X = morpho.reconstruction(Y, X2, selem =V4);

L, T8c = mes. label (X, structure =V4);

return T8, T8c, TT8

<figure><img src="image_776.png" /><figcaption>Figure 33.5: Extraction of the window centered in *x* = (1, 4), given as (row,col).</figcaption></figure>

### 33.3.2 Topological classification

The different types are given by the following code.
<figure><img src="image_777.png" /></figure>

1 def nc_type(X):
    # evaluates the connectivity numbers
    a, b, c = nc(X);
    if (a==0):
        y=1; # isolated point
    if ((a==1) and (b==1) and (c>1)):
        y=5; # border point
    if (b==0):
        y=7; # interior point
    if ((a==1) and (b==1) and (c==1)):
        y=6; # end point
    if (a==2):
        y=2; # 2-junction point
    if (a==3):
        y=3; # 3-junction point
    if (a==4):
        y=4; # 4-junction point
    return y;

In order to perform the classification of all pixels of an image, one has to loop
over all the pixels, except the ones at the sides. The results are presented in Fig.33.6

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

def classification (A):
    # for the whole image
    m, n = A.shape
    B = np(zeros((m,n)));
    for i in range(1, m-1):
        for j in range(1, n-1):
            if A[i, j] > 0:
                X = A[i-1:i+2, j-1:j+2];
                B[i, j] = nc_type(X);
    return B
<figure><img src="image_779.png" /><figcaption>Figure 33.6: Classification of all the pixels of the original image.</figcaption></figure>
<figure><img src="image_780.png" /></figure>

## 34 Image Characterization

This tutorial aims to characterize objects by geometrical and morphometrical measurements.

The different processes will be applied on synthetic images as well as images from the Kimia database [1, 39]:

<figure><img src="image_781.png" /><figcaption>(a) Bat.</figcaption></figure>

<figure><img src="image_782.png" /><figcaption>(b) Camel.</figcaption></figure>

<figure><img src="image_783.png" /><figcaption>(c) Ray.</figcaption></figure>

### 34.1 Perimeters

We are going to calculate the perimeter using the Crofton formula. This formula consists in integrating the intercept number of the object with lines of various orientation and positions. Its expression in the 2-D planar case is given by:

$$P(X) = \pi \int \chi(X \cap L)dL$$

where the Euler-poincaré characteristic $\chi$ is equal to the number of connected components of the intersection of $X$ with a line $L$.

In the discrete case, the Crofton formula can be estimated by considering the intercept numbers for the horizontal $i_0$, vertical $i_{\pi/2}$ and diagonal orientations $i_{\pi/4}$ and $i_{3\pi/4}$ as:

$$P(X) = \pi \times \frac{1}{4} \left( i_0 + \frac{1}{\sqrt{2}} i_{\pi/4} + i_{\pi/2} + \frac{1}{\sqrt{2}} i_{3\pi/4} \right)$$
<figure><img src="image_784.png" /></figure>

1. Calculate the intercept number of a binary object from the Kimia database with lines oriented in the following four directions: 0, $\pi/4$, $\pi/2$, $3\pi/4$.

2. Deduce the value of the Crofton perimeter.

3. Compare the result with the classical perimeter functions.

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

See perimeter from skimage.

## 34.2 Feret Diameter

The Feret diameter (a.k.a. the caliper diameter) is the length of the projection of an object in one specified direction Fig.34.1.

<figure><img src="image_786.png" /><figcaption>Figure 34.1: Feret diameter of the object in horizontal direction.</figcaption></figure>

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

1. Calculate the projections in different directions of a binary object from the Kimia database.

2. Deduce the minimum, maximum and mean value of the Feret diameters.

## 34.3 Circularity

We want to know if the object X is similar to a disk. For that, we define the following measurement (circularity criterion):

$$circ(X) = \frac{4\pi A(X)}{P(X)^2}$$

where A(X) and P(X) denote the area and perimeter of the object X.
<figure><img src="image_788.png" /></figure>

1. Show that the circularity of a disc is equal to 1.

2. Generate an array representing an object as a discrete disc.

3. Calculate its circularity and comment the results.

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

Use numpy.meshgrid.

## 34.4 Convexity

We want to know if the object X is convex. For that, we define the following measurement:

$$conv(X) = \frac{A(X)}{A(CH(X))}$$

where CH(X) denotes the filled convex hull of the object X.

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

1. Compute the convex hull of a pattern from the Kimia database.

2. Evaluate the area of the filled convex hull.

3. Deduce the convexity of the pattern.

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

See ConvexHull from spacy . spatial .
<figure><img src="image_792.png" /></figure>

## 34.5. Python correction

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

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

import matplotlib.pyplot as plt

2 import numpy as np
from scipy import ndimage

4 from scipy import signal
from scipy import misc

6 from scipy import spatial
from skimage import measure

### 34.5.1 Perimeters

The convolution is used here as an easy way of getting the borders of the object in one direction, i.e. counting the number of intercepts. This method is efficient because, for one given direction, a high number of lines are considered. Keep in mind that only 4 orientations are used in the proposed code.

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

1 def countIntercepts(I, h):
   B = np.abs( signal .convolve2d(I, h, mode='same'));

3 n = np.sum(B) / 2;
return n;

5
def perimCrofton(I):
    """

7
    Approximate the Crofton perimeter with 4 directions
9
    I is the input binary image
9
    return the perimeter, float value
    """

11
# defines an orientation
13 h = np.array([[ -1, 1]]);

15
n2 = countIntercepts(I, h `.transpose()`);

17
h = np.array([[1, 0], [0, -1]]);
19 n3 = countIntercepts(I, h `.transpose()`);

21
h = np.array([[0, 1], [-1, 0]]);
23 n4 = countIntercepts(I, h `.transpose()`);

25
perim_Crofton = np.pi/4 * (n1+n2+ 1/np .sqrt(2) *(n3+n4));
return perim_Crofton;
<figure><img src="image_796.png" /></figure>

1. Crofton perimeter: 1305.50015965
   Classical perimeter in N4: 1374.19509294

3. Classical perimeter in N8: 1642.82842712

### 34.5.2 Feret diameter

The code consists in rotating the image and evaluating the projected diameter. The rotation function can interpolate the pixel, the nearest method is thus required. The function `np.max` directly performs the projection in one direction.

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

```python
def feretDiameter(I):
    """
    I: input binary image
    Returns min, max and mean Feret diameter, which is the length of the
    projected object in one direction
    """
diameter = []

for angle in range(180):
    I2 = ndimage.interpolation Rotate (I, angle, mode='nearest');
    I3 = I2.max (axis=0);
    diameter.append (np.sum (I3 > 0));

return np.min (diameter), np.max (diameter), np.mean (diameter);
```

For the camel image [1], the Feret diameters are:

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

1. min, Max, average Feret diameters: 182 325 266.05

### 34.5.3 Circularity

For a disk, the perimeter is $\pi \cdot D$ and the surface is $\pi \cdot \frac{D^2}{4}$, which shows that the circularity criterion has value 1 for a disk.

In order to generate a binary image containing a disk, one simple way is to use the formula: $(x - x_0)^2 + (y - y_0)^2 \le R^2$. The efficient way to do this is to use `np.meshgrid`.
<figure><img src="image_799.png" /></figure>

```python
def disk(t, r):
    """
    Generates a binary array representing a disk, centered, of radius r
    an array of size [2t,2t] is generated
    """
    x = np.arange(-t, t, 1);
    X,Y = np.meshgrid(x, x);
    I = (X**2 + Y**2) <= r **2;
    return I;
```

2

4

6

8

The circularity uses the perimeters as evaluated earlier.

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

1 def circularity (I):
"""

3
    Circularity criterion
    4*pi*A/P**2

5
    returns crofton and classic
    """

7
    P = perimCrofton(I);
    print("Perimeter by crofton: ", P)

9
    A = np.sum(I);

C = np.pi * 4 * A / P**2;

11

p = measure.perimeter(I, neighbourhood=4);

13

print("Usual perimeter: ", p)

c = np.pi * 4 * A / p**2;

15

return C, c;

This gives for the camel image:

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

1 Perimeter by crofton: 2514.43300853

Usual perimeter: 2649.36074863

3 circularity by crofton: 0.999019146137

circularity usual 0.89985338478

### 34.5.4 Convexity

The convex hull is computed with the `scipy` spatial tools. The following function plots the result. Notice that the coordinates of an image and the coordinates of an array are different, and one has to flip the image array in order to display both data in the same figure.
<figure><img src="image_802.png" /><figcaption>Figure 34.2: Convex hull of the image Camel.</figcaption></figure>

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

def convexity(I):

2

### """
Evaluates the convexity criterion
I is a binary image (np array)
return convexity
### """

4

6

```python
# be careful that coordinates between images and arrays are inversed
# thus, image is flipped before extracting points coordinates
points = npurlatespose(npwhere(npflip(I,0)));
hull = spatiaIConvexHull(points);
A = npsum(I);
Ah = hull .volume;
```

10

12

14

```python
plt .figure ()
for simplex in hull .simplices :
    plt .plot (points [simplex, 1], points [simplex, 0], 'k-')
    plt .axis ( 'equal')
    plt .show()
return A/Ah;
```

16

18

This criterion is between 0 and 1. This is illustrated in Fig.34.2.

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

1

convexity of Camel: 0.644014426566
<figure><img src="image_805.png" /></figure>
<figure><img src="image_806.png" /></figure>

35 Shape Diagrams

The objective is to study some shape diagrams and the possibility to define
properties that may be useful in order to distinguish the different objects.

Shape diagrams are representations of single shapes (connected compact sets, see [30, 31, 32]) as points in the 2-D unit square plane. They are based on inequalities between 6 geometrical measurements: area A, perimeter P, radius of the inscribed circle r, radius of the circumscribed circle R, minimum Feret diameter ω and maximum Feret diameter d. In this way, the morphometrical functionals used in the different shape diagrams are normalized ratios of such geometrical functionals. The following table shows the morphometrical functionals for non-convex sets:

<table><caption>Table 35.1: Morphometrical functionals.</caption><thead><tr><th>Geometrical functionals</th><th>Inequalities</th><th>Morphological functionals</th></tr></thead><tbody><tr><td>r, R</td><td>r ≤ R</td><td>r/R</td></tr><tr><td>ω, R</td><td>ω ≤ 2R</td><td>ω/2R</td></tr><tr><td>A, R</td><td>A ≤ πR<sup>2</sup></td><td>A/πR<sup>2</sup></td></tr><tr><td>d, R</td><td>d ≤ 2R</td><td>d/2R</td></tr><tr><td>r, d</td><td>2r ≤ d</td><td>2r/d</td></tr><tr><td>ω, d</td><td>ω ≤ d</td><td>ω/d</td></tr><tr><td>A, d</td><td>4A ≤ πd<sup>2</sup></td><td>4A/πd<sup>2</sup></td></tr><tr><td>R, d</td><td>√3R ≤ d</td><td>√3R/d</td></tr><tr><td>r, P</td><td>2πr ≤ P</td><td>2πr/P</td></tr><tr><td>ω, P</td><td>πω ≤ P</td><td>πω/P</td></tr><tr><td>A, P</td><td>4πA ≤ P<sup>2</sup></td><td>4πA/P<sup>2</sup></td></tr><tr><td>d, P</td><td>2d ≤ P</td><td>2d/P</td></tr><tr><td>R, P</td><td>4R ≤ P</td><td>4R/P</td></tr><tr><td>r, A</td><td>πr<sup>2</sup> ≤ A</td><td>πr<sup>2</sup>/A</td></tr><tr><td>r, ω</td><td>2r ≤ ω</td><td>2r/ω</td></tr></tbody></table>
Figure 35.1: The different processes will be applied on images from the Kimia database [1, 39].

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

## 35.1 Geometrical functionals

?

Code functions in order to evaluate the different parameters:

* the area, Crofton perimeter and Feret diameters have been already presented in tutorial 32;

* the radius of the inscribed circle can be defined from the ultimate erosion of a set.

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

The function scipy.ndimage.morphology.distance_transform_cdt computes the distance map of a binary image (chamfer distance transform).

## 35.2 Morphometrical functionals

?

Code and evaluate some of the morphometrical functionals listed in the table 35.1. Note that each of them has a physical meaning, e.g. $\frac{4\pi A}{P^2}$ (circularity), $\frac{4A}{\pi d^2}$ (roundness), $2\omega/P$ (thinness).
35.3 Shape diagrams

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

* Visualize the different shape diagrams for all the images (from the Kimia database) within the three classes 'apple', 'bone' and 'camel'. The Fig.35.2 illustrates the result for the shape diagram ($x = 2r/d, y = P/\pi d$).

* Which shape diagram is the most appropriate for the discrimination of such objects?

<figure><img src="image_810.png" /><figcaption>Figure 35.2: Example of a shape diagram.</figcaption></figure>

35.4 Shape classification

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

* Use a K-means clustering method for automatic classification of such shapes.

* Propose a method to quantify the classification accuracy for each shape diagram.
<figure><img src="image_812.png" /></figure>

## 35.5. Python correction

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

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

import numpy as np

2

import spipy.ndimage

4

import imageio # imread and imwrite

import matplotlib.pyplot as plt

6

import skimage # some geometrical descriptors

8

# for reading files

import glob

10

from sklearn.cluster import KMeans

### 35.5.1 Geometrical functionals

The Crofton perimeter is defined by multiple projections, as well as the Feret diameter. Be careful while performing the rotation of the object (as it is a binary object, the interpolation method could introduce non integer values).

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

1

def crofton_perimeter(I):
    """ Computation of crofton perimeter
    """

3

inter = [];
h = np.array([[1, -1]]);

5

for i in range(4):
    II = np.copy(I);

7

I2 = spipy.misc.imrotate(II, 45*i, interp='nearest');

9

I3 = spipy.ndimage(convolve(I2, h));

11

inter.append(np.sum(I3>100));

13

crofton = np.pi/4. * (inter[0]+inter[2] + (inter[1]+inter[3])/np

15

return crofton
<figure><img src="image_816.png" /></figure>

def feret_diameter(I):

2

### Computation of the Feret diameter

4

minimum: d (meso-diameter)

maximum: D (exo-diameter)

6

Input: I binary image

8

### d = np.max(I.shape);

10

D = 0;

12

for a in np.arange(0, 180, 30):

I2 = scipy.misc.imrotate(I, a, interp='nearest');

14

F = np.max(I2, axis=0);

measure = np.sum(F>100);

16

if (measure<d):

18

d = measure;

if (measure>D):

20

D = measure;

return d,D;

The inscribed circle is just the maximum of the distance transform inside the object. The distance map for an image apple is shown in Fig.35.3a.

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

1

def inscribedRadius(I):

### computes the radius of the inscribed circle

3

### dm = scipy.ndimage.morphology.distance_transform_cdt(I>100);

5

radius = np.max(dm);

### return radius;

7
<caption>Figure 35.3: Illustration of the computation of two shape parameters.</caption>

(a) Distance map of an object apple. The inverse is actually displayed in order to see correctly the progression. The maximum of the distance map is the radius of the inscribed circle.

<figure><img src="image_818.png" /><figcaption>(b) Circumscribed circle.</figcaption></figure>

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

The smallest enclosing circle is computed by using the code from the Project Nayuki<fnref n="1" /> published under the GNU Lesser General Public License. The result is presented in Fig.35.3b.

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

1 def circumCircle(I):

3 this version uses a function provided by Project Nayuki
under GNU Lesser General Public License

5 points = np.argwhere(I > 100);

7 c = smallestenclosingcircle .make_circle(points);
return c;

35.5.2 Shape diagrams

The shape diagrams are constructed by reading all the images and computing the
shape descriptors.

<footnote n="1">¹https://www.nayuki.io/page/smallest-enclosing-circle</footnote>
<figure><img src="image_821.png" /></figure>

```python
def diagrams():
    name = ['apple-*.bmp', 'Bone-*.bmp', 'camel-*.bmp'];
    elongation = [];
    thinness = [];
    roundness = [];
    z = [];

    for pattern in name:
        namesList = glob.glob
        for fichier in namesList:
            I = imageio.imread(fichier);
            radius = inscribedRadius(I);
            d, D = feret_diameter(I);
            crofton = crofton_perimeter(I);

            elongation.append(d/D);
            thinness.append(2*radius / D);
            roundness.append(4*np.sum(I>100)/(np.pi * D**2));
            z.append(crofton / (np.pi * D));
```

The display of the different plots is just a use of the function `plt`. `plot`.

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

```python
plt plot(elongation [0:20], thinness [0:20], "o", label='Apple')
2

plt plot(elongation [20:40], thinness [20:40], "+", label='Bone')

3

plt plot(elongation [40:60], thinness [40:60], ".", label='Camel')

4

plt legend(name)

5

plt show

6
evaluateQuality(elongation, thinness);

8

plt figure ();

9

plot (z [0:20], roundness [0:20], "o", label='Apple')

10

plot (z [20:40], roundness [20:40], "+", label='Bone')

11

plot (z [40:60], roundness [40:60], ".", label='Camel')

12

plot legend(name)

13

plot show

14
evaluateQuality(z, roundness);

16

plot figure ();

17

plot (thinness [0:20], z [0:20], "o", label='Apple')

18

plot (thinness [20:40], z [20:40], "+", label='Bone')

19

plot (thinness [40:60], z [40:60], ".", label='Camel')

20

plot legend(name)

21

plot show

22
evaluateQuality(thinness, z);
```
35.5.3 Shape classification

The following code evaluates the quality by comparing the known class of the
shape with the segmented (via the kmeans method) class. The result is illustrated
in Fig.35.4 with an accuracy of 98.3%.

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

def evaluateQuality(x, y):

2

global i

n = 3;

4

k Means(init='k-means++', n clusters=n)

X = np .asarray (x);

6

Y = np .asarray (y);

pts = np .stack ((X, Y));

8

pts = pts .T;

#print (pts)

10

kmeans.fit (pts);

kmeans labels = kmeans.labels_;

12

kmeans_clusterCenters = kmeans.clusterCenters_;

# plot

14

fig = plt .figure ()

colors = ['#4EACC5', '#FF9C34', '#4E9A06']

16

# KMeans

for k, col in zip (range (n), colors):

18

my_members = kmeans labels == k

cluster_center = kmeans_clusterCenters [k]

20

plt .plot (pts [my_members, 0], pts [my_members, 1], 'o',

markerfacecolor =col, markersize =6)

22

plt .plot ( cluster_center [0], cluster_center [1], 'o', markerfacecolor =col,

markeredgecolor ='k', markersize =12)

24

plt .title ('KMeans')

plt.show()

26

fig .savefig ("kmeans"+str (i))+".pdf ");

i +=1;

28

###

Evaluation of the quality : count the number of shapes correctly detected

30

###

accuracy = np .sum (kmeans labels [0:20] ==scipy .stats .mode (kmeans labels
       ↪ → [0:20]));

32

accuracy +=np .sum (kmeans labels [20:40] ==scipy .stats .mode (kmeans labels
                                                    ↪ → [20:40]));

accuracy +=np .sum (kmeans labels [40:60] ==scipy .stats .mode (kmeans labels
                                                    ↪ → [40:60]));

34

accuracy = accuracy / 60 * 100;

print('Accuracy: {0:.2 f}%' .format(accuracy));
<figure><img src="image_824.png" /><figcaption>Figure 35.4: Illustration of the accuracy of the classification from a k-means method. The k-means method is not necessarily the adapted to these data. The measured accuracy if of 98.3%.</figcaption></figure>
<figure><img src="image_825.png" /></figure>
<figure><img src="image_826.png" /></figure>

## 36 Freeman Chain Code

This tutorial is focused on shape representation by the Freeman chain code. This code is an ordered sequence of connected segments (of specific sizes and directions) representing the contour of the shape to be analyzed. The direction of each segment is encoded by a number depending on the selected connectivity (Fig. 36.1). In this example, the contour of the shape and its Freeman code are given in 4-connectivity.

<figure><img src="image_827.png" /><figcaption>Figure 36.1: Freeman chain code examples.</figcaption></figure>

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

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

### Notations:

* $y$ is 4-adjacent to $x$ if $|y_1 - x_1| + |y_2 - x_2| \le 1$.

* $y$ is 8-adjacent to $x$ if $\max(|y_1 - x_1|, |y_2 - x_2|) \le 1$.

* $N_4(x) = \{y : y$ is 4-adjacent to $x\}$; $N_4^*(x) = N_4(x)\{x\}$.

* $N_8(x) = \{y : y$ is 8-adjacent to $x\}$; $N_8^*(x) = N_8(x)\{x\}$.
<figure><img src="image_830.png" /><figcaption>(a) $V_4(x)$</figcaption></figure>

<figure><img src="image_831.png" /><figcaption>(b) $V_8(x)$</figcaption></figure>

## 36.1 Shape contours

Before determining the Freeman chain code of the shape, it is necessary to extract its contour.

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

1. Generate or load a simple shape as a binary image A.

2. Extract its contour $C_4(A)$ or $C_8(A)$ according to the 4-connectivity or the 8-connectivity, respectively:

$$x \in C_4(A) \Leftrightarrow \exists y \in N_8(x) \quad y \in cA$$

$$x \in C_8(A) \Leftrightarrow \exists y \in N_4(x) \quad y \in cA$$

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

### Informations

You can erode the object and subtract this erosion to it, with a structuring element that corresponds to $N_4$ or $N_8$ if you want to have 8 or 4 connectivity, respectively.

## 36.2 Freeman chain code

From the shape contours, the Freeman chain code can be calculated.

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

1. From the binary array of pixels, extract the first point belonging to the shape (from left to right, top to bottom).
<figure><img src="image_835.png" /></figure>

Informations

You can use np.argwhere to locate the first point.

2. From this initial point, determine the Freeman chain code c (counter-clockwise direction) using the $N_4$ or $N_8$ connectivity.

## 36.3 Normalization

The Freeman chain code is depending on the initial point and is not invariant to shape rotation. It is then necessary to normalize this code.

1. The first step consists in defining a differential code $d$ from the code $c$:

$$d_k = c_k - c_{k-1} (\text{mod } 4 \text{ or } 8)$$

In this example, $d = 3003131331300133031130$;

2. The second step consists in normalizing the code $d$. We have to extract the lowest number $p$ from all the cyclic translations of $d$.

In this example, $p = 0013303113030031313313$.

This Freeman chain code $p$ is then independent from the initial point and invariant by rotations of angles $k * \pi / 2$ radians ($k \in \mathbb{Z}$).

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

* Code above steps 1 and 2. Prototypes of functions are given.

* Validate your code by rotating the shape of $3\pi / 2$ radians.

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

def codediff(fcc, connectivity=8):

2 # computes differential code

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

Use numpy roll for circularly shifting the Freeman code array. Code a small function that tests all elements of array, one by one, in order to get the minimum of two arrays.
# 36.4 Geometrical characterization

## 36.4.1 Perimeter

From the Freeman chain code, it is possible to make different geometrical measurement of the shape.

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

Calculate the perimeter of the shape (take care of the diagonals).

## 36.4.2 Area

The area is evaluated by the following algorithm:

* The area is initialized to 0 and the parameter $B$ is initialized to 0.

* For each iteration on the Freeman chain code $c$, the area and the parameter value $B$ are incremented with the following rules:

<table><thead><tr><td>8-code</td><td>area</td><td>$B$</td></tr></thead><tbody><tr><td>0</td><td>-B</td><td>0</td></tr><tr><td>1</td><td>-B-0.5</td><td>1</td></tr><tr><td>2</td><td>0</td><td>1</td></tr><tr><td>3</td><td>B+0.5</td><td>1</td></tr><tr><td>4</td><td>B</td><td>0</td></tr><tr><td>5</td><td>B-0.5</td><td>-1</td></tr><tr><td>6</td><td>0</td><td>-1</td></tr><tr><td>7</td><td>-B+0.5</td><td>-1</td></tr></tbody></table>

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

Calculate the area of the shape in 8-connectivity. Note that the area does not correspond to the total number of pixels belonging to the shape.
<figure><img src="image_841.png" /></figure>

## 36.5. Python correction

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

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

import matplotlib .pyplot as plt

2 import numpy as np

from skimage.morphology import binary_erosion, disk, rectangle

4 from skimage.measurement import perimeter

### 36.5.1 Shape contours

To generate a simple object, here is an example:

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

A = np(zeros((20,20))).astype('bool');

2 A[4:14, 9:17] = True;
A[1:18,11:16] = True;

The contours are computed in 4- or 8-connectivity, see Fig.36.2. This function uses the mathematical morphology in order to get the contour.

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

1 def bwperim(I, connectivity=8):
    """

3 Morphological inner contour, in 4 or 8 connectivity
I: binary image

5 return : binary image representing the contour
        """

7 if connectivity==8:
    SE = disk (1);

9 else :
    SE = rectangle (3,3);

11 E = binary_erosion(I, selem=SE);

13 return I^E;

15 \# compute the contours
17 contours8 = bwperim(A, 4);

17 contours4 = bwperim(A, 8);

### 36.5.2 Freeman chain code
<figure><img src="image_846.png" /><figcaption>(c) Sample object.</figcaption></figure>

<figure><img src="image_847.png" /><figcaption>(d) Contour in 4-connectivity.</figcaption></figure>

<figure><img src="image_848.png" /><figcaption>(e) Contour in 8-connectivity.</figcaption></figure>

<caption>Figure 36.2: Simple object and its contours in 4- or 8-connectivity.</caption>

## First point of the shape

The important thing is to find one point in the contour. The Freeman chain code is sensitive to this choice, but several methods can transform this code so that this first choice does not have any importance.

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

def firstPoint (C):

2

### find first point of contour

4

### returns point (as array)

### find first point of contour

6

p = np.argwhere(C);

return p[0];

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

1

>>[r0,c0]=firstPoint (A)

3

r0 =
5

5

c0 =
10

## Freeman chain code

The principle is to follow the contour, delete each pixel at each step, and find the direction of the next pixel.
<figure><img src="image_851.png" /></figure>

def freeman(C, p, connectivity=8):

2

def getIndex(contour, point, connectivity):

""" subfunction for getting the local direction

4

"""

if connectivity ==8:

6

lut = np.array ([ [ 1, 2, 3], [ 8, 0, 4], [ 7, 6, 5] ]) ;

else :

8

lut = np.array ([ [ 0, 2, 0], [ 8, 0, 4], [ 0, 6, 0] ]) ;

10

window =contour [ point[0]-1 :point [0]+2, point [1]-1 :point [1]+2] ;

window =window *lut;

12

index =np.max (window);

return index-1;

14

# Be careful that these LUTs consider coordinates from left to right, top to
  ↪ → bottom

16

lutx =np.array ([ -1, -1, -1, 0, 1, 1, 1, 0]) ;

18

luty =np.array ([ -1, 0, 1, 1, 1, 0, -1, -1]) ;

lutcode =np.array ([ [3, 2, 1, 0, 7, 6, 5, 4]] ) ;

20

nbrpoints =np.sum (C);

22

code =[];

point =p.copy() ;

24

C2 =C.copy() ;

for i in np.arange (nbrpoints):

26

C2[point [0], point [1]] =0;

28

index =getIndex (C2, point, connectivity);

30

if (index==0):

32

C2[p[0], p[1]] =1;

index =getIndex (C2, point, connectivity);

34

# new point

36

point [0] =point [0] +lutx[index];

point [1] =point [1] +luty[index];

38

# add code

40

code.append (lutcode[index]) ;

return code;

42

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

code =freeman(C8, p);
<figure><img src="image_853.png" /></figure>

1 code = array ([6, 6, 5, 4, 6, 6, 6, 6, 6, 6, 6, 6, 6, 0, 7, 6, 6, 6, 0, 0, 0, 0, 2,
↔ 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 4, 4, 4, 4])

36.5.3 Normalization

Differential code

This is the first step towards independence from the first point.

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

1 def codediff(fcc, connectivity=8):
    sr = npRoll(fcc, 1);

3 d = fcc - sr;
return d%connectivity;

Normalization

The differential code is then normalized, in order to get a rotation invariant code. All
the circular shifts are evaluated, and a criterion (the minimum value) is established
to be able to always find the same result, for every position of the first point.

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

def minmag(code):
    # high value for min computing
    codemin = np.max(code)* np.ones(code.shape);
    nb = len(code);
    for i in np.arange(len(code)):
        C = np. roll (code, i);
        for j in np.arange nb):
            if C[j] > codemin[j]:
                break;
            elif C[j] < codemin[j]:
                codemin = C;
                break;
        if j == nb:
            codemin = C;
    return codemin;

The differential code is evaluated in d8, the normalization gives shapenumber8:
<figure><img src="image_856.png" /></figure>

c = codediff (code, 8);

2 shapenumber8 = minmag(c);

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

c= array([2, 0, 7, 7, 2, 0, 0, 0, 0, 0, 0, 0, 2, 7, 7, 0, 0, 2, 0, 0, 0, 2, 0,
    → 0, 7, 1, 0, 0, 0, 0, 0, 0, 0, 1, 7, 0, 2, 0, 0, 0])

2
shapenumber8= array([0, 0, 0, 0, 0, 0, 0, 1, 7, 0, 2, 0, 0, 0, 2, 0, 7, 7, 2,
    → 0, 0, 0, 0, 0, 0, 0, 2, 7, 7, 0, 0, 2, 0, 0, 0, 2, 0, 0, 7, 1])

Validation

This validation shows the effect on a different starting point.

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

1 p = np.array([4, 9]);

Another test is to verify the result after a rotation. To prevent discretization
problems, we use 90 degrees and take the transpose of the matrix.

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

1 % check for rotation by 90 deg
contours8rot=contours8';

The same code should be found in both cases.

36.5.4 Geometrical characterization

Perimeter for 8-connectivity

We first need to extract the codes in the diagonal directions and apply a √2 factor,
then add the number of codes in vertical and horizontal directions.
<figure><img src="image_860.png" /></figure>

def Perimeter(fcode):

2

###

fcode: Freeman code

4

###

nb_diag = np.sum(np.array(fcode)%2);

6

perim = nb_diag* nplandsq(2) + len(fcode) - nb_diag;

return perim;

The perimeter is evaluated in the same way in skimage.

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

Perimeter: 43.65685424949238

2

skimagemeasure.perimeter: 43.65685424949238

Area for 8-connectivity

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

def Area(fcode):

2

###

###

4

area = 0;

B = 0;

6

lutB = np.array([0, 1, 1, 1, 0, -1, -1, -1]);

for i in np.matenp(len(fcode)):

8

lutArea = np.array([-B, -(B+0.5), 0, (B+0.5), B, (B-0.5), 0, -(B-0.5)]);

area = area + lutArea[fcode[i]];

10

B = B + lutB[fcode[i]];

return area;

Notice that the area evaluated by this way is different from the number of pixels.

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

1

Area: 93.0

Number of pixels (area): 115
<figure><img src="image_864.png" /></figure>

## 37 Machine Learning

The objective of this tutorial is to classify images by using machine learning techniques. Some images, as illustrated in Figure 37.1 of the Kimia database will be used [1, 39].

Figure 37.1: The different processes will be applied on images from the Kimia database. Here are some examples.

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

### 37.1

#### Feature extraction

The image database used in this tutorial is composed of 18 classes. Each class contains 12 images. All these 216 images come from the Kimia database. In order to classify these images, we first have to extract some features of each binary image.

?

1. For each image in the database, extract a number of different features, denoted *nbFeatures*.

2. Organize these features in an array of *nbFeatures* lines and 216 columns. In this way, each column *i* represents the different features of the image *i*.

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

You can use the Python function skimageMeasure.regionprops to extract the same geometrical features.
## 37.2 Image classification

In order to classify the images, we are going to use neural networks. Pattern recognition networks are feedforward networks that can be trained to classify inputs according to target classes. The inputs are the features of each image. The target data are the labels, indicating the class of each image.
?

1. Build the target data, representing the class of each image. The target data for pattern recognition networks should consist of vectors of all zero values except for a 1 in element i, where i is the class they are to represent. In our example, the target data will be an array of 18 lines and 216 columns.

2. The database will be divided into a training set (75%) and a test set (25%).

3. Run the training task and classify the test images.

4. Show the classification confusion matrix as well as the overall performance.

5. Try to run the same process. What can you conclude?

6. Try to change the parameters of the network to improve the classification performance.

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

Look at the Python module `sklearn.neural_network.MLPClassifier` to make the classification. Notice that your data may be normalized before performing the classification.
<figure><img src="image_868.png" /></figure>

## 37.3. Python correction

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

The following imports may be considered.

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

```
# image manipulation and features construction
```

```
import glob
```

```
from skimage import measure, io
```

```
import numpy as np
```

```
import matplotlib .pyplot as plt
```

```
# preprocessing data and normalization
```

```
from sklern .preprocessing import StandardScaler
```

```
from sklern .preprocessing .data import QuantileTransformer
```

```
# learning methods
```

```
from sklern .import svm
```

```
from sklern.neural_network import MLPClassifier
```

```
from sklern .model_selection import train_test_split
```

```
from sklern .metrics import classification_report ,confusion_matrix
```

```
# plot confusion matrix
```

```
import seaborn as sn
```

```
import pandas as pd
```

### 37.3.1 Feature extraction

We make a loop on the whole database to extract some features of each image. The 9 features used here are: area, convex area, eccentricity, equivalent diameter, extent, major axis length, minor axis length, perimeter and solidity.
<figure><img src="image_871.png" /></figure>

# Definitions of the database, classes and images

2. `rep = 'images_Kimia216/';`
3. `classes = ['bird', 'bone', 'brick', 'camel', 'car', 'children', 'classic', 'elephant', 'face', 'fork', 'fountain', 'glass', 'hammer', 'heart', 'key', 'misk', 'ray', 'turtle'];`
4.

6. `nbClasses = len classes;`
7. `nbImages = 12;`

8. # The features are manually computed
    properties = np(zeros((nbClasses*nbImages,9)));
10. `target = np zeros nbClasses * nbImages;`
    `index=0;`

12. `for ind_c, c in enumerate classes):`
    `filelist = glob.glob+c+'*');`

14. `for filename in filelist :`
    `I = io .imread(filename);`
16. `prop = measure.regionprops(I);`
    `properties [index, 0] = prop [0].area;`
18. `properties [index, 1] = prop [0].convex_area;`
    `properties [index, 2] = prop [0].eccentricity;`
20. `properties [index, 3] = prop [0].equivalent_diameter;`
    `properties [index, 4] = prop [0].extent;`
22. `properties [index, 5] = prop [0].major_axis_length;`
    `properties [index, 6] = prop [0].minor_axis_length;`
24. `properties [index, 7] = prop [0].perimeter;`
    `properties [index, 8] = prop [0].solidity;`
26. `target[index] = ind_c;`

28. `index = index +1;`

Note that in the same time, the target array (required in the following) is built within this loop. It represents the true classes of the objects.

## 37.3.2 Classification

We used a training set of 75% of the database and 25% for the test set. The scaler is used in order to rescale the data having different ranges and dimensions. Other scalers are proposed in `sklearn.preprocessing`.

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

# percentage of the data used for splitting into train / test

2. `percentTest = .25;`

# MLP Classifier (Multi-Layer Perceptron)

4. # the data are first scaled
    `propertiesMLP = StandardScaler() . fit_transform (properties);`

6. `prop_train, prop_test, target_train, target_test = train_test_split (propertiesMLP,`
    `→ target, test_size =percentTest, random_state=0);`
The network is created with 1 hidden layers of 10 neurons. One could also use SVM instead of MLP ( classifier = svm.SVC Misclasser='linear') ).

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

# feedforward neural network

2 # max_iter should be extended max_iter=100000 for adam or sgd solvers

mlp = MLPClassifier hiddenLayerSizes =(10,) , solver='lbfgs');

4 target_pred = mlp.fit (prop_train, target_train) . predict (prop_test)

print("Training set score: %f" % mlp.score (prop_train, target_train))

The training score should be around 1, according to the training dataset and
the different parameters employed.

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

1 Training set score: 1.000000

37.3.3 Performance

The confusion matrix gives an idea of the overall performance. The following
function uses the modules seaborn and pandas to generate a heatmap that displays the
confusion matrix.

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

def plot_cm(cm, classes, normalize=False, cmap= plt.cm.Blues ):

2

Plot confusion matrix

4

cm: confusion matrix, as output by sklearn.metrics.confusion_matrix

classes : labels to be used

6

normalize: display number (False by default) or fraction (True)

cmap: colormap

8

returns : figure that can be used for pdf export

###

10

if normalize:

cm = cm .astype('float') / cm .sum( axis =1) [ : , np .newaxis]

12

fmt ='.1f' if normalize else 'd'

df_cm =pd .Frame(cm, index=classes, columns=classes);

14

fig = plt .figure() ;

16

sn. set( font_scale =.8)

sn . heatmap(df_cm, annot =True, cmap =cmap, fmt =fmt);

18

plt .xlabel('Target label')

plt .ylabel('True label')

20

return fig ;
The previous function is then used, and illustrated in Fig.37.2.

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

1 cnf_matrix = confusion_matrix( target_test , target_pred ); fig =plot_cm(cnf_matrix, classes , normalize=True);

sklearn.metrics.classification_report can be used to display the performance.

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

<table><thead><tr><th></th><th></th><th>precision</th><th>recall</th><th>f1-score</th><th>support</th></tr></thead><tbody><tr><td>2</td><td>bird</td><td>0.75</td><td>1.00</td><td>0.86</td><td>3</td></tr><tr><td>4</td><td>bone</td><td>1.00</td><td>1.00</td><td>1.00</td><td>5</td></tr><tr><td></td><td>brick</td><td>1.00</td><td>1.00</td><td>1.00</td><td>1</td></tr><tr><td>6</td><td>camel</td><td>1.00</td><td>1.00</td><td>1.00</td><td>3</td></tr><tr><td></td><td>car</td><td>1.00</td><td>1.00</td><td>1.00</td><td>1</td></tr><tr><td>8</td><td>children</td><td>1.00</td><td>1.00</td><td>1.00</td><td>2</td></tr><tr><td></td><td>classic</td><td>1.00</td><td>1.00</td><td>1.00</td><td>5</td></tr><tr><td>10</td><td>elephant</td><td>1.00</td><td>1.00</td><td>1.00</td><td>3</td></tr><tr><td></td><td>face</td><td>1.00</td><td>1.00</td><td>1.00</td><td>4</td></tr><tr><td>12</td><td>fork</td><td>1.00</td><td>1.00</td><td>1.00</td><td>2</td></tr><tr><td></td><td>fountain</td><td>1.00</td><td>1.00</td><td>1.00</td><td>3</td></tr><tr><td>14</td><td>glass</td><td>1.00</td><td>1.00</td><td>1.00</td><td>5</td></tr><tr><td></td><td>hammer</td><td>1.00</td><td>0.75</td><td>0.86</td><td>4</td></tr><tr><td>16</td><td>heart</td><td>1.00</td><td>1.00</td><td>1.00</td><td>2</td></tr><tr><td></td><td>key</td><td>1.00</td><td>1.00</td><td>1.00</td><td>3</td></tr><tr><td>18</td><td>misk</td><td>1.00</td><td>1.00</td><td>1.00</td><td>3</td></tr><tr><td></td><td>ray</td><td>1.00</td><td>1.00</td><td>1.00</td><td>4</td></tr><tr><td>20</td><td>turtle</td><td>1.00</td><td>1.00</td><td>1.00</td><td>1</td></tr><tr><td>22</td><td>avg / total</td><td>0.99</td><td>0.98</td><td>0.98</td><td>54</td></tr></tbody></table>
<figure><img src="image_878.png" /><figcaption>Figure 37.2: Normalized confusion matrix of the classification result.</figcaption></figure>
<figure><img src="image_879.png" /></figure>

38 Harris corner detector

The aim of this tutorial is to develop a simple Harris corner detector. This is the first step in pattern matching, generally followed by a feature descriptor construction, and a matching process.

38.1 Corner detector and cornerness measure

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

Use the sobel and gaussian_filter from the sciPY ndimage module.

38.1.1 Gradient evaluation

The Harris corner detector is based on the gradients of the image, Iₓ and Iᵧ in x and y directions, respectively.

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

Apply a Sobel gradient in both directions in order to compute Iₓ and Iᵧ.

38.1.2 Structure tensor

The structure tensor is defined by the following matrix. The coefficients ω follow a gaussian law, and each summation represents a gaussian filtering process. W is an operating window.

$$
M = \begin{bmatrix}
\sum_{(u,v) \in W} \omega(u,v) I_x(u,v)^2 & \sum_{(u,v) \in W} \omega(u,v) I_x(u,v) I_y(u,v) \\
\sum_{(u,v) \in W} \omega(u,v) I_x(u,v) I_y(u,v) & \sum_{(u,v) \in W} \omega(u,v) I_y(u,v)^2
\end{bmatrix}
= \begin{bmatrix}
M_1 & M_2 \\
M_3 & M_4
\end{bmatrix}
$$

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

• Evaluate M₁ to M₄ for each pixel of the image.
### 38.1.3 Cornerness measure

The cornerness measure C, as proposed by Harris and Stephens, is defined as follows for every pixel of coordinates (x, y):

$$C(x, y) = \det(M) - K \operatorname{trace}(M)^2$$

with K between 0.04 and 0.15.

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

Compute C for all pixels and display it for several scales σ.

## 38.2 Corners detection

A so-called Harris corner is the result of keeping only local maxima above a certain threshold value. You can use the checkerboard image for testing, or load the sweden road sign image Fig.38.1.

Use the following function to generate a checkerboard pattern.

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

def checkerboard(nb_x=2, nb_y=2, s=10):

2 checkerboard generation

4 a grid of size 2*nb_x by 2*nb_y is generated

6 each square has s pixels.

8 C = 255*np.kron([1, 0] * nb_x, [0, 1] * nb_x] * nb_y, np.ones((s, s)))

<figure><img src="image_885.png" /><figcaption>Figure 38.1: Sweden road sign to be used for corner detection.</figcaption></figure>
?

* Evaluate the extended maxima of the image.

* Only the strongest values of the cornerness measure should be kept. Two strategies can be employed in conjunction:

- Use a threshold value $t$ on $C$: the choice of this value is not trivial, and it strongly depends on the considered image. An adaptive method would be preferred.

- Keep only the $n$ strongest values.

* The previous operations are affected by the borders of the image. Thus, eliminate the corner points near the borders.

* The detected corners may contain several pixels. Keep only the centroid of each cluster.
<figure><img src="image_886.png" /></figure>

## 38.3. Python correction

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

### 38.3.1 Cornerness measure

Figure 38.2: Cornerness measure for the checkerboard image. The next step is to locate the corners by thresholding the measure, extracting the local maxima, eliminating points near the edges...

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

The first step is to compute the gradient in both x and y directions.

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

1. `Ix = scipy.ndimage.sobel(I, axis=0);`

2. `Iy = scipy.ndimage.sobel(I, axis=1);`

Then, the coefficients of the matrix are computed.

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

1. `M1 = np multiply(Ix, Ix);`

2. `M2 = np multiply(Iy, Ix);`

3. `M4 = np multiply(Iy, Iy);`

In case of using a scale parameter, these coefficients should be filtered (for example via a gaussian filter).

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

1. `M1 = scipy.ndimage.gaussian_filter(M1, sigma);`

2. `M2 = scipy.ndimage.gaussian_filter(M2, sigma);`

3. `M4 = scipy.ndimage.gaussian_filter(M4, sigma);`
Finally, the cornerness measure is evaluated.

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

1. $C = (\text{np multiply}(M1, M4) - \text{np multiply}(M2, M2)) - K * \text{np multiply}(M1+M4, M1+M4);$

The cornerness measure is displayed in Fig.38.2.

### 38.3.2 Corners detection

In order to keep only the strongest corner points, a threshold value $t$ is applied on $C$. This value is really depending on the considered image, thus such a global threshold is not generally a good idea. One would probably prefer a h-maximum or equivalent operator. For the purpose of this tutorial, we will keep this strategy.

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

1. $C[C<t] = 0;$

The local maxima are then extracted.

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

1. $\text{corners} = \text{peak_local_max}(C, \text{indices}=\text{False}, \text{min_distance}=2);$

The result is a binary image, where some clusters of points are the corner points. To keep only one point per cluster, the centroid of each is detected. The final result is presented in Fig.38.2c.

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

1. $L = \text{measure.label}(\text{corners});$
   $\text{props} = \text{measure.regionprops}(L);$

2. $\text{centers} = [];$

3. for prop in props:

$\text{centers}.append(prop.centroid);$

### 38.3.3 Road sign image application

In this case, the values $t = 10^7$ and $\sigma = 3$ are used. The result is illustrated in Fig.38.3.
<figure><img src="image_896.png" /><figcaption>Figure 38.3: Harris corner detection with scale $σ = 3$ and threshold value $t = 10^7$.</figcaption></figure>
<figure><img src="image_897.png" /></figure>

39 Local Binary Patterns

This tutorial aims to study a texture descriptor named 'Local Binary Patterns'. The first objective is to implement this descriptor. Therefater, digital images of textures will be classified using this descriptor and the k-means algorithm.

The different processes will be applied on this kind of texture images:

<figure><img src="image_898.png" /><figcaption>(a) metal image</figcaption></figure>

<figure><img src="image_899.png" /><figcaption>(b) sand image</figcaption></figure>

<figure><img src="image_900.png" /><figcaption>(c) ground image</figcaption></figure>

39.1

Local Binary Patterns

The Local Binary Patterns (LBP) descriptor is a simple yet very efficient texture operator which labels the pixels of an image by thresholding the neighborhood of each pixel and considers the result as a binary number. Due to its discriminative power and computational simplicity, LBP texture operator has become a popular approach in various applications. It can be seen as a unifying approach to the traditionally divergent statistical and structural models of texture analysis. Perhaps the most important property of the LBP operator in real-world applications is its robustness to monotonic gray-scale changes caused, for example, by illumination variations. Another important property is its computational simplicity, which makes it possible to analyze images in challenging real-time settings.

The LBP feature vector, in its simplest form, is created in the following manner:

• For each pixel, compare the pixel to each of its 8 neighbors (on its left-top, left-middle, left-bottom, right-top, etc.). Follow the pixels along a circle, i.e. clockwise or counter-clockwise.

* Where the center pixel's value is greater than the neighbor's value, write "1". Otherwise, write "0". This gives an 8-digit binary number (which is usually converted to decimal for convenience).
* Compute the histogram of the frequency of each "number" occurring (i.e., each combination of which pixels are smaller and which are greater than the center).

* Normalize the histogram.

<figure><img src="image_901.png" /><figcaption>Figure 39.1: Local binary pattern. From wikipedia, author Xiawi, CC-By-SA.</figcaption></figure>

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

* Code a function for computing the Local Binary Patterns.

* Test this operator on a texture image from the given database.

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

Consider the function numpy.histogram for histogram computation.

## 39.2

# Classification of texture images

The objective is to classify the texture images from the given database by using the LBP descriptor.

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

1. Calculate the LBP descriptor for each image of the database.

2. Compare the descriptors for each class of images.
3. Compute the distance between each pair of images in order to get a dissimilarity matrix. Comment the result.

4. Use the k-means algorithm to classify the images of the database into k classes ($k = 3$).

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

See KMeans from sklearn cluster.
<figure><img src="image_906.png" /></figure>

## 39.3. Python correction

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

The following imports are used.

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

import numpy as np

2 from scipy import misc

import matplotlib .pyplot as plt

4 import glob

import seaborn as sn

6 import pandas as pd

import os

8 from sklearn.cluster import KMeans

### 39.3.1 LBP computation

Each pixel is given a specific 8 bits value according to a code as follows.

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

def LBP(I):

2 B = np(zeros(np.shape(I)));

code = np.array ([[1,2,4],[8,0,16],[32,64,128]]);

4

# loop over all pixels except border pixels

6

for i in np.arange(1, I.shape[0]-2):

for j in np.arange(1, I.shape[1]-2):

8

w = I[i-1:i+2, j-1:j+2];

w = w >= I[i, j];

10

w = w * code;

B[i, j] = np.sum(w);

Then, all values (except for border values) are summarized in the histogram.

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

1 h,edges = np.histogram(B[1:-1, 1:-1], density = True, bins = 256);

For the first image of sand, the histogram is shown in Fig.39.2.
<figure><img src="image_911.png" /><figcaption>Figure 39.2: Illustration of the Local Binary Pattern computed on an entire image.</figcaption></figure>

### 39.3.2 Classification

For all images of the same family, the LBP are computed and represented in the same graph. The histograms really look similar (see Fig.39.3). The following code is used for the “sand” family.

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

```python
1 classes = ['Terrain', 'Metal', 'Sand'];
names = [];
3 hh = [];
for c in classes:
    print(c);
    fig = plt.figure();
    for file in sorted(glob.glob(' ../ matlab/images/' + c + ' *.bmp')):
        names.append(os.path basename(file));
        I = imageio impossread(file);
        I = I [:-1];
    h, edges = LBP(I);
    plt plot(h);
    hh.append(h);
```
Figure 39.3: Illustration of the LBP of 4 images of each family. The histogram are almost equivalent, which shows that this descriptor can be employed to discriminate between the different families.

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

A distance criterion is used to compare the different histograms: the classical SAD (Sum of Absolute Differences) gives a numerical values. All pairs of distances are concatenated in a matrix, displayed as an image in Fig.39.4.

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

1. `# compute distance between LBP`
   `n = len(hh);`

3. `dists = np.zeros((n, n))`

for i in np.arange(n):
    for j in np.arange(n):
        dists[i, j] = np.sum(np.abs(hh[i] - hh[j]))

7. `fig = plot_dists(dists, names);`

In order to display this matrix, the module `seaborn` is used.
<figure><img src="image_915.png" /></figure>

```python
1 def plot_dists (dists, classes, cmap= plt.cm.Blues ):
    """
3 Plot matrix of distances
dists : all computed distances
5 classes : labels to be used
cmap: colormap
7 returns : figure that can be used for pdf export
        """
9 df_cm = pd.DataFrame (dists, index = classes, columns = classes ) ;
    fig = plt .figure () ;
11 sn. set ( font_scale =.8)
    sn. heatmap (df_cm, annot = True, cmap = cmap, fmt ='.2f') ;
13 return fig ;
```

Figure 39.4: Sum of Absolute Differences between the different LBP histograms of each image. 3 families of 4 textures are represented here, terrain images are in the first part, metal images in the second and sand images in the last. Black represents 0 distance and white is 1 (highest distance, the values are normalized).

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

The kmeans algorithm uses such a distance, and we can verify that the clustering
process works as expected. The result is presented in the next box.
<figure><img src="image_917.png" /></figure>

1. `# kmeans clustering`
   `n=3;`

3. `k Means(init='k-means++', n clusters=n, n init =10)`
   `kmeans.fit (hh);`

5. `print (kmeans.labels_)`

The result show that the kmeans algorithm perfectly performs the classification.

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

1. `[1 1 1 1 0 0 0 0 2 2 2 2]`
# Part VI Exams
<figure><img src="image_919.png" /></figure>
All printed documents allowed. Send your code via the campus website (zip files if several files are to be sent).

## 40.1 Segmentation (13 points)

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

* Segment the image.

<figure><img src="image_921.png" /><figcaption>Figure 40.1: Test image</figcaption></figure>

## 40.2 Integral image (7 points)

If the original image is denoted $f$, the integral image $I$ is defined by the following equation, for all pixels of coordinates $(x, y)$:

$$I(x, y) = f(x, y) + I(x - 1, y) + I(x, y - 1) - I(x - 1, y - 1)$$

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

* Code a function that computes this integral image. The prototype of this function must be `function II=int_image(I)`.
* Code a function that computes the local average of the image f. Notice that:

$$ \sum_{x=x_1}^{x_2} \sum_{y=y_1}^{y_2} f(x, y) = I(x_2, y_2) - I(x_2, y_1 - 1) \\ - I(x_1 - 1, y_2) + I(x_1 - 1, y_1 - 1) $$

* Compare it to a mean filter.
41.1

Basic notions

41.1.1 Image characterisation

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

* Cite the names of the major image file formats and their main differ-
  ences. What is DICOM?

• Define sampling in numerical images. Define the image resolution.

41.1.2 Sensors

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

• What is a Bayer filter?

• Explain the principles of demosaicing.

41.1.3 Image processing

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

• Define the operation of histogram equalization.

• What is the difference with histogram stretching?

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

* From the mathematical definition of the derivative, explain the construction of the gradient operator, and then of the Laplacian.

• Cite some method for contours detection, and list their pros and cons.
41.2 Open question

The Fig. 41.1 shows a human retina (eye fundus). Propose a method for segmenting:

• the vessels,

• the optical nerve (bright disk).

Justify your choices. Indicate all elements that seem important and the way to deal with them. Look carefully at the shapes and intensities of objects you need to segment.

For your record, ophthalmologists need to measure the diameter and the circularity of the optical nerve, as well as the length and the number of vessels (branches, tortuosity...).
<figure><img src="image_927.png" /><figcaption>Figure 41.1: Human retina.</figcaption></figure>
<figure><img src="image_928.png" /></figure>
42.1

Warm-up questions

42.1.1 Thresholding methods

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

Cite different methods in order to find a threshold value to binarize a grayscale image. Explain the Otsu's method. What are its limits?

42.1.2 Image filtering

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

* What is the difference between a rank filter and a linear filter?

* Cite the name of a filter of each type.

* In case of a salt and pepper noise present in the image, what type of filter would you require to restore it? Why?

* For the latter, precise its limits and propose a way to improve the result.

42.1.3 Retina images

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

As a project, you coded a method for retina vessels segmentation. Explain the principles of the algorithm.
## 42.2 Open question: QR code

A QR code is a binary code represented as a 2D matrix. A scanner is in charge of reading it (via a camera) and translating it into the actual code. The structure of a QR code is represented in Fig.42.1.

<figure><img src="image_932.png" /><figcaption>Figure 42.1: Conversion of a binary pattern into a binary matrix.</figcaption></figure>

Describe the different steps required to perform the acquisition and geometric transformation of the QRcode into a binary 2D matrix representing the code (1 value for 1 square, see Fig.42.1).

* List different situations that may complicate the image acquisition and analysis.
* Propose some (image processing) solutions in order to deal with these situations. Explain them, give the necessary details describing the methods.

Figure 42.2: Different QR codes. The image acquisition and analysis must be tolerant to different observation conditions... Images from Collectif Raspouteam http://raspouteam.org/. Other images from unknown sources.

<figure><img src="image_933.png" /></figure>
# Index

Aliasing, 14, 22

Characterization, 345, 356, 365, 373
  Circularity, 367
  Convexity, 367
  Diameter, 366
  Perimeter, 365

Classification, 374, 403

Color
  Chromaticity diagram, 119
  Color Matching Functions, 119
  Color Spaces, 117

Computational Geometry
  Convex Hull, 315, 367
  Delaunay Triangulation, 323
  Minimum Spanning Tree, 323
  Voronoi Diagram, 323

Features, 373
  Detection, 384
  Harris Corner Detector, 395
  Local Binary Patterns, 401

Filtering
  Choquet, 128
  Convolution, 14
  Degenerate Diffusion, 142
  High-pass, 15, 24
  Linear Diffusion, 140
  Low-pass, 14, 23
  Non Linear Diffusion, 140
  Partial Differential Equations, 139

Fourier Transform, 49, 79
  2D FFT, 40
  Application, 41
  Filtering, 40
  Inverse, 40, 81

## Frameworks
Color Logarithmic Image Processing, 117
General Adaptive Neighborhood Image Processing, 127
Logarithmic Image Processing, 107
Freeman Chain Code, 355

Geometry
  Crofton Perimeter, 334
  Integral Geometry, 331

Gradient
  Prewitt, 15, 24
  Sobel, 15, 24

Heat Equation, 140
Histogram, 12, 18
  Definition, 29
  Equalization, 29
  Matching, 30

Image
  Display, 12, 17
  Load, 12, 17

Laplacian, 15, 24

Mapping, 20
Mathematical Morphology, 257, 301
  Alternate Sequential Filters, 272
  Attribute Filters, 279
  Dilation, 258
  Erosion, 258
  General Adaptive Neighborhood, 130
  Geodesic Filtering, 271
  Granulometry, 307
  Hit-or-miss Transform, 285
Reconstruction, 258, 272

Skeleton, 285

Topological Skeleton, 285

Watershed, 293

Multiscale, 149

Filtering, 151

Kramer & Bruckner, 151

Pyramid, 149

Noise

Exponential, 65

Gaussian, 65

Salt & Pepper, 65

Uniform, 65

Point Spread Function, 79

Quantization, 13, 21

Random Fields, 227

Registration, 381

Restoration

Adaptive median filter, 66

Blind deconvolution, 84

Deconvolution, 79

Lucy-Richardson filter, 84

Median filter, 66

Van Cittert filter, 83

Wiener filter, 82

Segmentation, 301

Active Contours, 195

Distance Map, 293

Histogram, 169

Hough Transform, 187

K-means, 170, 347, 403

Otsu thresholding, 170

Region Growing, 181

Threshold, 169

Watershed, 293

Watershed by Markers, 294

Shape From Focus, 95

SML, 97

Tenengrad, 98

Variance, 97

Variance of Tenengrad, 98

Spatial Processes, 205

Boolean Models, 219

Characterization, 209, 219

Gibbs Point Process, 207

Marked Process, 209

Miles Formula, 221

Neyman-Scott Point Process, 206

Poisson Point Process, 205

Ripley Functions, 209

Stereology, 237

Random Chords, 239

Random Planes, 241

Tomography, 159

Backprojection, 160

Filtered Backprojection, 161

Topology

Classification, 341

Neighborhood, 331, 339, 355

Wavelets

Definition, 49
<figure><img src="image_934.png" /></figure>
Bibliography

[1] http://vision.lems.brown.edu/content/available-software-and-databases. 367, 371, 376, 395

[2] R. J. Adler. *The Geometry of Random Fields*. Wiley, Chichester, UK, June 1981. 307

[3] O. S. Ahmad. *Stochastic representation and analysis of rough surface topography by random fields and integral geometry – Application to the UHMWPE cup involved in total hip arthroplasty*. PhD thesis, École Nationale Supérieure des Mines de Saint-Etienne, 2013. 307

[4] D. Bradley and G. Roth. Adaptive thresholding using the integral image. *Journal of graphics tools*, 12(2):13–21, 2007. 419

[5] F. Catté, P.-L. Lions, J.-M. Morel, and T. Coll. Image Selective Smoothing and Edge Detection by Nonlinear Diffusion. *SIAM Journal on Numerical Analysis*, 29(1):182–193, 1992. 147

[6] M. Couprie and G. Bertrand. Discrete topological transformations for image processing. In *Digital Geometry Algorithms*, pages 73–107. Springer, 2012. 361

[7] M. Fernandes, Y. Gavet, and J.-C. Pinoli. Improving focus measurements using logarithmic image processing. *Journal of Microscopy*, 242(3):228–241, 2010. 104, 106

[8] M. Fernandes, Y. Gavet, and J.-C. Pinoli. Robust 3-D reconstruction of surfaces from image focus by local cross-sectional multivariate statistical analyses: Application to human *ex vivo* corneal endotheliums. *Medical Image Analysis*, 16(6):1293–1306, 2012. 104, 106

[9] Y. Gavet, J. Debayle, and J.-C. Pinoli. The Color Logarithmic Image Processing (CoLIP) Antagonist Space. In *Color Image and Video Enhancement*, pages 155–182. Springer International Publishing, 2015. 123

[10] Y. Gavet, J. Debayle, and J.-C. Pinoli. The Color Logarithmic Image Processing (CoLIP) Antagonist Space and Chromaticity Diagram. In *International Workshop on Computational Color Imaging*, pages 131–138. Springer International Publishing, 2015. 123
[11] R. C. Gonzalez and R. E. Woods. *Digital Image Processing*. Pearson, 2007. 75, 83

[12] H. Gouinaud. *Traitement logarithmique d'images couleur*. PhD thesis, Ecole Nationale Supérieure des Mines de Saint-Etienne, 2013. 123

[13] H. Gouinaud, Y. Gavet, J. Debayle, and J.-C. Pinoli. Color correction in the framework of color logarithmic image processing. In *Image and Signal Processing and Analysis (ISPA), 2011 7th International Symposium on*, pages 129–133. IEEE, 2011. 123

[14] E. Grisan, A. Paviotti, N. Laurenti, and A. Ruggeri. A lattice estimation approach for the automatic evaluation of corneal endothelium density. In *Engineering in Medicine and Biology Society, 2005. IEEE-EMBS 2005. 27th Annual International Conference of the*, pages 1700–1703, 2005. 54

[15] F. C. Groen, I. T. Young, and G. Ligthart. A comparison of different focus functions for use in autofocus algorithms. *Cytometry*, 6(2):81–91, 1985. 105

[16] M. Jourlin and J. C. Pinoli. Logarithmic image processing. *Acta Stereologica*, 6:651–656, 1987. 113

[17] M. Jourlin and J.-C. Pinoli. Image Dynamic Range Enhancement and Stabilization in the Context of the Logarithmic Image Processing Model. *Signal Process.*, 41(2):225–237, Jan. 1995. 115

[18] M. Kass, A. Witkin, and D. Terzopoulos. Snakes: Active contour models. *International journal of computer vision*, 1(4):321–331, 1988. 243

[19] J. Koenderink. The structure of images. *Biological Cybernetics*, 50(5):363–370, 1984. 146

[20] H. P. Kramer and J. B. Bruckner. Iterations of a non-linear transformation for enhancement of digital images. *Pattern Recognition*, 7(1):53 – 58, 1975. 158

[21] E. Krotkov. Focusing. *International Journal of Computer Vision*, 1(3):223–237, 1988. 105

[22] A. Lang and J. Potthoff. Fast simulation of Gaussian random fields. *Monte Carlo Methods and Applications*, 17(3):195–214, 2011. 307

[23] P. Maragos and M. Akmal Butt. Partial differential equations in image analysis: continuous modeling, discrete processing. In *Image Processing, 1996. Proceedings., International Conference on*, volume 3, pages 61–64 vol.3, Sep 1996. 148

[24] G. Marsaglia. Choosing a Point from the Surface of a Sphere. *Ann. Math. Statist.*, 43(2):645–646, 04 1972. 321
[25] S. K. Nayar and Y. Nakagawa. Shape from focus. *Pattern analysis and machine intelligence, IEEE Transactions on*, 16(8):824-831, 1994. 105

[26] N. Otsu. A Threshold Selection Method from Gray-Level Histograms. *Systems, Man and Cybernetics, IEEE Transactions on*, 9(1):62-66, Jan 1979. 218

[27] P. Perona and J. Malik. Scale-space and edge detection using anisotropic diffusion. *Pattern Analysis and Machine Intelligence, IEEE Transactions on*, 12(7):629-639, Jul 1990. 147

[28] J. C. Pinoli. *Contribution à la modélisation, au traitement et à l'analyse d'image*. PhD thesis, Département de Mathématiques, Université de Saint-Etienne, France, 1987. 113

[29] J. C. Pinoli. The logarithmic image processing model: Connections with human brightness perception and contrast estimators. *Journal of Mathematical Imaging and Vision*, 7(4):341-358, Oct. 1997. 113

[30] S. Rivollier, J. Debayle, and J.-C. Pinoli. Shape diagrams for 2d compact sets-part i: analytic convex sets. australian journal of. *The Australian Journal of Mathematical Analysis and applications*, 7(2), 2010. 375

[31] S. Rivollier, J. Debayle, and J.-C. Pinoli. Shape diagrams for 2d compact sets-part ii: analytic simply connected sets. *The Australian Journal of Mathematical Analysis and applications*, 7(2), 2010. 375

[32] S. Rivollier, J. Debayle, and J.-C. Pinoli. Shape diagrams for 2d compact sets-part iii: convexity discrimination for analytic and discretized simply connected sets. *The Australian Journal of Mathematical Analysis and applications*, 7(2), 2010. 375

[33] A. Ruggeri, E. Grisan, and J. Jaroszewski. A new system for the automatic estimation of endothelial cell density in donor corneas. *Br J Ophthalmol*, 89(3):306-311, 2005. 54

[34] A. Ruggeri, E. Grisan, and J. Schroeter. Evaluation of repeatability for the automatic estimation of endothelial cell density in donor corneas. *Br J Ophthalmol*, 0, 2007. 54

[35] J. C. Russ and R. T. Dehoff. *Practical Stereology*. Kluwer Academic, 2nd edition, 2000. 317

[36] D. Sage, L. Donati, F. Soulez, D. Fortun, G. Schmit, A. Seitz, R. Guiet, C. Vonesch, and M. Unser. Deconvolutionlab2: An open-source software for deconvolution microscopy. *Methods*, 115:28-41, 2017. Image Processing for Biologists. 87

[37] B. Selig, K. A. Vermeer, B. Rieger, T. Hillenaar, and C. L. Luengo Hendriks. Fully automatic evaluation of the corneal endothelium from in vivo confocal microscopy. *Biomed central*, 2015. 49, 54
[38] J. Serra. *Image Analysis and Mathematical Morphology*. London: academic press, 1982. 175

[39] D. Sharvit, J. Chan, H. Tek, and B. B. Kimia. Symmetry-based indexing of image databases. In *Content-Based Access of Image and Video Libraries, 1998. Proceedings. IEEE Workshop on*, pages 56–62. IEEE, 1998. 367, 376, 395

[40] P. Soille. *Morphological Image Analysis: Principles and Applications*. Springer-Verlag New York, Inc., Secaucus, NJ, USA, 2nd edition, 2003. 175, 209

[41] S. A. Sugimoto and Y. Ichioka. Digital composition of images with increased depth of focus considering depth information. *Applied optics*, 24(14):2076–2080, 1985. 105

[42] J. M. Tenenbaum. *Accommodation in computer vision*. PhD thesis, STANFORD UNIV CA DEPT OF COMPUTER SCIENCE, 1970. 105

[43] K. Worsley. The geometry of random fields. *Chance*, 9(1):27–40, 1997. 307
**Image processing tutorials with python** This book is a collection of tutorials and exercises given at MINES Saint-Étienne as part of the Master's Degree in Science and Executive Engineering ("Ingénieur Civil des Mines - ICM"). In recent years, project-based learning has been used to illustrate theoretical concepts with real and concrete applications.

Whether you are in the early years of your university studies, in preparatory classes for the French Grandes Ecoles or in an engineering school, or even as a teacher, this book is made for you. You will find a large number of tutorials, classified by field, to familiarize yourself with the theoretical concepts of image processing and analysis.

Go to http://iptutorials.science to download the complete codes in python.

**Yann GAVET** He graduated from Mines Saint-Étienne with a Master's Degree in Science and Executive Engineering ("Ingénieur Civil des Mines - ICM") in 2001, obtained a Master of Science in 2004 and defended his PhD thesis in 2008. He teaches signal-processing, image-processing and pattern-recognition as well as C programming at Master's level.

**Johan DEBAYLE** He received his master of science and PhD thesis in 2002 and 2005. He is the head of the master MISPA of MINES Saint-Étienne, and teaches signal-processing, image-processing and pattern-recognition at Master's level.

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

MINES
Saint-Étienne

Une école de l'IMT