scikit-image

1
import inspect
2

3
import numpy as np
4
import scipy.ndimage as ndi
5

6
from skimage.data import camera
7
from skimage import restoration, data, color
8
from skimage.morphology import binary_dilation
9

10
try:
11
    from skimage.morphology import disk
12
except ImportError:
13
    from skimage.morphology import circle as disk
14
from . import _channel_kwarg, _skip_slow
15

16
# inspect signature to automatically handle API changes across versions
17
if 'num_iter' in inspect.signature(restoration.richardson_lucy).parameters:
18
    rl_iter_kwarg = dict(num_iter=10)
19
else:
20
    rl_iter_kwarg = dict(iterations=10)
21

22

23
class RestorationSuite:
24
    """Benchmark for restoration routines in scikit image."""
25

26
    timeout = 120
27

28
    def setup(self):
29
        nz = 32
30
        self.volume_f64 = (
31
            np.stack(
32
                [
33
                    camera()[::2, ::2],
34
                ]
35
                * nz,
36
                axis=-1,
37
            ).astype(float)
38
            / 255
39
        )
40
        self.sigma = 0.05
41
        self.volume_f64 += self.sigma * np.random.randn(*self.volume_f64.shape)
42
        self.volume_f32 = self.volume_f64.astype(np.float32)
43

44
    def peakmem_setup(self):
45
        pass
46

47
    def time_denoise_nl_means_f64(self):
48
        restoration.denoise_nl_means(
49
            self.volume_f64,
50
            patch_size=3,
51
            patch_distance=2,
52
            sigma=self.sigma,
53
            h=0.7 * self.sigma,
54
            fast_mode=False,
55
            **_channel_kwarg(False),
56
        )
57

58
    def time_denoise_nl_means_f32(self):
59
        restoration.denoise_nl_means(
60
            self.volume_f32,
61
            patch_size=3,
62
            patch_distance=2,
63
            sigma=self.sigma,
64
            h=0.7 * self.sigma,
65
            fast_mode=False,
66
            **_channel_kwarg(False),
67
        )
68

69
    def time_denoise_nl_means_fast_f64(self):
70
        restoration.denoise_nl_means(
71
            self.volume_f64,
72
            patch_size=3,
73
            patch_distance=2,
74
            sigma=self.sigma,
75
            h=0.7 * self.sigma,
76
            fast_mode=True,
77
            **_channel_kwarg(False),
78
        )
79

80
    def time_denoise_nl_means_fast_f32(self):
81
        restoration.denoise_nl_means(
82
            self.volume_f32,
83
            patch_size=3,
84
            patch_distance=2,
85
            sigma=self.sigma,
86
            h=0.7 * self.sigma,
87
            fast_mode=True,
88
        )
89

90
    def peakmem_denoise_nl_means_f64(self):
91
        restoration.denoise_nl_means(
92
            self.volume_f64,
93
            patch_size=3,
94
            patch_distance=2,
95
            sigma=self.sigma,
96
            h=0.7 * self.sigma,
97
            fast_mode=False,
98
            **_channel_kwarg(False),
99
        )
100

101
    def peakmem_denoise_nl_means_f32(self):
102
        restoration.denoise_nl_means(
103
            self.volume_f32,
104
            patch_size=3,
105
            patch_distance=2,
106
            sigma=self.sigma,
107
            h=0.7 * self.sigma,
108
            fast_mode=False,
109
        )
110

111
    def peakmem_denoise_nl_means_fast_f64(self):
112
        restoration.denoise_nl_means(
113
            self.volume_f64,
114
            patch_size=3,
115
            patch_distance=2,
116
            sigma=self.sigma,
117
            h=0.7 * self.sigma,
118
            fast_mode=True,
119
            **_channel_kwarg(False),
120
        )
121

122
    def peakmem_denoise_nl_means_fast_f32(self):
123
        restoration.denoise_nl_means(
124
            self.volume_f32,
125
            patch_size=3,
126
            patch_distance=2,
127
            sigma=self.sigma,
128
            h=0.7 * self.sigma,
129
            fast_mode=True,
130
            **_channel_kwarg(False),
131
        )
132

133

134
class DeconvolutionSuite:
135
    """Benchmark for restoration routines in scikit image."""
136

137
    def setup(self):
138
        nz = 32
139
        self.volume_f64 = (
140
            np.stack(
141
                [
142
                    camera()[::2, ::2],
143
                ]
144
                * nz,
145
                axis=-1,
146
            ).astype(float)
147
            / 255
148
        )
149
        self.sigma = 0.02
150
        self.psf_f64 = np.ones((5, 5, 5)) / 125
151
        self.psf_f32 = self.psf_f64.astype(np.float32)
152
        self.volume_f64 = ndi.convolve(self.volume_f64, self.psf_f64)
153
        self.volume_f64 += self.sigma * np.random.randn(*self.volume_f64.shape)
154
        self.volume_f32 = self.volume_f64.astype(np.float32)
155

156
    def peakmem_setup(self):
157
        pass
158

159
    def time_richardson_lucy_f64(self):
160
        restoration.richardson_lucy(self.volume_f64, self.psf_f64, **rl_iter_kwarg)
161

162
    def time_richardson_lucy_f32(self):
163
        restoration.richardson_lucy(self.volume_f32, self.psf_f32, **rl_iter_kwarg)
164

165
    # use iterations=1 for peak-memory cases to save time
166
    def peakmem_richardson_lucy_f64(self):
167
        restoration.richardson_lucy(self.volume_f64, self.psf_f64, **rl_iter_kwarg)
168

169
    def peakmem_richardson_lucy_f32(self):
170
        restoration.richardson_lucy(self.volume_f32, self.psf_f32, **rl_iter_kwarg)
171

172

173
class RollingBall:
174
    """Benchmark Rolling Ball algorithm."""
175

176
    timeout = 120
177

178
    def time_rollingball(self, radius):
179
        restoration.rolling_ball(data.coins(), radius=radius)
180

181
    time_rollingball.params = [25, 50, 100, 200]
182
    time_rollingball.param_names = ["radius"]
183

184
    def peakmem_reference(self, *args):
185
        """Provide reference for memory measurement with empty benchmark.
186

187
        Peakmem benchmarks measure the maximum amount of RAM used by a
188
        function. However, this maximum also includes the memory used
189
        during the setup routine (as of asv 0.2.1; see [1]_).
190
        Measuring an empty peakmem function might allow us to disambiguate
191
        between the memory used by setup and the memory used by target (see
192
        other ``peakmem_`` functions below).
193

194
        References
195
        ----------
196
        .. [1]: https://asv.readthedocs.io/en/stable/writing_benchmarks.html#peak-memory
197
        """
198
        pass
199

200
    def peakmem_rollingball(self, radius):
201
        restoration.rolling_ball(data.coins(), radius=radius)
202

203
    peakmem_rollingball.params = [25, 50, 100, 200]
204
    peakmem_rollingball.param_names = ["radius"]
205

206
    def time_rollingball_nan(self, radius):
207
        image = data.coins().astype(float)
208
        pos = np.arange(np.min(image.shape))
209
        image[pos, pos] = np.nan
210
        restoration.rolling_ball(image, radius=radius, nansafe=True)
211

212
    time_rollingball_nan.params = [25, 50, 100, 200]
213
    time_rollingball_nan.param_names = ["radius"]
214

215
    def time_rollingball_ndim(self):
216
        from skimage.restoration._rolling_ball import ellipsoid_kernel
217

218
        image = data.cells3d()[:, 1, ...]
219
        kernel = ellipsoid_kernel((1, 100, 100), 100)
220
        restoration.rolling_ball(image, kernel=kernel)
221

222
    time_rollingball_ndim.setup = _skip_slow
223

224
    def time_rollingball_threads(self, threads):
225
        restoration.rolling_ball(data.coins(), radius=100, num_threads=threads)
226

227
    time_rollingball_threads.params = (0, 2, 4, 8)
228
    time_rollingball_threads.param_names = ["threads"]
229

230

231
class Inpaint:
232
    """Benchmark inpainting algorithm."""
233

234
    def setup(self):
235
        image = data.astronaut()
236

237
        # Create mask with six block defect regions
238
        mask = np.zeros(image.shape[:-1], dtype=bool)
239
        mask[20:60, :20] = 1
240
        mask[160:180, 70:155] = 1
241
        mask[30:60, 170:195] = 1
242
        mask[-60:-30, 170:195] = 1
243
        mask[-180:-160, 70:155] = 1
244
        mask[-60:-20, :20] = 1
245

246
        # add a few long, narrow defects
247
        mask[200:205, -200:] = 1
248
        mask[150:255, 20:23] = 1
249
        mask[365:368, 60:130] = 1
250

251
        # add randomly positioned small point-like defects
252
        rstate = np.random.RandomState(0)
253
        for radius in [0, 2, 4]:
254
            # larger defects are less common
255
            thresh = 2.75 + 0.25 * radius  # larger defects are less common
256
            tmp_mask = rstate.randn(*image.shape[:-1]) > thresh
257
            if radius > 0:
258
                tmp_mask = binary_dilation(tmp_mask, disk(radius, dtype=bool))
259
            mask[tmp_mask] = 1
260

261
        for layer in range(image.shape[-1]):
262
            image[np.where(mask)] = 0
263

264
        self.image_defect = image
265
        self.image_defect_gray = color.rgb2gray(image)
266
        self.mask = mask
267

268
    def time_inpaint_rgb(self):
269
        restoration.inpaint_biharmonic(
270
            self.image_defect, self.mask, **_channel_kwarg(True)
271
        )
272

273
    def time_inpaint_grey(self):
274
        restoration.inpaint_biharmonic(
275
            self.image_defect_gray, self.mask, **_channel_kwarg(False)
276
        )
277