1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
"""
============================================
Measure perimeters with different estimators
============================================
In this example, we show the error on measuring perimeters, comparing classic
approximations and Crofton ones. For that, we estimate the perimeter of an
object (either a square or a disk) and its rotated version, as we increase the
rotation angle.
"""
import matplotlib.pyplot as plt
import numpy as np
from skimage.measure import perimeter, perimeter_crofton
from skimage.transform import rotate
# scale parameter can be used to increase the grid size. The resulting curves
# should be smoothed with higher scales
scale = 10
# Construct two objects, a square and a disk
square = np.zeros((100 * scale, 100 * scale))
square[40 * scale : 60 * scale, 40 * scale : 60 * scale] = 1
[X, Y] = np.meshgrid(np.linspace(0, 100 * scale), np.linspace(0, 100 * scale))
R = 20 * scale
disk = (X - 50 * scale) ** 2 + (Y - 50 * scale) ** 2 <= R**2
fig, axes = plt.subplots(1, 2, figsize=(8, 5))
ax = axes.flatten()
dX = X[0, 1] - X[0, 0]
true_perimeters = [80 * scale, 2 * np.pi * R / dX]
# For each type of object, estimate its perimeter as the object is rotated,
# according to different approximations
for index, obj in enumerate([square, disk]):
# `neighborhood` value can be 4 or 8 for the classic perimeter estimator
for n in [4, 8]:
p = []
angles = range(90)
for i in angles:
rotated = rotate(obj, i, order=0)
p.append(perimeter(rotated, n))
ax[index].plot(angles, p)
# `directions` value can be 2 or 4 for the Crofton estimator
for d in [2, 4]:
p = []
angles = np.arange(0, 90, 2)
for i in angles:
rotated = rotate(obj, i, order=0)
p.append(perimeter_crofton(rotated, d))
ax[index].plot(angles, p)
ax[index].axhline(true_perimeters[index], linestyle='--', color='k')
ax[index].set_xlabel('Rotation angle')
ax[index].legend(
[
'N4 perimeter',
'N8 perimeter',
'Crofton 2 directions',
'Crofton 4 directions',
'Ground truth',
],
loc='best',
)
ax[index].set_ylabel('Perimeter of the rotated object')
ax[0].set_title('Square')
ax[1].set_title('Disk')
plt.tight_layout()
plt.show()
|
