Coverage for install/scipp/scipy/ndimage/__init__.py: 94%
77 statements
« prev ^ index » next coverage.py v7.6.1, created at 2024-11-17 01:51 +0000
1# SPDX-License-Identifier: BSD-3-Clause
2# Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
3# @author Simon Heybrock
4"""Sub-package for multidimensional image processing.
6This subpackage provides wrappers for a subset of functions from
7:py:mod:`scipy.ndimage`.
8"""
10from collections.abc import Callable, Mapping
11from functools import wraps
12from typing import Any, Protocol, TypeVar
14import scipy.ndimage
16from ...core import (
17 CoordError,
18 DataArray,
19 DimensionError,
20 Variable,
21 VariancesError,
22 empty_like,
23 islinspace,
24 ones,
25)
26from ...typing import VariableLike
28_T = TypeVar('_T', Variable, DataArray)
31def _ndfilter(
32 func: Callable[..., _T],
33) -> Callable[..., _T]:
34 @wraps(func)
35 def function(x: _T, **kwargs: Any) -> _T:
36 if 'output' in kwargs:
37 raise TypeError("The 'output' argument is not supported")
38 if x.variances is not None:
39 raise VariancesError(
40 "Filter cannot be applied to input array with variances."
41 )
42 if getattr(x, 'masks', None):
43 raise ValueError("Filter cannot be applied to input array with masks.")
44 return func(x, **kwargs)
46 return function
49def _delta_to_positional(
50 x: Any,
51 dim: str,
52 index: float | Variable | Mapping[str, float | Variable],
53 dtype: type,
54) -> Any:
55 if not isinstance(index, Variable):
56 return index
57 coord = x.coords[dim]
58 if not islinspace(coord, dim).value:
59 raise CoordError(
60 f"Data points not regularly spaced along {dim}. To ignore this, "
61 f"provide a plain value (convertible to {dtype.__name__}) instead of a "
62 "scalar variable. Note that this will correspond to plain positional "
63 "indices/offsets."
64 )
65 pos = (len(coord) - 1) * (index.to(unit=coord.unit) / (coord[-1] - coord[0])).value
66 return dtype(pos)
69def _require_matching_dims(
70 index: Mapping[str, float | Variable],
71 x: VariableLike,
72 name: str | None,
73) -> None:
74 if set(index) != set(x.dims):
75 raise KeyError(
76 f"Data has dims={x.dims} but input argument '{name}' provides "
77 f"values for {tuple(index)}"
78 )
81def _positional_index(
82 x: Any,
83 index: float | Variable | Mapping[str, float | Variable],
84 name: str | None = None,
85 dtype: type = int,
86) -> list[Any]:
87 if not isinstance(index, Mapping):
88 return [_delta_to_positional(x, dim, index, dtype=dtype) for dim in x.dims]
89 _require_matching_dims(index, x, name)
90 return [_delta_to_positional(x, dim, index[dim], dtype=dtype) for dim in x.dims]
93@_ndfilter
94def gaussian_filter(
95 x: _T,
96 /,
97 *,
98 sigma: float | Variable | Mapping[str, float | Variable],
99 order: int | Mapping[str, int] | None = 0,
100 **kwargs: Any,
101) -> _T:
102 """
103 Multidimensional Gaussian filter.
105 This is a wrapper around :py:func:`scipy.ndimage.gaussian_filter`. See there for
106 full argument description. There are two key differences:
108 - This wrapper uses explicit dimension labels in the ``sigma`` and ``order``
109 arguments. For example, instead of ``sigma=[4, 6]`` use
110 ``sigma={'time':4, 'space':6}``
111 (with appropriate dimension labels for the data).
112 - Coordinate values can be used (and should be preferred) for ``sigma``. For
113 example, instead of ``sigma=[4, 6]`` use
114 ``sigma={'time':sc.scalar(5.0, unit='ms'), 'space':sc.scalar(1.2, unit='mm')}``.
115 In this case it is required that the corresponding coordinates exist and form a
116 "linspace", i.e., are evenly spaced.
118 Warning
119 -------
120 If ``sigma`` is an integer or a mapping to integers then coordinate values are
121 ignored. That is, the filter is applied even if the data points are not evenly
122 spaced. The resulting filtered data may thus have no meaningful interpretation.
124 Parameters
125 ----------
126 x: scipp.typing.VariableLike
127 Input variable or data array.
128 sigma:
129 Standard deviation for Gaussian kernel. The standard deviations of the Gaussian
130 filter are given as a mapping from dimension labels to numbers or scalar
131 variables, or as a single number or scalar variable, in which case it is equal
132 for all axes.
133 order:
134 The order of the filter along each dimension, given as mapping from dimension
135 labels to integers, or as a single integer. An order of 0 corresponds to
136 convolution with a Gaussian kernel. A positive order corresponds to convolution
137 with that derivative of a Gaussian.
139 Returns
140 -------
141 : scipp.typing.VariableLike
142 Filtered variable or data array
144 Examples
145 --------
147 .. plot:: :context: close-figs
149 >>> from scipp.scipy.ndimage import gaussian_filter
150 >>> da = sc.data.data_xy()
151 >>> da.plot()
153 With sigma as integer:
155 .. plot:: :context: close-figs
157 >>> filtered = gaussian_filter(da, sigma=4)
158 >>> filtered.plot()
160 With sigma based on input coordinate values:
162 .. plot:: :context: close-figs
164 >>> filtered = gaussian_filter(da, sigma=sc.scalar(0.1, unit='mm'))
165 >>> filtered.plot()
167 With different sigma for different dimensions:
169 .. plot:: :context: close-figs
171 >>> filtered = gaussian_filter(da, sigma={'x':sc.scalar(0.01, unit='mm'),
172 ... 'y':sc.scalar(1.0, unit='mm')})
173 >>> filtered.plot()
174 """
175 sigma_values = _positional_index(x, sigma, name='sigma', dtype=float)
176 if isinstance(order, Mapping):
177 _require_matching_dims(order, x, 'order')
178 order = [order[dim] for dim in x.dims] # type: ignore[assignment]
179 out = empty_like(x)
180 scipy.ndimage.gaussian_filter(
181 x.values, sigma=sigma_values, order=order, output=out.values, **kwargs
182 )
183 return out
186def _make_footprint(
187 x: Variable | DataArray,
188 size: int | Variable | Mapping[str, int | Variable] | None,
189 footprint: Variable | None,
190) -> Variable:
191 if footprint is None:
192 if size is None:
193 raise ValueError("Provide either 'footprint' or 'size'.")
194 footprint = ones(
195 dims=x.dims, shape=_positional_index(x, size, name='size'), dtype='bool'
196 )
197 else:
198 if size is not None:
199 raise ValueError("Provide either 'size' or 'footprint', not both.")
200 if set(footprint.dims) != set(x.dims):
201 raise DimensionError(
202 f"Dimensions {footprint.dims} must match data dimensions {x.dim}"
203 )
204 return footprint
207class _FootprintFilter(Protocol):
208 def __call__(
209 self,
210 x: _T,
211 /,
212 *,
213 size: int | Variable | Mapping[str, int | Variable] | None = None,
214 footprint: Variable | None = None,
215 origin: int | Variable | Mapping[str, int | Variable] | None = 0,
216 **kwargs: Any,
217 ) -> _T: ...
220def _make_footprint_filter(
221 name: str, example: bool = True, extra_args: str = ''
222) -> _FootprintFilter:
223 def footprint_filter(
224 x: _T,
225 /,
226 *,
227 size: int | Variable | Mapping[str, int | Variable] | None = None,
228 footprint: Variable | None = None,
229 origin: int | Variable | Mapping[str, int | Variable] = 0,
230 **kwargs: Any,
231 ) -> _T:
232 footprint = _make_footprint(x, size=size, footprint=footprint)
233 out = empty_like(x)
234 scipy_filter = getattr(scipy.ndimage, name)
235 scipy_filter(
236 x.values,
237 footprint=footprint.values,
238 origin=_positional_index(x, origin, name='origin'),
239 output=out.values,
240 **kwargs,
241 )
242 return out
244 footprint_filter.__name__ = name
245 if extra_args:
246 extra_args = f', {extra_args}'
247 doc = f"""
248 Calculate a multidimensional {name.replace('_', ' ')}.
250 This is a wrapper around :py:func:`scipy.ndimage.{name}`. See there for full
251 argument description. There are two key differences:
253 - This wrapper uses explicit dimension labels in the ``size``, ``footprint``, and
254 ``origin`` arguments. For example, instead of ``size=[4, 6]`` use
255 ``size={{'time':4, 'space':6}}`` (with appropriate dimension labels for the data).
256 - Coordinate values can be used (and should be preferred) for ``size`` and
257 ``origin``. For example, instead of ``size=[4, 6]`` use
258 ``size={{'time':sc.scalar(5.0, unit='ms'), 'space':sc.scalar(1.2, unit='mm')}}``.
259 In this case it is required that the corresponding coordinates exist and form a
260 "linspace", i.e., are evenly spaced.
262 Warning
263 -------
264 When ``size`` is an integer or a mapping to integers or when ``footprint`` is
265 given, coordinate values are ignored. That is, the filter is applied even if the
266 data points are not evenly spaced. The resulting filtered data may thus have no
267 meaningful interpretation.
269 Parameters
270 ----------
271 x: scipp.typing.VariableLike
272 Input variable or data array.
273 size:
274 Integer or scalar variable or mapping from dimension labels to integers or
275 scalar variables. Defines the footprint (see below).
276 footprint:
277 Variable with same dimension labels (but different shape) as the input data.
278 The boolean values specify (implicitly) a shape, but also which of the elements
279 within this shape will get passed to the filter function.
280 origin:
281 Integer or scalar variable or mapping from dimension labels to integers or
282 scalar variables. Controls the placement of the filter on the input array.
284 Returns
285 -------
286 : scipp.typing.VariableLike
287 Filtered variable or data array
288 """
289 if example:
290 doc += f"""
291 Examples
292 --------
294 .. plot:: :context: close-figs
296 >>> from scipp.scipy.ndimage import {name}
297 >>> da = sc.data.data_xy()
298 >>> da.plot()
300 With size as integer:
302 .. plot:: :context: close-figs
304 >>> filtered = {name}(da, size=4{extra_args})
305 >>> filtered.plot()
307 With size based on input coordinate values:
309 .. plot:: :context: close-figs
311 >>> filtered = {name}(da, size=sc.scalar(0.2, unit='mm'){extra_args})
312 >>> filtered.plot()
314 With different size for different dimensions:
316 .. plot:: :context: close-figs
318 >>> filtered = {name}(da, size={{'x':sc.scalar(0.2, unit='mm'),
319 ... {' ' * len(name)}'y':sc.scalar(1.1, unit='mm')}}{extra_args})
320 >>> filtered.plot()
321 """ # noqa: E501
322 footprint_filter.__doc__ = doc
323 return _ndfilter(footprint_filter)
326generic_filter = _make_footprint_filter('generic_filter', example=False)
327maximum_filter = _make_footprint_filter('maximum_filter')
328median_filter = _make_footprint_filter('median_filter')
329minimum_filter = _make_footprint_filter('minimum_filter')
330percentile_filter = _make_footprint_filter(
331 'percentile_filter', extra_args='percentile=80'
332)
333rank_filter = _make_footprint_filter('rank_filter', extra_args='rank=3')
335__all__ = [
336 'gaussian_filter',
337 'generic_filter',
338 'maximum_filter',
339 'median_filter',
340 'minimum_filter',
341 'percentile_filter',
342 'rank_filter',
343]