Repository: pyxem/kikuchipy (language: python)

Issue:
Orientation correction after dictionary indexing (and PC correction) before refinement
We should provide dynamic orientation and projection center corrections following the procedure in [Singh et al. (2017)](https://doi.org/10.1107/S1600576717014200). Python versions were implemented by @friedkitteh in [his thesis](https://ntnuopen.ntnu.no/ntnu-xmlui/handle/11250/2785332). I think both corrections can be available as separate methods to the `EBSDDetector`.


Discussion / hints:
Duplicate of #352:

> We should add a projection/pattern center (PC) calculation method to kikuchipy.detectors.EBSDDetector which populates a map of navigation shape with PCs from the geometrical relations in Eqs. (8) and (9) in Appendix A in Singh et al. (2017).
> 
> This could then be used as a first approximation of PCs per point, instead of a constant PC, in orientation and/or PC refinement following dictionary indexing, given that the detector pixel size and navigation step size in both directons are known.
These resources from the xcdskd package could be of interest as well:
* https://github.com/wiai/xcdskd/blob/master/src/aloe/exp/calibpc.py
* https://github.com/wiai/xcdskd/blob/master/src/aloe/exp/pcfitter.py

Fix patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 92da7a44..a60064f2 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -18,6 +18,22 @@ Unreleased
 
 Added
 -----
+- Method ``EBSDDetector.fit_pc()`` to fit a plane using a projective or affine
+  transformation to projection centers following work by Winkelmann and co-workers
+  (2020). (`#587 <https://github.com/pyxem/kikuchipy/pull/587>`_)
+- Method ``EBSDDetector.extrapolate_pc()`` to return a new detector with a plane of
+  projection centers (PCs) extrapolated from a mean PC calculated from one or more PCs
+  following work by Singh et al. (2017).
+  (`#587 <https://github.com/pyxem/kikuchipy/pull/587>`_)
+- Methods ``EBSDDetector.estimate_xtilt()`` and ``EBSDDetector.estimate_xtilt_ztilt()``
+  to estimate the tilts about the detector X and Z axes which bring the detector plane
+  normal parallel to the sample plane normal, following work by Winkelmann and
+  co-workers (2020). (`#587 <https://github.com/pyxem/kikuchipy/pull/587>`_)
+- Method ``EBSDDetector.plot_pc()`` to plot projection centers (PCs) in maps, scatter
+  plots or in 3D. (`#587 <https://github.com/pyxem/kikuchipy/pull/587>`_)
+- Convenience function ``kikuchipy.draw.plot_pattern_positions_in_map()`` to plot
+  positions of selected patterns (typically calibration patterns) in a 2D map.
+  (`#587 <https://github.com/pyxem/kikuchipy/pull/587>`_)
 - EBSD signal returned from NORDIF calibration pattern reader tries to add the following
   new info to the original metadata: Shapes of area and region of interest (ROI), offset
   of ROI, calibration pattern indices and area overview image. All shapes and
@@ -84,6 +100,8 @@ Removed
 
 Fixed
 -----
+- ``EBSDDetector.pc_average`` no longer rounds the PC to three decimals.
+  (`#586 <https://github.com/pyxem/kikuchipy/pull/586>`_)
 - Microscope magnification is now read correctly from EDAX h5ebsd files.
   (`#586 <https://github.com/pyxem/kikuchipy/pull/586>`_)
 - kikuchipy h5ebsd reader can read a signal with an EBSD detector with a PC array of
diff --git a/examples/reference_frames/README.rst b/examples/reference_frames/README.rst
new file mode 100644
index 00000000..6dda0742
--- /dev/null
+++ b/examples/reference_frames/README.rst
@@ -0,0 +1,5 @@
+Reference frames
+================
+
+These examples cover use of tools like the :class:`~kikuchipy.detectors.EBSDDetector`
+for understanding reference frames relevant to EBSD.
\ No newline at end of file
diff --git a/examples/reference_frames/estimate_tilt_about_detector_xaxis.py b/examples/reference_frames/estimate_tilt_about_detector_xaxis.py
new file mode 100644
index 00000000..9b65ace7
--- /dev/null
+++ b/examples/reference_frames/estimate_tilt_about_detector_xaxis.py
@@ -0,0 +1,64 @@
+r"""
+Estimate tilt about the detector x axis
+=======================================
+
+This example shows how to (robustly) estimate the tilt about the detector :math:`X_d`
+axis which brings the sample plane normal into coincidence with the detector plane
+normal (but in the opposite direction) :cite:`winkelmann2020refined`.
+
+The estimate is found using :meth:`~kikuchipy.detectors.EBSDDetector.estimate_xtilt`
+which performs linear regression of :attr:`~kikuchipy.detectors.EBSDDetector.pcz` vs.
+:attr:`~kikuchipy.detectors.EBSDDetector.pcy`.
+
+To test the estimation, we add some noise to realistic projection center (PC) values
+(PCx, PCy, PCz). The realistic PCs are extrapolated from a PC in the upper left corner
+of a map, assuming a nominal sample tilt of 70 degrees, a detector tilt of 0 degrees, a
+detector pixel size of 70 microns and a sample step size of 50 microns.
+"""
+
+import kikuchipy as kp
+import numpy as np
+
+
+# Create an initial detector with one PC assumed to be for the upper
+# left corner of a map
+det0 = kp.detectors.EBSDDetector(
+    shape=(480, 480),
+    pc=(0.5, 0.3, 0.5),
+    sample_tilt=70,
+    tilt=0,
+    px_size=70,
+)
+
+# Extrapolate a map of PCs
+nav_shape = (15, 20)
+det = det0.extrapolate_pc(
+    pc_indices=[0, 0],
+    navigation_shape=nav_shape,
+    step_sizes=(50, 50),
+)
+
+# Add +/- 0.0025 as random noise to PCy and PCz
+rng = np.random.default_rng()
+det.pcy += rng.uniform(-0.001, 0.001, det.navigation_size).reshape(nav_shape)
+det.pcz += rng.uniform(-0.001, 0.001, det.navigation_size).reshape(nav_shape)
+
+# Add outliers by adding more noise
+outlier_idx1d = rng.choice(det.navigation_size, 10, replace=False)
+is_outlier = np.zeros(det.navigation_size, dtype=bool)
+is_outlier[outlier_idx1d] = True
+noise_outlier = rng.uniform(-0.01, 0.01, outlier_idx1d.size)
+outlier_idx2d = np.unravel_index(outlier_idx1d, shape=det.navigation_shape)
+det.pcz[outlier_idx2d] += noise_outlier
+
+# Robust estimation by detecting outliers
+xtilt, detected_outliers = det.estimate_xtilt(
+    detect_outliers=True, degrees=True, return_outliers=True
+)
+
+# Print true tilt and estimated tilt
+true_tilt = 90 - det.sample_tilt + det.tilt
+print(f"True/estimated tilt about detector x [deg]: {true_tilt:.2f}/{xtilt:.2f}")
+
+correct_outliers = np.isin(outlier_idx1d, np.where(detected_outliers)[0])
+print(f"{correct_outliers.sum()}/{outlier_idx1d.size} of added outliers detected")
diff --git a/examples/reference_frames/estimate_tilts_about_detector_xaxis_zaxis.py b/examples/reference_frames/estimate_tilts_about_detector_xaxis_zaxis.py
new file mode 100644
index 00000000..5a17a54d
--- /dev/null
+++ b/examples/reference_frames/estimate_tilts_about_detector_xaxis_zaxis.py
@@ -0,0 +1,72 @@
+r"""
+Estimate tilts about the detector x and z axis
+==============================================
+
+This example shows how to (robustly) estimate the tilts about the detector :math:`X_d`
+and :math:`Z_d` axes, which bring the sample plane normal into coincidence with the
+detector plane normal (but in the opposite direction) :cite:`winkelmann2020refined`.
+
+Estimates are found using
+:meth:`~kikuchipy.detectors.EBSDDetector.estimate_xtilt_ztilt`, which fits a hyperplane
+to :attr:`~kikuchipy.detectors.EBSDDetector.pc` using singular value decomposition.
+
+To test the estimations, we add some noise to realistic projection center (PC) values
+(PCx, PCy, PCz). The realistic PCs are extrapolated from a PC in the upper left corner
+of a map, assuming a nominal sample tilt of 70 degrees, a detector tilt of 0 degrees, a
+detector pixel size of 70 microns and a sample step size of 50 microns.
+"""
+
+import kikuchipy as kp
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+plt.rcParams["font.size"] = 18
+
+# Create an initial detector with one PC assumed to be for the upper
+# left corner of a map
+det0 = kp.detectors.EBSDDetector(
+    shape=(480, 480),
+    pc=(0.5, 0.3, 0.5),
+    sample_tilt=70,
+    tilt=0,
+    px_size=70,
+)
+
+# Extrapolate a map of PCs
+nav_shape = (15, 20)
+nav_size = np.prod(nav_shape)
+det = det0.extrapolate_pc(
+    pc_indices=[0, 0],
+    navigation_shape=nav_shape,
+    step_sizes=(50, 50),
+)
+
+# Add +/- 0.0025 as random noise
+dev_noise = 0.001
+rng = np.random.default_rng()
+det.pcx += rng.uniform(-dev_noise, dev_noise, nav_size).reshape(nav_shape)
+det.pcy += rng.uniform(-dev_noise, dev_noise, nav_size).reshape(nav_shape)
+det.pcz += rng.uniform(-dev_noise, dev_noise, nav_size).reshape(nav_shape)
+
+# Add outliers by adding more noise
+dev_outlier = 0.01
+n_outliers = 20
+outlier_idx1d = rng.choice(nav_size, n_outliers, replace=False)
+is_outlier = np.zeros(nav_size, dtype=bool)
+is_outlier[outlier_idx1d] = True
+outlier_idx2d = np.unravel_index(outlier_idx1d, shape=det.navigation_shape)
+det.pcx[outlier_idx2d] += rng.uniform(-dev_outlier, dev_outlier, n_outliers)
+det.pcy[outlier_idx2d] += rng.uniform(-dev_outlier, dev_outlier, n_outliers)
+det.pcz[outlier_idx2d] += rng.uniform(-dev_outlier, dev_outlier, n_outliers)
+
+# Plot PC values
+det.plot_pc("scatter")
+
+# Robust estimation by detecting outliers
+xtilt, ztilt = det.estimate_xtilt_ztilt(degrees=True)
+
+# Print true tilt and estimated tilt
+true_xtilt = 90 - det.sample_tilt + det.tilt
+print(f"True/estimated tilt about detector x [deg]: {true_xtilt:.2f}/{xtilt:.2f}")
+print(f"True/estimated tilt about detector z [deg]: {0}/{ztilt:.2f}")
diff --git a/examples/reference_frames/fit_plane_to_selected_pc.py b/examples/reference_frames/fit_plane_to_selected_pc.py
new file mode 100644
index 00000000..87d21f6d
--- /dev/null
+++ b/examples/reference_frames/fit_plane_to_selected_pc.py
@@ -0,0 +1,60 @@
+r"""
+Fit a plane to selected projection centers
+==========================================
+
+This example shows how to fit a plane to selected projection centers (PCs) using a
+projective transformation, following :cite:`winkelmann2020refined`.
+
+To test the fit, we add some noise to realistic projection center (PC) values
+(PCx, PCy, PCz) and fit a plane to a few of the PCs. The realistic PCs are extrapolated
+from a PC in the upper left corner of a map, assuming a nominal sample tilt of 70
+degrees, a detector tilt of 0 degrees, a detector pixel size of 70 microns and a sample
+step size of 50 microns.
+"""
+
+import kikuchipy as kp
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+plt.rcParams["font.size"] = 10
+
+# Create an initial detector with one PC assumed to be for the upper
+# left corner of a map
+det0 = kp.detectors.EBSDDetector(
+    shape=(480, 640),
+    pc=(0.5, 0.3, 0.5),
+    sample_tilt=70,
+    tilt=0,
+    px_size=70,
+)
+
+# Extrapolate a map of PCs
+nav_shape = (30, 45)
+det1 = det0.extrapolate_pc(
+    pc_indices=[0, 0],
+    navigation_shape=nav_shape,
+    step_sizes=(50, 50),
+)
+
+# Add random noise
+rng = np.random.default_rng()
+dev = 0.002
+det1.pcx += rng.uniform(-dev, dev, det1.navigation_size).reshape(nav_shape)
+det1.pcy += rng.uniform(-dev, dev, det1.navigation_size).reshape(nav_shape)
+det1.pcz += rng.uniform(-dev, dev, det1.navigation_size).reshape(nav_shape)
+
+# Extract a (7, 7) grid of PCs
+grid_shape = (7, 7)
+pc_indices = kp.signals.util.grid_indices(grid_shape, nav_shape=nav_shape)
+det2 = det1.deepcopy()
+det2.pc = det2.pc[tuple(pc_indices)].reshape(grid_shape + (3,))
+
+# Get a plane of PCs and plot the match at the same time
+map_indices = np.stack(np.indices(nav_shape))
+det_fit = det2.fit_pc(pc_indices=pc_indices, map_indices=map_indices)
+
+# Inspect the max. error and sample tilt
+max_err = abs(det_fit.pc_flattened - det1.pc_flattened).max(axis=0)
+print("Max error in (PCx, PCy, PCz):", max_err)
+print(f"Estimated sample tilt [deg]: {det_fit.sample_tilt:.2f}")
diff --git a/examples/selecting_data/extract_patterns_from_grid.py b/examples/selecting_data/extract_patterns_from_grid.py
index 7c585144..c235f627 100644
--- a/examples/selecting_data/extract_patterns_from_grid.py
+++ b/examples/selecting_data/extract_patterns_from_grid.py
@@ -9,9 +9,10 @@
 import hyperspy.api as hs
 import kikuchipy as kp
 import matplotlib.pyplot as plt
-import numpy as np
 
 
+plt.rcParams["font.size"] = 15
+
 # Silence progressbars
 hs.preferences.General.show_progressbar = False
 
@@ -31,11 +32,11 @@
 vbse_img.rescale_intensity(dtype_out="float32", percentiles=(0.5, 99.5))
 
 # Plot grid of extracted patterns
-fig, ax = plt.subplots()
-ax.imshow(vbse_img.data, cmap="gray")
-ax.scatter(*idx[::-1], c=np.arange(s2.axes_manager.navigation_size), s=400, ec="k")
-ax.axis("off")
-fig.tight_layout()
+kp.draw.plot_pattern_positions_in_map(
+    idx.reshape(2, -1).T,
+    roi_shape=s.axes_manager.navigation_shape[::-1],
+    roi_image=vbse_img.data,
+)
 
 # Plot extracted patterns
 s2.remove_static_background()
diff --git a/kikuchipy/detectors/ebsd_detector.py b/kikuchipy/detectors/ebsd_detector.py
index 4ec50912..7a8d236a 100644
--- a/kikuchipy/detectors/ebsd_detector.py
+++ b/kikuchipy/detectors/ebsd_detector.py
@@ -17,13 +17,31 @@
 
 from __future__ import annotations
 from copy import deepcopy
+import logging
 from typing import List, Optional, Tuple, Union
 
 from matplotlib.figure import Figure
 from matplotlib.markers import MarkerStyle
 import matplotlib.patches as mpatches
 import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
 import numpy as np
+from orix.quaternion import Rotation
+from orix.vector import Vector3d
+import scipy.stats as scs
+from skimage.transform import ProjectiveTransform
+from sklearn.linear_model import LinearRegression, RANSACRegressor
+
+
+_logger = logging.getLogger(__name__)
+
+CONVENTION_ALIAS = {
+    "bruker": ["bruker"],
+    "tsl": ["edax", "tsl", "amatek"],
+    "oxford": ["oxford", "aztec"],
+    "emsoft": ["emsoft", "emsoft4", "emsoft5"],
+}
+CONVENTION_ALIAS_ALL = list(np.concatenate(list(CONVENTION_ALIAS.values())))
 
 
 class EBSDDetector:
@@ -67,10 +85,11 @@ class EBSDDetector:
         ``[0.5, 0.5, 0.5]``.
     convention
         PC convention. If not given, Bruker's convention is assumed.
-        Options are ``"tsl"``/``"edax"``, ``"oxford"``, ``"bruker"``,
-        ``"emsoft"``, ``"emsoft4"``, and ``"emsoft5"``. ``"emsoft"`` and
-        ``"emsoft5"`` is the same convention. See *Notes* for
-        conversions between conventions.
+        Options are ``"tsl"``/``"edax"``/``"amatek"``,
+        ``"oxford"``/``"aztec"``, ``"bruker"``, ``"emsoft"``,
+        ``"emsoft4"``, and ``"emsoft5"``. ``"emsoft"`` and ``"emsoft5"``
+        is the same convention. See *Notes* for conversions between
+        conventions.
 
     Notes
     -----
@@ -159,9 +178,11 @@ def __init__(
         sample_tilt: float = 70,
         pc: Union[np.ndarray, list, tuple] = (0.5, 0.5, 0.5),
         convention: Optional[str] = None,
-    ):
-        """Create an EBSD detector with a shape, pixel size, binning,
-        and projection/pattern center(s) (PC(s)).
+    ) -> None:
+        """Create an EBSD detector with a shape, pixel size, binning
+        factor, sample and detector tilt about the detector X axis,
+        azimuthal tilt about the detector Y axis and one or more
+        projection/pattern centers (PCs).
         """
         self.shape = shape
         self.px_size = px_size
@@ -170,13 +191,14 @@ def __init__(
         self.azimuthal = azimuthal
         self.sample_tilt = sample_tilt
         self.pc = pc
-        self._set_pc_convention(convention)
+        self._set_pc_from_convention(convention)
 
     def __repr__(self) -> str:
+        pc_average = tuple(self.pc_average.round(3))
         return (
-            f"{self.__class__.__name__} {self.shape}, "
+            f"{type(self).__name__} {self.shape}, "
             f"px_size {self.px_size} um, binning {self.binning}, "
-            f"tilt {self.tilt}, azimuthal {self.azimuthal}, pc {tuple(self.pc_average)}"
+            f"tilt {self.tilt}, azimuthal {self.azimuthal}, pc {pc_average}"
         )
 
     @property
@@ -244,6 +266,13 @@ def pc(self, value: Union[np.ndarray, List, Tuple]):
         """Set all projection center coordinates."""
         self._pc = np.atleast_2d(value)
 
+    @property
+    def pc_flattened(self) -> np.ndarray:
+        """Return flattened array of projection center coordinates of
+        shape (:attr:`navigation_size`, 3).
+        """
+        return self.pc.reshape((-1, 3))
+
     @property
     def pcx(self) -> np.ndarray:
         """Return or set the projection center x coordinates.
@@ -307,7 +336,7 @@ def pc_average(self) -> np.ndarray:
             axis += (0,)
         elif ndim == 3:
             axis += (0, 1)
-        return np.nanmean(self.pc, axis=axis).round(3)
+        return np.nanmean(self.pc, axis=axis)
 
     @property
     def navigation_shape(self) -> tuple:
@@ -337,6 +366,11 @@ def navigation_dimension(self) -> int:
         """
         return len(self.navigation_shape)
 
+    @property
+    def navigation_size(self) -> int:
+        """Return the number of projection centers."""
+        return int(np.prod(self.navigation_shape))
+
     @property
     def bounds(self) -> np.ndarray:
         """Return the detector bounds ``[x0, x1, y0, y1]`` in pixel
@@ -440,7 +474,7 @@ def crop(self, extent: Union[Tuple[int, int, int, int], List[int]]) -> EBSDDetec
         >>> det.crop((1, 5, 2, 6))
         EBSDDetector (4, 4), px_size 1 um, binning 1, tilt 0, azimuthal 0, pc (0.25, 0.25, 0.75)
 
-        Plot a cropped detector with the PC on cropped a pattern
+        Plot a cropped detector with the PC on a cropped pattern
 
         >>> s = kp.data.nickel_ebsd_small()
         >>> s.remove_static_background(show_progressbar=False)
@@ -492,6 +526,493 @@ def deepcopy(self) -> EBSDDetector:
         """
         return deepcopy(self)
 
+    def estimate_xtilt(
+        self,
+        detect_outliers: bool = False,
+        plot: bool = True,
+        degrees: bool = False,
+        return_figure: bool = False,
+        return_outliers: bool = False,
+        figure_kwargs: Optional[dict] = None,
+    ) -> Union[
+        float,
+        Tuple[float, np.ndarray],
+        Tuple[float, plt.Figure],
+        Tuple[float, np.ndarray, plt.Figure],
+    ]:
+        r"""Estimate the tilt about the detector :math:`X_d` axis.
+
+        This tilt is assumed to bring the sample plane normal into
+        coincidence with the detector plane normal (but in the opposite
+        direction) :cite:`winkelmann2020refined`.
+
+        See the `reference frame tutorial
+        </doc/tutorials/reference_frames.ipynb>`_ for details on the
+        detector sample geometry.
+
+        An estimate is found by linear regression of :attr:`pcz` vs.
+        :attr:`pcy`.
+
+        Parameters
+        ----------
+        detect_outliers
+            Whether to attempt to detect outliers. If ``False``
+            (default), a linear fit to all points is performed. If
+            ``True``, a robust fit using the RANSAC algorithm is
+            performed instead, which also detects outliers.
+        plot
+            Whether to plot data points and the estimated line. Default
+            is ``True``.
+        degrees
+            Whether to return the estimated tilt in radians (``False``,
+            default) or degrees (``True``).
+        return_figure
+            Whether to return the plotted figure. Default is ``False``.
+        return_outliers
+            Whether to return a mask with ``True`` for PC values
+            considered outliers. Default is ``False``. If ``True``,
+            ``detect_outliers`` is assumed to be ``True`` and the value
+            passed is not considered.
+        figure_kwargs
+            Keyword arguments passed to
+            :func:`matplotlib.pyplot.Figure` if ``plot=True``.
+
+        Returns
+        -------
+        x_tilt
+            Estimated tilt about detector :math:`X_d` in radians
+            (``degrees=False``) or degrees (``degrees=True``).
+        outliers
+            Returned if ``return_outliers=True``.
+        fig
+            Returned if ``plot=True`` and ``return_figure=True``.
+
+        Notes
+        -----
+        This method is adapted from Aimo Winkelmann's function
+        ``fit_xtilt()`` in the *xcdskd* Python package. See
+        :cite:`winkelmann2020refined` for their use of related
+        functions.
+
+        See Also
+        --------
+        sklearn.linear_model.LinearRegression,
+        sklearn.linear_model.RANSACRegressor, estimate_xtilt_ztilt,
+        fit_pc
+        """
+        if self.navigation_size == 1:
+            raise ValueError("Estimation requires more than one projection center")
+
+        if return_outliers:
+            detect_outliers = True
+
+        # Get regressor
+        if detect_outliers:
+            regressor = RANSACRegressor()
+        else:
+            regressor = LinearRegression()
+
+        # Estimate slope
+        pcy = self.pcy.reshape((-1, 1))
+        pcz = self.pcz.reshape((-1, 1))
+        regressor.fit(pcz, pcy)
+        if detect_outliers:
+            slope = regressor.estimator_.coef_
+        else:
+            slope = regressor.coef_
+        slope = slope.squeeze()
+
+        # Get tilt from estimated coefficient
+        x_tilt = np.pi / 2 + np.arctan(slope)
+
+        # Get outliers (if possible)
+        if detect_outliers:
+            is_outlier = ~regressor.inlier_mask_
+        else:
+            is_outlier = np.zeros(self.navigation_size)
+
+        # Predict data of estimated model
+        pcz_fit = np.linspace(np.min(pcz), np.max(pcz), 2)
+        pcy_fit = regressor.predict(pcz_fit[:, np.newaxis])
+
+        if degrees:
+            x_tilt = np.rad2deg(x_tilt)
+
+        out = (x_tilt,)
+        if return_outliers:
+            out += (is_outlier,)
+
+        if plot:
+            if figure_kwargs is None:
+                figure_kwargs = {}
+            figure_kwargs.setdefault("layout", "tight")
+            fig = plt.figure(**figure_kwargs)
+            ax = fig.add_subplot()
+            ax.scatter(pcz, pcy, c="yellowgreen", ec="k", label="Data")
+            if detect_outliers:
+                ax.scatter(pcz[is_outlier], pcy[is_outlier], c="gold", label="Outliers")
+            fit_label = "Fit, tilt = " + f"{x_tilt:.2f}" + r"$^{\circ}$"
+            ax.plot(pcz_fit, pcy_fit, label=fit_label, c="C1")
+            ax.legend(bbox_to_anchor=(1.05, 1), loc="upper left", borderaxespad=0.0)
+            ax.set(aspect="equal", xlabel="PCz", ylabel="PCy")
+
+            if return_figure:
+                out += (fig,)
+
+        if len(out) == 1:
+            out = out[0]
+
+        return out
+
+    def estimate_xtilt_ztilt(
+        self,
+        degrees: bool = False,
+        is_outlier: Optional[Union[list, tuple, np.ndarray]] = None,
+    ) -> Union[float, Tuple[float, float]]:
+        r"""Estimate the tilts about the detector :math:`X_d` and
+        :math:`Z_d` axes.
+
+        These tilts bring the sample plane normal into coincidence with
+        the detector plane normal (but in the opposite direction)
+        :cite:`winkelmann2020refined`.
+
+        See the `reference frame tutorial
+        </doc/tutorials/reference_frames.ipynb>`_ for details on the
+        detector sample geometry.
+
+        Estimates are found by fitting a hyperplane to :attr:`pc` using
+        singular value decomposition.
+
+        Parameters
+        ----------
+        degrees
+            Whether to return the estimated tilts in radians (``False``,
+            default) or degrees (``True``).
+        is_outlier
+            Boolean array with ``True`` for PCs to not include in the
+            fit. If not given, all PCs are used. Must be of
+            :attr:`navigation_shape`.
+
+        Returns
+        -------
+        x_tilt
+            Estimated tilt about detector :math:`X_d` in radians
+            (``degrees=False``) or degrees (``degrees=True``).
+        z_tilt
+            Estimated tilt about detector :math:`Z_d` in radians
+            (``degrees=False``) or degrees (``degrees=True``).
+
+        Notes
+        -----
+        This method is adapted from Aimo Winkelmann's function
+        ``fit_plane()`` in the *xcdskd* Python package. Its use is
+        described in :cite:`winkelmann2020refined`.
+
+        Winkelmann refers to Gander & Hrebicek, "Solving Problems in
+        Scientific Computing", 3rd Ed., Chapter 6, p. 97 for the
+        implementation of the hyperplane fitting.
+
+        See Also
+        --------
+        estimate_xtilt, fit_pc
+        """
+        if self.navigation_size == 1:
+            raise ValueError("Estimation requires more than one projection center")
+
+        pc = self.pc
+        if isinstance(is_outlier, np.ndarray):
+            pc = pc[~is_outlier]
+        pc = pc.reshape((-1, 3))
+
+        pc_centered = pc - pc.mean(axis=0)
+        x_tilt, z_tilt, *_ = _fit_hyperplane(pc_centered)
+
+        if degrees:
+            x_tilt = np.rad2deg(x_tilt)
+            z_tilt = np.rad2deg(z_tilt)
+
+        return x_tilt, z_tilt
+
+    def extrapolate_pc(
+        self,
+        pc_indices: Union[tuple, list, np.ndarray],
+        navigation_shape: tuple,
+        step_sizes: tuple,
+        shape: Optional[tuple] = None,
+        px_size: float = None,
+        binning: int = None,
+        is_outlier: Optional[Union[tuple, list, np.ndarray]] = None,
+    ):
+        r"""Return a new detector with projection centers (PCs)
+        extrapolated from an average PC.
+
+        The average PC :math:`\bar{PC}` is calculated from :attr:`pc`,
+        possibly excluding some PCs based on the ``is_outlier`` mask.
+        The sample position having this PC, :math:`(\bar{x}, \bar{y})`,
+        is assumed to be the one obtained by averaging ``pc_indices``.
+        All other PCs :math:`(PC_x, PC_y, PC_z)` in positions
+        :math:`(x, y)` are then extrapolated based on the following
+        equations given in appendix A in :cite:`singh2017application`:
+
+        .. math::
+            PC_x &= \bar{PC_x} + (\bar{x} - x) \cdot \Delta x / (\delta \cdot N_x \cdot b),\\
+            PC_y &= \bar{PC_y} + (\bar{y} - y) \cdot \Delta y \cdot \cos{\alpha} / (\delta \cdot N_y \cdot b),\\
+            PC_z &= \bar{PC_z} - (\bar{y} - y) \cdot \Delta y \cdot \sin{\alpha} / (\delta \cdot N_y \cdot b),\\
+
+        where :math:`(\Delta y, \Delta x)` are the vertical and
+        horizontal step sizes, respectively, :math:`(N_y, N_x)` are the
+        number of binned detector rows and columns, respectively, the
+        angle :math:`\alpha = 90^{\circ} - \sigma + \theta`, where
+        :math:`\sigma` is the sample tilt and :math:`\theta` is the
+        detector tilt, :math:`\delta` is the unbinned detector pixel
+        size and :math:`b` is the binning factor.
+
+        Parameters
+        ----------
+        pc_indices
+            2D map pixel coordinates (row, column) of each :attr:`pc`,
+            possibly outside ``navigation_shape``. Must be a flattened
+            array of shape (2,) + :attr:`navigation_size`.
+        navigation_shape
+            Shape of the output PC array (n rows, n columns).
+        step_sizes
+            Vertical and horizontal step sizes (dy, dx).
+        shape
+            Detector (signal) shape (n rows, n columns). If not given,
+            (:attr:`nrows`, :attr:`ncols`) is used.
+        px_size
+            Unbinned detector pixel size. If not given, :attr:`px_size`
+            is used.
+        binning
+            Detector binning factor. If not given, :attr:`binning` is
+            used.
+        is_outlier
+            Boolean array with ``True`` for PCs to not include in the
+            fit. If not given, all PCs are used. Must be of
+            :attr:`navigation_shape`.
+
+        Returns
+        -------
+        new_detector
+            Detector with :attr:`navigation_shape` given by
+            input ``navigation_shape``.
+        """
+        # Parse input parameters
+        dy, dx = step_sizes
+        pc_indices = np.atleast_2d(pc_indices).T
+        ny, nx = navigation_shape
+        if not shape:
+            shape = self.shape
+        nrows, ncols = shape
+        if not px_size:
+            px_size = self.px_size
+        if not binning:
+            binning = self.binning
+
+        pc = self.pc_flattened
+
+        if is_outlier is not None:
+            is_outlier = np.asarray(is_outlier)
+            pc = pc[~is_outlier]
+            pc_indices = pc_indices[:, ~is_outlier]
+
+        # Calculate mean PC and position
+        pc_mean = pc.mean(axis=0)
+        pc_indices_mean = pc_indices.mean(axis=1)
+        pc_indices_mean = np.round(pc_indices_mean).astype(int)
+
+        # Make PC plane
+        alpha = np.deg2rad(90 - self.sample_tilt + self.tilt)
+        y, x = np.indices((ny, nx), dtype=float)
+        factor = px_size * binning
+        d_pcx = -(pc_indices_mean[1] - x) * dx / (factor * ncols)
+        d_pcy = -(pc_indices_mean[0] - y) * dy * np.cos(alpha) / (factor * nrows)
+        d_pcz = +(pc_indices_mean[0] - y) * dy * np.sin(alpha) / (factor * nrows)
+        pcx = pc_mean[0] - d_pcx
+        pcy = pc_mean[1] - d_pcy
+        pcz = pc_mean[2] - d_pcz
+
+        new_detector = self.deepcopy()
+        new_detector.pc = np.stack((pcx, pcy, pcz), axis=2)
+        new_detector.shape = shape
+        new_detector.px_size = px_size
+        new_detector.binning = binning
+
+        return new_detector
+
+    def fit_pc(
+        self,
+        pc_indices: Union[list, tuple, np.ndarray],
+        map_indices: Union[list, tuple, np.ndarray],
+        transformation: str = "projective",
+        is_outlier: Optional[np.ndarray] = None,
+        plot: bool = True,
+        return_figure: bool = False,
+        figure_kwargs: Optional[dict] = None,
+    ) -> Union[EBSDDetector, Tuple[EBSDDetector, plt.Figure]]:
+        """Return a new detector with interpolated projection centers
+        (PCs) for all points in a map by fitting a plane to :attr:`pc`
+        :cite:`winkelmann2020refined`.
+
+        Parameters
+        ----------
+        pc_indices
+            2D coordinates (row, column) of each :attr:`pc` in
+            ``map_coordinates``. Must be a flattened array of shape
+            (2,) + :attr:`navigation_shape`.
+        map_indices
+            2D coordinates (row, column) of all map points in a regular
+            grid to interpolate PCs for. Must be a flattened array of
+            shape ``(2,) + map_shape``.
+        transformation
+            Which transformation function to use when fitting PCs,
+            either ``"projective"`` (default) or ``"affine"``. Both
+            transformations perserve co-planarity of map points, while
+            the projective transformation allows parallel lines in the
+            map point grid to become non-parallel within the sample
+            plane.
+        is_outlier
+            Boolean array with ``True`` for PCs to not include in the
+            fit. If not given, all PCs are used. Must be of
+            :attr:`navigation_shape`.
+        plot
+            Whether to plot the experimental and estimated PCs (default
+            is ``True``).
+        return_figure
+            Whether to return the figure if ``plot=True`` (default is
+            ``False``).
+        figure_kwargs
+            Keyword arguments passed to
+            :func:`matplotlib.pyplot.Figure` if ``plot=True``.
+
+        Returns
+        -------
+        new_detector
+            New detector with as many interpolated PCs as indices given
+            in ``map_indices`` and an estimated sample tilt. The
+            detector tilt is assumed to be constant.
+        fig
+            Figure of experimental and estimated PCs, returned if
+            ``plot=True`` and ``return_figure=True``.
+
+        Raises
+        ------
+        ValueError
+            If :attr:`navigation_size` is 1 or if the ``pc_indices`` or
+            ``map_indices`` arrays have the incorrect shape.
+
+        See Also
+        --------
+        estimate_xtilt, estimate_xtilt_ztilt, extrapolate_pc
+
+        Notes
+        -----
+        This method is adapted from Aimo Winkelmann's functions
+        ``fit_affine()`` and ``fit_projective()`` in the *xcdskd* Python
+        package. Their uses are described in
+        :cite:`winkelmann2020refined`. Winkelmann refers to a code
+        example from StackOverflow
+        (https://stackoverflow.com/a/20555267/3228100) for the affine
+        transformation.
+        """
+        n_pc = self.navigation_size
+        if n_pc == 1:
+            raise ValueError("Fitting requires multiple projection centers (PCs)")
+
+        pc_indices = np.asarray(pc_indices)
+        map_indices = np.asarray(map_indices)
+
+        nav_shape = self.navigation_shape
+        if pc_indices.shape != (2,) + nav_shape:
+            raise ValueError(
+                f"`pc_indices` array shape {pc_indices.shape} must be equal to "
+                f"{(2,) + nav_shape}"
+            )
+
+        if map_indices.ndim not in [2, 3] or map_indices.shape[0] != 2:
+            raise ValueError(
+                f"`map_indices` array shape {map_indices.shape} must be (2, m columns) "
+                "or (2, n rows, m columns)"
+            )
+
+        if is_outlier is not None and (
+            not isinstance(is_outlier, np.ndarray)
+            or not np.issubdtype(is_outlier.dtype, np.bool_)
+            or is_outlier.size != n_pc
+        ):
+            raise ValueError(
+                "`is_outlier` must be a boolean array of a size equal to the number of "
+                "PCs"
+            )
+
+        # Prepare PCs and PC map indices for fitting
+        pc_flat = self.pc_flattened
+        pc_indices_flat = pc_indices.reshape((2, -1)).T
+        pc_indices_flat = np.column_stack((pc_indices_flat, np.ones(n_pc)))
+
+        # Prepare all map indices for projection
+        map_indices_flat = map_indices.reshape((2, -1)).T
+        map_indices_flat = np.column_stack(
+            (map_indices_flat, np.ones(map_indices_flat.shape[0]))
+        )
+
+        if isinstance(is_outlier, np.ndarray):
+            is_outlier = is_outlier.ravel()
+            nav_shape = (np.sum(~is_outlier),)
+            pc_flat = pc_flat[~is_outlier]
+            pc_indices_flat = pc_indices_flat[~is_outlier]
+
+        if transformation == "projective":
+            pc_average = np.mean(pc_flat, axis=0)
+            pc_flat_centered = pc_flat - pc_average
+            pc_fit, pc_fit_map = _fit_pc_projective(
+                pc_flat_centered,
+                pc_indices_flat,
+                map_indices_flat,
+            )
+            pc_fit += pc_average
+            pc_fit_map += pc_average
+        else:
+            pc_fit, pc_fit_map = _fit_pc_affine(
+                pc_flat, pc_indices_flat, map_indices_flat
+            )
+
+        max_err = np.max(np.abs(pc_flat - pc_fit), axis=0)
+        _logger.debug(f"Max. error for (PCx, PCy, PCz): {max_err}")
+
+        # Linear fit to YZ-plane to estimate tilt angle about detector X
+        slope, intercept, _, _, _ = scs.linregress(pc_fit[..., 2], pc_fit[..., 1])
+        x_tilt = np.pi / 2 + np.arctan(slope)
+        x_tilt_deg = np.rad2deg(x_tilt)
+        _logger.debug(f"Estimated detector X tilt: {x_tilt_deg:.5f} deg")
+
+        # Reshape new fitted PC array to desired map shape
+        pc_fit_map = pc_fit_map.reshape(map_indices.shape[1:] + (3,))
+
+        new_detector = self.deepcopy()
+        new_detector.pc = pc_fit_map
+        new_detector.sample_tilt = 90 - x_tilt_deg - new_detector.tilt
+
+        if plot:
+            if figure_kwargs is None:
+                figure_kwargs = {}
+            pc_fit_2d = pc_fit.reshape(nav_shape + (3,))
+            fig = _plot_pc_fit(
+                pc_flat,
+                pc_fit_2d,
+                intercept,
+                slope,
+                return_figure=return_figure,
+                figure_kwargs=figure_kwargs,
+            )
+        else:
+            fig = None
+
+        if return_figure and fig:
+            return new_detector, fig
+        else:
+            return new_detector
+
     def pc_emsoft(self, version: int = 5) -> np.ndarray:
         r"""Return PC in the EMsoft convention.
 
@@ -669,6 +1190,7 @@ def plot(
 
         Examples
         --------
+        >>> import matplotlib.pyplot as plt
         >>> import kikuchipy as kp
         >>> det = kp.detectors.EBSDDetector(
         ...     shape=(60, 60),
@@ -677,6 +1199,7 @@ def plot(
         ...     sample_tilt=70,
         ... )
         >>> det.plot()
+        >>> plt.show()
 
         Plot with gnomonic coordinates and circles
 
@@ -685,12 +1208,12 @@ def plot(
         ...     draw_gnomonic_circles=True,
         ...     gnomonic_circles_kwargs={"edgecolor": "b", "alpha": 0.3}
         ... )
+        >>> plt.show()
 
-        Plot a pattern on the detector and save it
+        Plot a pattern on the detector and return it for saving etc.
 
         >>> s = kp.data.nickel_ebsd_small()
         >>> fig = det.plot(pattern=s.inav[0, 0].data, return_figure=True)
-        >>> # fig.savefig("detector.png")
         """
         sy, sx = self.shape
         pcx, pcy = self.pc_average[:2]
@@ -710,9 +1233,7 @@ def plot(
 
         fig, ax = plt.subplots()
         ax.axis(zoom * bounds)
-        ax.set_aspect("equal")
-        ax.set_xlabel(x_label)
-        ax.set_ylabel(y_label)
+        ax.set(xlabel=x_label, ylabel=y_label, aspect="equal")
 
         # Plot a pattern on the detector
         if isinstance(pattern, np.ndarray):
@@ -773,35 +1294,192 @@ def plot(
         if return_figure:
             return fig
 
+    def plot_pc(
+        self,
+        mode: str = "map",
+        return_figure: bool = False,
+        orientation: str = "horizontal",
+        annotate: bool = False,
+        figure_kwargs: Optional[dict] = None,
+        **kwargs,
+    ) -> Union[None, plt.Figure]:
+        """Plot all projection centers (PCs).
+
+        Parameters
+        ----------
+        mode
+            String describing how to plot PCs. Options are ``"map"``
+            (default), ``"scatter"`` and ``"3d"``. If ``mode="map"``,
+            :attr:`navigation_dimension` must be 2.
+        return_figure
+            Whether to return the figure (default is ``False``).
+        orientation
+            Whether to align the plots in a ``"horizontal"`` (default)
+            or ``"vertical"`` orientation.
+        annotate
+            Whether to label each pattern with its 1D index into
+            :attr:`pc_flattened` when ``mode="scatter"``. Default is
+            ``False``.
+        figure_kwargs
+            Keyword arguments to pass to
+            :func:`matplotlib.pyplot.figure` upon figure creation. Note
+            that ``layout="tight"`` is used by default unless another
+            layout is passed.
+        **kwargs
+            Keyword arguments passed to the plotting function, which is
+            :meth:`~matplotlib.axes.Axes.imshow` if ``mode="map"``,
+            :meth:`~matplotlib.axes.Axes.scatter` if ``mode="scatter"``
+            and :meth:`~mpl_toolkits.mplot3d.axes3d.Axes3D.scatter`
+            if ``mode="3d"``.
+
+        Returns
+        -------
+        fig
+            Figure is returned if ``return_figure=True``.
+
+        Examples
+        --------
+        Create a detector with smoothly changing PC values, extrapolated
+        from a single PC (assumed to be in the upper left corner of a
+        map)
+
+        >>> import matplotlib.pyplot as plt
+        >>> import kikuchipy as kp
+        >>> det0 = kp.detectors.EBSDDetector(
+        ... shape=(480, 640), pc=(0.4, 0.3, 0.5), px_size=70, sample_tilt=70
+        ... )
+        >>> det0
+        EBSDDetector (480, 640), px_size 70 um, binning 1, tilt 0, azimuthal 0, pc (0.4, 0.3, 0.5)
+        >>> det = det0.extrapolate_pc(
+        ... pc_indices=[0, 0], navigation_shape=(5, 10), step_sizes=(20, 20)
+        ... )
+        >>> det
+        EBSDDetector (480, 640), px_size 70 um, binning 1, tilt 0, azimuthal 0, pc (0.398, 0.299, 0.5)
+
+        Plot PC values in maps
+
+        >>> det.plot_pc()
+        >>> plt.show()
+
+        Plot in scatter plots in vertical orientation
+
+        >>> det.plot_pc("scatter", orientation="vertical", annotate=True)
+        >>> plt.show()
+
+        Plot in a 3D scatter plot, returning the figure for saving etc.
+
+        >>> fig = det.plot_pc("3d", return_figure=True)
+        """
+        # Ensure there are PCs to plot
+        if self.navigation_size == 1:
+            raise ValueError(
+                "Detector must have more than one projection center value to plot"
+            )
+        if mode == "map" and self.navigation_dimension != 2:
+            raise ValueError("Detector's `navigation_dimension` must be 2D")
+
+        # Ensure mode is OK
+        modes = ["map", "scatter", "3d"]
+        if not isinstance(mode, str) or mode.lower() not in modes:
+            raise ValueError(
+                f"Plot mode '{mode}' must be one of the following strings {modes}"
+            )
+        mode = mode.lower()
+
+        if figure_kwargs is None:
+            figure_kwargs = {}
+        figure_kwargs.setdefault("layout", "tight")
+
+        # Prepare keyword arguments common to at least two modes
+        if mode in ["map", "scatter"]:
+            w, h = plt.rcParams["figure.figsize"]
+            k = max(w, h) / 3
+            if orientation == "horizontal":
+                figure_kwargs.setdefault("figsize", (6 * k, 2 * k))
+                subplots_kw = dict(ncols=3)
+            else:
+                figure_kwargs.setdefault("figsize", (2 * k, 6 * k))
+                subplots_kw = dict(nrows=3)
+
+        if mode in ["scatter", "3d"]:
+            kwargs.setdefault("c", np.arange(self.navigation_size))
+            kwargs.setdefault("ec", "k")
+            kwargs.setdefault("clip_on", False)
+
+        fig = plt.figure(**figure_kwargs)
+
+        labels = ["PCx", "PCy", "PCz"]
+        if mode == "map":
+            axes = fig.subplots(**subplots_kw)
+            for i, ax in enumerate(axes):
+                ax.set(xlabel="Column", ylabel="Row", aspect="equal")
+                im = ax.imshow(self.pc[..., i], **kwargs)
+                divider = make_axes_locatable(ax)
+                cax = divider.append_axes(position="right", size="5%", pad=0.1)
+                fig.colorbar(im, cax=cax, label=labels[i])
+        elif mode == "scatter":
+            pc_flat = self.pc_flattened
+            axes = fig.subplots(**subplots_kw)
+            for i, (j, k) in enumerate([[0, 1], [0, 2], [2, 1]]):
+                x_coord = pc_flat[:, j]
+                y_coord = pc_flat[:, k]
+                axes[i].scatter(x_coord, y_coord, **kwargs)
+                axes[i].set(xlabel=labels[j], ylabel=labels[k], aspect="equal")
+                if annotate:
+                    for l, (x, y) in enumerate(zip(x_coord, y_coord)):
+                        axes[i].text(x, y, l, ha="left", va="bottom")
+            axes[0].invert_xaxis()
+            axes[1].invert_xaxis()
+            axes[1].invert_yaxis()
+        else:
+            ax = fig.add_subplot(projection="3d")
+
+            pcx, pcy, pcz = self.pc_flattened.T
+            ax.scatter(pcx, pcz, pcy, **kwargs)
+            nav_axes = tuple(np.arange(len(self.pc.shape))[: self.navigation_dimension])
+            extent_min = np.min(self.pc, axis=nav_axes)
+            extent_max = np.max(self.pc, axis=nav_axes)
+            ax.set(
+                xlabel=labels[0],
+                ylabel=labels[2],
+                zlabel=labels[1],
+                xlim=[extent_min[0], extent_max[0]],
+                ylim=[extent_min[2], extent_max[2]],
+                zlim=[extent_min[1], extent_max[1]],
+            )
+            ax.invert_zaxis()
+
+            if annotate:
+                for i, (x, z, y) in enumerate(zip(pcx, pcz, pcy)):
+                    ax.text(x, z, y, i)
+
+        if return_figure:
+            return fig
+
     # ------------------------ Private methods ----------------------- #
 
-    def _set_pc_convention(self, convention: Optional[str] = None):
-        if convention is None or convention.lower() == "bruker":
-            pass
-        elif convention.lower() in ["tsl", "edax", "amatek"]:
-            self.pc = self._pc_tsl2bruker()
-        elif convention.lower() == "oxford":
-            self.pc = self._pc_tsl2bruker()
-        elif convention.lower() in ["emsoft", "emsoft4", "emsoft5"]:
+    def _get_pc_from_convention(self, convention: Optional[str] = None) -> np.ndarray:
+        if convention is None or convention.lower() in CONVENTION_ALIAS["bruker"]:
+            return self.pc
+
+        conv = convention.lower()
+        if conv in CONVENTION_ALIAS["tsl"] + CONVENTION_ALIAS["oxford"]:
+            return self._pc_tsl2bruker()
+        elif conv in CONVENTION_ALIAS["emsoft"]:
             try:
                 version = int(convention[-1])
             except ValueError:
                 version = 5
-            self.pc = self._pc_emsoft2bruker(version=version)
+            return self._pc_emsoft2bruker(version)
         else:
-            conventions = [
-                "bruker",
-                "emsoft",
-                "emsoft4",
-                "emsoft5",
-                "oxford",
-                "tsl",
-            ]
             raise ValueError(
                 f"Projection center convention '{convention}' not among the "
-                f"recognised conventions {conventions}."
+                f"recognised conventions {CONVENTION_ALIAS_ALL}."
             )
 
+    def _set_pc_from_convention(self, convention: Optional[str] = None):
+        self.pc = self._get_pc_from_convention(convention)
+
     def _pc_emsoft2bruker(self, version: int = 5) -> np.ndarray:
         new_pc = np.zeros_like(self.pc, dtype=float)
         pcx = self.pcx
@@ -832,3 +1510,158 @@ def _pc_bruker2tsl(self) -> np.ndarray:
         new_pc[..., 1] = (1 - self.pcy) / self.aspect_ratio
         new_pc[..., 2] /= self.aspect_ratio
         return new_pc
+
+
+def _fit_hyperplane(
+    pc_centered: np.ndarray,
+) -> Tuple[float, float, Rotation, Rotation, np.ndarray]:
+    # Hyperplane fit
+    pc_trim_mean = scs.trim_mean(pc_centered, proportiontocut=0.1)
+    pc_trim_centered = pc_centered - pc_trim_mean[np.newaxis, :]
+    # u @ np.diag(s) @ vh = (u * s) @ vh
+    u, s, vh = np.linalg.svd(pc_trim_centered, full_matrices=False)
+
+    # Check handedness of the coordinate system spanned by the plane
+    # normals. The determinant should be 1.
+    determinant = np.linalg.det(vh)
+    _logger.debug(f"Determinant of SVD: {determinant:.7f}")
+
+    # Extract estimated sample plane unit normal vector
+    sample_normal = Vector3d(vh[2])
+    sample_normal = sample_normal.unit
+    if sample_normal.z < 0:
+        # Make normal point towards detector screen
+        sample_normal = -sample_normal
+    _logger.debug(f"Sample plane normal from SVD: {sample_normal.data.squeeze()}")
+
+    # Tilt about the detector x and z axes
+    vx, vy, vz = sample_normal.data.squeeze()
+    x_tilt = np.arccos(vz)
+    z_tilt = np.pi / 2 - np.arctan2(vy, vx)
+
+    # Check that rotation of [001] gives the surface normal in the
+    # detector system
+    rot_xtilt = Rotation.from_axes_angles([1, 0, 0], -x_tilt)
+    rot_ztilt = Rotation.from_axes_angles([0, 0, 1], -z_tilt)
+    sample_normal_tilts = rot_ztilt * rot_xtilt * Vector3d.zvector()
+    _logger.debug(
+        f"Sample plane normal from tilts: {sample_normal_tilts.data.squeeze()}"
+    )
+
+    return x_tilt, z_tilt, rot_xtilt, rot_ztilt, pc_trim_mean
+
+
+def _fit_pc_projective(
+    pc_centered_flat: np.ndarray,
+    pc_indices_flat: np.ndarray,
+    map_indices_flat: np.ndarray,
+) -> Tuple[np.ndarray, np.ndarray]:
+    *_, rot_xtilt, rot_ztilt, pc_trim_mean = _fit_hyperplane(pc_centered_flat)
+
+    v_pc_centered = Vector3d(pc_centered_flat)
+    v_pc_trim_mean = Vector3d(pc_trim_mean)
+    v_pc_plane = ~(rot_ztilt * rot_xtilt) * (v_pc_centered - v_pc_trim_mean)
+
+    # Get transformation matrix
+    tform = ProjectiveTransform()
+    status = tform.estimate(pc_indices_flat[:, :2], v_pc_plane.data[:, :2])
+    _logger.debug(f"Status of projective transformation: {status}")
+    matrix = tform.params.T
+
+    # PC coordinates projected from beam indices and tilt parameters
+    pc_fit = np.dot(pc_indices_flat, matrix)
+    pc_fit /= pc_fit[:, 2, None]
+    pc_fit[:, 2] = 0
+    v_pc_fit = rot_ztilt * rot_xtilt * Vector3d(pc_fit)
+    v_pc_fit += v_pc_trim_mean
+
+    # Do the same for interpolated PC coordinates
+    pc_fit_map = np.dot(map_indices_flat, matrix)
+    pc_fit_map /= pc_fit_map[:, 2, None]
+    pc_fit_map[:, 2] = 0
+    v_pc_fit_map = rot_ztilt * rot_xtilt * Vector3d(pc_fit_map)
+    v_pc_fit_map += v_pc_trim_mean
+
+    return v_pc_fit.data, v_pc_fit_map.data
+
+
+def _fit_pc_affine(
+    pc_flat: np.ndarray, pc_indices_flat: np.ndarray, map_indices_flat: np.ndarray
+) -> Tuple[np.array, np.ndarray]:
+    # Solve the least squares problem X * A = Y
+    # Source: https://stackoverflow.com/a/20555267/3228100
+    matrix, res, *_ = np.linalg.lstsq(pc_indices_flat, pc_flat, rcond=None)
+    _logger.debug(f"Residuals of least squares fit: {res}")
+
+    pc_fit = np.dot(pc_indices_flat, matrix)
+    pc_fit_map = np.dot(map_indices_flat, matrix)
+
+    return pc_fit, pc_fit_map
+
+
+def _plot_pc_fit(
+    pc: np.ndarray,
+    pc_fit: np.ndarray,
+    fit_intercept: float,
+    fit_slope: float,
+    figure_kwargs: dict,
+    return_figure: bool = False,
+) -> Union[None, plt.Figure]:
+    pcx, pcy, pcz = pc.T
+    pcx_fit_2d, pcy_fit_2d, pcz_fit_2d = pc_fit.T
+    pcx_fit = pcx_fit_2d.ravel()
+    pcy_fit = pcy_fit_2d.ravel()
+    pcz_fit = pcz_fit_2d.ravel()
+
+    pcy_fit_line = fit_intercept + fit_slope * pcz_fit
+
+    data_kw = dict(s=25, c="k")
+    fit_kw = dict(s=50, fc="gray", alpha=0.5, ec="k")
+
+    w, h = plt.rcParams["figure.figsize"]
+    figure_kwargs.setdefault("layout", "compressed")
+    figure_kwargs.setdefault("figsize", (w, h))
+
+    fig = plt.figure()
+    # PCx v PCy
+    ax0 = fig.add_subplot(221)
+    ax0.set(xlabel="PCx", ylabel="PCy", aspect="equal")
+    ax0.scatter(pcx_fit, pcy_fit, **fit_kw)
+    ax0.scatter(pcx, pcy, **data_kw)
+    ax0.invert_xaxis()
+    # PCx v PCz
+    ax1 = fig.add_subplot(222)
+    ax1.set(xlabel="PCx", ylabel="PCz", aspect="equal")
+    ax1.scatter(pcx, pcz, **data_kw)
+    ax1.scatter(pcx_fit, pcz_fit, **fit_kw)
+    ax1.invert_xaxis()
+    ax1.invert_yaxis()
+    # PCz v PCy
+    ax2 = fig.add_subplot(223)
+    ax2.set(xlabel="PCz", ylabel="PCy", aspect="equal")
+    ax2.scatter(pcz, pcy, label="Data", **data_kw)
+    ax2.scatter(pcz_fit, pcy_fit, label="Fit", **fit_kw)
+    ax2.plot(pcz_fit, pcy_fit_line, "r-", label="Linear fit")
+    ax2.legend(bbox_to_anchor=(1.05, 1), loc="upper left", borderaxespad=0.0)
+    # Plot PC values in 3D
+    ax3 = fig.add_subplot(224, projection="3d")
+    ax3.scatter(pcx, pcz, pcy, **data_kw)
+    ax3.scatter(pcx_fit, pcz_fit, pcy_fit, **fit_kw)
+    ax3.set(
+        xlabel="PCx",
+        ylabel="PCz",
+        zlabel="PCy",
+        xlim=[np.min([pcx, pcx_fit]), np.max([pcx, pcx_fit])],
+        ylim=[np.min([pcz, pcz_fit]), np.max([pcz, pcz_fit])],
+        zlim=[np.min([pcy, pcy_fit]), np.max([pcy, pcy_fit])],
+    )
+    ax3.invert_zaxis()
+
+    # Add 3D plane
+    if pcx_fit_2d.ndim == 2:
+        ax3.plot_surface(pcx_fit_2d, pcz_fit_2d, pcy_fit_2d, color="r", alpha=0.5)
+
+    fig.tight_layout()
+
+    if return_figure:
+        return fig
diff --git a/kikuchipy/draw/__init__.py b/kikuchipy/draw/__init__.py
index e7c854ee..d1380fb0 100644
--- a/kikuchipy/draw/__init__.py
+++ b/kikuchipy/draw/__init__.py
@@ -17,7 +17,10 @@
 
 """Tools for use in plotting of signals."""
 
-__all__ = ["get_rgb_navigator"]
+__all__ = [
+    "get_rgb_navigator",
+    "plot_pattern_positions_in_map",
+]
 
 
 def __dir__():
@@ -27,6 +30,7 @@ def __dir__():
 def __getattr__(name):
     _import_mapping = {
         "get_rgb_navigator": "_navigators",
+        "plot_pattern_positions_in_map": "_plot_pattern_positions_in_map",
     }
     if name in __all__:
         import importlib
diff --git a/kikuchipy/draw/_plot_pattern_positions_in_map.py b/kikuchipy/draw/_plot_pattern_positions_in_map.py
new file mode 100644
index 00000000..8289d401
--- /dev/null
+++ b/kikuchipy/draw/_plot_pattern_positions_in_map.py
@@ -0,0 +1,138 @@
+# Copyright 2019-2023 The kikuchipy developers
+#
+# This file is part of kikuchipy.
+#
+# kikuchipy is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# kikuchipy is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# 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 plotting positions of patterns in a 2D map.
+
+These functions were created specifically to plot positions of NORDIF
+calibration patterns. These are acquired manually within a region of
+interest within a larger area.
+"""
+
+from typing import Optional
+
+import matplotlib.patches as mpatches
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+def plot_pattern_positions_in_map(
+    rc: np.ndarray,
+    roi_shape: tuple,
+    roi_origin: tuple = (0, 0),
+    area_shape: Optional[tuple] = None,
+    roi_image: Optional[np.ndarray] = None,
+    area_image: Optional[np.ndarray] = None,
+    axis: Optional[plt.Axes] = None,
+    return_figure: bool = False,
+    color: Optional[str] = "k",
+) -> Optional[plt.Figure]:
+    """Plot pattern positions in a 2D map within a region of interest
+    (ROI), the ROI potentially within a larger area.
+
+    Parameters
+    ----------
+    rc
+        Position coordinates (row, column) in an array of shape (n, 2).
+        If ``area_shape`` is passed, coordinates are assumed to be given
+        with respect to the area origin, and so if ``roi_origin`` is
+        passed, the origin is subtracted from the coordinates.
+    roi_shape
+        Shape of the ROI as (n rows, n columns).
+    roi_origin
+        Origin (row, column) of the ROI with respect to the area. If
+        this and ``area_shape`` is passed, the origin is subtracted
+        from the ``rc`` coordinates.
+    area_shape
+        Shape of the area including the ROI as (n rows, n columns). If
+        this and ``roi_origin`` is passed, the origin is subtracted from
+        the ``rc`` coordinates.
+    roi_image
+        Image to plot within the ROI, of the same aspect ratio.
+    area_image
+        Image to plot within the area, of the same aspect ratio.
+    axis
+        Existing Matplotlib axis to add the positions to. If not passed,
+        a new figure will be created. If passed, only the coordinate
+        markers and labels are added to the axis. E.g. ``roi_image`` or
+        ``area_image`` will not be used.
+    return_figure
+        Whether to return the created figure. Default is ``False``.
+    color
+        Color of position markers and labels. Default is ``k``. Must be
+        a valid Matplotlib color identifier.
+
+    Returns
+    -------
+    fig
+        Created figure, returned if ``return_figure=True``.
+    """
+    # Set different styling of rectangle(s), markers and text. Default
+    # values assume no area_shape or images are passed.
+    roi_rect_kw = dict(fc="none", lw=2, clip_on=False, zorder=5)
+    roi_ny, roi_nx = roi_shape
+
+    if isinstance(axis, plt.Axes):
+        new_axis = False
+        ax = axis
+        fig = ax.figure
+    else:
+        new_axis = True
+        fig, ax = plt.subplots()
+
+    if new_axis and area_shape:
+        area_ny, area_nx = area_shape
+        roi_rect_kw.update(dict(ec="r", fc="none"))
+        if isinstance(area_image, np.ndarray):
+            area_image_flipped = area_image[::-1, :]  # Inverse y-axis
+            ax.imshow(area_image_flipped, extent=(0, area_nx, 0, area_ny))
+        ax.text(0, 0, s="Area", va="bottom")
+        xlim = (0, area_nx)
+        ylim = (0, area_ny)
+    else:
+        if roi_origin != (0, 0):
+            rc = rc.copy()
+            rc -= np.array(roi_origin)
+            roi_origin = (0, 0)
+        xlim = (0, roi_nx)
+        ylim = (0, roi_ny)
+
+    if new_axis:
+        if isinstance(roi_image, np.ndarray):
+            extent = (
+                roi_origin[1],
+                roi_origin[1] + roi_nx,
+                roi_origin[0] + roi_ny,
+                roi_origin[0],
+            )
+            ax.imshow(roi_image, extent=extent, zorder=1)
+            roi_rect_kw.update(dict(fc="none"))
+
+        rect_roi = mpatches.Rectangle(roi_origin[::-1], roi_nx, roi_ny, **roi_rect_kw)
+        if area_shape:
+            ax.text(*rect_roi.xy, s="Region of interest", va="bottom")
+            ax.add_artist(rect_roi)
+
+    for i, (r, c) in enumerate(rc):
+        ax.scatter(c, r, s=50, marker="+", c=color)
+        ax.text(c, r, str(i), va="bottom", ha="left", c=color)
+
+    if new_axis:
+        ax.set(xlim=xlim, ylim=ylim, xlabel="Column", ylabel="Row", aspect="equal")
+        ax.invert_yaxis()
+
+    if return_figure:
+        return fig
diff --git a/kikuchipy/indexing/_refinement/_refinement.py b/kikuchipy/indexing/_refinement/_refinement.py
index b2f1e8f9..faa9c846 100644
--- a/kikuchipy/indexing/_refinement/_refinement.py
+++ b/kikuchipy/indexing/_refinement/_refinement.py
@@ -75,8 +75,10 @@ def compute_refine_orientation_results(
         time_start = time()
         computed_results = results.compute().reshape((-1, 4))
         total_time = time() - time_start
-        patterns_per_second = int(np.floor(n_patterns / total_time))
-        print(f"Refinement speed: {patterns_per_second} patterns/s", file=sys.stdout)
+        patterns_per_second = n_patterns / total_time
+        print(
+            f"Refinement speed: {patterns_per_second:.5f} patterns/s", file=sys.stdout
+        )
         # (n, score, phi1, Phi, phi2)
         computed_results = np.array(computed_results)
         xmap_refined = CrystalMap(
@@ -123,8 +125,10 @@ def compute_refine_projection_center_results(
         time_start = time()
         computed_results = results.compute().reshape((-1, 4))
         total_time = time() - time_start
-        patterns_per_second = int(np.floor(n_patterns / total_time))
-        print(f"Refinement speed: {patterns_per_second} patterns/s", file=sys.stdout)
+        patterns_per_second = n_patterns / total_time
+        print(
+            f"Refinement speed: {patterns_per_second:.5f} patterns/s", file=sys.stdout
+        )
         # (n, score, PCx, PCy, PCz)
         computed_results = np.array(computed_results)
         new_detector = detector.deepcopy()
@@ -180,8 +184,10 @@ def compute_refine_orientation_projection_center_results(
         time_start = time()
         computed_results = results.compute().reshape((-1, 7))
         total_time = time() - time_start
-        patterns_per_second = int(np.floor(n_patterns / total_time))
-        print(f"Refinement speed: {patterns_per_second} patterns/s", file=sys.stdout)
+        patterns_per_second = n_patterns / total_time
+        print(
+            f"Refinement speed: {patterns_per_second:.5f} patterns/s", file=sys.stdout
+        )
         computed_results = np.array(computed_results)
         # (n, score, phi1, Phi, phi2, PCx, PCy, PCz)
         xmap_refined = CrystalMap(
@@ -827,7 +833,7 @@ def __init__(
         self.rotations_array = da.from_array(rot, chunks=self.chunks)
 
         # Relevant data from the detector
-        self.unique_pc = np.prod(detector.navigation_shape) != 1 and self.nav_size > 1
+        self.unique_pc = detector.navigation_size != 1 and self.nav_size > 1
         dtype = np.float64
         pc_shape = self.nav_shape + (3,)
         if self.unique_pc:
@@ -841,7 +847,10 @@ def __init__(
         self.pc_array = da.from_array(pc, chunks=self.chunks)
 
         if mode == "ori_pc":
-            self.rotations_pc_array = da.dstack((self.rotations_array, self.pc_array))
+            nav_ndim = len(self.nav_shape)
+            self.rotations_pc_array = da.concatenate(
+                [self.rotations_array, self.pc_array], axis=nav_ndim
+            )
 
         # Keyword arguments passed to Dask when iterating over chunks
         self.map_blocks_kwargs.update(
@@ -1045,8 +1054,17 @@ def get_bound_constraints(
             upper_bounds = lower_bounds.copy()
             return lower_bounds, upper_bounds
 
-        eu_lower = 3 * [0]
-        eu_upper = [2 * np.pi, np.pi, 2 * np.pi]
+        # Absolute constraints for Euler angles and PCs. Note that
+        # constraints on the angles need some leeway since NLopt can
+        # throw a RuntimeError when a starting Euler angle is too close
+        # to the constraint.
+        angle_leeway = np.deg2rad(5)
+        eu_lower = 3 * [-angle_leeway]
+        eu_upper = [
+            2 * np.pi + angle_leeway,
+            np.pi + angle_leeway,
+            2 * np.pi + angle_leeway,
+        ]
         pc_lower = 3 * [-2]
         pc_upper = 3 * [2]
 


Test patch:
diff --git a/kikuchipy/_util/tests/test_import.py b/kikuchipy/_util/tests/test_import.py
index d8da5183..e300bd50 100644
--- a/kikuchipy/_util/tests/test_import.py
+++ b/kikuchipy/_util/tests/test_import.py
@@ -224,7 +224,10 @@ def test_dir_detectors(self):
     def test_dir_draw(self):
         import kikuchipy.draw
 
-        assert dir(kikuchipy.draw) == ["get_rgb_navigator"]
+        assert dir(kikuchipy.draw) == [
+            "get_rgb_navigator",
+            "plot_pattern_positions_in_map",
+        ]
 
     def test_dir_filters(self):
         import kikuchipy.filters
diff --git a/kikuchipy/conftest.py b/kikuchipy/conftest.py
index 2ef6cae8..2b18d009 100644
--- a/kikuchipy/conftest.py
+++ b/kikuchipy/conftest.py
@@ -146,8 +146,8 @@ def dummy_background():
 def ebsd_with_axes_and_random_data(request):
     """EBSD signal with minimally defined axes and random data.
 
-    Parameters
-    ----------
+    Parameters expected in `request`
+    -------------------------------
     navigation_shape : tuple
     signal_shape : tuple
     lazy : bool
@@ -239,7 +239,7 @@ def get_single_phase_xmap(rotations):
     def _get_single_phase_xmap(
         nav_shape,
         rotations_per_point=5,
-        prop_names=["scores", "simulation_indices"],
+        prop_names=("scores", "simulation_indices"),
         name="a",
         phase_id=0,
         step_sizes=None,
@@ -448,7 +448,6 @@ def ebsd_directory(tmpdir, request):
 
 @pytest.fixture(autouse=True)
 def doctest_setup_teardown(request):
-    # Setup
     # Temporarily turn off interactive plotting with Matplotlib
     plt.ioff()
 
diff --git a/kikuchipy/detectors/tests/test_ebsd_detector.py b/kikuchipy/detectors/tests/test_ebsd_detector.py
index c1643bea..d0b84843 100644
--- a/kikuchipy/detectors/tests/test_ebsd_detector.py
+++ b/kikuchipy/detectors/tests/test_ebsd_detector.py
@@ -18,6 +18,7 @@
 from copy import deepcopy
 
 import matplotlib
+import matplotlib.collections as mcollections
 import matplotlib.colors as mcolors
 import matplotlib.pyplot as plt
 import numpy as np
@@ -151,7 +152,9 @@ def test_set_pc_coordinates(self, pc1):
     def test_pc_average(self, pc, desired_pc_average):
         """Calculation of PC average."""
         assert np.allclose(
-            kp.detectors.EBSDDetector(pc=pc).pc_average, desired_pc_average
+            kp.detectors.EBSDDetector(pc=pc).pc_average,
+            desired_pc_average,
+            atol=1e-3,
         )
 
     @pytest.mark.parametrize(
@@ -585,3 +588,415 @@ def test_crop_simulated(self):
         sim2 = mp.get_patterns(rot, det2, **kwds)
 
         assert np.allclose(sim2.data, sim1.isig[20:60, :50].data)
+
+
+class TestPlotPC:
+    det = kp.detectors.EBSDDetector(
+        shape=(60, 60),
+        pc=np.stack(
+            (
+                np.repeat(np.linspace(0.55, 0.45, 30), 20).reshape(30, 20).T,
+                np.repeat(np.linspace(0.75, 0.70, 20), 30).reshape(20, 30),
+                np.repeat(np.linspace(0.50, 0.55, 20), 30).reshape(20, 30),
+            ),
+            axis=2,
+        ),
+        sample_tilt=70,
+    )
+
+    def test_plot_pc_raises(self):
+        det = self.det.deepcopy()
+        det.pc = det.pc_average
+        with pytest.raises(ValueError, match="Detector must have more than one "):
+            det.plot_pc()
+
+        det2 = self.det.deepcopy()
+        det2.pc = det2.pc[0]
+        with pytest.raises(ValueError, match="Detector's `navigation_dimension` must "):
+            det2.plot_pc()
+
+        with pytest.raises(ValueError, match="Plot mode 'stereographic' must be one "):
+            self.det.plot_pc("stereographic")
+
+    def test_plot_pc_map_horizontal(self):
+        fig = self.det.plot_pc(return_figure=True)
+
+        figsize = fig.get_size_inches()
+        assert (figsize[0] / figsize[1]) > 1
+
+        ax = fig.axes
+        assert len(ax) == 6
+        assert all([a.get_xlabel() == "Column" for a in ax[:3]])
+        assert all(
+            [a.get_ylabel() == f"PC{l}" for a, l in zip(ax[3:], ["x", "y", "z"])]
+        )
+
+        plt.close(fig)
+
+    def test_plot_pc_map_vertical(self):
+        fig = self.det.plot_pc(return_figure=True, orientation="vertical")
+
+        figsize = fig.get_size_inches()
+        assert (figsize[0] / figsize[1]) < 1
+
+        ax = fig.axes
+        assert len(ax) == 6
+        assert all([a.get_xlabel() == "Column" for a in ax[:3]])
+        assert all(
+            [a.get_ylabel() == f"PC{l}" for a, l in zip(ax[3:], ["x", "y", "z"])]
+        )
+
+        plt.close(fig)
+
+    def test_plot_pc_scatter_horizontal(self):
+        fig = self.det.plot_pc("scatter", return_figure=True, annotate=True)
+
+        figsize = fig.get_size_inches()
+        assert (figsize[0] / figsize[1]) > 1
+
+        ax = fig.axes
+        assert len(ax) == 3
+        assert all(
+            [a.get_xlabel() == f"PC{l}" for a, l in zip(ax[3:], ["x", "x", "z"])]
+        )
+        assert all(
+            [a.get_ylabel() == f"PC{l}" for a, l in zip(ax[3:], ["y", "z", "y"])]
+        )
+
+        texts = ax[0].texts
+        assert len(texts) == self.det.navigation_size
+        assert texts[0].get_text() == "0"
+        assert texts[-1].get_text() == "599"
+
+        plt.close(fig)
+
+    def test_plot_pc_scatter_vertical(self):
+        fig = self.det.plot_pc("scatter", return_figure=True, orientation="vertical")
+
+        figsize = fig.get_size_inches()
+        assert (figsize[0] / figsize[1]) < 1
+
+        ax = fig.axes
+        assert len(ax) == 3
+        assert all(
+            [a.get_xlabel() == f"PC{l}" for a, l in zip(ax[3:], ["x", "x", "z"])]
+        )
+        assert all(
+            [a.get_ylabel() == f"PC{l}" for a, l in zip(ax[3:], ["y", "z", "y"])]
+        )
+
+        plt.close(fig)
+
+    def test_plot_pc_3d(self):
+        fig = self.det.plot_pc("3d", return_figure=True, annotate=True)
+
+        texts = fig.axes[0].texts
+        assert len(texts) == self.det.navigation_size
+        assert texts[0].get_text() == "0"
+        assert texts[-1].get_text() == "599"
+
+        plt.close(fig)
+
+    def test_plot_pc_figure(self):
+        fig1 = self.det.plot_pc(figure_kwargs=dict(figsize=(9, 3)), return_figure=True)
+        assert fig1.get_tight_layout()
+
+        fig2 = self.det.plot_pc(
+            figure_kwargs=dict(figsize=(6, 3), layout="constrained"), return_figure=True
+        )
+        assert fig2.get_constrained_layout()
+        assert not np.allclose(fig1.get_size_inches(), fig2.get_size_inches())
+
+        plt.close("all")
+
+
+class TestEstimateTilts:
+    det0 = kp.detectors.EBSDDetector(
+        shape=(480, 480),
+        pc=(0.5, 0.3, 0.5),
+        sample_tilt=70,
+        tilt=0,
+        px_size=70,
+    )
+
+    def test_estimate_xtilt_raises(self):
+        with pytest.raises(ValueError, match="Estimation requires more than one "):
+            _ = self.det0.estimate_xtilt()
+
+    def test_estimate_xtilt(self):
+        det = self.det0.extrapolate_pc(
+            pc_indices=[0, 0],
+            navigation_shape=(15, 20),
+            step_sizes=(50, 50),
+        )
+
+        xtilt = det.estimate_xtilt(degrees=True)
+        assert np.isclose(xtilt, 90 - self.det0.sample_tilt + self.det0.tilt)
+        assert plt.get_fignums() == [1]
+        ax = plt.gca()
+        assert not any(["Outliers" in t.get_text() for t in ax.get_legend().texts])
+
+        xtilt, is_outliers, fig = det.estimate_xtilt(
+            return_outliers=True, return_figure=True
+        )
+        assert isinstance(is_outliers, np.ndarray)
+        assert is_outliers.sum() == 0
+
+        plt.close("all")
+
+    def test_estimate_xtilt_outliers(self):
+        det = self.det0.extrapolate_pc(
+            pc_indices=[0, 0],
+            navigation_shape=(15, 20),
+            step_sizes=(50, 50),
+        )
+        det.pc[0, 0] = (0.5, 0.5, 0.5)
+
+        xtilt, is_outliers, fig = det.estimate_xtilt(
+            return_outliers=True, return_figure=True
+        )
+        assert isinstance(is_outliers, np.ndarray)
+        assert is_outliers.sum() == 1
+        assert np.allclose(np.where(is_outliers)[0], [0, 0])
+        assert any(["Outliers" in t.get_text() for t in fig.axes[0].get_legend().texts])
+
+    def test_estimate_xtilt_ztilt(self):
+        det1 = self.det0.extrapolate_pc(
+            pc_indices=[0, 0],
+            navigation_shape=(15, 20),
+            step_sizes=(1, 1),
+        )
+        xtilt, ztilt = det1.estimate_xtilt_ztilt(degrees=True)
+        assert np.isclose(xtilt, 20)
+        assert np.isclose(ztilt, 0)
+
+        # Add outliers
+        det2 = det1.deepcopy()
+        outlier_idx = [[0, 0], [0, 10]]
+        det2.pc[tuple(outlier_idx)] = (0.5, 0.4, 0.5)
+
+        np.random.seed(42)
+
+        xtilt2, ztilt2 = det2.estimate_xtilt_ztilt(degrees=True)
+        assert np.isclose(xtilt2, 0.169, atol=1e-3)
+        assert np.isclose(ztilt2, -74.339, atol=1e-3)
+
+        is_outlier = np.ravel_multi_index(outlier_idx, det1.navigation_shape)
+        xtilt3, ztilt3 = det2.estimate_xtilt_ztilt(degrees=True, is_outlier=is_outlier)
+        assert np.isclose(xtilt3, 20)
+        assert np.isclose(ztilt3, 0)
+
+    def test_estimate_xtilt_ztilt_raises(self):
+        with pytest.raises(ValueError, match="Estimation requires more than one "):
+            _ = self.det0.estimate_xtilt_ztilt()
+
+
+class TestExtrapolatePC:
+    det0 = kp.detectors.EBSDDetector(
+        shape=(240, 240),
+        pc=(0.5, 0.3, 0.5),
+        sample_tilt=70,
+        tilt=0,
+        px_size=70,
+        binning=2,
+    )
+
+    def test_extrapolate_pc(self):
+        det = self.det0.extrapolate_pc(
+            pc_indices=[7, 15],
+            navigation_shape=(15, 31),
+            step_sizes=(50, 50),
+        )
+        assert det.navigation_shape == (15, 31)
+        assert np.allclose(
+            [
+                self.det0.nrows,
+                self.det0.ncols,
+                self.det0.sample_tilt,
+                self.det0.tilt,
+                self.det0.px_size,
+                self.det0.binning,
+                self.det0.azimuthal,
+            ],
+            [
+                det.nrows,
+                det.ncols,
+                det.sample_tilt,
+                det.tilt,
+                det.px_size,
+                det.binning,
+                det.azimuthal,
+            ],
+        )
+        assert np.allclose(det.pc_average, self.det0.pc)
+        assert np.allclose(det.pc_flattened.min(0), [0.4777, 0.2902, 0.4964], atol=1e-4)
+        assert np.allclose(det.pc_flattened.max(0), [0.5223, 0.3098, 0.5036], atol=1e-4)
+
+    def test_extrapolate_pc_multiple_indices(self):
+        det1 = self.det0.deepcopy()
+        # Specify PC values in four corners of the map, visualized here
+        # as they would show up in the (PCx, PCy) scatter plot
+        # fmt: off
+        det1.pc = [
+            [0.5, 0.3, 0.5],                [0.3, 0.3, 0.5],
+
+
+
+
+            [0.5, 0.2, 0.6],                [0.3, 0.2, 0.6],
+        ]
+        # fmt: on
+
+        det2 = det1.extrapolate_pc(
+            pc_indices=[[0, 0], [0, 10], [20, 0], [20, 10]],
+            navigation_shape=(11, 21),
+            step_sizes=(11, 11),
+        )
+        assert np.allclose(det1.pc_average, det2.pc_average, atol=1e-2)
+
+        det3 = det1.extrapolate_pc(
+            pc_indices=[[0, 0], [0, 10], [20, 0], [20, 10]],
+            navigation_shape=(11, 21),
+            step_sizes=(5, 5),
+            shape=(60, 60),
+            binning=8,
+            px_size=70 * 8,
+        )
+        assert np.allclose(det1.pc_average, det3.pc_average, atol=1e-2)
+        assert det3.shape == (60, 60)
+        assert det3.binning == 8
+        assert det3.px_size == 70 * 8
+
+    def test_extrapolate_pc_outliers(self):
+        det1 = self.det0.deepcopy()
+        det1.pc = [[0.5, 0.3, 0.5], [0.3, 0.3, 0.5], [0.5, 0.2, 0.6], [0.3, 0.2, 0.6]]
+        det2 = det1.extrapolate_pc(
+            pc_indices=[[0, 0], [0, 10], [20, 0], [20, 10]],
+            navigation_shape=(11, 21),
+            step_sizes=(11, 11),
+            is_outlier=[True, False, False, False],
+        )
+        assert np.allclose(det2.pc_average, [0.366, 0.236, 0.566], atol=1e-3)
+
+
+class TestFitPC:
+    def setup_method(self):
+        """Create a plane of PCs with a known mean, add some noise,
+        extract selected patterns, and try to 'reconstruct' the plane
+        by fitting.
+        """
+        det0 = kp.detectors.EBSDDetector(
+            shape=(240, 240),
+            pc=(0.5, 0.3, 0.5),
+            sample_tilt=70,
+        )
+
+        det = det0.extrapolate_pc(
+            pc_indices=[7, 15], navigation_shape=(15, 31), step_sizes=(50, 50)
+        )
+
+        # Add noise
+        rng = np.random.default_rng(42)
+        v = 0.005
+        det.pcy += rng.uniform(-v, v, det.navigation_size).reshape(det.navigation_shape)
+        det.pcz += rng.uniform(-v, v, det.navigation_size).reshape(det.navigation_shape)
+
+        self.det0 = det0
+        self.det = det
+        self.map_indices = np.stack(np.indices(det.navigation_shape))
+
+    def test_fit_pc_corner_patterns(self):
+        """Test projective fit."""
+        pc_indices = [[0, 0, 14, 14], [0, 30, 0, 30]]
+        det2 = self.det.deepcopy()
+        det2.pc = det2.pc[tuple(pc_indices)].reshape((2, 2, 3))
+        pc_indices = np.array(pc_indices).reshape((2, 2, 2))
+        det_fit, fig = det2.fit_pc(
+            pc_indices=pc_indices, map_indices=self.map_indices, return_figure=True
+        )
+        assert np.all(abs(det_fit.pc_flattened - self.det.pc_flattened).max(0) < 0.009)
+
+        # We have a plane in the 3D plot
+        assert isinstance(fig.axes[3].collections[2], mcollections.PolyCollection)
+
+        plt.close("all")
+
+    @pytest.mark.parametrize(
+        "grid_shape, max_error", [((3, 3), 0.009), ((5, 5), 0.0091), ((7, 7), 0.0064)]
+    )
+    def test_fit_pc_grid_patterns_33(self, grid_shape, max_error):
+        """Test projective fit."""
+        pc_indices = kp.signals.util.grid_indices(grid_shape, self.det.navigation_shape)
+        pc_indices = pc_indices.reshape(2, -1)
+        det2 = self.det.deepcopy()
+        det2.pc = det2.pc[tuple(pc_indices)]
+        det_fit = det2.fit_pc(
+            pc_indices=pc_indices, map_indices=self.map_indices, plot=False
+        )
+        assert np.all(
+            abs(det_fit.pc_flattened - self.det.pc_flattened).max(0) < max_error
+        )
+
+    def test_fit_pc_affine_outliers(self):
+        grid_shape = (7, 7)
+        pc_indices = kp.signals.util.grid_indices(grid_shape, self.det.navigation_shape)
+        pc_indices = pc_indices.reshape(2, -1)
+        det2 = self.det.deepcopy()
+        det2.pc = det2.pc[tuple(pc_indices)]
+
+        # Add outliers to extracted PCs
+        det2.pc = np.append(det2.pc, [[0.55, 0.15, 0.55], [0.6, 0.10, 0.6]], axis=0)
+        is_outlier = np.zeros(det2.navigation_size, dtype=bool)
+        is_outlier[[-2, -1]] = True
+        pc_indices = np.append(pc_indices, [[1, 1], [1, 2]], axis=1)
+
+        # Bad fit
+        det_fit1 = det2.fit_pc(
+            pc_indices=pc_indices,
+            map_indices=self.map_indices,
+            transformation="affine",
+        )
+        assert np.allclose(
+            abs(det_fit1.pc_flattened - self.det.pc_flattened).max(0),
+            [0.70, 0.35, 0.13],
+            atol=1e-2,
+        )
+
+        # Good fit
+        det_fit2, fig = det2.fit_pc(
+            pc_indices=pc_indices,
+            map_indices=self.map_indices,
+            transformation="affine",
+            is_outlier=is_outlier,
+            return_figure=True,
+        )
+        assert np.all(abs(det_fit2.pc_flattened - self.det.pc_flattened).max(0) < 0.009)
+        assert isinstance(fig, plt.Figure)
+        assert len(fig.axes) == 4
+
+        plt.close("all")
+
+    def test_fit_pc_raises(self):
+        grid_shape = (7, 5)
+        pc_indices = kp.signals.util.grid_indices(grid_shape, self.det.navigation_shape)
+        pc_indices = pc_indices.reshape(2, -1)
+        det2 = self.det.deepcopy()
+        det2.pc = det2.pc[tuple(pc_indices)]
+
+        with pytest.raises(ValueError, match="Fitting requires multiple projection "):
+            _ = self.det0.fit_pc(pc_indices, self.map_indices)
+
+        with pytest.raises(ValueError, match=r"`pc_indices` array shape \(2, 7, 5\) "):
+            _ = det2.fit_pc(pc_indices.reshape((2,) + grid_shape), self.map_indices)
+
+        with pytest.raises(ValueError, match=r"`map_indices` array shape \(930,\) "):
+            _ = det2.fit_pc(pc_indices, self.map_indices.ravel())
+
+        det2.pc = det2.pc.reshape(grid_shape + (3,))
+        with pytest.raises(ValueError, match=r"`pc_indices` array shape \(2, 35\) "):
+            _ = det2.fit_pc(pc_indices, self.map_indices)
+        det2.pc = det2.pc.reshape(-1, 3)
+
+        is_outlier = np.zeros(det2.navigation_size - 1, dtype=bool)
+        with pytest.raises(ValueError, match="`is_outlier` must be a boolean array of"):
+            _ = det2.fit_pc(pc_indices, self.map_indices, is_outlier=is_outlier)
diff --git a/kikuchipy/draw/tests/test_plot_pattern_positions_in_map.py b/kikuchipy/draw/tests/test_plot_pattern_positions_in_map.py
new file mode 100644
index 00000000..3ae8b6a1
--- /dev/null
+++ b/kikuchipy/draw/tests/test_plot_pattern_positions_in_map.py
@@ -0,0 +1,153 @@
+# Copyright 2019-2023 The kikuchipy developers
+#
+# This file is part of kikuchipy.
+#
+# kikuchipy is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# kikuchipy is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with kikuchipy. If not, see <http://www.gnu.org/licenses/>.
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import kikuchipy as kp
+
+
+class TestPlotPatternPositionsInMap:
+    # Original metadata similar to the original metadata returned with
+    # the Ni gain 0 calibration pattern dataset
+    omd = {
+        "roi": {
+            "origin": (185, 73),
+            "origin_scaled": (49, 19),
+            "shape": (558, 750),
+            "shape_scaled": (149, 200),
+        },
+        "area": {"shape": (1000, 1000), "shape_scaled": (267, 267)},
+        "area_image": np.ones((1000, 1000), dtype="uint8"),
+        "calibration_patterns": {
+            "indices": np.array(
+                [
+                    [447, 425],
+                    [532, 294],
+                    [543, 573],
+                    [378, 596],
+                    [369, 308],
+                    [632, 171],
+                    [668, 704],
+                    [269, 696],
+                    [247, 152],
+                ]
+            ),
+            "indices_scaled": np.array(
+                [
+                    [119, 113],
+                    [142, 78],
+                    [145, 153],
+                    [101, 159],
+                    [98, 82],
+                    [169, 46],
+                    [178, 188],
+                    [72, 186],
+                    [66, 41],
+                ]
+            ),
+        },
+    }
+
+    def test_plot_pos_no_area(self):
+        rc = self.omd["calibration_patterns"]["indices_scaled"]
+        rc_shifted = rc - np.array(self.omd["roi"]["origin_scaled"])
+
+        fig = kp.draw.plot_pattern_positions_in_map(
+            rc=rc,
+            roi_shape=self.omd["roi"]["shape_scaled"],
+            roi_origin=self.omd["roi"]["origin_scaled"],
+            return_figure=True,
+        )
+        ax = fig.axes[0]
+
+        assert ax.get_xlabel() == "Column"
+        assert ax.get_ylabel() == "Row"
+        assert np.allclose(ax.get_xlim(), (0, 200))
+        assert np.allclose(ax.get_ylim(), (149, 0))  # Checks inversion
+
+        coll = ax.collections
+        assert len(coll) == len(rc)
+        # Ensure it's a cross
+        assert np.allclose(
+            coll[0].get_paths()[0].vertices, [[-0.5, 0], [0.5, 0], [0, -0.5], [0, 0.5]]
+        )
+
+        texts = ax.texts
+        assert len(texts) == 9
+        for i, t in enumerate(texts):
+            assert int(t.get_text()) == i
+            assert np.allclose(t.get_position(), rc_shifted[i][::-1])
+            assert t.get_color() == "k"
+
+        plt.close(fig)
+
+    def test_plot_pos_area(self):
+        rc = self.omd["calibration_patterns"]["indices"]
+
+        area_shape = self.omd["area"]["shape"]
+        area_size = int(np.prod(area_shape))
+
+        fig = kp.draw.plot_pattern_positions_in_map(
+            rc=rc,
+            roi_shape=self.omd["roi"]["shape"],
+            roi_origin=self.omd["roi"]["origin"],
+            area_shape=area_shape,
+            area_image=np.random.random(area_size).reshape(area_shape),
+            roi_image=np.zeros(self.omd["roi"]["shape"]),
+            return_figure=True,
+            color="w",
+        )
+        ax = fig.axes[0]
+
+        assert np.allclose(ax.get_xlim(), (0, 1000))
+        assert np.allclose(ax.get_ylim(), (1000, 0))
+
+        images = ax.images
+        assert images[0].get_array().shape == area_shape
+        assert images[1].get_array().shape == self.omd["roi"]["shape"]
+
+        # Marker labels plus "Area" and "Region of interest"
+        texts = ax.texts
+        assert len(texts) == 11
+        assert texts[0].get_color() == "black"
+        assert texts[2].get_color() == "w"
+        assert texts[0].get_text() == "Area"
+        assert texts[1].get_text() == "Region of interest"
+
+        # Red rectangle for the ROI
+        rectangle = ax.patches[0]
+        assert np.allclose(rectangle.get_edgecolor(), [1, 0, 0, 1])
+        assert rectangle.get_xy() == self.omd["roi"]["origin"][::-1]
+
+        plt.close(fig)
+
+    def test_plot_pos_existing_axis(self):
+        fig = plt.figure()
+        ax = fig.add_subplot()
+
+        kp.draw.plot_pattern_positions_in_map(
+            rc=self.omd["calibration_patterns"]["indices_scaled"],
+            roi_shape=self.omd["roi"]["shape_scaled"],
+            roi_origin=self.omd["roi"]["origin_scaled"],
+            axis=ax,
+        )
+
+        assert len(ax.collections) == len(self.omd["calibration_patterns"]["indices"])
+        assert len(ax.texts) == 9
+
+        plt.close(fig)