1import math
2
3import numpy as np
4from scipy import ndimage as ndi
5
6from .._shared import utils
7from .._shared.filters import gaussian
8from .._shared.utils import convert_to_float
9from ..transform import resize
10
11
12def _smooth(image, sigma, mode, cval, multichannel=None):
13    """Return image with each channel smoothed by the Gaussian filter."""
14    smoothed = np.empty_like(image)
15
16    # apply Gaussian filter to all channels independently
17    if multichannel:
18        sigma = (sigma, ) * (image.ndim - 1) + (0, )
19        channel_axis = -1
20    else:
21        channel_axis = None
22    gaussian(image, sigma, output=smoothed, mode=mode, cval=cval,
23             channel_axis=channel_axis)
24    return smoothed
25
26
27def _check_factor(factor):
28    if factor <= 1:
29        raise ValueError('scale factor must be greater than 1')
30
31
32@utils.channel_as_last_axis()
33@utils.deprecate_multichannel_kwarg(multichannel_position=6)
34def pyramid_reduce(image, downscale=2, sigma=None, order=1,
35                   mode='reflect', cval=0, multichannel=False,
36                   preserve_range=False, *, channel_axis=-1):
37    """Smooth and then downsample image.
38
39    Parameters
40    ----------
41    image : ndarray
42        Input image.
43    downscale : float, optional
44        Downscale factor.
45    sigma : float, optional
46        Sigma for Gaussian filter. Default is `2 * downscale / 6.0` which
47        corresponds to a filter mask twice the size of the scale factor that
48        covers more than 99% of the Gaussian distribution.
49    order : int, optional
50        Order of splines used in interpolation of downsampling. See
51        `skimage.transform.warp` for detail.
52    mode : {'reflect', 'constant', 'edge', 'symmetric', 'wrap'}, optional
53        The mode parameter determines how the array borders are handled, where
54        cval is the value when mode is equal to 'constant'.
55    cval : float, optional
56        Value to fill past edges of input if mode is 'constant'.
57    multichannel : bool, optional
58        Whether the last axis of the image is to be interpreted as multiple
59        channels or another spatial dimension. This argument is deprecated:
60        specify `channel_axis` instead.
61    preserve_range : bool, optional
62        Whether to keep the original range of values. Otherwise, the input
63        image is converted according to the conventions of `img_as_float`.
64        Also see https://scikit-image.org/docs/dev/user_guide/data_types.html
65    channel_axis : int or None, optional
66        If None, the image is assumed to be a grayscale (single channel) image.
67        Otherwise, this parameter indicates which axis of the array corresponds
68        to channels.
69
70        .. versionadded:: 0.19
71           ``channel_axis`` was added in 0.19.
72
73    Returns
74    -------
75    out : array
76        Smoothed and downsampled float image.
77
78    References
79    ----------
80    .. [1] http://persci.mit.edu/pub_pdfs/pyramid83.pdf
81
82    """
83    _check_factor(downscale)
84    multichannel = channel_axis is not None
85
86    image = convert_to_float(image, preserve_range)
87
88    out_shape = tuple([math.ceil(d / float(downscale)) for d in image.shape])
89    if multichannel:
90        out_shape = out_shape[:-1]
91
92    if sigma is None:
93        # automatically determine sigma which covers > 99% of distribution
94        sigma = 2 * downscale / 6.0
95
96    smoothed = _smooth(image, sigma, mode, cval, multichannel)
97    out = resize(smoothed, out_shape, order=order, mode=mode, cval=cval,
98                 anti_aliasing=False)
99
100    return out
101
102
103@utils.channel_as_last_axis()
104@utils.deprecate_multichannel_kwarg(multichannel_position=6)
105def pyramid_expand(image, upscale=2, sigma=None, order=1,
106                   mode='reflect', cval=0, multichannel=False,
107                   preserve_range=False, *, channel_axis=-1):
108    """Upsample and then smooth image.
109
110    Parameters
111    ----------
112    image : ndarray
113        Input image.
114    upscale : float, optional
115        Upscale factor.
116    sigma : float, optional
117        Sigma for Gaussian filter. Default is `2 * upscale / 6.0` which
118        corresponds to a filter mask twice the size of the scale factor that
119        covers more than 99% of the Gaussian distribution.
120    order : int, optional
121        Order of splines used in interpolation of upsampling. See
122        `skimage.transform.warp` for detail.
123    mode : {'reflect', 'constant', 'edge', 'symmetric', 'wrap'}, optional
124        The mode parameter determines how the array borders are handled, where
125        cval is the value when mode is equal to 'constant'.
126    cval : float, optional
127        Value to fill past edges of input if mode is 'constant'.
128    multichannel : bool, optional
129        Whether the last axis of the image is to be interpreted as multiple
130        channels or another spatial dimension. This argument is deprecated:
131        specify `channel_axis` instead.
132    preserve_range : bool, optional
133        Whether to keep the original range of values. Otherwise, the input
134        image is converted according to the conventions of `img_as_float`.
135        Also see https://scikit-image.org/docs/dev/user_guide/data_types.html
136    channel_axis : int or None, optional
137        If None, the image is assumed to be a grayscale (single channel) image.
138        Otherwise, this parameter indicates which axis of the array corresponds
139        to channels.
140
141        .. versionadded:: 0.19
142           ``channel_axis`` was added in 0.19.
143
144    Returns
145    -------
146    out : array
147        Upsampled and smoothed float image.
148
149    References
150    ----------
151    .. [1] http://persci.mit.edu/pub_pdfs/pyramid83.pdf
152
153    """
154    _check_factor(upscale)
155    multichannel = channel_axis is not None
156
157    image = convert_to_float(image, preserve_range)
158
159    out_shape = tuple([math.ceil(upscale * d) for d in image.shape])
160    if multichannel:
161        out_shape = out_shape[:-1]
162
163    if sigma is None:
164        # automatically determine sigma which covers > 99% of distribution
165        sigma = 2 * upscale / 6.0
166
167    resized = resize(image, out_shape, order=order,
168                     mode=mode, cval=cval, anti_aliasing=False)
169    out = _smooth(resized, sigma, mode, cval, multichannel)
170
171    return out
172
173
174@utils.channel_as_last_axis()
175@utils.deprecate_multichannel_kwarg(multichannel_position=7)
176def pyramid_gaussian(image, max_layer=-1, downscale=2, sigma=None, order=1,
177                     mode='reflect', cval=0, multichannel=False,
178                     preserve_range=False, *, channel_axis=-1):
179    """Yield images of the Gaussian pyramid formed by the input image.
180
181    Recursively applies the `pyramid_reduce` function to the image, and yields
182    the downscaled images.
183
184    Note that the first image of the pyramid will be the original, unscaled
185    image. The total number of images is `max_layer + 1`. In case all layers
186    are computed, the last image is either a one-pixel image or the image where
187    the reduction does not change its shape.
188
189    Parameters
190    ----------
191    image : ndarray
192        Input image.
193    max_layer : int, optional
194        Number of layers for the pyramid. 0th layer is the original image.
195        Default is -1 which builds all possible layers.
196    downscale : float, optional
197        Downscale factor.
198    sigma : float, optional
199        Sigma for Gaussian filter. Default is `2 * downscale / 6.0` which
200        corresponds to a filter mask twice the size of the scale factor that
201        covers more than 99% of the Gaussian distribution.
202    order : int, optional
203        Order of splines used in interpolation of downsampling. See
204        `skimage.transform.warp` for detail.
205    mode : {'reflect', 'constant', 'edge', 'symmetric', 'wrap'}, optional
206        The mode parameter determines how the array borders are handled, where
207        cval is the value when mode is equal to 'constant'.
208    cval : float, optional
209        Value to fill past edges of input if mode is 'constant'.
210    multichannel : bool, optional
211        Whether the last axis of the image is to be interpreted as multiple
212        channels or another spatial dimension. This argument is deprecated:
213        specify `channel_axis` instead.
214    preserve_range : bool, optional
215        Whether to keep the original range of values. Otherwise, the input
216        image is converted according to the conventions of `img_as_float`.
217        Also see https://scikit-image.org/docs/dev/user_guide/data_types.html
218    channel_axis : int or None, optional
219        If None, the image is assumed to be a grayscale (single channel) image.
220        Otherwise, this parameter indicates which axis of the array corresponds
221        to channels.
222
223        .. versionadded:: 0.19
224           ``channel_axis`` was added in 0.19.
225
226    Returns
227    -------
228    pyramid : generator
229        Generator yielding pyramid layers as float images.
230
231    References
232    ----------
233    .. [1] http://persci.mit.edu/pub_pdfs/pyramid83.pdf
234
235    """
236    _check_factor(downscale)
237
238    # cast to float for consistent data type in pyramid
239    image = convert_to_float(image, preserve_range)
240
241    layer = 0
242    current_shape = image.shape
243
244    prev_layer_image = image
245    yield image
246
247    # build downsampled images until max_layer is reached or downscale process
248    # does not change image size
249    while layer != max_layer:
250        layer += 1
251
252        layer_image = pyramid_reduce(prev_layer_image, downscale, sigma, order,
253                                     mode, cval, channel_axis=channel_axis)
254
255        prev_shape = np.asarray(current_shape)
256        prev_layer_image = layer_image
257        current_shape = np.asarray(layer_image.shape)
258
259        # no change to previous pyramid layer
260        if np.all(current_shape == prev_shape):
261            break
262
263        yield layer_image
264
265
266@utils.channel_as_last_axis()
267@utils.deprecate_multichannel_kwarg(multichannel_position=7)
268def pyramid_laplacian(image, max_layer=-1, downscale=2, sigma=None, order=1,
269                      mode='reflect', cval=0, multichannel=False,
270                      preserve_range=False, *, channel_axis=-1):
271    """Yield images of the laplacian pyramid formed by the input image.
272
273    Each layer contains the difference between the downsampled and the
274    downsampled, smoothed image::
275
276        layer = resize(prev_layer) - smooth(resize(prev_layer))
277
278    Note that the first image of the pyramid will be the difference between the
279    original, unscaled image and its smoothed version. The total number of
280    images is `max_layer + 1`. In case all layers are computed, the last image
281    is either a one-pixel image or the image where the reduction does not
282    change its shape.
283
284    Parameters
285    ----------
286    image : ndarray
287        Input image.
288    max_layer : int, optional
289        Number of layers for the pyramid. 0th layer is the original image.
290        Default is -1 which builds all possible layers.
291    downscale : float, optional
292        Downscale factor.
293    sigma : float, optional
294        Sigma for Gaussian filter. Default is `2 * downscale / 6.0` which
295        corresponds to a filter mask twice the size of the scale factor that
296        covers more than 99% of the Gaussian distribution.
297    order : int, optional
298        Order of splines used in interpolation of downsampling. See
299        `skimage.transform.warp` for detail.
300    mode : {'reflect', 'constant', 'edge', 'symmetric', 'wrap'}, optional
301        The mode parameter determines how the array borders are handled, where
302        cval is the value when mode is equal to 'constant'.
303    cval : float, optional
304        Value to fill past edges of input if mode is 'constant'.
305    multichannel : bool, optional
306        Whether the last axis of the image is to be interpreted as multiple
307        channels or another spatial dimension. This argument is deprecated:
308        specify `channel_axis` instead.
309    preserve_range : bool, optional
310        Whether to keep the original range of values. Otherwise, the input
311        image is converted according to the conventions of `img_as_float`.
312        Also see https://scikit-image.org/docs/dev/user_guide/data_types.html
313    channel_axis : int or None, optional
314        If None, the image is assumed to be a grayscale (single channel) image.
315        Otherwise, this parameter indicates which axis of the array corresponds
316        to channels.
317
318        .. versionadded:: 0.19
319           ``channel_axis`` was added in 0.19.
320
321    Returns
322    -------
323    pyramid : generator
324        Generator yielding pyramid layers as float images.
325
326    References
327    ----------
328    .. [1] http://persci.mit.edu/pub_pdfs/pyramid83.pdf
329    .. [2] http://sepwww.stanford.edu/data/media/public/sep/morgan/texturematch/paper_html/node3.html
330
331    """
332    _check_factor(downscale)
333    multichannel = channel_axis is not None
334
335    # cast to float for consistent data type in pyramid
336    image = convert_to_float(image, preserve_range)
337
338    if sigma is None:
339        # automatically determine sigma which covers > 99% of distribution
340        sigma = 2 * downscale / 6.0
341
342    current_shape = image.shape
343
344    smoothed_image = _smooth(image, sigma, mode, cval, multichannel)
345    yield image - smoothed_image
346
347    # build downsampled images until max_layer is reached or downscale process
348    # does not change image size
349    if max_layer == -1:
350        max_layer = int(np.ceil(math.log(np.max(current_shape), downscale)))
351
352    for layer in range(max_layer):
353
354        out_shape = tuple(
355            [math.ceil(d / float(downscale)) for d in current_shape])
356
357        if multichannel:
358            out_shape = out_shape[:-1]
359
360        resized_image = resize(smoothed_image, out_shape, order=order,
361                               mode=mode, cval=cval, anti_aliasing=False)
362        smoothed_image = _smooth(resized_image, sigma, mode, cval,
363                                 multichannel)
364        current_shape = np.asarray(resized_image.shape)
365
366        yield resized_image - smoothed_image
367