Repository: pyxem/kikuchipy (language: python)

Issue:
Alternative downsampling to the one from HyperSpy
We could provide an alternative downsampling to `EBSD.rebin()` from HyperSpy, e.g. based on https://github.com/wiai/xcdskd/blob/master/src/aloe/image/downsample.py.

This should be as simple as possible. The user can provide a binning factor, and we can crop if the pattern shape is not a multiple of the factor. Can be called `EBSD.downsample()`.


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 1387a0a8..93f34698 100644
--- a/CHANGELOG.rst
+++ b/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>`_)
diff --git a/examples/pattern_processing/pattern_binning.py b/examples/pattern_processing/pattern_binning.py
index 0734a1bc..e7429b6e 100644
--- a/examples/pattern_processing/pattern_binning.py
+++ b/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)
diff --git a/kikuchipy/pattern/__init__.py b/kikuchipy/pattern/__init__.py
index dac76384..b4aec23a 100644
--- a/kikuchipy/pattern/__init__.py
+++ b/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",
diff --git a/kikuchipy/pattern/_pattern.py b/kikuchipy/pattern/_pattern.py
index 9a829673..708ebc42 100644
--- a/kikuchipy/pattern/_pattern.py
+++ b/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)
diff --git a/kikuchipy/pattern/chunk.py b/kikuchipy/pattern/chunk.py
index 35d8b369..9699bcf3 100644
--- a/kikuchipy/pattern/chunk.py
+++ b/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,
diff --git a/kikuchipy/signals/ebsd.py b/kikuchipy/signals/ebsd.py
index b547225b..22799a1c 100644
--- a/kikuchipy/signals/ebsd.py
+++ b/kikuchipy/signals/ebsd.py
@@ -65,6 +65,7 @@
 from kikuchipy.pattern._pattern import (
     fft_frequency_vectors,
     fft_filter,
+    _downsample2d,
     _dynamic_background_frequency_space_setup,
     _get_image_quality,
     _remove_static_background_subtract,
@@ -834,1534 +835,1615 @@ def adaptive_histogram_equalization(
         else:
             self.data = equalized_patterns
 
-    def get_image_quality(
+    def fft_filter(
         self,
-        normalize: bool = True,
+        transfer_function: Union[np.ndarray, Window],
+        function_domain: str,
+        shift: bool = False,
         show_progressbar: Optional[bool] = None,
-    ) -> Union[np.ndarray, da.Array]:
-        """Compute the image quality map of patterns in an EBSD scan.
+    ) -> None:
+        """Filter an EBSD scan inplace in the frequency domain.
 
-        The image quality :math:`Q` is calculated based on the procedure
-        defined by Krieger Lassen :cite:`lassen1994automated`.
+        Patterns are transformed via the Fast Fourier Transform (FFT) to
+        the frequency domain, where their spectrum is multiplied by the
+        ``transfer_function``, and the filtered spectrum is subsequently
+        transformed to the spatial domain via the inverse FFT (IFFT).
+        Filtered patterns are rescaled to input data type range.
+
+        Note that if ``function_domain`` is ``"spatial"``, only real
+        valued FFT and IFFT is used.
 
         Parameters
         ----------
-        normalize
-            Whether to normalize patterns to a mean of zero and standard
-            deviation of 1 before calculating the image quality. Default
-            is ``True``.
+        transfer_function
+            Filter to apply to patterns. This can either be a transfer
+            function in the frequency domain of pattern shape or a
+            kernel in the spatial domain. What is passed is determined
+            from ``function_domain``.
+        function_domain
+            Options are ``"frequency"`` and ``"spatial"``, indicating,
+            respectively, whether the filter function passed to
+            ``filter_function`` is a transfer function in the frequency
+            domain or a kernel in the spatial domain.
+        shift
+            Whether to shift the zero-frequency component to the center.
+            Default is ``False``. This is only used when
+            ``function_domain="frequency"``.
         show_progressbar
             Whether to show a progressbar. If not given, the value of
             :obj:`hyperspy.api.preferences.General.show_progressbar`
             is used.
 
-        Returns
-        -------
-        image_quality_map
-            Image quality map of same shape as navigation axes. This is
-            a Dask array if the signal is lazy.
-
         See Also
         --------
-        kikuchipy.pattern.get_image_quality
+        kikuchipy.filters.Window
 
         Examples
         --------
-        Load an example dataset, remove the static and dynamic
-        background and compute :math:`Q`
+        Applying a Gaussian low pass filter with a cutoff frequency of
+        20:
 
         >>> import kikuchipy as kp
         >>> s = kp.data.nickel_ebsd_small()
-        >>> s
-        <EBSD, title: patterns Scan 1, dimensions: (3, 3|60, 60)>
-        >>> s.remove_static_background()
-        >>> s.remove_dynamic_background()
-        >>> iq = s.get_image_quality()
-        >>> iq
-        array([[0.19935645, 0.16657268, 0.18803978],
-               [0.19040637, 0.1616931 , 0.17834103],
-               [0.19411428, 0.16031407, 0.18413563]], dtype=float32)
+        >>> pattern_shape = s.axes_manager.signal_shape[::-1]
+        >>> w = kp.filters.Window(
+        ...     "lowpass", cutoff=20, shape=pattern_shape
+        ... )
+        >>> s.fft_filter(
+        ...     transfer_function=w,
+        ...     function_domain="frequency",
+        ...     shift=True,
+        ... )
         """
-        # Calculate frequency vectors
-        sx, sy = self.axes_manager.signal_shape
-        frequency_vectors = fft_frequency_vectors((sy, sx))
-        inertia_max = np.sum(frequency_vectors) / (sy * sx)
+        dtype_out = self.data.dtype.type
 
-        image_quality_map = self.map(
-            _get_image_quality,
-            show_progressbar=show_progressbar,
-            parallel=True,
-            inplace=False,
-            output_dtype=np.float32,
-            normalize=normalize,
-            frequency_vectors=frequency_vectors,
-            inertia_max=inertia_max,
+        dtype = np.float32
+        dask_array = get_dask_array(signal=self, dtype=dtype)
+
+        kwargs = {}
+        if function_domain == "frequency":
+            filter_func = fft_filter
+            kwargs["shift"] = shift
+        elif function_domain == "spatial":
+            filter_func = _fft_filter  # Barnes
+            kwargs["window_shape"] = transfer_function.shape
+
+            # FFT filter setup
+            (
+                kwargs["fft_shape"],
+                transfer_function,  # Padded
+                kwargs["offset_before_fft"],
+                kwargs["offset_after_ifft"],
+            ) = _fft_filter_setup(
+                image_shape=self.axes_manager.signal_shape[::-1],
+                window=transfer_function,
+            )
+        else:
+            function_domains = ["frequency", "spatial"]
+            raise ValueError(f"{function_domain} must be either of {function_domains}.")
+
+        filtered_patterns = dask_array.map_blocks(
+            func=chunk.fft_filter,
+            filter_func=filter_func,
+            transfer_function=transfer_function,
+            dtype_out=dtype_out,
+            dtype=dtype_out,
+            **kwargs,
         )
 
-        return image_quality_map.data
+        # Overwrite signal patterns
+        if not self._lazy:
+            pbar = ProgressBar()
+            if show_progressbar or (
+                show_progressbar is None and hs.preferences.General.show_progressbar
+            ):
+                pbar.register()
 
-    def hough_indexing(
+            filtered_patterns.store(self.data, compute=True)
+
+            try:
+                pbar.unregister()
+            except KeyError:
+                pass
+        else:
+            self.data = filtered_patterns
+
+    def average_neighbour_patterns(
         self,
-        phase_list: PhaseList,
-        indexer: "EBSDIndexer",
-        chunksize: int = 528,
-        verbose: int = 1,
-        return_index_data: bool = False,
-        return_band_data: bool = False,
-    ) -> Union[
-        CrystalMap,
-        Tuple[CrystalMap, np.ndarray],
-        Tuple[CrystalMap, np.ndarray, np.ndarray],
-    ]:
-        """Index patterns by Hough indexing using :mod:`pyebsdindex`.
+        window: Union[str, np.ndarray, da.Array, Window] = "circular",
+        window_shape: Tuple[int, ...] = (3, 3),
+        show_progressbar: Optional[bool] = None,
+        **kwargs,
+    ) -> None:
+        """Average patterns inplace with its neighbours within a window.
 
-        See :meth:`~pyebsdindex.ebsd_index.EBSDIndexer` and
-        :meth:`~pyebsdindex.ebsd_index.EBSDIndexer.index_pats` for
-        details.
+        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.
 
-        Currently, PyEBSDIndex only supports indexing with a single
-        mean projection center (PC).
+        Averaging is accomplished by correlating the window with the
+        extended array of patterns using
+        :func:`scipy.ndimage.correlate`.
 
         Parameters
         ----------
-        phase_list
-            List of phases. The list can only contain one face-centered
-            cubic (FCC) phase, one body-centered cubic (BCC) phase or
-            both types.
-        indexer
-            PyEBSDIndex EBSD indexer instance of which the
-            :meth:`~pyebsdindex.ebsd_index.EBSDIndexer.index_pats`
-            method is called. Its `phaselist` must be compatible with
-            the given ``phase_list``, and the ``indexer.vendor`` must be
-            ``"KIKUCHIPY"``. An indexer can be obtained with
-            :meth:`~kikuchipy.detectors.EBSDDetector.get_indexer`.
-        chunksize
-            Number of patterns to index at a time. Default is the
-            minimum of 528 or the number of patterns in the signal.
-            Increasing the chunksize may give faster indexing but
-            increases memory use.
-        verbose
-            Which information to print from PyEBSDIndex. Options are
-            0 - no output, 1 - timings (default), 2 - timings and the
-            Hough transform of the first pattern with detected bands
-            highlighted.
-        return_index_data
-            Whether to return the index data array returned from
-            ``EBSDIndexer.index_pats()`` in addition to the resulting
-            crystal map. Default is ``False``.
-        return_band_data
-            Whether to return the band data array returned from
-            ``EBSDIndexer.index_pats()``. Default is ``False``.
-
-        Returns
-        -------
-        xmap
-            Crystal map with indexing results.
-        index_data
-            Array returned from ``EBSDIndexer.index_pats()``, returned
-            if ``return_index_data=True``.
-        band_data
-            Array returned from ``EBSDIndexer.index_pats()``, returned
-            if ``return_band_data=True``.
-
-        Notes
-        -----
-        Requires :mod:`pyebsdindex` to be installed. See
-        :ref:`optional-dependencies` for further details.
+        window
+            Name of averaging window or an array. Available types are
+            listed in :func:`scipy.signal.windows.get_window`, in
+            addition to a ``"circular"`` window (default) filled with
+            ones in which corner coefficients are set to zero. A window
+            element is considered to be in a corner if its radial
+            distance to the origin (window center) is shorter or equal
+            to the half width of the window's longest axis. A 1D or 2D
+            :class:`~numpy.ndarray`, :class:`~dask.array.Array` or
+            :class:`~kikuchipy.filters.Window` can also be passed.
+        window_shape
+            Shape of averaging window. Not used if a custom window or
+            :class:`~kikuchipy.filters.Window` is passed to ``window``.
+            This can be either 1D or 2D, and can be asymmetrical.
+            Default is ``(3, 3)``.
+        show_progressbar
+            Whether to show a progressbar. If not given, the value of
+            :obj:`hyperspy.api.preferences.General.show_progressbar`
+            is used.
+        **kwargs
+            Keyword arguments passed to the available window type listed
+            in :func:`~scipy.signal.windows.get_window`. If not given,
+            the default values of that particular window are used.
 
-        This wrapper of PyEBSDIndex is meant for convenience more than
-        speed. It uses the GPU if :mod:`pyopencl` is installed, but only
-        uses a single thread. If you need the fastest indexing, refer to
-        the PyEBSDIndex documentation for multi-threading and more.
+        See Also
+        --------
+        kikuchipy.filters.Window, scipy.signal.windows.get_window,
+        scipy.ndimage.correlate
         """
-        if not _pyebsdindex_installed:  # pragma: no cover
-            raise ValueError(
-                "Hough indexing requires pyebsdindex to be installed. Install it with "
-                "pip install pyebsdindex. See "
-                "https://kikuchipy.org/en/stable/user/installation.html for details."
-            )
-        if self._lazy:
-            from pyebsdindex import _pyopencl_installed
+        if isinstance(window, Window) and window.is_valid:
+            averaging_window = copy.copy(window)
+        else:
+            averaging_window = Window(window=window, shape=window_shape, **kwargs)
 
-            if not _pyopencl_installed:  # pragma: no cover
-                raise ValueError(
-                    "Hough indexing of lazy signals must use the GPU, which requires "
-                    "pyopencl to be installed. See "
-                    "https://documen.tician.de/pyopencl/misc.html for installation "
-                    "instructions."
-                )
+        nav_shape = self.axes_manager.navigation_shape[::-1]
+        window_shape = averaging_window.shape
+        if window_shape in [(1,), (1, 1)]:
+            # Do nothing if a window of shape (1,) or (1, 1) is passed
+            return warnings.warn(
+                f"A window of shape {window_shape} was passed, no averaging is "
+                "therefore performed."
+            )
+        elif len(nav_shape) > len(window_shape):
+            averaging_window = averaging_window.reshape(window_shape + (1,))
 
-        am = self.axes_manager
-        nav_shape = am.navigation_shape[::-1]
-        nav_size = int(np.prod(nav_shape))
-        sig_shape = am.signal_shape[::-1]
-        step_sizes = tuple([a.scale for a in am.navigation_axes[::-1]])
+        # Get sum of window data for each pattern, to normalize with
+        # after correlation
+        window_sums = correlate(
+            input=np.ones(nav_shape, dtype=int),
+            weights=averaging_window,
+            mode="constant",
+        )
 
-        # Check indexer
-        phase_list_pei = _get_pyebsdindex_phaselist(phase_list)
-        _ = _indexer_is_compatible_with_kikuchipy(
-            indexer, sig_shape, nav_size, raise_if_not=True
+        # Add signal dimensions to window array to enable its use with
+        # Dask's map_overlap()
+        sig_dim = self.axes_manager.signal_dimension
+        averaging_window = averaging_window.reshape(
+            averaging_window.shape + (1,) * sig_dim
         )
-        if indexer.phaselist != phase_list_pei:
-            raise ValueError(
-                f"`indexer.phaselist` {indexer.phaselist} and the list determined from"
-                f" `phase_list` {phase_list_pei} must be the same"
-            )
 
-        # Prepare patterns
-        chunksize = min(chunksize, max(am.navigation_size, 1))
-        patterns = self.data.reshape((-1,) + sig_shape)
-        if self._lazy:  # pragma: no cover
-            patterns = patterns.rechunk({0: chunksize, 1: -1, 2: -1})
+        # Create dask array of signal patterns and do processing on this
+        if self._lazy:
+            old_chunks = self.data.chunks
+        dask_array = get_dask_array(signal=self, chunk_bytes=8e6, rechunk=True)
 
-        xmap, index_data, band_data = _hough_indexing(
-            patterns=patterns,
-            phase_list=phase_list,
-            nav_shape=nav_shape,
-            step_sizes=step_sizes,
-            indexer=indexer,
-            chunksize=chunksize,
-            verbose=verbose,
+        # Add signal dimensions to array be able to use with Dask's
+        # map_overlap()
+        nav_dim = self.axes_manager.navigation_dimension
+        for i in range(sig_dim):
+            window_sums = np.expand_dims(window_sums, axis=window_sums.ndim)
+        window_sums = da.from_array(
+            window_sums, chunks=dask_array.chunks[:nav_dim] + (1,) * sig_dim
         )
 
-        # Set scan unit
-        if len(nav_shape) > 0:  # Navigation shape can be (1,)
-            scan_unit = str(am.navigation_axes[0].units)
-            if scan_unit != "<undefined>":
-                xmap.scan_unit = scan_unit
+        # Create overlap between chunks to enable correlation with the
+        # window using Dask's map_overlap()
+        window_dim = averaging_window.ndim
+        overlap_depth = {}
+        for i in range(nav_dim):
+            if i < window_dim and dask_array.chunks[i][0] < dask_array.shape[i]:
+                overlap_depth[i] = (window_shape[i] // 2) + 1
+            else:
+                overlap_depth[i] = 1
+        overlap_depth.update(
+            {i: 0 for i in self.axes_manager.signal_indices_in_array[::-1]}
+        )
 
-        if return_index_data and return_band_data:
-            return xmap, index_data, band_data
-        elif return_index_data and not return_band_data:
-            return xmap, index_data
-        elif not return_index_data and return_band_data:
-            return xmap, band_data
+        dtype_out = self.data.dtype
+        omin, omax = dtype_range[dtype_out.type]
+        averaged_patterns = da.overlap.map_overlap(
+            _average_neighbour_patterns,
+            dask_array,
+            window_sums,
+            window=averaging_window,
+            dtype_out=dtype_out,
+            omin=omin,
+            omax=omax,
+            dtype=dtype_out,
+            depth=overlap_depth,
+            boundary="none",
+        )
+
+        # Overwrite signal patterns
+        if not self._lazy:
+            pbar = ProgressBar()
+            if show_progressbar or (
+                show_progressbar is None and hs.preferences.General.show_progressbar
+            ):
+                pbar.register()
+
+            averaged_patterns.store(self.data, compute=True)
+
+            try:
+                pbar.unregister()
+            except KeyError:
+                pass
         else:
-            return xmap
+            # Revert original chunks
+            averaged_patterns = averaged_patterns.rechunk(old_chunks)
 
-    def hough_indexing_optimize_pc(
-        self,
-        pc0: Union[list, tuple, np.ndarray],
-        indexer: "EBSDIndexer",
-        batch: bool = False,
-        method: str = "Nelder-Mead",
-    ) -> "EBSDDetector":
-        """Return a detector with one projection center (PC) per
-        pattern optimized using Hough indexing from :mod:`pyebsdindex`.
+            self.data = averaged_patterns
 
-        See :mod:`~pyebsdindex.pcopt.optimize` for details.
+        # Don't sink
+        gc.collect()
+
+    def downsample(
+        self,
+        factor: int,
+        dtype_out: Optional[str] = None,
+        show_progressbar: Optional[bool] = None,
+    ) -> None:
+        r"""Downsample the pattern shape by an integer factor and
+        rescale intensities to fill the data type range inplace.
 
         Parameters
         ----------
-        pc0
-            A single initial guess of PC for all patterns in Brukers
-            convention, (PCx, PCy, PCz).
-        indexer
-            PyEBSDIndex EBSD indexer instance to pass on to the
-            optimization function. Its `phaselist` must be compatible
-            with ``phase_list``, if given, and the ``indexer.vendor``
-            must be ``"KIKUCHIPY"``. An indexer can be obtained with
-            :meth:`~kikuchipy.detectors.EBSDDetector.get_indexer`.
-        batch
-            Whether the fit for the patterns should be optimized using
-            the cumulative fit for all patterns (``False``, default), or
-            if an optimization is run for each pattern individually.
-        method
-            Which optimization method to use, either ``"Nelder-Mead"``
-            from SciPy (default) or ``"PSO"`` (particle swarm). The
-            latter method does not support ``batch=True``.
+        factor
+            Integer binning factor to downsample by. Must be a divisor
+            of the initial pattern shape :math:`(s_y, s_x)`. The new
+            pattern shape given by the ``factor`` :math:`k` is
+            :math:`(s_y / k, s_x / k)`.
+        dtype_out
+            Name of the data type of the new patterns overwriting
+            :attr:`data`. Contrast between patterns is lost. If not
+            given, patterns maintain their data type and. Patterns are
+            rescaled to fill the data type range.
+        show_progressbar
+            Whether to show a progressbar. If not given, the value of
+            :obj:`hyperspy.api.preferences.General.show_progressbar`
+            is used.
 
-        Returns
-        -------
-        new_detector
-            EBSD detector with one PC if ``batch=False`` or one PC per
-            pattern if ``batch=True``. The detector attributes are
-            extracted from ``indexer.sampleTilt`` etc.
+        See Also
+        --------
+        rebin, crop
 
         Notes
         -----
-        Requires :mod:`pyebsdindex` to be installed. See
-        :ref:`optional-dependencies` for further details.
+        This method differs from :meth:`rebin` in that intensities are
+        rescaled after binning in order to maintain the data type. If
+        rescaling is undesirable, use :meth:`rebin` instead.
         """
-        if not _pyebsdindex_installed:  # pragma: no cover
-            raise ValueError(
-                "Hough indexing requires pyebsdindex to be installed. Install it with "
-                "pip install pyebsdindex. See "
-                "https://kikuchipy.org/en/stable/user/installation.html for details."
-            )
-        if self._lazy:
-            from pyebsdindex import _pyopencl_installed
-
-            if not _pyopencl_installed:  # pragma: no cover
-                raise ValueError(
-                    "Hough indexing of lazy signals must use the GPU, which requires "
-                    "pyopencl to be installed. See "
-                    "https://documen.tician.de/pyopencl/misc.html for installation "
-                    "instructions."
-                )
-
-        pc0 = np.asarray(pc0)
-        if pc0.size != 3:
-            raise ValueError("`pc0` must be of size 3")
-        pc0 = list(pc0)
+        if not isinstance(factor, int) or factor <= 1:
+            raise ValueError(f"Binning `factor` {factor} must be an integer > 1.")
+        else:
+            factor = np.int64(factor)
 
-        supported_methods = ["nelder-mead", "pso"]
-        method = method.lower()
-        if method not in supported_methods:
+        sig_shape_old = self.axes_manager.signal_shape
+        rest = np.mod(sig_shape_old, factor)
+        if not all(rest == 0):
             raise ValueError(
-                f"`method` '{method}' must be one of the supported methods "
-                f"{supported_methods}."
+                f"Binning `factor` {factor} must be a divisor of the initial pattern "
+                f"shape {sig_shape_old}, but {tuple(rest)} pixels remain.\n"
+                "You might try to crop away these pixels first using EBSD.crop()."
             )
-        elif batch and method == "pso":
-            raise ValueError("PSO optimization method does not support `batch=True`.")
-
-        am = self.axes_manager
-        nav_shape = am.navigation_shape[::-1]
-        nav_size = int(np.prod(nav_shape))
-        sig_shape = am.signal_shape[::-1]
+        sig_shape_new = tuple(np.array(sig_shape_old) // factor)
 
-        # Check indexer
-        _ = _indexer_is_compatible_with_kikuchipy(
-            indexer, sig_shape, nav_size, check_pc=False, raise_if_not=True
-        )
+        if dtype_out is not None:
+            dtype_out = np.dtype(dtype_out).type
+        else:
+            dtype_out = self.data.dtype.type
+        omin, omax = dtype_range[dtype_out]
 
-        # Prepare patterns
-        patterns = self.data.reshape((-1,) + sig_shape)
-        if self._lazy:  # pragma: no cover
-            patterns = patterns.rechunk({0: "auto", 1: -1, 2: -1})
+        attrs = self._get_custom_attributes()
 
-        pc = _optimize_pc(
-            pc0=pc0, patterns=patterns, indexer=indexer, batch=batch, method=method
-        )
+        # Update static background
+        static_bg = attrs["static_background"]
+        if static_bg is not None:
+            static_bg_new = _downsample2d(static_bg, factor, omin, omax, dtype_out)
+            attrs["static_background"] = static_bg_new
 
-        if batch:
-            pc = pc.reshape(nav_shape + (3,))
+        # Update detector shape and binning factor
+        attrs["detector"].shape = sig_shape_new
+        attrs["detector"].binning *= factor
 
-        new_detector = EBSDDetector(
-            shape=sig_shape,
-            pc=pc,
-            sample_tilt=indexer.sampleTilt,
-            tilt=indexer.camElev,
+        self.map(
+            _downsample2d,
+            show_progressbar=show_progressbar,
+            parallel=True,
+            output_dtype=dtype_out,
+            factor=factor,
+            omin=omin,
+            omax=omax,
+            dtype_out=dtype_out,
         )
+        self._set_custom_attributes(attrs)
 
-        return new_detector
-
-    def dictionary_indexing(
+    def get_neighbour_dot_product_matrices(
         self,
-        dictionary: EBSD,
-        metric: Union[SimilarityMetric, str] = "ncc",
-        keep_n: int = 20,
-        n_per_iteration: Optional[int] = None,
-        signal_mask: Optional[np.ndarray] = None,
-        rechunk: bool = False,
-        dtype: Union[str, np.dtype, type, None] = None,
-    ) -> CrystalMap:
-        """Index patterns by matching each pattern to a dictionary of
-        simulated patterns of known orientations
-        :cite:`chen2015dictionary,jackson2019dictionary`.
+        window: Optional[Window] = None,
+        zero_mean: bool = True,
+        normalize: bool = True,
+        dtype_out: Union[str, np.dtype, type] = "float32",
+        show_progressbar: Optional[bool] = None,
+    ) -> Union[np.ndarray, da.Array]:
+        """Get an array with dot products of a pattern and its
+        neighbours within a window.
 
         Parameters
         ----------
-        dictionary
-            EBSD signal with dictionary patterns. The signal must have a
-            1D navigation axis, an :attr:`xmap` property with crystal
-            orientations set, and equal detector shape.
-        metric
-            Similarity metric, by default ``"ncc"`` (normalized
-            cross-correlation). ``"ndp"`` (normalized dot product) is
-            also available. A valid user-defined similarity metric
-            may be used instead. The metric must be a class implementing
-            the :class:`~kikuchipy.indexing.SimilarityMetric` abstract
-            class methods. See
-            :class:`~kikuchipy.indexing.NormalizedCrossCorrelationMetric`
-            and :class:`~kikuchipy.indexing.NormalizedDotProductMetric`
-            for examples.
-        keep_n
-            Number of best matches to keep, by default 20 or the number
-            of dictionary patterns if fewer than 20 are available.
-        n_per_iteration
-            Number of dictionary patterns to compare to all experimental
-            patterns in each indexing iteration. If not given, and the
-            dictionary is a ``LazyEBSD`` signal, it is equal to the
-            chunk size of the first pattern array axis, while if if is
-            an ``EBSD`` signal, it is set equal to the number of
-            dictionary patterns, yielding only one iteration. This
-            parameter can be increased to use less memory during
-            indexing, but this will increase the computation time.
-
-            .. versionadded:: 0.5
-        signal_mask
-            A boolean mask equal to the experimental patterns' detector
-            shape ``(n rows, n columns)``, where only pixels equal to
-            ``False`` are matched. If not given, all pixels are used.
-
-            .. versionadded:: 0.5
-        rechunk
-            Whether ``metric`` is allowed to rechunk experimental and
-            dictionary patterns before matching. Default is ``False``.
-            Rechunking usually makes indexing faster, but uses more
-            memory. If a custom ``metric`` is passed, whatever
-            :attr:`~kikuchipy.indexing.SimilarityMetric.rechunk` is set
-            to will be used.
-
-            .. versionadded:: 0.5
-        dtype
-            Which data type ``metric`` shall cast the patterns to before
-            matching. If not given, ``"float32"`` will be used unless a
-            custom ``metric`` is passed and it has set the
-            :attr:`~kikuchipy.indexing.SimilarityMetric.dtype`, which
-            will then be used instead. ``"float32"`` and ``"float64"``
-            are allowed for the available ``"ncc"`` and ``"ndp"``
-            metrics.
-
-            .. versionadded:: 0.5
+        window
+            Window with integer coefficients defining the neighbours to
+            calculate the dot products with. If not given, the four
+            nearest neighbours are used. Must have the same number of
+            dimensions as signal navigation dimensions.
+        zero_mean
+            Whether to subtract the mean of each pattern individually to
+            center the intensities about zero before calculating the
+            dot products. Default is ``True``.
+        normalize
+            Whether to normalize the pattern intensities to a standard
+            deviation of 1 before calculating the dot products. This
+            operation is performed after centering the intensities if
+            ``zero_mean=True``. Default is ``True``.
+        dtype_out
+            Data type of the output map. Default is ``"float32"``.
+        show_progressbar
+            Whether to show a progressbar. If not given, the value of
+            :obj:`hyperspy.api.preferences.General.show_progressbar`
+            is used.
 
         Returns
         -------
-        xmap
-            A crystal map with ``keep_n`` rotations per point with the
-            sorted best matching orientations in the dictionary. The
-            corresponding best scores and indices into the dictionary
-            are stored in the ``xmap.prop`` dictionary as ``"scores"``
-            and ``"simulation_indices"``.
+        dp_matrices
+            Dot products between a pattern and its nearest neighbours.
+        """
+        if self.axes_manager.navigation_dimension == 0:
+            raise ValueError("Signal must have at least one navigation dimension")
 
-        See Also
-        --------
-        kikuchipy.indexing.SimilarityMetric
-        kikuchipy.indexing.NormalizedCrossCorrelationMetric
-        kikuchipy.indexing.NormalizedDotProductMetric
-        kikuchipy.indexing.merge_crystal_maps :
-            Merge multiple single phase crystal maps into one multi
-            phase map.
-        kikuchipy.indexing.orientation_similarity_map :
-            Calculate an orientation similarity map.
+        # Create dask array of signal patterns and do processing on this
+        dask_array = get_dask_array(signal=self)
 
-        Notes
-        -----
-        .. versionchanged:: 0.5
-           Only one dictionary can be passed and the
-           ``return_merged_crystal_map`` and
-           ``get_orientation_similarity_map`` parameters were removed.
-        """
-        exp_am = self.axes_manager
-        dict_am = dictionary.axes_manager
-        dict_size = dict_am.navigation_size
+        # Default to the nearest neighbours
+        nav_dim = self.axes_manager.navigation_dimension
+        if window is None:
+            window = Window(window="circular", shape=(3, 3)[:nav_dim])
 
-        if n_per_iteration is None:
-            if isinstance(dictionary.data, da.Array):
-                n_per_iteration = dictionary.data.chunksize[0]
-            else:
-                n_per_iteration = dict_size
+        # Set overlap depth between navigation chunks equal to the max.
+        # number of nearest neighbours in each navigation axis
+        overlap_depth = _get_chunk_overlap_depth(
+            window=window,
+            axes_manager=self.axes_manager,
+            chunksize=dask_array.chunksize,
+        )
 
-        exp_sig_shape = exp_am.signal_shape[::-1]
-        dict_sig_shape = dict_am.signal_shape[::-1]
-        if exp_sig_shape != dict_sig_shape:
-            raise ValueError(
-                f"Experimental {exp_sig_shape} and dictionary {dict_sig_shape} signal "
-                "shapes must be identical"
-            )
+        dtype_out = np.dtype(dtype_out)
 
-        dict_xmap = dictionary.xmap
-        if dict_xmap is None or dict_xmap.shape != (dict_size,):
-            raise ValueError(
-                "Dictionary signal must have a non-empty `EBSD.xmap` property of equal "
-                "size as the number of dictionary patterns, and both the signal and"
-                "crystal map must have only one navigation dimension"
-            )
+        dp_matrices = dask_array.map_overlap(
+            _get_neighbour_dot_product_matrices,
+            window=window,
+            sig_dim=self.axes_manager.signal_dimension,
+            sig_size=self.axes_manager.signal_size,
+            zero_mean=zero_mean,
+            normalize=normalize,
+            dtype_out=dtype_out,
+            drop_axis=self.axes_manager.signal_indices_in_array[::-1],
+            new_axis=tuple(np.arange(window.ndim) + nav_dim),
+            dtype=dtype_out,
+            depth=overlap_depth,
+            boundary="none",
+        )
 
-        metric = self._prepare_metric(metric, signal_mask, dtype, rechunk, dict_size)
+        if not self._lazy:
+            pbar = ProgressBar()
+            if show_progressbar or (
+                show_progressbar is None and hs.preferences.General.show_progressbar
+            ):
+                pbar.register()
 
-        with dask.config.set(**{"array.slicing.split_large_chunks": False}):
-            return _dictionary_indexing(
-                experimental=self.data,
-                experimental_nav_shape=exp_am.navigation_shape[::-1],
-                dictionary=dictionary.data,
-                step_sizes=tuple(a.scale for a in exp_am.navigation_axes[::-1]),
-                dictionary_xmap=dictionary.xmap,
-                keep_n=keep_n,
-                n_per_iteration=n_per_iteration,
-                metric=metric,
-            )
+            dp_matrices = dp_matrices.compute()
 
-    def refine_orientation(
-        self,
-        xmap: CrystalMap,
-        detector: EBSDDetector,
-        master_pattern: "EBSDMasterPattern",
-        energy: Union[int, float],
-        signal_mask: Optional[np.ndarray] = None,
-        method: Optional[str] = "minimize",
-        method_kwargs: Optional[dict] = None,
-        trust_region: Union[tuple, list, np.ndarray, None] = None,
-        initial_step: Union[float] = None,
-        rtol: float = 1e-4,
-        maxeval: Optional[int] = None,
-        compute: bool = True,
-        rechunk: bool = True,
-        chunk_kwargs: Optional[dict] = None,
-    ) -> Union[CrystalMap, da.Array]:
-        r"""Refine orientations by searching orientation space around
-        the best indexed solution using fixed projection centers.
+            try:
+                pbar.unregister()
+            except KeyError:
+                pass
 
-        Refinement attempts to maximize the similarity between patterns
-        in this signal and simulated patterns projected from a master
-        pattern. The similarity metric used is the normalized
-        cross-correlation (NCC). The orientation, represented by a
-        Euler angle triplet (:math:`\phi_1`, :math:`\Phi`,
-        :math:`\phi_2`) relative to the EDAX TSL sample reference frame
-        RD-TD-ND, is optimized during refinement, while the
-        sample-detector geometry, represented by the three projection
-        center (PC) parameters (PCx, PCy, PCz) in the Bruker convention,
-        is fixed.
-        A subset of the optimization methods in *SciPy* and *NLopt* are
-        available:
+        return dp_matrices
 
-        - Local optimization:
-            - :func:`~scipy.optimize.minimize` (includes Nelder-Mead,
-              Powell etc.).
-            - Nelder-Mead via `nlopt.LN_NELDERMEAD
-              <https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#nelder-mead-simplex>`_
-        - Global optimization:
-            - :func:`~scipy.optimize.differential_evolution`
-            - :func:`~scipy.optimize.dual_annealing`
-            - :func:`~scipy.optimize.basinhopping`
-            - :func:`~scipy.optimize.shgo`
+    def get_image_quality(
+        self,
+        normalize: bool = True,
+        show_progressbar: Optional[bool] = None,
+    ) -> Union[np.ndarray, da.Array]:
+        """Compute the image quality map of patterns in an EBSD scan.
+
+        The image quality :math:`Q` is calculated based on the procedure
+        defined by Krieger Lassen :cite:`lassen1994automated`.
 
         Parameters
         ----------
-        xmap
-            Single phase crystal map with at least one orientation per
-            point.
-        detector
-            Detector describing the detector-sample geometry with either
-            one PC to be used for all map points or one for each point.
-        master_pattern
-            Master pattern in the square Lambert projection of the same
-            phase as the one in the crystal map.
-        energy
-            Accelerating voltage of the electron beam in kV specifying
-            which master pattern energy to use during projection of
-            simulated patterns.
-        signal_mask
-            A boolean mask equal to the experimental patterns' detector
-            shape ``(n rows, n columns)``, where only pixels equal to
-            ``False`` are matched. If not given, all pixels are used.
-        method
-            Name of the :mod:`scipy.optimize` or *NLopt* optimization
-            method, among ``"minimize"``, ``"differential_evolution"``,
-            ``"dual_annealing"``, ``"basinhopping"``, ``"shgo"`` and
-            ``"ln_neldermead"`` (from *NLopt*). Default is
-            ``"minimize"``, which by default performs local optimization
-            with the Nelder-Mead method, unless another ``"minimize"``
-            method is passed to ``method_kwargs``.
-        method_kwargs
-            Keyword arguments passed to the :mod:`scipy.optimize`
-            ``method``. For example, to perform refinement with the
-            modified Powell algorithm from *SciPy*, pass
-            ``method="minimize"`` and
-            ``method_kwargs=dict(method="Powell")``. Not used if
-            ``method="LN_NELDERMEAD"``.
-        trust_region
-            List of three +/- angular deviations in degrees used to
-            determine the bound constraints on the three Euler angles
-            per navigation point, e.g. ``[2, 2, 2]``. Not passed to
-            *SciPy* ``method`` if it does not support bounds. The
-            definition ranges of the Euler angles are
-            :math:`\phi_1 \in [0, 360]`, :math:`\Phi \in [0, 180]` and
-            :math:`\phi_2 \in [0, 360]` in radians.
-        initial_step
-            A single initial step size for all Euler angle, in degrees.
-            Only used if ``method="LN_NELDERMEAD"``. If not given, this
-            is not set for the *NLopt* optimizer.
-        rtol
-            Stop optimization of a pattern when the difference in NCC
-            score between two iterations is below this value (relative
-            tolerance). Default is ``1e-4``. Only used if
-            ``method="LN_NELDERMEAD"``.
-        maxeval
-            Stop optimization of a pattern when the number of function
-            evaluations exceeds this value, e.g. ``100``. Only used if
-            ``method="LN_NELDERMEAD"``.
-        compute
-            Whether to refine now (``True``) or later (``False``).
-            Default is ``True``. See :meth:`~dask.array.Array.compute`
-            for more details.
-        rechunk
-            If ``True`` (default), rechunk the dask array with patterns
-            used in refinement (not the signal data inplace) if it is
-            returned from :func:`~kikuchipy.signals.util.get_dask_array`
-            in a single chunk. This ensures small data sets are
-            rechunked so as to utilize multiple CPUs.
-        chunk_kwargs
-            Keyword arguments passed to
-            :func:`~kikuchipy.signals.util.get_chunking` if
-            ``rechunk=True`` and the dask array with patterns used in
-            refinement is returned from
-            :func:`~kikuchipy.signals.util.get_dask_array` in a single
-            chunk.
+        normalize
+            Whether to normalize patterns to a mean of zero and standard
+            deviation of 1 before calculating the image quality. Default
+            is ``True``.
+        show_progressbar
+            Whether to show a progressbar. If not given, the value of
+            :obj:`hyperspy.api.preferences.General.show_progressbar`
+            is used.
 
         Returns
         -------
-        out
-            Crystal map with refined orientations and similarity metrics
-            in a ``"scores"`` property if ``compute=True``. If
-            ``compute=False``, a dask array of navigation shape + (4,)
-            is returned, to be computed later. See
-            :func:`~kikuchipy.indexing.compute_refine_orientation_results`.
-            Each navigation point has the optimized score and the three
-            Euler angles in radians in element 0, 1, 2, and 3,
-            respectively.
+        image_quality_map
+            Image quality map of same shape as navigation axes. This is
+            a Dask array if the signal is lazy.
 
         See Also
         --------
-        scipy.optimize, refine_projection_center,
-        refine_orientation_projection_center
+        kikuchipy.pattern.get_image_quality
 
-        Notes
-        -----
-        *NLopt* is for now an optional dependency, see
-        :ref:`optional-dependencies` for details. Be aware that *NLopt*
-        does not fail gracefully. If continued use of *NLopt* proves
-        stable enough, its implementation of the Nelder-Mead algorithm
-        might become the default.
+        Examples
+        --------
+        Load an example dataset, remove the static and dynamic
+        background and compute :math:`Q`
+
+        >>> import kikuchipy as kp
+        >>> s = kp.data.nickel_ebsd_small()
+        >>> s
+        <EBSD, title: patterns Scan 1, dimensions: (3, 3|60, 60)>
+        >>> s.remove_static_background()
+        >>> s.remove_dynamic_background()
+        >>> iq = s.get_image_quality()
+        >>> iq
+        array([[0.19935645, 0.16657268, 0.18803978],
+               [0.19040637, 0.1616931 , 0.17834103],
+               [0.19411428, 0.16031407, 0.18413563]], dtype=float32)
         """
-        self._check_refinement_parameters(
-            xmap=xmap,
-            detector=detector,
-            master_pattern=master_pattern,
-            signal_mask=signal_mask,
-        )
-        patterns, signal_mask = self._prepare_patterns_for_refinement(
-            signal_mask=signal_mask, rechunk=rechunk, chunk_kwargs=chunk_kwargs
-        )
-        return _refine_orientation(
-            xmap=xmap,
-            detector=detector,
-            master_pattern=master_pattern,
-            energy=energy,
-            patterns=patterns,
-            signal_mask=signal_mask,
-            method=method,
-            method_kwargs=method_kwargs,
-            trust_region=trust_region,
-            initial_step=initial_step,
-            rtol=rtol,
-            maxeval=maxeval,
-            compute=compute,
+        # Calculate frequency vectors
+        sx, sy = self.axes_manager.signal_shape
+        frequency_vectors = fft_frequency_vectors((sy, sx))
+        inertia_max = np.sum(frequency_vectors) / (sy * sx)
+
+        image_quality_map = self.map(
+            _get_image_quality,
+            show_progressbar=show_progressbar,
+            parallel=True,
+            inplace=False,
+            output_dtype=np.float32,
+            normalize=normalize,
+            frequency_vectors=frequency_vectors,
+            inertia_max=inertia_max,
         )
 
-    def refine_projection_center(
+        return image_quality_map.data
+
+    def get_average_neighbour_dot_product_map(
         self,
-        xmap: CrystalMap,
-        detector: EBSDDetector,
-        master_pattern: "EBSDMasterPattern",
-        energy: Union[int, float],
-        signal_mask: Optional[np.ndarray] = None,
-        method: Optional[str] = "minimize",
-        method_kwargs: Optional[dict] = None,
-        trust_region: Union[tuple, list, np.ndarray, None] = None,
-        initial_step: Union[float] = None,
-        rtol: float = 1e-4,
-        maxeval: Optional[int] = None,
-        compute: bool = True,
-        rechunk: bool = True,
-        chunk_kwargs: Optional[dict] = None,
-    ) -> Union[Tuple[np.ndarray, EBSDDetector], da.Array]:
-        """Refine projection centers by searching the parameter space
-        using fixed orientations.
+        window: Optional[Window] = None,
+        zero_mean: bool = True,
+        normalize: bool = True,
+        dtype_out: Union[str, np.dtype, type] = "float32",
+        dp_matrices: Optional[np.ndarray] = None,
+        show_progressbar: Optional[bool] = None,
+    ) -> Union[np.ndarray, da.Array]:
+        """Get a map of the average dot product between patterns and
+        their neighbours within an averaging window.
 
-        Refinement attempts to maximize the similarity between patterns
-        in this signal and simulated patterns projected from a master
-        pattern. The similarity metric used is the normalized
-        cross-correlation (NCC). The sample-detector geometry,
-        represented by the three projection center (PC) parameters
-        (PCx, PCy, PCz) in the Bruker convention, is updated during
-        refinement, while the orientations, defined relative to the EDAX
-        TSL sample reference frame RD-TD-ND, are fixed.
+        Parameters
+        ----------
+        window
+            Window with integer coefficients defining the neighbours to
+            calculate the average with. If not given, the four nearest
+            neighbours are used. Must have the same number of dimensions
+            as signal navigation dimensions.
+        zero_mean
+            Whether to subtract the mean of each pattern individually to
+            center the intensities about zero before calculating the
+            dot products. Default is ``True``.
+        normalize
+            Whether to normalize the pattern intensities to a standard
+            deviation of 1 before calculating the dot products. This
+            operation is performed after centering the intensities if
+            ``zero_mean=True``. Default is ``True``.
+        dtype_out
+            Data type of the output map. Default is ``"float32"``.
+        dp_matrices
+            Optional pre-calculated dot product matrices, by default
+            ``None``. If an array is passed, the average dot product map
+            is calculated from this array. The ``dp_matrices`` array can
+            be obtained from :meth:`get_neighbour_dot_product_matrices`.
+            Its shape must correspond to the signal's navigation shape
+            and the window's shape.
+        show_progressbar
+            Whether to show a progressbar. If not given, the value of
+            :obj:`hyperspy.api.preferences.General.show_progressbar`
+            is used.
 
-        A subset of the optimization methods in *SciPy* and *NLopt* are
-        available:
+        Returns
+        -------
+        adp
+            Average dot product map.
+        """
+        if self.axes_manager.navigation_dimension == 0:
+            raise ValueError("Signal must have at least one navigation dimension")
 
-        - Local optimization:
-            - :func:`~scipy.optimize.minimize` (includes Nelder-Mead,
-              Powell etc.).
-            - Nelder-Mead via `nlopt.LN_NELDERMEAD
-              <https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#nelder-mead-simplex>`_
-        - Global optimization:
-            - :func:`~scipy.optimize.differential_evolution`
-            - :func:`~scipy.optimize.dual_annealing`
-            - :func:`~scipy.optimize.basinhopping`
-            - :func:`~scipy.optimize.shgo`
+        # Default to the nearest neighbours
+        nav_dim = self.axes_manager.navigation_dimension
+        if window is None:
+            window = Window(window="circular", shape=(3, 3)[:nav_dim])
 
-        Parameters
-        ----------
-        xmap
-            Single phase crystal map with at least one orientation per
-            point.
-        detector
-            Detector describing the detector-sample geometry with either
-            one PC to be used for all map points or one for each point.
-        master_pattern
-            Master pattern in the square Lambert projection of the same
-            phase as the one in the crystal map.
-        energy
-            Accelerating voltage of the electron beam in kV specifying
-            which master pattern energy to use during projection of
-            simulated patterns.
-        signal_mask
-            A boolean mask equal to the experimental patterns' detector
-            shape ``(n rows, n columns)``, where only pixels equal to
-            ``False`` are matched. If not given, all pixels are used.
-        method
-            Name of the :mod:`scipy.optimize` or *NLopt* optimization
-            method, among ``"minimize"``, ``"differential_evolution"``,
-            ``"dual_annealing"``, ``"basinhopping"``, ``"shgo"`` and
-            ``"ln_neldermead"`` (from *NLopt*). Default is
-            ``"minimize"``, which by default performs local optimization
-            with the Nelder-Mead method, unless another ``"minimize"``
-            method is passed to ``method_kwargs``.
-        method_kwargs
-            Keyword arguments passed to the :mod:`scipy.optimize`
-            ``method``. For example, to perform refinement with the
-            modified Powell algorithm from *SciPy*, pass
-            ``method="minimize"`` and
-            ``method_kwargs=dict(method="Powell")``. Not used if
-            ``method="LN_NELDERMEAD"``.
-        trust_region
-            List of three +/- deviations in the range [0, 1] used to
-            determine the bounds constraints on the PC parameters per
-            navigation point, e.g. ``[0.05, 0.05, 0.05]``. Not passed to
-            *SciPy* ``method`` if it does not support bounds. The
-            definition range of the PC parameters are assumed to be
-            [-2, 2].
-        initial_step
-            A single initial step size for all PC parameters in the
-            range [0, 1]. Only used if ``method="LN_NELDERMEAD"``.
-        rtol
-            Stop optimization of a pattern when the difference in NCC
-            score between two iterations is below this value (relative
-            tolerance). Default is ``1e-4``. Only used if
-            ``method="LN_NELDERMEAD"``.
-        maxeval
-            Stop optimization of a pattern when the number of function
-            evaluations exceeds this value, e.g. ``100``. Only used if
-            ``method="LN_NELDERMEAD"``.
-        compute
-            Whether to refine now (``True``) or later (``False``).
-            Default is ``True``. See :meth:`~dask.array.Array.compute`
-            for more details.
-        rechunk
-            If ``True`` (default), rechunk the dask array with patterns
-            used in refinement (not the signal data inplace) if it is
-            returned from :func:`~kikuchipy.signals.util.get_dask_array`
-            in a single chunk. This ensures small data sets are
-            rechunked so as to utilize multiple CPUs.
-        chunk_kwargs
-            Keyword arguments passed to
-            :func:`~kikuchipy.signals.util.get_chunking` if
-            ``rechunk=True`` and the dask array with patterns used in
-            refinement is returned from
-            :func:`~kikuchipy.signals.util.get_dask_array` in a single
-            chunk.
-
-        Returns
-        -------
-        out
-            New similarity metrics and a new EBSD detector instance with
-            the refined PCs if ``compute=True``. If ``compute=False``,
-            a dask array of navigation shape + (4,) is returned, to be
-            computed later. See
-            :func:`~kikuchipy.indexing.compute_refine_projection_center_results`.
-            Each navigation point has the optimized score and the three
-            PC parameters in the Bruker convention in element 0, 1, 2,
-            and 3, respectively.
+        if dp_matrices is not None:
+            nan_slices = [slice(None) for _ in range(nav_dim)]
+            nan_slices += [slice(i, i + 1) for i in window.origin]
+            dp_matrices2 = dp_matrices.copy()
+            dp_matrices2[tuple(nan_slices)] = np.nan
+            if nav_dim == 1:
+                mean_axis = 1
+            else:  # == 2
+                mean_axis = (2, 3)
+            return np.nanmean(dp_matrices2, axis=mean_axis)
 
-        See Also
-        --------
-        scipy.optimize, refine_orientation,
-        refine_orientation_projection_center
+        # Create dask array of signal data and do processing on this
+        dask_array = get_dask_array(signal=self)
 
-        Notes
-        -----
-        *NLopt* is for now an optional dependency, see
-        :ref:`optional-dependencies` for details. Be aware that *NLopt*
-        does not fail gracefully. If continued use of *NLopt* proves
-        stable enough, its implementation of the Nelder-Mead algorithm
-        might become the default.
-        """
-        self._check_refinement_parameters(
-            xmap=xmap,
-            detector=detector,
-            master_pattern=master_pattern,
-            signal_mask=signal_mask,
+        # Set overlap depth between navigation chunks equal to the max.
+        # number of nearest neighbours in each navigation axis
+        overlap_depth = _get_chunk_overlap_depth(
+            window=window,
+            axes_manager=self.axes_manager,
+            chunksize=dask_array.chunksize,
         )
-        patterns, signal_mask = self._prepare_patterns_for_refinement(
-            signal_mask=signal_mask, rechunk=rechunk, chunk_kwargs=chunk_kwargs
+
+        dtype_out = np.dtype(dtype_out)
+
+        adp = dask_array.map_overlap(
+            _get_average_dot_product_map,
+            window=window,
+            sig_dim=self.axes_manager.signal_dimension,
+            sig_size=self.axes_manager.signal_size,
+            zero_mean=zero_mean,
+            normalize=normalize,
+            dtype_out=dtype_out,
+            drop_axis=self.axes_manager.signal_indices_in_array,
+            dtype=dtype_out,
+            depth=overlap_depth,
+            boundary="none",
         )
-        return _refine_pc(
-            xmap=xmap,
-            detector=detector,
-            master_pattern=master_pattern,
-            energy=energy,
-            patterns=patterns,
-            signal_mask=signal_mask,
-            method=method,
-            method_kwargs=method_kwargs,
-            trust_region=trust_region,
-            initial_step=initial_step,
-            rtol=rtol,
-            maxeval=maxeval,
-            compute=compute,
+        chunks = get_chunking(
+            data_shape=self.axes_manager.navigation_shape[::-1],
+            nav_dim=nav_dim,
+            sig_dim=0,
+            dtype=dtype_out,
         )
+        adp = adp.rechunk(chunks=chunks)
 
-    def refine_orientation_projection_center(
+        if not self._lazy:
+            pbar = ProgressBar()
+            if show_progressbar or (
+                show_progressbar is None and hs.preferences.General.show_progressbar
+            ):
+                pbar.register()
+
+            adp = adp.compute()
+
+            try:
+                pbar.unregister()
+            except KeyError:
+                pass
+
+        return adp
+
+    def plot_virtual_bse_intensity(
         self,
-        xmap: CrystalMap,
-        detector: EBSDDetector,
-        master_pattern: "EBSDMasterPattern",
-        energy: Union[int, float],
-        signal_mask: Optional[np.ndarray] = None,
-        method: Optional[str] = "minimize",
-        method_kwargs: Optional[dict] = None,
-        trust_region: Union[tuple, list, np.ndarray, None] = None,
-        initial_step: Union[tuple, list, np.ndarray, None] = None,
-        rtol: Optional[float] = 1e-4,
-        maxeval: Optional[int] = None,
-        compute: bool = True,
-        rechunk: bool = True,
-        chunk_kwargs: Optional[dict] = None,
-    ) -> Union[Tuple[CrystalMap, EBSDDetector], da.Array]:
-        r"""Refine orientations and projection centers simultaneously by
-        searching the orientation and PC parameter space.
+        roi: BaseInteractiveROI,
+        out_signal_axes: Union[Iterable[int], Iterable[str], None] = None,
+        **kwargs,
+    ) -> None:
+        """Plot an interactive virtual backscatter electron (VBSE)
+        image formed from intensities within a specified and adjustable
+        region of interest (ROI) on the detector.
 
-        Refinement attempts to maximize the similarity between patterns
-        in this signal and simulated patterns projected from a master
-        pattern. The only supported similarity metric is the normalized
-        cross-correlation (NCC). The orientation, represented by a
-        Euler angle triplet (:math:`\phi_1`, :math:`\Phi`,
-        :math:`\phi_2`) relative to the EDAX TSL sample reference frame
-        RD-TD-ND, is optimized during refinement, while the
-        sample-detector geometry, represented by the three projection
-        center (PC) parameters (PCx, PCy, PCz) in the Bruker convention,
-        is fixed.
+        Adapted from
+        :meth:`pyxem.signals.common_diffraction.CommonDiffraction.plot_integrated_intensity`.
 
-        A subset of the optimization methods in *SciPy* and *NLopt* are
-        available:
+        Parameters
+        ----------
+        roi
+            Any interactive ROI detailed in HyperSpy.
+        out_signal_axes
+            Which navigation axes to use as signal axes in the virtual
+            image. If not given, the first two navigation axes are used.
+        **kwargs:
+            Keyword arguments passed to the ``plot()`` method of the
+            virtual image.
 
-        - Local optimization:
-            - :func:`~scipy.optimize.minimize` (includes Nelder-Mead,
-              Powell etc.).
-            - Nelder-Mead via `nlopt.LN_NELDERMEAD
-              <https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#nelder-mead-simplex>`_
-        - Global optimization:
-            - :func:`~scipy.optimize.differential_evolution`
-            - :func:`~scipy.optimize.dual_annealing`
-            - :func:`~scipy.optimize.basinhopping`
-            - :func:`~scipy.optimize.shgo`
+        See Also
+        --------
+        get_virtual_bse_intensity
+
+        Examples
+        --------
+        >>> import hyperspy.api as hs
+        >>> import kikuchipy as kp
+        >>> s = kp.data.nickel_ebsd_small()
+        >>> rect_roi = hs.roi.RectangularROI(
+        ...     left=0, right=5, top=0, bottom=5
+        ... )
+        >>> s.plot_virtual_bse_intensity(rect_roi)
+        """
+        # Plot signal if necessary
+        if self._plot is None or not self._plot.is_active:
+            self.plot()
+
+        # Get the sliced signal from the ROI
+        sliced_signal = roi.interactive(self, axes=self.axes_manager.signal_axes)
+
+        # Create an output signal for the virtual backscatter electron
+        # calculation
+        out = self._get_sum_signal(self, out_signal_axes)
+        out.metadata.General.title = "Virtual backscatter electron intensity"
+
+        # Create the interactive signal
+        hs.interactive(
+            f=sliced_signal.nansum,
+            axis=sliced_signal.axes_manager.signal_axes,
+            event=roi.events.changed,
+            recompute_out_event=None,
+            out=out,
+        )
+
+        # Plot the result
+        out.plot(**kwargs)
+
+    def get_virtual_bse_intensity(
+        self,
+        roi: BaseInteractiveROI,
+        out_signal_axes: Union[Iterable[int], Iterable[str], None] = None,
+    ) -> VirtualBSEImage:
+        """Get a virtual backscatter electron (VBSE) image formed from
+        intensities within a region of interest (ROI) on the detector.
+
+        Adapted from
+        :meth:`pyxem.signals.common_diffraction.CommonDiffraction.get_integrated_intensity`.
 
         Parameters
         ----------
-        xmap
-            Single phase crystal map with at least one orientation per
-            point.
-        detector
-            Detector describing the detector-sample geometry with either
-            one PC to be used for all map points or one for each point.
-        master_pattern
-            Master pattern in the square Lambert projection of the same
-            phase as the one in the crystal map.
-        energy
-            Accelerating voltage of the electron beam in kV specifying
-            which master pattern energy to use during projection of
-            simulated patterns.
-        signal_mask
-            A boolean mask equal to the experimental patterns' detector
-            shape ``(n rows, n columns)``, where only pixels equal to
-            ``False`` are matched. If not given, all pixels are used.
-        method
-            Name of the :mod:`scipy.optimize` or *NLopt* optimization
-            method, among ``"minimize"``, ``"differential_evolution"``,
-            ``"dual_annealing"``, ``"basinhopping"``, ``"shgo"`` and
-            ``"ln_neldermead"`` (from *NLopt*). Default is
-            ``"minimize"``, which by default performs local optimization
-            with the Nelder-Mead method, unless another ``"minimize"``
-            method is passed to ``method_kwargs``.
-        method_kwargs
-            Keyword arguments passed to the :mod:`scipy.optimize`
-            ``method``. For example, to perform refinement with the
-            modified Powell algorithm from *SciPy*, pass
-            ``method="minimize"`` and
-            ``method_kwargs=dict(method="Powell")``. Not used if
-            ``method="LN_NELDERMEAD"``.
-        trust_region
-            List of three +/- angular deviations in degrees as bound
-            constraints on the three Euler angles and three +/-
-            deviations in the range [0, 1] as bound constraints on the
-            PC parameters, e.g. ``[2, 2, 2, 0.05, 0.05, 0.05]``. Not
-            passed to *SciPy* ``method`` if it does not support bounds.
-            The definition ranges of the Euler angles are
-            :math:`\phi_1 \in [0, 360]`, :math:`\Phi \in [0, 180]` and
-            :math:`\phi_2 \in [0, 360]` in radians, while the definition
-            range of the PC parameters are assumed to be [-2, 2].
-        initial_step
-            A list of two initial step sizes to use, one in degrees for
-            all Euler angles and one in the range [0, 1] for all PC
-            parameters. Only used if ``method="LN_NELDERMEAD"``.
-        rtol
-            Stop optimization of a pattern when the difference in NCC
-            score between two iterations is below this value (relative
-            tolerance). Only used if ``method="LN_NELDERMEAD"``. If not
-            given, this is set to ``1e-4``.
-        maxeval
-            Stop optimization of a pattern when the number of function
-            evaluations exceeds this value, e.g. ``100``. Only used if
-            ``method="LN_NELDERMEAD"``.
-        compute
-            Whether to refine now (``True``) or later (``False``).
-            Default is ``True``. See :meth:`~dask.array.Array.compute`
-            for more details.
-        rechunk
-            If ``True`` (default), rechunk the dask array with patterns
-            used in refinement (not the signal data inplace) if it is
-            returned from :func:`~kikuchipy.signals.util.get_dask_array`
-            in a single chunk. This ensures small data sets are
-            rechunked so as to utilize multiple CPUs.
-        chunk_kwargs
-            Keyword arguments passed to
-            :func:`~kikuchipy.signals.util.get_chunking` if
-            ``rechunk=True`` and the dask array with patterns used in
-            refinement is returned from
-            :func:`~kikuchipy.signals.util.get_dask_array` in a single
-            chunk.
+        roi
+            Any interactive ROI detailed in HyperSpy.
+        out_signal_axes
+            Which navigation axes to use as signal axes in the virtual
+            image. If not given, the first two navigation axes are used.
 
         Returns
         -------
-        out
-            Crystal map with refined orientations and a new EBSD
-            detector instance with the refined PCs, if ``compute=True``.
-            If ``compute=False``, a dask array of navigation shape +
-            (7,) is returned, to be computed later. See
-            :func:`~kikuchipy.indexing.compute_refine_orientation_projection_center_results`.
-            Each navigation point has the optimized score, the three
-            Euler angles in radians, and the three PC parameters in
-            element 0, 1, 2, 3, 4, 5, and 6, respectively.
+        virtual_image
+            VBSE image formed from detector intensities within an ROI
+            on the detector.
 
         See Also
         --------
-        scipy.optimize, refine_orientation, refine_projection_center
-
-        Notes
-        -----
-        The method attempts to refine the orientations and projection
-        center at the same time for each map point. The optimization
-        landscape is sloppy :cite:`pang2020global`, where the
-        orientation and PC can make up for each other. Thus, it is
-        possible that the parameters that yield the highest similarity
-        are incorrect. As always, it is left to the user to ensure that
-        the output is reasonable.
+        plot_virtual_bse_intensity
 
-        *NLopt* is for now an optional dependency, see
-        :ref:`optional-dependencies` for details. Be aware that *NLopt*
-        does not fail gracefully. If continued use of *NLopt* proves
-        stable enough, its implementation of the Nelder-Mead algorithm
-        might become the default.
+        Examples
+        --------
+        >>> import hyperspy.api as hs
+        >>> import kikuchipy as kp
+        >>> rect_roi = hs.roi.RectangularROI(
+        ...     left=0, right=5, top=0, bottom=5
+        ... )
+        >>> s = kp.data.nickel_ebsd_small()
+        >>> vbse_image = s.get_virtual_bse_intensity(rect_roi)
         """
-        self._check_refinement_parameters(
-            xmap=xmap,
-            detector=detector,
-            master_pattern=master_pattern,
-            signal_mask=signal_mask,
-        )
-        patterns, signal_mask = self._prepare_patterns_for_refinement(
-            signal_mask=signal_mask, rechunk=rechunk, chunk_kwargs=chunk_kwargs
-        )
-        return _refine_orientation_pc(
-            xmap=xmap,
-            detector=detector,
-            master_pattern=master_pattern,
-            energy=energy,
-            patterns=patterns,
-            signal_mask=signal_mask,
-            method=method,
-            method_kwargs=method_kwargs,
-            trust_region=trust_region,
-            initial_step=initial_step,
-            rtol=rtol,
-            maxeval=maxeval,
-            compute=compute,
-        )
+        vbse = roi(self, axes=self.axes_manager.signal_axes)
+        vbse_sum = self._get_sum_signal(vbse, out_signal_axes)
+        vbse_sum.metadata.General.title = "Virtual backscatter electron image"
+        vbse_sum.set_signal_type("VirtualBSEImage")
+        return vbse_sum
 
-    def fft_filter(
+    def hough_indexing(
         self,
-        transfer_function: Union[np.ndarray, Window],
-        function_domain: str,
-        shift: bool = False,
-        show_progressbar: Optional[bool] = None,
-    ) -> None:
-        """Filter an EBSD scan inplace in the frequency domain.
+        phase_list: PhaseList,
+        indexer: "EBSDIndexer",
+        chunksize: int = 528,
+        verbose: int = 1,
+        return_index_data: bool = False,
+        return_band_data: bool = False,
+    ) -> Union[
+        CrystalMap,
+        Tuple[CrystalMap, np.ndarray],
+        Tuple[CrystalMap, np.ndarray, np.ndarray],
+    ]:
+        """Index patterns by Hough indexing using :mod:`pyebsdindex`.
 
-        Patterns are transformed via the Fast Fourier Transform (FFT) to
-        the frequency domain, where their spectrum is multiplied by the
-        ``transfer_function``, and the filtered spectrum is subsequently
-        transformed to the spatial domain via the inverse FFT (IFFT).
-        Filtered patterns are rescaled to input data type range.
+        See :meth:`~pyebsdindex.ebsd_index.EBSDIndexer` and
+        :meth:`~pyebsdindex.ebsd_index.EBSDIndexer.index_pats` for
+        details.
 
-        Note that if ``function_domain`` is ``"spatial"``, only real
-        valued FFT and IFFT is used.
+        Currently, PyEBSDIndex only supports indexing with a single
+        mean projection center (PC).
 
         Parameters
         ----------
-        transfer_function
-            Filter to apply to patterns. This can either be a transfer
-            function in the frequency domain of pattern shape or a
-            kernel in the spatial domain. What is passed is determined
-            from ``function_domain``.
-        function_domain
-            Options are ``"frequency"`` and ``"spatial"``, indicating,
-            respectively, whether the filter function passed to
-            ``filter_function`` is a transfer function in the frequency
-            domain or a kernel in the spatial domain.
-        shift
-            Whether to shift the zero-frequency component to the center.
-            Default is ``False``. This is only used when
-            ``function_domain="frequency"``.
-        show_progressbar
-            Whether to show a progressbar. If not given, the value of
-            :obj:`hyperspy.api.preferences.General.show_progressbar`
-            is used.
+        phase_list
+            List of phases. The list can only contain one face-centered
+            cubic (FCC) phase, one body-centered cubic (BCC) phase or
+            both types.
+        indexer
+            PyEBSDIndex EBSD indexer instance of which the
+            :meth:`~pyebsdindex.ebsd_index.EBSDIndexer.index_pats`
+            method is called. Its `phaselist` must be compatible with
+            the given ``phase_list``, and the ``indexer.vendor`` must be
+            ``"KIKUCHIPY"``. An indexer can be obtained with
+            :meth:`~kikuchipy.detectors.EBSDDetector.get_indexer`.
+        chunksize
+            Number of patterns to index at a time. Default is the
+            minimum of 528 or the number of patterns in the signal.
+            Increasing the chunksize may give faster indexing but
+            increases memory use.
+        verbose
+            Which information to print from PyEBSDIndex. Options are
+            0 - no output, 1 - timings (default), 2 - timings and the
+            Hough transform of the first pattern with detected bands
+            highlighted.
+        return_index_data
+            Whether to return the index data array returned from
+            ``EBSDIndexer.index_pats()`` in addition to the resulting
+            crystal map. Default is ``False``.
+        return_band_data
+            Whether to return the band data array returned from
+            ``EBSDIndexer.index_pats()``. Default is ``False``.
 
-        See Also
-        --------
-        kikuchipy.filters.Window
+        Returns
+        -------
+        xmap
+            Crystal map with indexing results.
+        index_data
+            Array returned from ``EBSDIndexer.index_pats()``, returned
+            if ``return_index_data=True``.
+        band_data
+            Array returned from ``EBSDIndexer.index_pats()``, returned
+            if ``return_band_data=True``.
 
-        Examples
-        --------
-        Applying a Gaussian low pass filter with a cutoff frequency of
-        20:
+        Notes
+        -----
+        Requires :mod:`pyebsdindex` to be installed. See
+        :ref:`optional-dependencies` for further details.
 
-        >>> import kikuchipy as kp
-        >>> s = kp.data.nickel_ebsd_small()
-        >>> pattern_shape = s.axes_manager.signal_shape[::-1]
-        >>> w = kp.filters.Window(
-        ...     "lowpass", cutoff=20, shape=pattern_shape
-        ... )
-        >>> s.fft_filter(
-        ...     transfer_function=w,
-        ...     function_domain="frequency",
-        ...     shift=True,
-        ... )
+        This wrapper of PyEBSDIndex is meant for convenience more than
+        speed. It uses the GPU if :mod:`pyopencl` is installed, but only
+        uses a single thread. If you need the fastest indexing, refer to
+        the PyEBSDIndex documentation for multi-threading and more.
         """
-        dtype_out = self.data.dtype.type
+        if not _pyebsdindex_installed:  # pragma: no cover
+            raise ValueError(
+                "Hough indexing requires pyebsdindex to be installed. Install it with "
+                "pip install pyebsdindex. See "
+                "https://kikuchipy.org/en/stable/user/installation.html for details."
+            )
+        if self._lazy:
+            from pyebsdindex import _pyopencl_installed
 
-        dtype = np.float32
-        dask_array = get_dask_array(signal=self, dtype=dtype)
+            if not _pyopencl_installed:  # pragma: no cover
+                raise ValueError(
+                    "Hough indexing of lazy signals must use the GPU, which requires "
+                    "pyopencl to be installed. See "
+                    "https://documen.tician.de/pyopencl/misc.html for installation "
+                    "instructions."
+                )
 
-        kwargs = {}
-        if function_domain == "frequency":
-            filter_func = fft_filter
-            kwargs["shift"] = shift
-        elif function_domain == "spatial":
-            filter_func = _fft_filter  # Barnes
-            kwargs["window_shape"] = transfer_function.shape
+        am = self.axes_manager
+        nav_shape = am.navigation_shape[::-1]
+        nav_size = int(np.prod(nav_shape))
+        sig_shape = am.signal_shape[::-1]
+        step_sizes = tuple([a.scale for a in am.navigation_axes[::-1]])
 
-            # FFT filter setup
-            (
-                kwargs["fft_shape"],
-                transfer_function,  # Padded
-                kwargs["offset_before_fft"],
-                kwargs["offset_after_ifft"],
-            ) = _fft_filter_setup(
-                image_shape=self.axes_manager.signal_shape[::-1],
-                window=transfer_function,
+        # Check indexer
+        phase_list_pei = _get_pyebsdindex_phaselist(phase_list)
+        _ = _indexer_is_compatible_with_kikuchipy(
+            indexer, sig_shape, nav_size, raise_if_not=True
+        )
+        if indexer.phaselist != phase_list_pei:
+            raise ValueError(
+                f"`indexer.phaselist` {indexer.phaselist} and the list determined from"
+                f" `phase_list` {phase_list_pei} must be the same"
             )
-        else:
-            function_domains = ["frequency", "spatial"]
-            raise ValueError(f"{function_domain} must be either of {function_domains}.")
 
-        filtered_patterns = dask_array.map_blocks(
-            func=chunk.fft_filter,
-            filter_func=filter_func,
-            transfer_function=transfer_function,
-            dtype_out=dtype_out,
-            dtype=dtype_out,
-            **kwargs,
-        )
+        # Prepare patterns
+        chunksize = min(chunksize, max(am.navigation_size, 1))
+        patterns = self.data.reshape((-1,) + sig_shape)
+        if self._lazy:  # pragma: no cover
+            patterns = patterns.rechunk({0: chunksize, 1: -1, 2: -1})
 
-        # Overwrite signal patterns
-        if not self._lazy:
-            pbar = ProgressBar()
-            if show_progressbar or (
-                show_progressbar is None and hs.preferences.General.show_progressbar
-            ):
-                pbar.register()
+        xmap, index_data, band_data = _hough_indexing(
+            patterns=patterns,
+            phase_list=phase_list,
+            nav_shape=nav_shape,
+            step_sizes=step_sizes,
+            indexer=indexer,
+            chunksize=chunksize,
+            verbose=verbose,
+        )
 
-            filtered_patterns.store(self.data, compute=True)
+        # Set scan unit
+        if len(nav_shape) > 0:  # Navigation shape can be (1,)
+            scan_unit = str(am.navigation_axes[0].units)
+            if scan_unit != "<undefined>":
+                xmap.scan_unit = scan_unit
 
-            try:
-                pbar.unregister()
-            except KeyError:
-                pass
+        if return_index_data and return_band_data:
+            return xmap, index_data, band_data
+        elif return_index_data and not return_band_data:
+            return xmap, index_data
+        elif not return_index_data and return_band_data:
+            return xmap, band_data
         else:
-            self.data = filtered_patterns
+            return xmap
 
-    def get_neighbour_dot_product_matrices(
+    def hough_indexing_optimize_pc(
         self,
-        window: Optional[Window] = None,
-        zero_mean: bool = True,
-        normalize: bool = True,
-        dtype_out: Union[str, np.dtype, type] = "float32",
-        show_progressbar: Optional[bool] = None,
-    ) -> Union[np.ndarray, da.Array]:
-        """Get an array with dot products of a pattern and its
-        neighbours within a window.
+        pc0: Union[list, tuple, np.ndarray],
+        indexer: "EBSDIndexer",
+        batch: bool = False,
+        method: str = "Nelder-Mead",
+    ) -> "EBSDDetector":
+        """Return a detector with one projection center (PC) per
+        pattern optimized using Hough indexing from :mod:`pyebsdindex`.
+
+        See :mod:`~pyebsdindex.pcopt.optimize` for details.
 
         Parameters
         ----------
-        window
-            Window with integer coefficients defining the neighbours to
-            calculate the dot products with. If not given, the four
-            nearest neighbours are used. Must have the same number of
-            dimensions as signal navigation dimensions.
-        zero_mean
-            Whether to subtract the mean of each pattern individually to
-            center the intensities about zero before calculating the
-            dot products. Default is ``True``.
-        normalize
-            Whether to normalize the pattern intensities to a standard
-            deviation of 1 before calculating the dot products. This
-            operation is performed after centering the intensities if
-            ``zero_mean=True``. Default is ``True``.
-        dtype_out
-            Data type of the output map. Default is ``"float32"``.
-        show_progressbar
-            Whether to show a progressbar. If not given, the value of
-            :obj:`hyperspy.api.preferences.General.show_progressbar`
-            is used.
+        pc0
+            A single initial guess of PC for all patterns in Brukers
+            convention, (PCx, PCy, PCz).
+        indexer
+            PyEBSDIndex EBSD indexer instance to pass on to the
+            optimization function. Its `phaselist` must be compatible
+            with ``phase_list``, if given, and the ``indexer.vendor``
+            must be ``"KIKUCHIPY"``. An indexer can be obtained with
+            :meth:`~kikuchipy.detectors.EBSDDetector.get_indexer`.
+        batch
+            Whether the fit for the patterns should be optimized using
+            the cumulative fit for all patterns (``False``, default), or
+            if an optimization is run for each pattern individually.
+        method
+            Which optimization method to use, either ``"Nelder-Mead"``
+            from SciPy (default) or ``"PSO"`` (particle swarm). The
+            latter method does not support ``batch=True``.
 
         Returns
         -------
-        dp_matrices
-            Dot products between a pattern and its nearest neighbours.
+        new_detector
+            EBSD detector with one PC if ``batch=False`` or one PC per
+            pattern if ``batch=True``. The detector attributes are
+            extracted from ``indexer.sampleTilt`` etc.
+
+        Notes
+        -----
+        Requires :mod:`pyebsdindex` to be installed. See
+        :ref:`optional-dependencies` for further details.
         """
-        if self.axes_manager.navigation_dimension == 0:
-            raise ValueError("Signal must have at least one navigation dimension")
+        if not _pyebsdindex_installed:  # pragma: no cover
+            raise ValueError(
+                "Hough indexing requires pyebsdindex to be installed. Install it with "
+                "pip install pyebsdindex. See "
+                "https://kikuchipy.org/en/stable/user/installation.html for details."
+            )
+        if self._lazy:
+            from pyebsdindex import _pyopencl_installed
 
-        # Create dask array of signal patterns and do processing on this
-        dask_array = get_dask_array(signal=self)
+            if not _pyopencl_installed:  # pragma: no cover
+                raise ValueError(
+                    "Hough indexing of lazy signals must use the GPU, which requires "
+                    "pyopencl to be installed. See "
+                    "https://documen.tician.de/pyopencl/misc.html for installation "
+                    "instructions."
+                )
 
-        # Default to the nearest neighbours
-        nav_dim = self.axes_manager.navigation_dimension
-        if window is None:
-            window = Window(window="circular", shape=(3, 3)[:nav_dim])
+        pc0 = np.asarray(pc0)
+        if pc0.size != 3:
+            raise ValueError("`pc0` must be of size 3")
+        pc0 = list(pc0)
 
-        # Set overlap depth between navigation chunks equal to the max.
-        # number of nearest neighbours in each navigation axis
-        overlap_depth = _get_chunk_overlap_depth(
-            window=window,
-            axes_manager=self.axes_manager,
-            chunksize=dask_array.chunksize,
-        )
+        supported_methods = ["nelder-mead", "pso"]
+        method = method.lower()
+        if method not in supported_methods:
+            raise ValueError(
+                f"`method` '{method}' must be one of the supported methods "
+                f"{supported_methods}."
+            )
+        elif batch and method == "pso":
+            raise ValueError("PSO optimization method does not support `batch=True`.")
 
-        dtype_out = np.dtype(dtype_out)
+        am = self.axes_manager
+        nav_shape = am.navigation_shape[::-1]
+        nav_size = int(np.prod(nav_shape))
+        sig_shape = am.signal_shape[::-1]
 
-        dp_matrices = dask_array.map_overlap(
-            _get_neighbour_dot_product_matrices,
-            window=window,
-            sig_dim=self.axes_manager.signal_dimension,
-            sig_size=self.axes_manager.signal_size,
-            zero_mean=zero_mean,
-            normalize=normalize,
-            dtype_out=dtype_out,
-            drop_axis=self.axes_manager.signal_indices_in_array[::-1],
-            new_axis=tuple(np.arange(window.ndim) + nav_dim),
-            dtype=dtype_out,
-            depth=overlap_depth,
-            boundary="none",
+        # Check indexer
+        _ = _indexer_is_compatible_with_kikuchipy(
+            indexer, sig_shape, nav_size, check_pc=False, raise_if_not=True
         )
 
-        if not self._lazy:
-            pbar = ProgressBar()
-            if show_progressbar or (
-                show_progressbar is None and hs.preferences.General.show_progressbar
-            ):
-                pbar.register()
+        # Prepare patterns
+        patterns = self.data.reshape((-1,) + sig_shape)
+        if self._lazy:  # pragma: no cover
+            patterns = patterns.rechunk({0: "auto", 1: -1, 2: -1})
 
-            dp_matrices = dp_matrices.compute()
+        pc = _optimize_pc(
+            pc0=pc0, patterns=patterns, indexer=indexer, batch=batch, method=method
+        )
 
-            try:
-                pbar.unregister()
-            except KeyError:
-                pass
+        if batch:
+            pc = pc.reshape(nav_shape + (3,))
 
-        return dp_matrices
+        new_detector = EBSDDetector(
+            shape=sig_shape,
+            pc=pc,
+            sample_tilt=indexer.sampleTilt,
+            tilt=indexer.camElev,
+        )
 
-    def get_average_neighbour_dot_product_map(
+        return new_detector
+
+    def dictionary_indexing(
         self,
-        window: Optional[Window] = None,
-        zero_mean: bool = True,
-        normalize: bool = True,
-        dtype_out: Union[str, np.dtype, type] = "float32",
-        dp_matrices: Optional[np.ndarray] = None,
-        show_progressbar: Optional[bool] = None,
-    ) -> Union[np.ndarray, da.Array]:
-        """Get a map of the average dot product between patterns and
-        their neighbours within an averaging window.
+        dictionary: EBSD,
+        metric: Union[SimilarityMetric, str] = "ncc",
+        keep_n: int = 20,
+        n_per_iteration: Optional[int] = None,
+        signal_mask: Optional[np.ndarray] = None,
+        rechunk: bool = False,
+        dtype: Union[str, np.dtype, type, None] = None,
+    ) -> CrystalMap:
+        """Index patterns by matching each pattern to a dictionary of
+        simulated patterns of known orientations
+        :cite:`chen2015dictionary,jackson2019dictionary`.
 
         Parameters
         ----------
-        window
-            Window with integer coefficients defining the neighbours to
-            calculate the average with. If not given, the four nearest
-            neighbours are used. Must have the same number of dimensions
-            as signal navigation dimensions.
-        zero_mean
-            Whether to subtract the mean of each pattern individually to
-            center the intensities about zero before calculating the
-            dot products. Default is ``True``.
-        normalize
-            Whether to normalize the pattern intensities to a standard
-            deviation of 1 before calculating the dot products. This
-            operation is performed after centering the intensities if
-            ``zero_mean=True``. Default is ``True``.
-        dtype_out
-            Data type of the output map. Default is ``"float32"``.
-        dp_matrices
-            Optional pre-calculated dot product matrices, by default
-            ``None``. If an array is passed, the average dot product map
-            is calculated from this array. The ``dp_matrices`` array can
-            be obtained from :meth:`get_neighbour_dot_product_matrices`.
-            Its shape must correspond to the signal's navigation shape
-            and the window's shape.
-        show_progressbar
-            Whether to show a progressbar. If not given, the value of
-            :obj:`hyperspy.api.preferences.General.show_progressbar`
-            is used.
+        dictionary
+            EBSD signal with dictionary patterns. The signal must have a
+            1D navigation axis, an :attr:`xmap` property with crystal
+            orientations set, and equal detector shape.
+        metric
+            Similarity metric, by default ``"ncc"`` (normalized
+            cross-correlation). ``"ndp"`` (normalized dot product) is
+            also available. A valid user-defined similarity metric
+            may be used instead. The metric must be a class implementing
+            the :class:`~kikuchipy.indexing.SimilarityMetric` abstract
+            class methods. See
+            :class:`~kikuchipy.indexing.NormalizedCrossCorrelationMetric`
+            and :class:`~kikuchipy.indexing.NormalizedDotProductMetric`
+            for examples.
+        keep_n
+            Number of best matches to keep, by default 20 or the number
+            of dictionary patterns if fewer than 20 are available.
+        n_per_iteration
+            Number of dictionary patterns to compare to all experimental
+            patterns in each indexing iteration. If not given, and the
+            dictionary is a ``LazyEBSD`` signal, it is equal to the
+            chunk size of the first pattern array axis, while if if is
+            an ``EBSD`` signal, it is set equal to the number of
+            dictionary patterns, yielding only one iteration. This
+            parameter can be increased to use less memory during
+            indexing, but this will increase the computation time.
 
-        Returns
-        -------
-        adp
-            Average dot product map.
-        """
-        if self.axes_manager.navigation_dimension == 0:
-            raise ValueError("Signal must have at least one navigation dimension")
+            .. versionadded:: 0.5
+        signal_mask
+            A boolean mask equal to the experimental patterns' detector
+            shape ``(n rows, n columns)``, where only pixels equal to
+            ``False`` are matched. If not given, all pixels are used.
 
-        # Default to the nearest neighbours
-        nav_dim = self.axes_manager.navigation_dimension
-        if window is None:
-            window = Window(window="circular", shape=(3, 3)[:nav_dim])
+            .. versionadded:: 0.5
+        rechunk
+            Whether ``metric`` is allowed to rechunk experimental and
+            dictionary patterns before matching. Default is ``False``.
+            Rechunking usually makes indexing faster, but uses more
+            memory. If a custom ``metric`` is passed, whatever
+            :attr:`~kikuchipy.indexing.SimilarityMetric.rechunk` is set
+            to will be used.
 
-        if dp_matrices is not None:
-            nan_slices = [slice(None) for _ in range(nav_dim)]
-            nan_slices += [slice(i, i + 1) for i in window.origin]
-            dp_matrices2 = dp_matrices.copy()
-            dp_matrices2[tuple(nan_slices)] = np.nan
-            if nav_dim == 1:
-                mean_axis = 1
-            else:  # == 2
-                mean_axis = (2, 3)
-            return np.nanmean(dp_matrices2, axis=mean_axis)
+            .. versionadded:: 0.5
+        dtype
+            Which data type ``metric`` shall cast the patterns to before
+            matching. If not given, ``"float32"`` will be used unless a
+            custom ``metric`` is passed and it has set the
+            :attr:`~kikuchipy.indexing.SimilarityMetric.dtype`, which
+            will then be used instead. ``"float32"`` and ``"float64"``
+            are allowed for the available ``"ncc"`` and ``"ndp"``
+            metrics.
 
-        # Create dask array of signal data and do processing on this
-        dask_array = get_dask_array(signal=self)
+            .. versionadded:: 0.5
 
-        # Set overlap depth between navigation chunks equal to the max.
-        # number of nearest neighbours in each navigation axis
-        overlap_depth = _get_chunk_overlap_depth(
-            window=window,
-            axes_manager=self.axes_manager,
-            chunksize=dask_array.chunksize,
-        )
+        Returns
+        -------
+        xmap
+            A crystal map with ``keep_n`` rotations per point with the
+            sorted best matching orientations in the dictionary. The
+            corresponding best scores and indices into the dictionary
+            are stored in the ``xmap.prop`` dictionary as ``"scores"``
+            and ``"simulation_indices"``.
+
+        See Also
+        --------
+        kikuchipy.indexing.SimilarityMetric
+        kikuchipy.indexing.NormalizedCrossCorrelationMetric
+        kikuchipy.indexing.NormalizedDotProductMetric
+        kikuchipy.indexing.merge_crystal_maps :
+            Merge multiple single phase crystal maps into one multi
+            phase map.
+        kikuchipy.indexing.orientation_similarity_map :
+            Calculate an orientation similarity map.
 
-        dtype_out = np.dtype(dtype_out)
+        Notes
+        -----
+        .. versionchanged:: 0.5
+           Only one dictionary can be passed and the
+           ``return_merged_crystal_map`` and
+           ``get_orientation_similarity_map`` parameters were removed.
+        """
+        exp_am = self.axes_manager
+        dict_am = dictionary.axes_manager
+        dict_size = dict_am.navigation_size
 
-        adp = dask_array.map_overlap(
-            _get_average_dot_product_map,
-            window=window,
-            sig_dim=self.axes_manager.signal_dimension,
-            sig_size=self.axes_manager.signal_size,
-            zero_mean=zero_mean,
-            normalize=normalize,
-            dtype_out=dtype_out,
-            drop_axis=self.axes_manager.signal_indices_in_array,
-            dtype=dtype_out,
-            depth=overlap_depth,
-            boundary="none",
-        )
-        chunks = get_chunking(
-            data_shape=self.axes_manager.navigation_shape[::-1],
-            nav_dim=nav_dim,
-            sig_dim=0,
-            dtype=dtype_out,
-        )
-        adp = adp.rechunk(chunks=chunks)
+        if n_per_iteration is None:
+            if isinstance(dictionary.data, da.Array):
+                n_per_iteration = dictionary.data.chunksize[0]
+            else:
+                n_per_iteration = dict_size
 
-        if not self._lazy:
-            pbar = ProgressBar()
-            if show_progressbar or (
-                show_progressbar is None and hs.preferences.General.show_progressbar
-            ):
-                pbar.register()
+        exp_sig_shape = exp_am.signal_shape[::-1]
+        dict_sig_shape = dict_am.signal_shape[::-1]
+        if exp_sig_shape != dict_sig_shape:
+            raise ValueError(
+                f"Experimental {exp_sig_shape} and dictionary {dict_sig_shape} signal "
+                "shapes must be identical"
+            )
 
-            adp = adp.compute()
+        dict_xmap = dictionary.xmap
+        if dict_xmap is None or dict_xmap.shape != (dict_size,):
+            raise ValueError(
+                "Dictionary signal must have a non-empty `EBSD.xmap` property of equal "
+                "size as the number of dictionary patterns, and both the signal and"
+                "crystal map must have only one navigation dimension"
+            )
 
-            try:
-                pbar.unregister()
-            except KeyError:
-                pass
+        metric = self._prepare_metric(metric, signal_mask, dtype, rechunk, dict_size)
 
-        return adp
+        with dask.config.set(**{"array.slicing.split_large_chunks": False}):
+            return _dictionary_indexing(
+                experimental=self.data,
+                experimental_nav_shape=exp_am.navigation_shape[::-1],
+                dictionary=dictionary.data,
+                step_sizes=tuple(a.scale for a in exp_am.navigation_axes[::-1]),
+                dictionary_xmap=dictionary.xmap,
+                keep_n=keep_n,
+                n_per_iteration=n_per_iteration,
+                metric=metric,
+            )
 
-    def average_neighbour_patterns(
+    def refine_orientation(
         self,
-        window: Union[str, np.ndarray, da.Array, Window] = "circular",
-        window_shape: Tuple[int, ...] = (3, 3),
-        show_progressbar: Optional[bool] = None,
-        **kwargs,
-    ) -> None:
-        """Average patterns inplace with its neighbours within a window.
+        xmap: CrystalMap,
+        detector: EBSDDetector,
+        master_pattern: "EBSDMasterPattern",
+        energy: Union[int, float],
+        signal_mask: Optional[np.ndarray] = None,
+        method: Optional[str] = "minimize",
+        method_kwargs: Optional[dict] = None,
+        trust_region: Union[tuple, list, np.ndarray, None] = None,
+        initial_step: Union[float] = None,
+        rtol: float = 1e-4,
+        maxeval: Optional[int] = None,
+        compute: bool = True,
+        rechunk: bool = True,
+        chunk_kwargs: Optional[dict] = None,
+    ) -> Union[CrystalMap, da.Array]:
+        r"""Refine orientations by searching orientation space around
+        the best indexed solution using fixed projection centers.
 
-        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.
+        Refinement attempts to maximize the similarity between patterns
+        in this signal and simulated patterns projected from a master
+        pattern. The similarity metric used is the normalized
+        cross-correlation (NCC). The orientation, represented by a
+        Euler angle triplet (:math:`\phi_1`, :math:`\Phi`,
+        :math:`\phi_2`) relative to the EDAX TSL sample reference frame
+        RD-TD-ND, is optimized during refinement, while the
+        sample-detector geometry, represented by the three projection
+        center (PC) parameters (PCx, PCy, PCz) in the Bruker convention,
+        is fixed.
+        A subset of the optimization methods in *SciPy* and *NLopt* are
+        available:
 
-        Averaging is accomplished by correlating the window with the
-        extended array of patterns using
-        :func:`scipy.ndimage.correlate`.
+        - Local optimization:
+            - :func:`~scipy.optimize.minimize` (includes Nelder-Mead,
+              Powell etc.).
+            - Nelder-Mead via `nlopt.LN_NELDERMEAD
+              <https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#nelder-mead-simplex>`_
+        - Global optimization:
+            - :func:`~scipy.optimize.differential_evolution`
+            - :func:`~scipy.optimize.dual_annealing`
+            - :func:`~scipy.optimize.basinhopping`
+            - :func:`~scipy.optimize.shgo`
 
         Parameters
         ----------
-        window
-            Name of averaging window or an array. Available types are
-            listed in :func:`scipy.signal.windows.get_window`, in
-            addition to a ``"circular"`` window (default) filled with
-            ones in which corner coefficients are set to zero. A window
-            element is considered to be in a corner if its radial
-            distance to the origin (window center) is shorter or equal
-            to the half width of the window's longest axis. A 1D or 2D
-            :class:`~numpy.ndarray`, :class:`~dask.array.Array` or
-            :class:`~kikuchipy.filters.Window` can also be passed.
-        window_shape
-            Shape of averaging window. Not used if a custom window or
-            :class:`~kikuchipy.filters.Window` is passed to ``window``.
-            This can be either 1D or 2D, and can be asymmetrical.
-            Default is ``(3, 3)``.
-        show_progressbar
-            Whether to show a progressbar. If not given, the value of
-            :obj:`hyperspy.api.preferences.General.show_progressbar`
-            is used.
-        **kwargs
-            Keyword arguments passed to the available window type listed
-            in :func:`~scipy.signal.windows.get_window`. If not given,
-            the default values of that particular window are used.
-
-        See Also
-        --------
-        kikuchipy.filters.Window, scipy.signal.windows.get_window,
-        scipy.ndimage.correlate
-        """
-        if isinstance(window, Window) and window.is_valid:
-            averaging_window = copy.copy(window)
-        else:
-            averaging_window = Window(window=window, shape=window_shape, **kwargs)
-
-        nav_shape = self.axes_manager.navigation_shape[::-1]
-        window_shape = averaging_window.shape
-        if window_shape in [(1,), (1, 1)]:
-            # Do nothing if a window of shape (1,) or (1, 1) is passed
-            return warnings.warn(
-                f"A window of shape {window_shape} was passed, no averaging is "
-                "therefore performed."
-            )
-        elif len(nav_shape) > len(window_shape):
-            averaging_window = averaging_window.reshape(window_shape + (1,))
-
-        # Get sum of window data for each pattern, to normalize with
-        # after correlation
-        window_sums = correlate(
-            input=np.ones(nav_shape, dtype=int),
-            weights=averaging_window,
-            mode="constant",
-        )
+        xmap
+            Single phase crystal map with at least one orientation per
+            point.
+        detector
+            Detector describing the detector-sample geometry with either
+            one PC to be used for all map points or one for each point.
+        master_pattern
+            Master pattern in the square Lambert projection of the same
+            phase as the one in the crystal map.
+        energy
+            Accelerating voltage of the electron beam in kV specifying
+            which master pattern energy to use during projection of
+            simulated patterns.
+        signal_mask
+            A boolean mask equal to the experimental patterns' detector
+            shape ``(n rows, n columns)``, where only pixels equal to
+            ``False`` are matched. If not given, all pixels are used.
+        method
+            Name of the :mod:`scipy.optimize` or *NLopt* optimization
+            method, among ``"minimize"``, ``"differential_evolution"``,
+            ``"dual_annealing"``, ``"basinhopping"``, ``"shgo"`` and
+            ``"ln_neldermead"`` (from *NLopt*). Default is
+            ``"minimize"``, which by default performs local optimization
+            with the Nelder-Mead method, unless another ``"minimize"``
+            method is passed to ``method_kwargs``.
+        method_kwargs
+            Keyword arguments passed to the :mod:`scipy.optimize`
+            ``method``. For example, to perform refinement with the
+            modified Powell algorithm from *SciPy*, pass
+            ``method="minimize"`` and
+            ``method_kwargs=dict(method="Powell")``. Not used if
+            ``method="LN_NELDERMEAD"``.
+        trust_region
+            List of three +/- angular deviations in degrees used to
+            determine the bound constraints on the three Euler angles
+            per navigation point, e.g. ``[2, 2, 2]``. Not passed to
+            *SciPy* ``method`` if it does not support bounds. The
+            definition ranges of the Euler angles are
+            :math:`\phi_1 \in [0, 360]`, :math:`\Phi \in [0, 180]` and
+            :math:`\phi_2 \in [0, 360]` in radians.
+        initial_step
+            A single initial step size for all Euler angle, in degrees.
+            Only used if ``method="LN_NELDERMEAD"``. If not given, this
+            is not set for the *NLopt* optimizer.
+        rtol
+            Stop optimization of a pattern when the difference in NCC
+            score between two iterations is below this value (relative
+            tolerance). Default is ``1e-4``. Only used if
+            ``method="LN_NELDERMEAD"``.
+        maxeval
+            Stop optimization of a pattern when the number of function
+            evaluations exceeds this value, e.g. ``100``. Only used if
+            ``method="LN_NELDERMEAD"``.
+        compute
+            Whether to refine now (``True``) or later (``False``).
+            Default is ``True``. See :meth:`~dask.array.Array.compute`
+            for more details.
+        rechunk
+            If ``True`` (default), rechunk the dask array with patterns
+            used in refinement (not the signal data inplace) if it is
+            returned from :func:`~kikuchipy.signals.util.get_dask_array`
+            in a single chunk. This ensures small data sets are
+            rechunked so as to utilize multiple CPUs.
+        chunk_kwargs
+            Keyword arguments passed to
+            :func:`~kikuchipy.signals.util.get_chunking` if
+            ``rechunk=True`` and the dask array with patterns used in
+            refinement is returned from
+            :func:`~kikuchipy.signals.util.get_dask_array` in a single
+            chunk.
 
-        # Add signal dimensions to window array to enable its use with
-        # Dask's map_overlap()
-        sig_dim = self.axes_manager.signal_dimension
-        averaging_window = averaging_window.reshape(
-            averaging_window.shape + (1,) * sig_dim
-        )
+        Returns
+        -------
+        out
+            Crystal map with refined orientations and similarity metrics
+            in a ``"scores"`` property if ``compute=True``. If
+            ``compute=False``, a dask array of navigation shape + (4,)
+            is returned, to be computed later. See
+            :func:`~kikuchipy.indexing.compute_refine_orientation_results`.
+            Each navigation point has the optimized score and the three
+            Euler angles in radians in element 0, 1, 2, and 3,
+            respectively.
 
-        # Create dask array of signal patterns and do processing on this
-        if self._lazy:
-            old_chunks = self.data.chunks
-        dask_array = get_dask_array(signal=self, chunk_bytes=8e6, rechunk=True)
+        See Also
+        --------
+        scipy.optimize, refine_projection_center,
+        refine_orientation_projection_center
 
-        # Add signal dimensions to array be able to use with Dask's
-        # map_overlap()
-        nav_dim = self.axes_manager.navigation_dimension
-        for i in range(sig_dim):
-            window_sums = np.expand_dims(window_sums, axis=window_sums.ndim)
-        window_sums = da.from_array(
-            window_sums, chunks=dask_array.chunks[:nav_dim] + (1,) * sig_dim
+        Notes
+        -----
+        *NLopt* is for now an optional dependency, see
+        :ref:`optional-dependencies` for details. Be aware that *NLopt*
+        does not fail gracefully. If continued use of *NLopt* proves
+        stable enough, its implementation of the Nelder-Mead algorithm
+        might become the default.
+        """
+        self._check_refinement_parameters(
+            xmap=xmap,
+            detector=detector,
+            master_pattern=master_pattern,
+            signal_mask=signal_mask,
         )
-
-        # Create overlap between chunks to enable correlation with the
-        # window using Dask's map_overlap()
-        window_dim = averaging_window.ndim
-        overlap_depth = {}
-        for i in range(nav_dim):
-            if i < window_dim and dask_array.chunks[i][0] < dask_array.shape[i]:
-                overlap_depth[i] = (window_shape[i] // 2) + 1
-            else:
-                overlap_depth[i] = 1
-        overlap_depth.update(
-            {i: 0 for i in self.axes_manager.signal_indices_in_array[::-1]}
+        patterns, signal_mask = self._prepare_patterns_for_refinement(
+            signal_mask=signal_mask, rechunk=rechunk, chunk_kwargs=chunk_kwargs
         )
-
-        dtype_out = self.data.dtype
-        omin, omax = dtype_range[dtype_out.type]
-        averaged_patterns = da.overlap.map_overlap(
-            _average_neighbour_patterns,
-            dask_array,
-            window_sums,
-            window=averaging_window,
-            dtype_out=dtype_out,
-            omin=omin,
-            omax=omax,
-            dtype=dtype_out,
-            depth=overlap_depth,
-            boundary="none",
+        return _refine_orientation(
+            xmap=xmap,
+            detector=detector,
+            master_pattern=master_pattern,
+            energy=energy,
+            patterns=patterns,
+            signal_mask=signal_mask,
+            method=method,
+            method_kwargs=method_kwargs,
+            trust_region=trust_region,
+            initial_step=initial_step,
+            rtol=rtol,
+            maxeval=maxeval,
+            compute=compute,
         )
 
-        # Overwrite signal patterns
-        if not self._lazy:
-            pbar = ProgressBar()
-            if show_progressbar or (
-                show_progressbar is None and hs.preferences.General.show_progressbar
-            ):
-                pbar.register()
-
-            averaged_patterns.store(self.data, compute=True)
-
-            try:
-                pbar.unregister()
-            except KeyError:
-                pass
-        else:
-            # Revert original chunks
-            averaged_patterns = averaged_patterns.rechunk(old_chunks)
-
-            self.data = averaged_patterns
+    def refine_projection_center(
+        self,
+        xmap: CrystalMap,
+        detector: EBSDDetector,
+        master_pattern: "EBSDMasterPattern",
+        energy: Union[int, float],
+        signal_mask: Optional[np.ndarray] = None,
+        method: Optional[str] = "minimize",
+        method_kwargs: Optional[dict] = None,
+        trust_region: Union[tuple, list, np.ndarray, None] = None,
+        initial_step: Union[float] = None,
+        rtol: float = 1e-4,
+        maxeval: Optional[int] = None,
+        compute: bool = True,
+        rechunk: bool = True,
+        chunk_kwargs: Optional[dict] = None,
+    ) -> Union[Tuple[np.ndarray, EBSDDetector], da.Array]:
+        """Refine projection centers by searching the parameter space
+        using fixed orientations.
 
-        # Don't sink
-        gc.collect()
+        Refinement attempts to maximize the similarity between patterns
+        in this signal and simulated patterns projected from a master
+        pattern. The similarity metric used is the normalized
+        cross-correlation (NCC). The sample-detector geometry,
+        represented by the three projection center (PC) parameters
+        (PCx, PCy, PCz) in the Bruker convention, is updated during
+        refinement, while the orientations, defined relative to the EDAX
+        TSL sample reference frame RD-TD-ND, are fixed.
 
-    def plot_virtual_bse_intensity(
-        self,
-        roi: BaseInteractiveROI,
-        out_signal_axes: Union[Iterable[int], Iterable[str], None] = None,
-        **kwargs,
-    ) -> None:
-        """Plot an interactive virtual backscatter electron (VBSE)
-        image formed from intensities within a specified and adjustable
-        region of interest (ROI) on the detector.
+        A subset of the optimization methods in *SciPy* and *NLopt* are
+        available:
 
-        Adapted from
-        :meth:`pyxem.signals.common_diffraction.CommonDiffraction.plot_integrated_intensity`.
+        - Local optimization:
+            - :func:`~scipy.optimize.minimize` (includes Nelder-Mead,
+              Powell etc.).
+            - Nelder-Mead via `nlopt.LN_NELDERMEAD
+              <https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#nelder-mead-simplex>`_
+        - Global optimization:
+            - :func:`~scipy.optimize.differential_evolution`
+            - :func:`~scipy.optimize.dual_annealing`
+            - :func:`~scipy.optimize.basinhopping`
+            - :func:`~scipy.optimize.shgo`
 
         Parameters
         ----------
-        roi
-            Any interactive ROI detailed in HyperSpy.
-        out_signal_axes
-            Which navigation axes to use as signal axes in the virtual
-            image. If not given, the first two navigation axes are used.
-        **kwargs:
-            Keyword arguments passed to the ``plot()`` method of the
-            virtual image.
-
-        See Also
-        --------
-        get_virtual_bse_intensity
-
-        Examples
-        --------
-        >>> import hyperspy.api as hs
-        >>> import kikuchipy as kp
-        >>> s = kp.data.nickel_ebsd_small()
-        >>> rect_roi = hs.roi.RectangularROI(
-        ...     left=0, right=5, top=0, bottom=5
-        ... )
-        >>> s.plot_virtual_bse_intensity(rect_roi)
-        """
-        # Plot signal if necessary
-        if self._plot is None or not self._plot.is_active:
-            self.plot()
-
-        # Get the sliced signal from the ROI
-        sliced_signal = roi.interactive(self, axes=self.axes_manager.signal_axes)
+        xmap
+            Single phase crystal map with at least one orientation per
+            point.
+        detector
+            Detector describing the detector-sample geometry with either
+            one PC to be used for all map points or one for each point.
+        master_pattern
+            Master pattern in the square Lambert projection of the same
+            phase as the one in the crystal map.
+        energy
+            Accelerating voltage of the electron beam in kV specifying
+            which master pattern energy to use during projection of
+            simulated patterns.
+        signal_mask
+            A boolean mask equal to the experimental patterns' detector
+            shape ``(n rows, n columns)``, where only pixels equal to
+            ``False`` are matched. If not given, all pixels are used.
+        method
+            Name of the :mod:`scipy.optimize` or *NLopt* optimization
+            method, among ``"minimize"``, ``"differential_evolution"``,
+            ``"dual_annealing"``, ``"basinhopping"``, ``"shgo"`` and
+            ``"ln_neldermead"`` (from *NLopt*). Default is
+            ``"minimize"``, which by default performs local optimization
+            with the Nelder-Mead method, unless another ``"minimize"``
+            method is passed to ``method_kwargs``.
+        method_kwargs
+            Keyword arguments passed to the :mod:`scipy.optimize`
+            ``method``. For example, to perform refinement with the
+            modified Powell algorithm from *SciPy*, pass
+            ``method="minimize"`` and
+            ``method_kwargs=dict(method="Powell")``. Not used if
+            ``method="LN_NELDERMEAD"``.
+        trust_region
+            List of three +/- deviations in the range [0, 1] used to
+            determine the bounds constraints on the PC parameters per
+            navigation point, e.g. ``[0.05, 0.05, 0.05]``. Not passed to
+            *SciPy* ``method`` if it does not support bounds. The
+            definition range of the PC parameters are assumed to be
+            [-2, 2].
+        initial_step
+            A single initial step size for all PC parameters in the
+            range [0, 1]. Only used if ``method="LN_NELDERMEAD"``.
+        rtol
+            Stop optimization of a pattern when the difference in NCC
+            score between two iterations is below this value (relative
+            tolerance). Default is ``1e-4``. Only used if
+            ``method="LN_NELDERMEAD"``.
+        maxeval
+            Stop optimization of a pattern when the number of function
+            evaluations exceeds this value, e.g. ``100``. Only used if
+            ``method="LN_NELDERMEAD"``.
+        compute
+            Whether to refine now (``True``) or later (``False``).
+            Default is ``True``. See :meth:`~dask.array.Array.compute`
+            for more details.
+        rechunk
+            If ``True`` (default), rechunk the dask array with patterns
+            used in refinement (not the signal data inplace) if it is
+            returned from :func:`~kikuchipy.signals.util.get_dask_array`
+            in a single chunk. This ensures small data sets are
+            rechunked so as to utilize multiple CPUs.
+        chunk_kwargs
+            Keyword arguments passed to
+            :func:`~kikuchipy.signals.util.get_chunking` if
+            ``rechunk=True`` and the dask array with patterns used in
+            refinement is returned from
+            :func:`~kikuchipy.signals.util.get_dask_array` in a single
+            chunk.
 
-        # Create an output signal for the virtual backscatter electron
-        # calculation
-        out = self._get_sum_signal(self, out_signal_axes)
-        out.metadata.General.title = "Virtual backscatter electron intensity"
+        Returns
+        -------
+        out
+            New similarity metrics and a new EBSD detector instance with
+            the refined PCs if ``compute=True``. If ``compute=False``,
+            a dask array of navigation shape + (4,) is returned, to be
+            computed later. See
+            :func:`~kikuchipy.indexing.compute_refine_projection_center_results`.
+            Each navigation point has the optimized score and the three
+            PC parameters in the Bruker convention in element 0, 1, 2,
+            and 3, respectively.
 
-        # Create the interactive signal
-        hs.interactive(
-            f=sliced_signal.nansum,
-            axis=sliced_signal.axes_manager.signal_axes,
-            event=roi.events.changed,
-            recompute_out_event=None,
-            out=out,
-        )
+        See Also
+        --------
+        scipy.optimize, refine_orientation,
+        refine_orientation_projection_center
 
-        # Plot the result
-        out.plot(**kwargs)
+        Notes
+        -----
+        *NLopt* is for now an optional dependency, see
+        :ref:`optional-dependencies` for details. Be aware that *NLopt*
+        does not fail gracefully. If continued use of *NLopt* proves
+        stable enough, its implementation of the Nelder-Mead algorithm
+        might become the default.
+        """
+        self._check_refinement_parameters(
+            xmap=xmap,
+            detector=detector,
+            master_pattern=master_pattern,
+            signal_mask=signal_mask,
+        )
+        patterns, signal_mask = self._prepare_patterns_for_refinement(
+            signal_mask=signal_mask, rechunk=rechunk, chunk_kwargs=chunk_kwargs
+        )
+        return _refine_pc(
+            xmap=xmap,
+            detector=detector,
+            master_pattern=master_pattern,
+            energy=energy,
+            patterns=patterns,
+            signal_mask=signal_mask,
+            method=method,
+            method_kwargs=method_kwargs,
+            trust_region=trust_region,
+            initial_step=initial_step,
+            rtol=rtol,
+            maxeval=maxeval,
+            compute=compute,
+        )
 
-    def get_virtual_bse_intensity(
+    def refine_orientation_projection_center(
         self,
-        roi: BaseInteractiveROI,
-        out_signal_axes: Union[Iterable[int], Iterable[str], None] = None,
-    ) -> VirtualBSEImage:
-        """Get a virtual backscatter electron (VBSE) image formed from
-        intensities within a region of interest (ROI) on the detector.
+        xmap: CrystalMap,
+        detector: EBSDDetector,
+        master_pattern: "EBSDMasterPattern",
+        energy: Union[int, float],
+        signal_mask: Optional[np.ndarray] = None,
+        method: Optional[str] = "minimize",
+        method_kwargs: Optional[dict] = None,
+        trust_region: Union[tuple, list, np.ndarray, None] = None,
+        initial_step: Union[tuple, list, np.ndarray, None] = None,
+        rtol: Optional[float] = 1e-4,
+        maxeval: Optional[int] = None,
+        compute: bool = True,
+        rechunk: bool = True,
+        chunk_kwargs: Optional[dict] = None,
+    ) -> Union[Tuple[CrystalMap, EBSDDetector], da.Array]:
+        r"""Refine orientations and projection centers simultaneously by
+        searching the orientation and PC parameter space.
 
-        Adapted from
-        :meth:`pyxem.signals.common_diffraction.CommonDiffraction.get_integrated_intensity`.
+        Refinement attempts to maximize the similarity between patterns
+        in this signal and simulated patterns projected from a master
+        pattern. The only supported similarity metric is the normalized
+        cross-correlation (NCC). The orientation, represented by a
+        Euler angle triplet (:math:`\phi_1`, :math:`\Phi`,
+        :math:`\phi_2`) relative to the EDAX TSL sample reference frame
+        RD-TD-ND, is optimized during refinement, while the
+        sample-detector geometry, represented by the three projection
+        center (PC) parameters (PCx, PCy, PCz) in the Bruker convention,
+        is fixed.
+
+        A subset of the optimization methods in *SciPy* and *NLopt* are
+        available:
+
+        - Local optimization:
+            - :func:`~scipy.optimize.minimize` (includes Nelder-Mead,
+              Powell etc.).
+            - Nelder-Mead via `nlopt.LN_NELDERMEAD
+              <https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#nelder-mead-simplex>`_
+        - Global optimization:
+            - :func:`~scipy.optimize.differential_evolution`
+            - :func:`~scipy.optimize.dual_annealing`
+            - :func:`~scipy.optimize.basinhopping`
+            - :func:`~scipy.optimize.shgo`
 
         Parameters
         ----------
-        roi
-            Any interactive ROI detailed in HyperSpy.
-        out_signal_axes
-            Which navigation axes to use as signal axes in the virtual
-            image. If not given, the first two navigation axes are used.
+        xmap
+            Single phase crystal map with at least one orientation per
+            point.
+        detector
+            Detector describing the detector-sample geometry with either
+            one PC to be used for all map points or one for each point.
+        master_pattern
+            Master pattern in the square Lambert projection of the same
+            phase as the one in the crystal map.
+        energy
+            Accelerating voltage of the electron beam in kV specifying
+            which master pattern energy to use during projection of
+            simulated patterns.
+        signal_mask
+            A boolean mask equal to the experimental patterns' detector
+            shape ``(n rows, n columns)``, where only pixels equal to
+            ``False`` are matched. If not given, all pixels are used.
+        method
+            Name of the :mod:`scipy.optimize` or *NLopt* optimization
+            method, among ``"minimize"``, ``"differential_evolution"``,
+            ``"dual_annealing"``, ``"basinhopping"``, ``"shgo"`` and
+            ``"ln_neldermead"`` (from *NLopt*). Default is
+            ``"minimize"``, which by default performs local optimization
+            with the Nelder-Mead method, unless another ``"minimize"``
+            method is passed to ``method_kwargs``.
+        method_kwargs
+            Keyword arguments passed to the :mod:`scipy.optimize`
+            ``method``. For example, to perform refinement with the
+            modified Powell algorithm from *SciPy*, pass
+            ``method="minimize"`` and
+            ``method_kwargs=dict(method="Powell")``. Not used if
+            ``method="LN_NELDERMEAD"``.
+        trust_region
+            List of three +/- angular deviations in degrees as bound
+            constraints on the three Euler angles and three +/-
+            deviations in the range [0, 1] as bound constraints on the
+            PC parameters, e.g. ``[2, 2, 2, 0.05, 0.05, 0.05]``. Not
+            passed to *SciPy* ``method`` if it does not support bounds.
+            The definition ranges of the Euler angles are
+            :math:`\phi_1 \in [0, 360]`, :math:`\Phi \in [0, 180]` and
+            :math:`\phi_2 \in [0, 360]` in radians, while the definition
+            range of the PC parameters are assumed to be [-2, 2].
+        initial_step
+            A list of two initial step sizes to use, one in degrees for
+            all Euler angles and one in the range [0, 1] for all PC
+            parameters. Only used if ``method="LN_NELDERMEAD"``.
+        rtol
+            Stop optimization of a pattern when the difference in NCC
+            score between two iterations is below this value (relative
+            tolerance). Only used if ``method="LN_NELDERMEAD"``. If not
+            given, this is set to ``1e-4``.
+        maxeval
+            Stop optimization of a pattern when the number of function
+            evaluations exceeds this value, e.g. ``100``. Only used if
+            ``method="LN_NELDERMEAD"``.
+        compute
+            Whether to refine now (``True``) or later (``False``).
+            Default is ``True``. See :meth:`~dask.array.Array.compute`
+            for more details.
+        rechunk
+            If ``True`` (default), rechunk the dask array with patterns
+            used in refinement (not the signal data inplace) if it is
+            returned from :func:`~kikuchipy.signals.util.get_dask_array`
+            in a single chunk. This ensures small data sets are
+            rechunked so as to utilize multiple CPUs.
+        chunk_kwargs
+            Keyword arguments passed to
+            :func:`~kikuchipy.signals.util.get_chunking` if
+            ``rechunk=True`` and the dask array with patterns used in
+            refinement is returned from
+            :func:`~kikuchipy.signals.util.get_dask_array` in a single
+            chunk.
 
         Returns
         -------
-        virtual_image
-            VBSE image formed from detector intensities within an ROI
-            on the detector.
+        out
+            Crystal map with refined orientations and a new EBSD
+            detector instance with the refined PCs, if ``compute=True``.
+            If ``compute=False``, a dask array of navigation shape +
+            (7,) is returned, to be computed later. See
+            :func:`~kikuchipy.indexing.compute_refine_orientation_projection_center_results`.
+            Each navigation point has the optimized score, the three
+            Euler angles in radians, and the three PC parameters in
+            element 0, 1, 2, 3, 4, 5, and 6, respectively.
 
         See Also
         --------
-        plot_virtual_bse_intensity
+        scipy.optimize, refine_orientation, refine_projection_center
 
-        Examples
-        --------
-        >>> import hyperspy.api as hs
-        >>> import kikuchipy as kp
-        >>> rect_roi = hs.roi.RectangularROI(
-        ...     left=0, right=5, top=0, bottom=5
-        ... )
-        >>> s = kp.data.nickel_ebsd_small()
-        >>> vbse_image = s.get_virtual_bse_intensity(rect_roi)
+        Notes
+        -----
+        The method attempts to refine the orientations and projection
+        center at the same time for each map point. The optimization
+        landscape is sloppy :cite:`pang2020global`, where the
+        orientation and PC can make up for each other. Thus, it is
+        possible that the parameters that yield the highest similarity
+        are incorrect. As always, it is left to the user to ensure that
+        the output is reasonable.
+
+        *NLopt* is for now an optional dependency, see
+        :ref:`optional-dependencies` for details. Be aware that *NLopt*
+        does not fail gracefully. If continued use of *NLopt* proves
+        stable enough, its implementation of the Nelder-Mead algorithm
+        might become the default.
         """
-        vbse = roi(self, axes=self.axes_manager.signal_axes)
-        vbse_sum = self._get_sum_signal(vbse, out_signal_axes)
-        vbse_sum.metadata.General.title = "Virtual backscatter electron image"
-        vbse_sum.set_signal_type("VirtualBSEImage")
-        return vbse_sum
+        self._check_refinement_parameters(
+            xmap=xmap,
+            detector=detector,
+            master_pattern=master_pattern,
+            signal_mask=signal_mask,
+        )
+        patterns, signal_mask = self._prepare_patterns_for_refinement(
+            signal_mask=signal_mask, rechunk=rechunk, chunk_kwargs=chunk_kwargs
+        )
+        return _refine_orientation_pc(
+            xmap=xmap,
+            detector=detector,
+            master_pattern=master_pattern,
+            energy=energy,
+            patterns=patterns,
+            signal_mask=signal_mask,
+            method=method,
+            method_kwargs=method_kwargs,
+            trust_region=trust_region,
+            initial_step=initial_step,
+            rtol=rtol,
+            maxeval=maxeval,
+            compute=compute,
+        )
 
     # ------ Methods overwritten from hyperspy.signals.Signal2D ------ #
 
@@ -2555,10 +2637,10 @@ def rebin(self, *args, **kwargs):
         nav_shape_new = am_new.navigation_shape[::-1]
         sig_shape_new = am_new.signal_shape[::-1]
 
-        props = self._get_custom_attributes(make_deepcopy=True)
+        attrs = self._get_custom_attributes(make_deepcopy=True)
 
         # Update static background
-        static_bg = props["static_background"]
+        static_bg = attrs["static_background"]
         if sig_shape_new != sig_shape_old and static_bg is not None:
             sig_idx = am_old.signal_indices_in_array[::-1]
             if params["new_shape"] is not None:
@@ -2569,23 +2651,23 @@ def rebin(self, *args, **kwargs):
             s_static_bg2 = s_static_bg.rebin(**params)
             static_bg2 = s_static_bg2.data
             static_bg2 = static_bg2.astype(new.data.dtype)
-            props["static_background"] = static_bg2
+            attrs["static_background"] = static_bg2
 
         # Update detector shape and binning factor
-        props["detector"].shape = sig_shape_new
+        attrs["detector"].shape = sig_shape_new
         factors = np.array(sig_shape_old) / np.array(sig_shape_new)
-        binning = props["detector"].binning * factors
+        binning = attrs["detector"].binning * factors
         if binning[0] == binning[1] and np.allclose(binning, binning.round(0)):
-            props["detector"].binning = int(binning[0])
+            attrs["detector"].binning = int(binning[0])
         else:
-            props["detector"].binning = 1
+            attrs["detector"].binning = 1
 
         if nav_shape_new != nav_shape_old:
-            props["xmap"] = None
-            props["detector"].pc = np.full(nav_shape_new + (3,), 0.5)
-            props["static_background"] = None
+            attrs["xmap"] = None
+            attrs["detector"].pc = np.full(nav_shape_new + (3,), 0.5)
+            attrs["static_background"] = None
 
-        new._set_custom_attributes(props)
+        new._set_custom_attributes(attrs)
 
         return new
 


Test patch:
diff --git a/kikuchipy/pattern/tests/test_chunk.py b/kikuchipy/pattern/tests/test_chunk.py
index 6f272473..547df228 100644
--- a/kikuchipy/pattern/tests/test_chunk.py
+++ b/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",
diff --git a/kikuchipy/pattern/tests/test_pattern.py b/kikuchipy/pattern/tests/test_pattern.py
index be3d0c04..f415894a 100644
--- a/kikuchipy/pattern/tests/test_pattern.py
+++ b/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)
diff --git a/kikuchipy/signals/tests/test_ebsd.py b/kikuchipy/signals/tests/test_ebsd.py
index e222b87f..ffe8884e 100644
--- a/kikuchipy/signals/tests/test_ebsd.py
+++ b/kikuchipy/signals/tests/test_ebsd.py
@@ -2115,3 +2115,57 @@ def test_extract_grid_1d(self):
         s5, idx5 = s4.extract_grid((6,), return_indices=True)
         assert s5.data.shape == (6, 60, 60)
         assert np.allclose(idx5, [10, 21, 32, 43, 54, 65])
+
+
+class TestDownsample:
+    def test_downsample(self):
+        s = kp.data.nickel_ebsd_small()
+        s2 = s.deepcopy()
+
+        s2.downsample(factor=2)
+
+        # Initial signal unaffected
+        assert s.data.shape == (3, 3, 60, 60)
+        assert s.detector.shape == (60, 60)
+        assert s.static_background.shape == (60, 60)
+
+        # Correct shape and data type
+        assert s2.data.shape == (3, 3, 30, 30)
+        assert s2.data.dtype == s.data.dtype
+        assert s2.detector.shape == (30, 30)
+        assert s2.static_background.shape == (30, 30)
+        assert s2.static_background.dtype == s.static_background.dtype
+
+    def test_downsample_lazy(self):
+        s = kp.data.nickel_ebsd_small(lazy=True)
+
+        s.downsample(factor=3)
+
+        assert isinstance(s, kp.signals.LazyEBSD)
+        assert isinstance(s.data, da.Array)
+
+        s.compute()
+
+        assert isinstance(s, kp.signals.EBSD)
+        assert isinstance(s.data, np.ndarray)
+
+    def test_downsample_dtype(self):
+        s = kp.data.nickel_ebsd_small()
+
+        s2 = s.deepcopy()
+        s2.downsample(factor=2, dtype_out="uint16")
+        assert s2.data.dtype.name == "uint16"
+
+    def test_downsample_raises(self):
+        s = kp.data.nickel_ebsd_small(lazy=True)
+
+        with pytest.raises(ValueError, match="Binning `factor` 2.5 must be an integer"):
+            s.downsample(2.5)
+
+        with pytest.raises(ValueError, match="Binning `factor` 1 must be an integer >"):
+            s.downsample(1)
+
+        with pytest.raises(
+            ValueError, match="Binning `factor` 7 must be a divisor of the initial "
+        ):
+            s.downsample(7)
diff --git a/kikuchipy/signals/tests/test_ebsd_hough_indexing.py b/kikuchipy/signals/tests/test_ebsd_hough_indexing.py
index d07f2a47..dcb91f83 100644
--- a/kikuchipy/signals/tests/test_ebsd_hough_indexing.py
+++ b/kikuchipy/signals/tests/test_ebsd_hough_indexing.py
@@ -282,11 +282,11 @@ def test_optimize_pc(self):
         not kp._pyebsdindex_installed, reason="pyebsdindex is not installed"
     )
     def test_optimize_pc_pso(self):
-        np.random.seed(42)  # Make PSO deterministic
         det = self.signal.hough_indexing_optimize_pc(
             self.detector.pc_average, self.indexer, method="PSO"
         )
-        assert np.allclose(det.pc_average, [0.422, 0.218, 0.500], atol=1e-3)
+        # Results are not deterministic, so we give a wide range here...
+        assert abs(self.detector.pc_average - det.pc_average).max() < 0.05
 
     @pytest.mark.skipif(
         not kp._pyebsdindex_installed, reason="pyebsdindex is not installed"