1import collections.abc as cabc
2import itertools
3
4import pyxu.abc as pxa
5import pyxu.info.deps as pxd
6import pyxu.info.ptype as pxt
7import pyxu.operator.linop.base as pxlb
8import pyxu.operator.linop.diff as pxld
9import pyxu.operator.linop.pad as pxlp
10import pyxu.operator.linop.stencil.stencil as pxls
11import pyxu.runtime as pxrt
12import pyxu.util as pxu
13
14try:
15 import scipy.ndimage._filters as scif
16except ImportError:
17 import scipy.ndimage.filters as scif
18
19import functools
20import typing as typ
21
22import numpy as np
23
24IndexSpec = cabc.Sequence[pxt.Integer]
25KernelSpec = pxls.Stencil.KernelSpec
26ModeSpec = pxlp.Pad.ModeSpec
27
28
29__all__ = [
30 "MovingAverage",
31 "GaussianFilter",
32 "DifferenceOfGaussians",
33 "DoG",
34 "Laplace",
35 "Sobel",
36 "Prewitt",
37 "Scharr",
38 "StructureTensor",
39]
40
41
42def _to_canonical_form(_, dim_shape):
43 if not isinstance(_, cabc.Sequence):
44 _ = (_,) * len(dim_shape)
45 else:
46 assert len(_) == len(dim_shape)
47 return _
48
49
50def _get_axes(axis, ndim):
51 if axis is None:
52 axes = list(range(ndim))
53 elif np.isscalar(axis):
54 axes = [axis]
55 else:
56 axes = axis
57 return axes
58
59
60def _sanitize_inputs(dim_shape, dtype, gpu):
61 ndim = len(dim_shape)
62
63 if dtype is None:
64 dtype = pxrt.Width.DOUBLE.value
65
66 if gpu:
67 assert pxd.CUPY_ENABLED
68 import cupy as xp
69 else:
70 import numpy as xp
71 return ndim, dtype, xp
72
73
[docs]
74def MovingAverage(
75 dim_shape: pxt.NDArrayShape,
76 size: typ.Union[typ.Tuple, pxt.Integer],
77 center: typ.Optional[IndexSpec] = None,
78 mode: ModeSpec = "constant",
79 gpu: bool = False,
80 dtype: typ.Optional[pxt.DType] = None,
81):
82 r"""
83 Multi-dimensional moving average or uniform filter.
84
85 Notes
86 -----
87 This operator performs a convolution between the input :math:`D`-dimensional NDArray :math:`\mathbf{x} \in
88 \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}` and a uniform :math:`D`-dimensional filter :math:`\mathbf{h} \in
89 \mathbb{R}^{\text{size} \times \cdots \times \text{size}}` that computes the `size`-point local mean values using
90 separable kernels for improved performance.
91
92 .. math::
93
94 y_{i} = \frac{1}{|\mathcal{N}_{i}|}\sum_{j \in \mathcal{N}_{i}} x_{j}
95
96 Where :math:`\mathcal{N}_{i}` is the set of elements neighbouring the :math:`i`-th element of the input array, and
97 :math:`\mathcal{N}_{i}` denotes its cardinality, i.e. the total number of neighbors.
98
99
100 Parameters
101 ----------
102 dim_shape: NDArrayShape
103 Shape of the input array.
104 size: int, tuple
105 Size of the moving average kernel.
106
107 If a single integer value is provided, then the moving average filter will have as many dimensions as the input
108 array. If a tuple is provided, it should contain as many elements as `dim_shape`. For example, the ``size=(1,
109 3)`` will convolve the input image with the filter ``[[1, 1, 1]] / 3``.
110
111 center: ~pyxu.operator.linop.filter.IndexSpec
112 (i_1, ..., i_D) index of the kernel's center.
113
114 `center` defines how a kernel is overlaid on inputs to produce outputs.
115 For odd `size`, it defaults to the central element (``center=size//2``).
116 For even `size` the desired center indices must be provided.
117
118 mode: str, list[str]
119 Boundary conditions.
120 Multiple forms are accepted:
121
122 * str: unique mode shared amongst dimensions.
123 Must be one of:
124
125 * 'constant' (default): zero-padding
126 * 'wrap'
127 * 'reflect'
128 * 'symmetric'
129 * 'edge'
130 * tuple[str, ...]: the `d`-th dimension uses `mode[d]` as boundary condition.
131
132 (See :py:func:`numpy.pad` for details.)
133 gpu: bool
134 Input NDArray type (`True` for GPU, `False` for CPU). Defaults to `False`.
135 dtype: DType
136 Working precision of the linear operator.
137
138 Returns
139 -------
140 op: OpT
141 MovingAverage
142
143 Example
144 -------
145
146 .. plot::
147
148 import matplotlib.pyplot as plt
149 import numpy as np
150 from pyxu.operator import MovingAverage
151
152 dim_shape = (11, 11)
153 image = np.zeros(dim_shape)
154 image[5, 5] = 1.
155
156 ma = MovingAverage(dim_shape, size=5)
157 out = ma(image)
158 plt.figure(figsize=(10, 5))
159 plt.subplot(121)
160 plt.imshow(image)
161 plt.colorbar()
162 plt.subplot(122)
163 plt.imshow(out)
164 plt.colorbar()
165
166 See Also
167 --------
168 :py:class:`~pyxu.operator.GaussianFilter`
169 """
170 size = _to_canonical_form(size, dim_shape)
171 if center is None:
172 assert all([s % 2 == 1 for s in size]), (
173 "Can only infer center for odd `size`s. For even `size`s, please " "provide the desired `center`s."
174 )
175 center = [s // 2 for s in size]
176
177 ndim, dtype, xp = _sanitize_inputs(dim_shape, dtype, gpu)
178
179 kernel = [xp.ones(s, dtype=dtype) for s in size] # use separable filters
180 scale = 1 / np.prod(size)
181
182 op = scale * pxls.Stencil(dim_shape=dim_shape, kernel=kernel, center=center, mode=mode)
183 op._name = "MovingAverage"
184 return op
185
186
187def GaussianFilter(
188 dim_shape: pxt.NDArrayShape,
189 sigma: typ.Union[typ.Tuple[pxt.Real], pxt.Real] = 1.0,
190 truncate: typ.Union[typ.Tuple[pxt.Real], pxt.Real] = 3.0,
191 order: typ.Union[typ.Tuple[pxt.Integer], pxt.Integer] = 0,
192 mode: ModeSpec = "constant",
193 sampling: typ.Union[pxt.Real, cabc.Sequence[pxt.Real, ...]] = 1,
194 gpu: bool = False,
195 dtype: typ.Optional[pxt.DType] = None,
196):
197 r"""
198 Multi-dimensional Gaussian filter.
199
200 Notes
201 -----
202 This operator performs a convolution between the input :math:`D`-dimensional NDArray :math:`\mathbf{x} \in
203 \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}` and a Gaussian :math:`D`-dimensional filter :math:`\mathbf{h} \in
204 \mathbb{R}^{\text{size} \times \cdots \times \text{size}}` using separable kernels for improved performance.
205
206 .. math::
207
208 y_{i} = \sum_{j \in \mathcal{N}_{i}} a_{j} x_{j} \exp(\frac{d_{ij}^{2}}{\sigma^{2}})
209
210 Where :math:`\mathcal{N}_{i}` is the set of elements neighbouring the :math:`i`-th element of the input array, and
211 :math:`a_{j} = \sum_{j \in \mathcal{N}_{i}} a_{j} \exp(\frac{d_{ij}^{2}}{\sigma^{2}})` normalizes the kernel to sum
212 to one.
213
214
215 Parameters
216 ----------
217 dim_shape: NDArrayShape
218 Shape of the input array.
219 sigma: float, tuple
220 Standard deviation of the Gaussian kernel.
221
222 If a scalar value is provided, then the Gaussian filter will have as many dimensions as the input array.
223 If a tuple is provided, it should contain as many elements as `dim_shape`.
224 Use ``0`` to prevent filtering in a given dimension.
225 For example, the ``sigma=(0, 3)`` will convolve the input image in its last dimension.
226 truncate: float, tuple
227 Truncate the filter at this many standard deviations.
228 Defaults to 3.0.
229 order: int, tuple
230 Gaussian derivative order.
231
232 Use ``0`` for the standard Gaussian kernel.
233 mode: str, list[str]
234 Boundary conditions.
235 Multiple forms are accepted:
236
237 * str: unique mode shared amongst dimensions.
238 Must be one of:
239
240 * 'constant' (default): zero-padding
241 * 'wrap'
242 * 'reflect'
243 * 'symmetric'
244 * 'edge'
245 * tuple[str, ...]: the `d`-th dimension uses `mode[d]` as boundary condition.
246
247 (See :py:func:`numpy.pad` for details.)
248 sampling: Real, list[Real]
249 Sampling step (i.e. distance between two consecutive elements of an array).
250 Defaults to 1.
251 gpu: bool
252 Input NDArray type (`True` for GPU, `False` for CPU).
253 Defaults to `False`.
254 dtype: DType
255 Working precision of the linear operator.
256
257 Returns
258 -------
259 op: OpT
260 GaussianFilter
261
262 Example
263 -------
264
265 .. plot::
266
267 import matplotlib.pyplot as plt
268 import numpy as np
269 from pyxu.operator import GaussianFilter
270
271 dim_shape = (11, 11)
272 image = np.zeros(dim_shape)
273 image[5, 5] = 1.
274
275 gaussian = GaussianFilter(dim_shape, sigma=3)
276 out = gaussian(image)
277 plt.figure(figsize=(10, 5))
278 plt.subplot(121)
279 plt.imshow(image)
280 plt.colorbar()
281 plt.subplot(122)
282 plt.imshow(out)
283 plt.colorbar()
284
285 See Also
286 --------
287 :py:class:`~pyxu.operator.MovingAverage`,
288 :py:class:`~pyxu.operator.DifferenceOfGaussians`
289 """
290
291 ndim, dtype, xp = _sanitize_inputs(dim_shape, dtype, gpu)
292 sigma = _to_canonical_form(sigma, dim_shape)
293 truncate = _to_canonical_form(truncate, dim_shape)
294 order = _to_canonical_form(order, dim_shape)
295 sampling = _to_canonical_form(sampling, dim_shape)
296
297 kernel = [
298 xp.array([1], dtype=dtype),
299 ] * len(dim_shape)
300 center = [0 for _ in range(len(dim_shape))]
301
302 for i, (sigma_, truncate_, order_, sampling_) in enumerate(zip(sigma, truncate, order, sampling)):
303 if sigma_:
304 sigma_pix = sigma_ / sampling_
305 radius = int(truncate_ * float(sigma_pix) + 0.5)
306 kernel[i] = xp.asarray(np.flip(scif._gaussian_kernel1d(sigma_pix, order_, radius)), dtype=dtype)
307 kernel[i] /= sampling_**order_
308 center[i] = radius
309
310 op = pxls.Stencil(dim_shape=dim_shape, kernel=kernel, center=center, mode=mode)
311 op._name = "GaussianFilter"
312 return op
313
314
[docs]
315def DifferenceOfGaussians(
316 dim_shape: pxt.NDArrayShape,
317 low_sigma=1.0,
318 high_sigma=None,
319 low_truncate=3.0,
320 high_truncate=3.0,
321 mode: ModeSpec = "constant",
322 sampling: typ.Union[pxt.Real, cabc.Sequence[pxt.Real, ...]] = 1,
323 gpu: bool = False,
324 dtype: typ.Optional[pxt.DType] = None,
325):
326 r"""
327 Multi-dimensional Difference of Gaussians filter.
328
329 Notes
330 -----
331
332 This operator uses the Difference of Gaussians (DoG) method to a :math:`D`-dimensional NDArray :math:`\mathbf{x} \in
333 \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}` using separable kernels for improved performance.
334 The DoG method blurs the input image with two Gaussian kernels with different sigma, and subtracts the more-blurred
335 signal from the less-blurred image.
336 This creates an output signal containing only the information from the original signal at the spatial scale
337 indicated by the two sigmas.
338
339 Parameters
340 ----------
341 dim_shape: NDArrayShape
342 Shape of the input array.
343 low_sigma: float, tuple
344 Standard deviation of the Gaussian kernel with smaller sigmas across all axes.
345
346 If a scalar value is provided, then the Gaussian filter will have as many dimensions as the input array.
347 If a tuple is provided, it should contain as many elements as `dim_shape`.
348 Use ``0`` to prevent filtering in a given dimension.
349 For example, the ``low_sigma=(0, 3)`` will convolve the input image in its last dimension.
350 high_sigma: float, tuple, None
351 Standard deviation of the Gaussian kernel with larger sigmas across all axes.
352 If ``None`` is given (default), sigmas for all axes are calculated as ``1.6 * low_sigma``.
353 low_truncate: float, tuple
354 Truncate the filter at this many standard deviations.
355 Defaults to 3.0.
356 high_truncate: float, tuple
357 Truncate the filter at this many standard deviations.
358 Defaults to 3.0.
359 mode: str, list[str]
360 Boundary conditions.
361 Multiple forms are accepted:
362
363 * str: unique mode shared amongst dimensions.
364 Must be one of:
365
366 * 'constant' (default): zero-padding
367 * 'wrap'
368 * 'reflect'
369 * 'symmetric'
370 * 'edge'
371 * tuple[str, ...]: the `d`-th dimension uses `mode[d]` as boundary condition.
372
373 (See :py:func:`numpy.pad` for details.)
374 sampling: Real, list[Real]
375 Sampling step (i.e. distance between two consecutive elements of an array).
376 Defaults to 1.
377 gpu: bool
378 Input NDArray type (`True` for GPU, `False` for CPU). Defaults to `False`.
379 dtype: DType
380 Working precision of the linear operator.
381
382 Returns
383 -------
384 op: OpT
385 DifferenceOfGaussians
386
387 Example
388 -------
389
390 .. plot::
391
392 import matplotlib.pyplot as plt
393 import numpy as np
394 from pyxu.operator import DoG
395
396 dim_shape = (11, 11)
397 image = np.zeros(dim_shape)
398 image[5, 5] = 1.
399
400 dog = DoG(dim_shape, low_sigma=3)
401 out = dog(image)
402 plt.figure(figsize=(10, 5))
403 plt.subplot(121)
404 plt.imshow(image)
405 plt.colorbar()
406 plt.subplot(122)
407 plt.imshow(out)
408 plt.colorbar()
409
410 See Also
411 --------
412 :py:class:`~pyxu.operator.Gaussian`,
413 :py:class:`~pyxu.operator.Sobel`,
414 :py:class:`~pyxu.operator.Prewitt`,
415 :py:class:`~pyxu.operator.Scharr`,
416 :py:class:`~pyxu.operator.StructureTensor`
417 """
418
419 low_sigma = _to_canonical_form(low_sigma, dim_shape)
420 if high_sigma is None:
421 high_sigma = tuple(s * 1.6 for s in low_sigma)
422
423 high_sigma = _to_canonical_form(high_sigma, dim_shape)
424 low_truncate = _to_canonical_form(low_truncate, dim_shape)
425 high_truncate = _to_canonical_form(high_truncate, dim_shape)
426
427 kwargs = {
428 "dim_shape": dim_shape,
429 "order": 0,
430 "mode": mode,
431 "gpu": gpu,
432 "dtype": dtype,
433 "sampling": sampling,
434 }
435 op_low = GaussianFilter(sigma=low_sigma, truncate=low_truncate, **kwargs)
436 op_high = GaussianFilter(sigma=high_sigma, truncate=high_truncate, **kwargs)
437 op = op_low - op_high
438 op._name = "DifferenceOfGaussians"
439 return op
440
441
442DoG = DifferenceOfGaussians #: Alias of :py:func:`~pyxu.operator.DifferenceOfGaussians`.
443
444
[docs]
445def Laplace(
446 dim_shape: pxt.NDArrayShape,
447 mode: ModeSpec = "constant",
448 sampling: typ.Union[pxt.Real, cabc.Sequence[pxt.Real, ...]] = 1,
449 gpu: bool = False,
450 dtype: typ.Optional[pxt.DType] = None,
451):
452 r"""
453 Multi-dimensional Laplace filter.
454
455 The implementation is based on second derivatives approximated via finite differences.
456
457 Notes
458 -----
459 This operator uses the applies the Laplace kernel :math:`[1 -2 1]` to a :math:`D`-dimensional NDArray
460 :math:`\mathbf{x} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}` using separable kernels for improved
461 performance. The Laplace filter is commonly used to find high-frequency components in the signal, such as for
462 example, the edges in an image.
463
464 Parameters
465 ----------
466 dim_shape: NDArrayShape
467 Shape of the input array.
468 mode: str, list[str]
469 Boundary conditions.
470 Multiple forms are accepted:
471
472 * str: unique mode shared amongst dimensions.
473 Must be one of:
474
475 * 'constant' (default): zero-padding
476 * 'wrap'
477 * 'reflect'
478 * 'symmetric'
479 * 'edge'
480 * tuple[str, ...]: the `d`-th dimension uses `mode[d]` as boundary condition.
481
482 (See :py:func:`numpy.pad` for details.)
483 sampling: Real, list[Real]
484 Sampling step (i.e. distance between two consecutive elements of an array).
485 Defaults to 1.
486 gpu: bool
487 Input NDArray type (`True` for GPU, `False` for CPU). Defaults to `False`.
488 dtype: DType
489 Working precision of the linear operator.
490
491 Returns
492 -------
493 op: OpT
494 DifferenceOfGaussians
495
496 Example
497 -------
498
499 .. plot::
500
501 import matplotlib.pyplot as plt
502 import numpy as np
503 from pyxu.operator import Laplace
504
505 dim_shape = (11, 11)
506 image = np.zeros(dim_shape)
507 image[5, 5] = 1.
508
509 laplace = Laplace(dim_shape)
510 out = laplace(image)
511 plt.figure(figsize=(10, 5))
512 plt.subplot(121)
513 plt.imshow(image)
514 plt.colorbar()
515 plt.subplot(122)
516 plt.imshow(out)
517 plt.colorbar()
518
519 See Also
520 --------
521 :py:class:`~pyxu.operator.Sobel`,
522 :py:class:`~pyxu.operator.Prewitt`,
523 :py:class:`~pyxu.operator.Scharr`,
524 """
525
526 ndim, dtype, xp = _sanitize_inputs(dim_shape, dtype, gpu)
527 sampling = _to_canonical_form(sampling, dim_shape)
528 centers = [[1 if i == dim else 0 for i in range(ndim)] for dim in range(ndim)]
529 kernels = [
530 xp.array([1.0, -2.0, 1.0], dtype=dtype).reshape([-1 if i == dim else 1 for i in range(ndim)]) / sampling[dim]
531 for dim in range(ndim)
532 ]
533 ops = [pxls.Stencil(dim_shape=dim_shape, kernel=k, center=c, mode=mode) for (k, c) in zip(kernels, centers)]
534 op = functools.reduce(lambda x, y: x + y, ops)
535 op._name = "Laplace"
536 return op
537
538
[docs]
539def Sobel(
540 dim_shape: pxt.NDArrayShape,
541 axis: typ.Optional[typ.Tuple] = None,
542 mode: ModeSpec = "constant",
543 sampling: typ.Union[pxt.Real, cabc.Sequence[pxt.Real, ...]] = 1,
544 gpu: bool = False,
545 dtype: typ.Optional[pxt.DType] = None,
546):
547 r"""
548 Multi-dimensional Sobel filter.
549
550 Notes
551 -----
552
553 This operator uses the applies the multi-dimensional Sobel filter to a :math:`D`-dimensional NDArray
554 :math:`\mathbf{x} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}` using separable kernels for improved
555 performance. The Sobel filter applies the following edge filter in the dimensions of interest: ``[1, 0, -1]`` and
556 the smoothing filter on the rest of dimensions: ``[1, 2, 1] / 4``. The Sobel filter is commonly used to find
557 high-frequency components in the signal, such as for example, the edges in an image.
558
559 Parameters
560 ----------
561 dim_shape: NDArrayShape
562 Shape of the input array.
563 axis: int, tuple
564 Compute the edge filter along this axis. If not provided, the edge magnitude is computed.
565
566 This is defined as: ``np.sqrt(sum([sobel(array, axis=i)**2 for i in range(array.ndim)]) / array.ndim)`` The
567 magnitude is also computed if axis is a sequence.
568
569 mode: str, list[str]
570 Boundary conditions.
571 Multiple forms are accepted:
572
573 * str: unique mode shared amongst dimensions.
574 Must be one of:
575
576 * 'constant' (default): zero-padding
577 * 'wrap'
578 * 'reflect'
579 * 'symmetric'
580 * 'edge'
581 * tuple[str, ...]: the `d`-th dimension uses `mode[d]` as boundary condition.
582
583 (See :py:func:`numpy.pad` for details.)
584 sampling: Real, list[Real]
585 Sampling step (i.e. distance between two consecutive elements of an array).
586 Defaults to 1.
587 gpu: bool
588 Input NDArray type (`True` for GPU, `False` for CPU). Defaults to `False`.
589 dtype: DType
590 Working precision of the linear operator.
591
592 Returns
593 -------
594 op: OpT
595 Sobel
596
597 Example
598 -------
599
600 .. plot::
601
602 import matplotlib.pyplot as plt
603 import numpy as np
604 from pyxu.operator import Sobel
605
606 dim_shape = (11, 11)
607 image = np.zeros(dim_shape)
608 image[5, 5] = 1.
609
610 sobel = Sobel(dim_shape)
611 out = sobel(image)
612 plt.figure(figsize=(10, 5))
613 plt.subplot(121)
614 plt.imshow(image)
615 plt.colorbar()
616 plt.subplot(122)
617 plt.imshow(out)
618 plt.colorbar()
619
620 See Also
621 --------
622 :py:class:`~pyxu.operator.Prewitt`,
623 :py:class:`~pyxu.operator.Scharr`,
624 """
625 smooth_kernel = np.array([1, 2, 1]) / 4
626 return _EdgeFilter(
627 dim_shape=dim_shape,
628 smooth_kernel=smooth_kernel,
629 filter_name="SobelFilter",
630 axis=axis,
631 mode=mode,
632 sampling=sampling,
633 gpu=gpu,
634 dtype=dtype,
635 )
636
637
[docs]
638def Prewitt(
639 dim_shape: pxt.NDArrayShape,
640 axis: typ.Optional[typ.Tuple] = None,
641 mode: ModeSpec = "constant",
642 sampling: typ.Union[pxt.Real, cabc.Sequence[pxt.Real, ...]] = 1,
643 gpu: bool = False,
644 dtype: typ.Optional[pxt.DType] = None,
645):
646 r"""
647 Multi-dimensional Prewitt filter.
648
649 Notes
650 -----
651
652 This operator uses the applies the multi-dimensional Prewitt filter to a :math:`D`-dimensional NDArray
653 :math:`\mathbf{x} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}` using separable kernels for improved
654 performance. The Prewitt filter applies the following edge filter in the dimensions of interest: ``[1, 0, -1]``,
655 and the smoothing filter on the rest of dimensions: ``[1, 1, 1] / 3``. The Prewitt filter is commonly used to find
656 high-frequency components in the signal, such as for example, the edges in an image.
657
658 Parameters
659 ----------
660 dim_shape: NDArrayShape
661 Shape of the input array.
662 axis: int, tuple
663 Compute the edge filter along this axis. If not provided, the edge magnitude is computed. This is defined as:
664
665 ``np.sqrt(sum([prewitt(array, axis=i)**2 for i in range(array.ndim)]) / array.ndim)`` The magnitude is also
666 computed if axis is a sequence.
667
668 mode: str, list[str]
669 Boundary conditions.
670 Multiple forms are accepted:
671
672 * str: unique mode shared amongst dimensions.
673 Must be one of:
674
675 * 'constant' (default): zero-padding
676 * 'wrap'
677 * 'reflect'
678 * 'symmetric'
679 * 'edge'
680 * tuple[str, ...]: the `d`-th dimension uses `mode[d]` as boundary condition.
681
682 (See :py:func:`numpy.pad` for details.)
683 sampling: Real, list[Real]
684 Sampling step (i.e. distance between two consecutive elements of an array).
685 Defaults to 1.
686 gpu: bool
687 Input NDArray type (`True` for GPU, `False` for CPU). Defaults to `False`.
688 dtype: DType
689 Working precision of the linear operator.
690
691 Returns
692 -------
693 op: OpT
694 Prewitt
695
696 Example
697 -------
698
699 .. plot::
700
701 import matplotlib.pyplot as plt
702 import numpy as np
703 from pyxu.operator import Prewitt
704
705 dim_shape = (11, 11)
706 image = np.zeros(dim_shape)
707 image[5, 5] = 1.
708
709 prewitt = Prewitt(dim_shape)
710 out = prewitt(image)
711 plt.figure(figsize=(10, 5))
712 plt.subplot(121)
713 plt.imshow(image)
714 plt.colorbar()
715 plt.subplot(122)
716 plt.imshow(out)
717 plt.colorbar()
718
719 See Also
720 --------
721 :py:class:`~pyxu.operator.Sobel`,
722 :py:class:`~pyxu.operator.Scharr`,
723 """
724 smooth_kernel = np.full((3,), 1 / 3)
725 return _EdgeFilter(
726 dim_shape=dim_shape,
727 smooth_kernel=smooth_kernel,
728 filter_name="Prewitt",
729 axis=axis,
730 mode=mode,
731 sampling=sampling,
732 gpu=gpu,
733 dtype=dtype,
734 )
735
736
[docs]
737def Scharr(
738 dim_shape: pxt.NDArrayShape,
739 axis: typ.Optional[typ.Tuple] = None,
740 mode: ModeSpec = "constant",
741 sampling: typ.Union[pxt.Real, cabc.Sequence[pxt.Real, ...]] = 1,
742 gpu: bool = False,
743 dtype: typ.Optional[pxt.DType] = None,
744):
745 r"""
746 Multi-dimensional Scharr filter.
747
748 Notes
749 -----
750
751 This operator uses the applies the multi-dimensional Scharr filter to a :math:`D`-dimensional NDArray
752 :math:`\mathbf{x} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}` using separable kernels for improved
753 performance. The Scharr filter applies the following edge filter in the dimensions of interest: ``[1, 0, -1]``, and
754 the smoothing filter on the rest of dimensions: ``[3, 10, 3] / 16``. The Scharr filter is commonly used to find
755 high-frequency components in the signal, such as for example, the edges in an image.
756
757 Parameters
758 ----------
759 dim_shape: NDArrayShape
760 Shape of the input array.
761 axis: int, tuple
762 Compute the edge filter along this axis. If not provided, the edge magnitude is computed. This is defined as:
763
764 ``np.sqrt(sum([scharr(array, axis=i)**2 for i in range(array.ndim)]) / array.ndim)`` The magnitude is also
765 computed if axis is a sequence.
766 mode: str, list[str]
767 Boundary conditions.
768 Multiple forms are accepted:
769
770 * str: unique mode shared amongst dimensions.
771 Must be one of:
772
773 * 'constant' (default): zero-padding
774 * 'wrap'
775 * 'reflect'
776 * 'symmetric'
777 * 'edge'
778 * tuple[str, ...]: the `d`-th dimension uses `mode[d]` as boundary condition.
779
780 (See :py:func:`numpy.pad` for details.)
781 sampling: Real, list[Real]
782 Sampling step (i.e. distance between two consecutive elements of an array).
783 Defaults to 1.
784 gpu: bool
785 Input NDArray type (`True` for GPU, `False` for CPU). Defaults to `False`.
786 dtype: DType
787 Working precision of the linear operator.
788
789 Returns
790 -------
791 op: OpT
792 Scharr
793
794 Example
795 -------
796
797 .. plot::
798
799 import matplotlib.pyplot as plt
800 import numpy as np
801 from pyxu.operator import Scharr
802
803 dim_shape = (11, 11)
804 image = np.zeros(dim_shape)
805 image[5, 5] = 1.
806
807 scharr = Scharr(dim_shape)
808 out = scharr(image)
809 plt.figure(figsize=(10, 5))
810 plt.subplot(121)
811 plt.imshow(image)
812 plt.colorbar()
813 plt.subplot(122)
814 plt.imshow(out)
815 plt.colorbar()
816
817 See Also
818 --------
819 :py:class:`~pyxu.operator.Sobel`,
820 :py:class:`~pyxu.operator.Prewitt`,
821 """
822 smooth_kernel = np.array([3, 10, 3]) / 16
823 return _EdgeFilter(
824 dim_shape=dim_shape,
825 smooth_kernel=smooth_kernel,
826 filter_name="Scharr",
827 axis=axis,
828 mode=mode,
829 sampling=sampling,
830 gpu=gpu,
831 dtype=dtype,
832 )
833
834
835def _EdgeFilter(
836 dim_shape: pxt.NDArrayShape,
837 smooth_kernel: KernelSpec,
838 filter_name: str,
839 axis: typ.Optional[typ.Tuple] = None,
840 mode: ModeSpec = "constant",
841 sampling: typ.Union[pxt.Real, cabc.Sequence[pxt.Real, ...]] = 1,
842 gpu: bool = False,
843 dtype: typ.Optional[pxt.DType] = None,
844):
845 from pyxu.operator import Sqrt, Square
846
847 square = Square(dim_shape=dim_shape)
848 sqrt = Sqrt(dim_shape=dim_shape)
849
850 ndim, dtype, xp = _sanitize_inputs(dim_shape, dtype, gpu)
851 sampling = _to_canonical_form(sampling, dim_shape)
852
853 axes = _get_axes(axis, ndim)
854
855 return_magnitude = len(axes) > 1
856
857 op_list = []
858 for edge_dim in axes:
859 kernel = [xp.array(1, dtype=dtype)] * len(dim_shape)
860 center = np.ones(len(dim_shape), dtype=int)
861 # We define the kernel reversed compared to Scipy or Skimage because we use correlation instead of convolution
862 kernel[edge_dim] = xp.array([-1, 0, 1], dtype=dtype) / sampling[edge_dim]
863 smooth_axes = list(set(range(ndim)) - {edge_dim})
864 for smooth_dim in smooth_axes:
865 kernel[smooth_dim] = xp.asarray(smooth_kernel, dtype=dtype) / sampling[smooth_dim]
866
867 if return_magnitude:
868 op_list.append(square * pxls.Stencil(dim_shape=dim_shape, kernel=kernel, center=center, mode=mode))
869 else:
870 op_list.append(pxls.Stencil(dim_shape=dim_shape, kernel=kernel, center=center, mode=mode))
871
872 op = functools.reduce(lambda x, y: x + y, op_list)
873 if return_magnitude:
874 op = (1 / np.sqrt(ndim)) * (sqrt * op)
875
876 op._name = filter_name
877 return op
878
879
[docs]
880class StructureTensor(pxa.DiffMap):
881 r"""
882 Structure tensor operator.
883
884 Notes
885 -----
886 The Structure Tensor, also known as the second-order moment tensor or the inertia tensor, is a matrix derived from
887 the gradient of a function. It describes the distribution of the gradient (i.e., its prominent directions) in a
888 specified neighbourhood around a point, and the degree to which those directions are coherent. The structure tensor
889 of a :math:`D`-dimensional signal :math:`\mathbf{f} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}` can be
890 written as:
891
892 .. math::
893
894 \mathbf{S}_\sigma \mathbf{f} = \mathbf{g}_{\sigma} * \nabla\mathbf{f} (\nabla\mathbf{f})^{\top} = \mathbf{g}_{\sigma} *
895 \begin{bmatrix}
896 \left( \dfrac{ \partial\mathbf{f} }{ \partial x_{0} } \right)^2 & \dfrac{ \partial^{2}\mathbf{f} }{ \partial x_{0}\,\partial x_{1} } & \cdots & \dfrac{ \partial\mathbf{f} }{ \partial x_{0} } \dfrac{ \partial\mathbf{f} }{ \partial x_{D-1} } \\
897 \dfrac{ \partial\mathbf{f} }{ \partial x_{1} } \dfrac{ \partial\mathbf{f} }{ \partial x_{0} } & \left( \dfrac{ \partial\mathbf{f} }{ \partial x_{1} }\right)^2 & \cdots & \dfrac{ \partial\mathbf{f} }{ \partial x_{1} } \dfrac{ \partial\mathbf{f} }{ \partial x_{D-1} } \\
898 \vdots & \vdots & \ddots & \vdots \\
899 \dfrac{ \partial\mathbf{f} }{ \partial x_{D-1} } \dfrac{ \partial\mathbf{f} }{ \partial x_{0} } & \dfrac{ \partial\mathbf{f} }{ \partial x_{D-1} } \dfrac{ \partial\mathbf{f} }{ \partial x_{1} } & \cdots & \left( \dfrac{ \partial\mathbf{f} }{ \partial x_{D-1}} \right)^2
900 \end{bmatrix},
901
902 where :math:`\mathbf{g}_{\sigma} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}` is a discrete Gaussian filter
903 with standard variation :math:`\sigma` with which a convolution is performed elementwise.
904
905 However, due to the symmetry of the structure tensor, only the upper triangular part is computed in practice:
906
907 .. math::
908
909 \mathbf{H}_{\mathbf{v}_1, \ldots ,\mathbf{v}_m} \mathbf{f} = \mathbf{g}_{\sigma} * \begin{bmatrix}
910 \left( \dfrac{ \partial\mathbf{f} }{ \partial x_{0} } \right)^2 \\
911 \dfrac{ \partial^{2}\mathbf{f} }{ \partial x_{0}\,\partial x_{1} } \\
912 \vdots \\
913 \left( \dfrac{ \partial\mathbf{f} }{ \partial x_{D-1}} \right)^2
914 \end{bmatrix} \mathbf{f} \in \mathbb{R}^{\frac{D (D-1)}{2} \times N_0 \times \cdots \times N_{D-1}}
915
916 Remark
917 ------
918 In case of using the finite differences (`diff_type="fd"`), the finite difference scheme defaults to `central` (see
919 :py:class:`~pyxu.operator.PartialDerivative`).
920
921 Example
922 -------
923
924 .. plot::
925
926 import numpy as np
927 import matplotlib.pyplot as plt
928 from pyxu.operator import StructureTensor
929 from pyxu.util.misc import peaks
930
931 # Define input image
932 n = 1000
933 x = np.linspace(-3, 3, n)
934 xx, yy = np.meshgrid(x, x)
935 image = peaks(xx, yy)
936 nsamples = 2
937 dim_shape = image.shape # (1000, 1000)
938 images = np.tile(image, (nsamples, 1, 1))
939 # Instantiate structure tensor operator
940 structure_tensor = StructureTensor(dim_shape=dim_shape)
941
942 outputs = structure_tensor(images)
943 print(outputs.shape) # (2, 3, 1000, 1000)
944 # Plot
945 plt.figure()
946 plt.imshow(images[0])
947 plt.colorbar()
948 plt.title("Image")
949 plt.axis("off")
950
951 plt.figure()
952 plt.imshow(outputs[0][0])
953 plt.colorbar()
954 plt.title(r"$\hat{S}_{xx}$")
955 plt.axis("off")
956
957 plt.figure()
958 plt.imshow(outputs[0][1])
959 plt.colorbar()
960 plt.title(r"$\hat{S}_{xy}$")
961 plt.axis("off")
962
963 plt.figure()
964 plt.imshow(outputs[0][2])
965 plt.colorbar()
966 plt.title(r"$\hat{S}_{yy}$")
967 plt.axis("off")
968
969 See Also
970 --------
971 :py:class:`~pyxu.operator.PartialDerivative`,
972 :py:class:`~pyxu.operator.Gradient`,
973 :py:class:`~pyxu.operator.Hessian`
974 """
975
976 def __init__(
977 self,
978 dim_shape: pxt.NDArrayShape,
979 diff_method="fd",
980 smooth_sigma: typ.Union[pxt.Real, tuple[pxt.Real, ...]] = 1.0,
981 smooth_truncate: typ.Union[pxt.Real, tuple[pxt.Real, ...]] = 3.0,
982 mode: ModeSpec = "constant",
983 sampling: typ.Union[pxt.Real, tuple[pxt.Real, ...]] = 1,
984 gpu: bool = False,
985 dtype: typ.Optional[pxt.DType] = None,
986 parallel: bool = False,
987 **diff_kwargs,
988 ):
989 ndim = len(dim_shape)
990 ntriu = (ndim * (ndim + 1)) // 2
991 super().__init__(dim_shape=dim_shape, codim_shape=(ntriu,) + dim_shape)
992 self.directions = tuple(
993 list(_) for _ in itertools.combinations_with_replacement(np.arange(len(dim_shape)).astype(int), 2)
994 )
995
996 if diff_method == "fd":
997 diff_kwargs.update({"scheme": diff_kwargs.pop("scheme", "central")})
998
999 self.grad = pxld.Gradient(
1000 dim_shape=dim_shape,
1001 directions=None,
1002 diff_method=diff_method,
1003 mode=mode,
1004 gpu=gpu,
1005 dtype=dtype,
1006 sampling=sampling,
1007 parallel=parallel,
1008 **diff_kwargs,
1009 )
1010
1011 if smooth_sigma:
1012 self.smooth = GaussianFilter(
1013 dim_shape=dim_shape,
1014 sigma=smooth_sigma,
1015 truncate=smooth_truncate,
1016 order=0,
1017 mode=mode,
1018 sampling=sampling,
1019 gpu=gpu,
1020 dtype=dtype,
1021 )
1022 else:
1023 self.smooth = pxlb.IdentityOp(dim_shape=dim_shape)
1024
[docs]
1025 def apply(self, arr):
1026 xp = pxu.get_array_module(arr)
1027 sh = arr.shape[: -self.dim_rank]
1028 grad = self.grad(arr)
1029
1030 def slice_(ax):
1031 return (slice(None),) * len(sh) + (slice(ax, ax + 1),)
1032
1033 return xp.concatenate(
1034 [self.smooth((grad[slice_(i)] * grad[slice_(j)])) for i, j in self.directions],
1035 axis=len(sh),
1036 ).reshape(*sh, len(self.directions), *self.dim_shape)
