Merge pull request #592 from hakonanes/add-downsample-method

Add method to downsample EBSD patterns while maintaining data type range

CHANGELOG.rst

@@ -18,6 +18,8 @@ Unreleased
 
 Added
 -----
+- Downsampling of EBSD patterns which maintain the data type by also rescaling to the
+  data type range. (`#592 <https://github.com/pyxem/kikuchipy/pull/592>`_)
 - Method to get a PyEBSDIndex ``EBSDIndexer`` instance from an ``EBSDDetector``,
   convenient for either indexing with PyEBSDIndex or for use with kikuchipy.
   (`#590 <https://github.com/pyxem/kikuchipy/pull/590>`_)

examples/pattern_processing/pattern_binning.py

@@ -2,8 +2,14 @@
 Pattern binning
 ===============
 
-This example shows how to bin :class:`~kikuchipy.signals.EBSD` patterns using HyperSpy's
+This example shows how to bin :class:`~kikuchipy.signals.EBSD` patterns using
+:meth:`~kikuchipy.signals.EBSD.downsample` or HyperSpy's
 :meth:`~kikuchipy.signals.EBSD.rebin` (see :ref:`hyperspy:rebin-label` for details).
+
+.. note::
+
+    In general, better contrast is obtained by removing the static (and dynamic)
+    background prior to binning instead of after it.
 """
 
 import hyperspy.api as hs
@@ -19,34 +25,29 @@
 print(s.detector)
 
 ########################################################################################
-# Rebin by passing the new shape (use ``(1, 1, 60, 60)`` if binning a 2D map). Note how
-# the :attr:`~kikuchipy.signals.EBSD.static_background` and
+# Downsample by a factor of 8 while maintaining the data type (achieved by rescaling the
+# pattern intensity). Note how the :attr:`~kikuchipy.signals.EBSD.static_background` and
 # :attr:`~kikuchipy.signals.EBSD.detector` attributes are updated.
 
-s2 = s.rebin(new_shape=(60, 60))
-
-_ = hs.plot.plot_images(
-    [s, s2],
-    axes_decor="off",
-    tight_layout=True,
-    label=None,
-    colorbar=False,
-)
+s2 = s.deepcopy()
+s2.downsample(8)
+_ = hs.plot.plot_images([s, s2], axes_decor="off", tight_layout=True, label=None)
 print(s2.static_background.shape)
 print(s2.detector)
 
 ########################################################################################
-# Rebin by passing the new:old pixel ratio (again, use ``(1, 1, 8, 8)`` for a 2D map)
+# Rebin by passing the new shape (use ``(1, 1, 60, 60)`` if binning a 2D map). Note how
+# the pattern is not rescaled and the data type is cast to either ``int64`` or
+# ``float64`` depending on the initial data type.
 
-s3 = s.rebin(scale=(8, 8))
-print(np.allclose(s2.data, s3.data))
-print(s3.static_background.shape)
-print(s3.detector)
+s3 = s.rebin(new_shape=(60, 60))
+print(s3.data.dtype)
+print(s3.data.min(), s3.data.max())
 
 ########################################################################################
-# Notice how ``rebin()`` casts the data to ``uint64``, increasing the memory use by a
-# factor of eight (so be careful...). Rescale intensities with
-# :meth:`~kikuchipy.signals.EBSD.rescale_intensity` if desirable.
+# The latter method is more flexible in that it allows for different binning factors in
+# each axis, the factors are not restricted to being integers and the factors do not
+# have to be divisors of the initial signal shape.
 
-print(s.data.dtype)
-print(s3.data.dtype)
+s4 = s.rebin(scale=(8, 9))
+print(s4.data.shape)

kikuchipy/pattern/__init__.py

@@ -15,7 +15,7 @@
 # You should have received a copy of the GNU General Public License
 # along with kikuchipy. If not, see <http://www.gnu.org/licenses/>.
 
-"""Single and chunk pattern processing used by signals."""
+"""Single pattern processing (used by signals)."""
 
 __all__ = [
     "chunk",

kikuchipy/pattern/_pattern.py

@@ -143,7 +143,7 @@ def _normalize(patterns: np.ndarray, axis: Union[int, tuple]) -> np.ndarray:
     return patterns / patterns_norm_squared
 
 
-@njit("float32[:](float32[:],bool_[:])", cache=True, nogil=True, fastmath=True)
+@njit("float32[:](float32[:], bool_[:])", cache=True, nogil=True, fastmath=True)
 def _mask_pattern(pattern: np.ndarray, mask: np.ndarray) -> np.ndarray:
     # Used in refinement solvers
     return pattern[mask].reshape(-1)
@@ -751,3 +751,37 @@ def _get_image_quality_numba(
     inertia = np.sum(spectrum * frequency_vectors) / np.sum(spectrum)
 
     return 1 - (inertia / inertia_max)
+
+
+@njit("float32[:, :](float32[:, :], int64)", cache=True, fastmath=True, nogil=True)
+def _bin2d(pattern: np.ndarray, factor: int) -> np.ndarray:
+    n_rows_new = pattern.shape[0] // factor
+    n_cols_new = pattern.shape[1] // factor
+
+    new_pattern = np.zeros((n_rows_new, n_cols_new), dtype=pattern.dtype)
+
+    for r in range(n_rows_new):
+        for rr in range(r * factor, (r + 1) * factor):
+            for c in range(n_cols_new):
+                value = new_pattern[r, c]
+                for cc in range(c * factor, (c + 1) * factor):
+                    value += pattern[rr, cc]
+                new_pattern[r, c] = value
+
+    return new_pattern
+
+
+@njit(cache=True, fastmath=True, nogil=True)
+def _downsample2d(
+    pattern: np.ndarray,
+    factor: int,
+    omin: Union[int, float],
+    omax: Union[int, float],
+    dtype_out: np.dtype,
+) -> np.ndarray:
+    pattern = pattern.astype(np.dtype("float32"))
+    binned_pattern = _bin2d(pattern, factor)
+    imin = binned_pattern.min()
+    imax = binned_pattern.max()
+    rescaled_pattern = _rescale_with_min_max(binned_pattern, imin, imax, omin, omax)
+    return rescaled_pattern.astype(dtype_out)

kikuchipy/pattern/chunk.py

@@ -15,14 +15,8 @@
 # You should have received a copy of the GNU General Public License
 # along with kikuchipy. If not, see <http://www.gnu.org/licenses/>.
 
-"""Functions for operating on :class:`numpy.ndarray` or
+"""Private functions for operating on :class:`numpy.ndarray` or
 :class:`dask.array.Array` chunks of EBSD patterns.
-
-.. warning::
-
-    This module will be become private for internal use only in v0.7.
-    If you need to process multiple EBSD patterns at once, please do
-    this using :class:`~kikuchipy.signals.EBSD`.
 """
 
 from typing import Union, Tuple, List
@@ -297,55 +291,6 @@ def normalize_intensity(
     return normalized_patterns
 
 
-def average_neighbour_patterns(
-    patterns: np.ndarray,
-    window_sums: np.ndarray,
-    window: Union[np.ndarray, Window],
-    dtype_out: Union[str, np.dtype, type, None] = None,
-) -> np.ndarray:
-    """Average a chunk of patterns with its neighbours within a window.
-
-    The amount of averaging is specified by the window coefficients.
-    All patterns are averaged with the same window. Map borders are
-    extended with zeros. Resulting pattern intensities are rescaled
-    to fill the input patterns' data type range individually.
-
-    Parameters
-    ----------
-    patterns
-        Patterns to average, with some overlap with surrounding chunks.
-    window_sums
-        Sum of window data for each image.
-    window
-        Averaging window.
-    dtype_out
-        Data type of averaged patterns. If None (default), it is set to
-        the same data type as the input patterns.
-
-    Returns
-    -------
-    averaged_patterns
-        Averaged patterns.
-    """
-    if dtype_out is None:
-        dtype_out = patterns.dtype
-    else:
-        dtype_out = np.dtype(dtype_out)
-
-    # Correlate patterns with window
-    correlated_patterns = correlate(patterns, weights=window, mode="constant", cval=0)
-
-    # Divide convolved patterns by number of neighbours averaged with
-    averaged_patterns = np.empty_like(correlated_patterns, dtype=dtype_out)
-    for nav_idx in np.ndindex(patterns.shape[:-2]):
-        averaged_patterns[nav_idx] = pattern_processing.rescale_intensity(
-            pattern=correlated_patterns[nav_idx] / window_sums[nav_idx],
-            dtype_out=dtype_out,
-        )
-
-    return averaged_patterns
-
-
 def _average_neighbour_patterns(
     patterns: np.ndarray,
     window_sums: np.ndarray,

kikuchipy/pattern/tests/test_chunk.py

@@ -294,52 +294,6 @@ def test_adaptive_histogram_equalization_chunk(self, dummy_signal):
         assert np.allclose(equalized_patterns[0, 0].compute(), ADAPT_EQ_UINT8)
 
 
-class TestAverageNeighbourPatternsChunk:
-    @pytest.mark.parametrize("dtype_in", [None, np.uint8])
-    def test_average_neighbour_patterns_chunk(self, dummy_signal, dtype_in):
-        w = Window()
-
-        # Get array to operate on
-        dask_array = get_dask_array(dummy_signal)
-        dtype_out = dask_array.dtype
-
-        # Get sum of window data for each image
-        nav_shape = dummy_signal.axes_manager.navigation_shape
-        w_sums = convolve(
-            input=np.ones(nav_shape[::-1], dtype=int),
-            weights=w.data,
-            mode="constant",
-            cval=0,
-        )
-
-        # Add signal dimensions to arrays to enable their use with
-        # Dask's map_blocks()
-        sig_dim = dummy_signal.axes_manager.signal_dimension
-        nav_dim = dummy_signal.axes_manager.navigation_dimension
-        for _ in range(sig_dim):
-            w_sums = np.expand_dims(w_sums, axis=w_sums.ndim)
-            w = np.expand_dims(w, axis=w.ndim)
-        w_sums = da.from_array(
-            w_sums, chunks=dask_array.chunks[:nav_dim] + (1,) * sig_dim
-        )
-
-        averaged_patterns = dask_array.map_blocks(
-            func=chunk.average_neighbour_patterns,
-            window_sums=w_sums,
-            window=w,
-            dtype_out=dtype_in,
-            dtype=dtype_out,
-        )
-
-        answer = np.array(
-            [255, 109, 218, 218, 36, 236, 255, 36, 0], dtype=np.uint8
-        ).reshape((3, 3))
-
-        # Check for correct data type and gives expected output intensities
-        assert averaged_patterns.dtype == dtype_out
-        assert np.allclose(averaged_patterns[0, 0].compute(), answer)
-
-
 class TestFFTFilterChunk:
     @pytest.mark.parametrize(
         "shift, transfer_function, kwargs, dtype_out, expected_spectrum_sum",

kikuchipy/pattern/tests/test_pattern.py

@@ -15,6 +15,7 @@
 # You should have received a copy of the GNU General Public License
 # along with kikuchipy. If not, see <http://www.gnu.org/licenses/>.
 
+import dask.array as da
 import numpy as np
 import pytest
 from scipy.fft import fft2
@@ -30,6 +31,8 @@
     normalize_intensity,
     rescale_intensity,
     remove_dynamic_background,
+    _bin2d,
+    _downsample2d,
     _dynamic_background_frequency_space_setup,
     _get_image_quality_numba,
     _mask_pattern,
@@ -524,3 +527,24 @@ def test_mask_pattern_numba(self):
         assert p_masked.size == 90
         assert np.isclose(p_masked.mean(), 54.5)
         assert np.allclose(p_masked, p_masked2)
+
+
+class TestDownsample:
+    def test_downsample_numba(self):
+        data = np.arange(1000, dtype="float32").reshape((20, 50))
+        data_da = da.from_array(data)
+
+        data_binned_da = da.coarsen(np.sum, data_da, {0: 2, 1: 2})
+        data_binned_kp = _bin2d(data, 2)
+        data_binned_kp2 = _bin2d.py_func(data, 2)
+
+        assert np.allclose(data_binned_da, data_binned_kp)
+        assert np.allclose(data_binned_kp, data_binned_kp2)
+
+        data_downsampled = _rescale_with_min_max(
+            data_binned_kp, data_binned_kp.min(), data_binned_kp.max(), omin=-1, omax=1
+        )
+        data_downsampled2 = _downsample2d(data, 2, -1, 1, np.float32)
+        data_downsampled3 = _downsample2d.py_func(data, 2, -1, 1, np.float32)
+        assert np.allclose(data_downsampled, data_downsampled2)
+        assert np.allclose(data_downsampled, data_downsampled3)