Orientational sites

src/gemdat/orientations.py

@@ -6,6 +6,7 @@
 from pymatgen.symmetry.groups import PointGroup
 
 from gemdat.trajectory import Trajectory
+from gemdat.volume import OrientationalVolume
 from gemdat.utils import fft_autocorrelation, cartesian_to_spherical
 
 
@@ -280,37 +281,28 @@ def plot_autocorrelation(self, **kwargs):
 
         return plots.autocorrelation(orientations=self, **kwargs)
 
+    def to_volume(
+        self,
+        shape: tuple[int, int] = (90, 360),
+        normalize_area: bool = False,
+    ) -> OrientationalVolume:
+        """Calculate density volume from orientation.
 
-def calculate_spherical_areas(shape: tuple[int, int], radius: float = 1) -> np.ndarray:
-    """Calculate the areas of a section of a sphere, defined in spherical
-    coordinates. Useful for normalization purposes.
+        The orientations are converted in spherical coordinates,
+        and the azimutal and elevation angles are binned into a 2d histogram.
 
-    Parameters
-    ----------
-    shape : tuple
-        Shape of the grid in integer degrees
-    radius : float
-        Radius of the sphere
-
-    Returns
-    -------
-    areas : np.ndarray
-        Areas of the section
-    """
-    elevation_angles = np.linspace(0, 180, shape[0])
-
-    areas = np.zeros(shape, dtype=float)
-    azimuthal_increment = np.deg2rad(1)
-    elevation_increment = np.deg2rad(1)
-
-    for i in range(shape[1]):
-        for j in range(shape[0]):
-            areas[j, i] = (
-                (radius**2)
-                * azimuthal_increment
-                * np.sin(np.deg2rad(elevation_angles[j]))
-                * elevation_increment
-            )
-            # hacky way to get rid of singularity on poles
-            areas[0, :] = areas[-1, 0]
-    return areas
+        Parameters
+        ----------
+        shape : tuple
+            The shape of the spherical sector in which the trajectory is
+        normalize_area : bool
+            If True, normalize the histogram by the area of the bins
+
+        Returns
+        -------
+        vol : OrientationalVolume
+            Output volume
+        """
+        from gemdat.volume import orientations_to_volume
+
+        return orientations_to_volume(self, shape=shape, normalize_area=normalize_area)

src/gemdat/plots/matplotlib/_rectilinear.py

@@ -3,17 +3,16 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from gemdat.orientations import (
-    Orientations,
-    calculate_spherical_areas,
-)
+from gemdat.orientations import Orientations
 
 
 def rectilinear(
     *,
     orientations: Orientations,
     shape: tuple[int, int] = (90, 360),
     normalize_histo: bool = True,
+    add_peaks: bool = False,
+    **kwargs,
 ) -> plt.Figure:
     """Plot a rectilinear projection of a spherical function. This function
     uses the transformed trajectory.
@@ -26,34 +25,21 @@ def rectilinear(
         The shape of the spherical sector in which the trajectory is plotted
     normalize_histo : bool, optional
         If True, normalize the histogram by the area of the bins, by default True
+    add_peaks : bool, optional
+        If True, plot the peaks of the histogram
+    **kwargs : dict
+        Additional keyword arguments for the `orientational_peaks` method
 
     Returns
     -------
     fig : matplotlib.figure.Figure
         Output figure
     """
-    az, el, _ = orientations.vectors_spherical.T
-    az = az.flatten()
-    el = el.flatten()
-
-    hist, *_ = np.histogram2d(el, az, shape)
-
-    if normalize_histo:
-        areas = calculate_spherical_areas(shape)
-        hist = hist / areas
-        # Drop the bins at the poles where normalization is not possible
-        hist = hist[1:-1, :]
-
-    axis_theta, axis_phi = hist.shape
-
-    phi = np.linspace(0, 360, axis_phi)
-    theta = np.linspace(0, 180, axis_theta)
-
-    x, y = np.meshgrid(phi, theta)
+    ov = orientations.to_volume(shape=shape, normalize_area=normalize_histo)
+    theta, phi, values = ov.data.T
 
     fig, ax = plt.subplots()
-
-    cs = ax.contourf(x, y, hist)
+    cs = ax.contourf(phi, theta, values)
 
     ax.set_xticks(np.linspace(0, 360, 9))
     ax.set_yticks(np.linspace(0, 180, 5))
@@ -63,4 +49,9 @@ def rectilinear(
 
     fig.colorbar(cs, label='Areal probability', format='')
 
+    if add_peaks:
+        peaks = ov.orientational_peaks(**kwargs)
+        xp, yp = zip(*peaks)
+        ax.plot(xp, yp, 'ro', markersize=5)
+
     return fig

src/gemdat/plots/plotly/_rectilinear.py

@@ -4,10 +4,8 @@
 
 import numpy as np
 
-from gemdat.orientations import (
-    Orientations,
-    calculate_spherical_areas,
-)
+from gemdat.orientations import Orientations
+from gemdat.utils import calculate_spherical_areas
 
 
 def rectilinear(

src/gemdat/utils.py

@@ -309,3 +309,38 @@ def fft_autocorrelation(coords: np.ndarray) -> np.ndarray:
     autocorrelation = autocorrelation / autocorrelation[:, 0, np.newaxis]
 
     return autocorrelation
+
+
+def calculate_spherical_areas(shape: tuple[int, int], radius: float = 1) -> np.ndarray:
+    """Calculate the areas of a section of a sphere, defined in spherical
+    coordinates. Useful for normalization purposes.
+
+    Parameters
+    ----------
+    shape : tuple
+        Shape of the grid in integer degrees
+    radius : float
+        Radius of the sphere
+
+    Returns
+    -------
+    areas : np.ndarray
+        Areas of the section
+    """
+    elevation_angles = np.linspace(0, 180, shape[0])
+
+    areas = np.zeros(shape, dtype=float)
+    azimuthal_increment = np.deg2rad(1)
+    elevation_increment = np.deg2rad(1)
+
+    for i in range(shape[1]):
+        for j in range(shape[0]):
+            areas[j, i] = (
+                (radius**2)
+                * azimuthal_increment
+                * np.sin(np.deg2rad(elevation_angles[j]))
+                * elevation_increment
+            )
+            # hacky way to get rid of singularity on poles
+            areas[0, :] = areas[-1, 0]
+    return areas

src/gemdat/volume.py

@@ -14,6 +14,7 @@
 from scipy.constants import physical_constants
 from skimage.feature import blob_dog
 from skimage.measure import regionprops
+from gemdat.utils import calculate_spherical_areas
 
 from .segmentation import watershed_pbc
 
@@ -24,6 +25,7 @@
 
     from gemdat.path import Pathway
     from gemdat.trajectory import Trajectory
+    from gemdat.orientations import Orientations
 
 
 @dataclass
@@ -448,6 +450,51 @@ def optimal_percolating_path(self, **kwargs) -> Pathway | None:
         return optimal_percolating_path(self, **kwargs)
 
 
+class OrientationalVolume(Volume):
+    """Container for orientational volume data, that are in spherical
+    coordinates."""
+
+    shape: tuple[int, int]
+    radii: np.ndarray
+
+    def __init__(self, shape, radii, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.shape = shape
+        self.radii = radii
+
+    def orientational_peaks(
+        self,
+        **kwargs,
+    ) -> list:
+        """Find peaks using the [Difference of
+        Gaussian][skimage.feature.blob_dog] function in [scikit-
+        image][skimage].
+
+        Parameters
+        ----------
+        **kwargs
+            Additional keyword arguments are passed to [skimage.feature.blob_dog][]
+
+        Returns
+        -------
+        peaks : list
+            List of coordinates of the peaks
+        """
+        kwargs.setdefault('threshold', 0.1)
+
+        theta, phi, values = self.data.T
+
+        x = theta
+        y = phi
+        z = values / values.max()
+
+        blobs = blob_dog(z, **kwargs)
+
+        peaks = [(y[int(i), int(j)], x[int(i), int(j)]) for i, j, sigma in blobs]
+
+        return peaks
+
+
 def trajectory_to_volume(
     trajectory: Trajectory,
     resolution: float = 0.2,
@@ -510,3 +557,60 @@ def trajectory_to_volume(
         label='trajectory',
         units=None,
     )
+
+
+def orientations_to_volume(
+    orientations: Orientations,
+    shape: tuple[int, int],
+    normalize_area: bool,
+) -> OrientationalVolume:
+    """Calculate density volume from orientation.
+
+    The orientations are converted in spherical coordinates,
+    and the azimutal and elevation angles are binned into a 2d histogram.
+
+    Parameters
+    ----------
+    orientations : Orientations
+        Input orientations
+    shape : tuple
+        The shape of the spherical sector in which the trajectory is
+    normalize_area : bool
+        If True, normalize the histogram by the area of the bins
+
+    Returns
+    -------
+    vol : OrientationalVolume
+        Output volume
+    """
+    az, el, r = orientations.vectors_spherical.T
+    az = az.flatten()
+    el = el.flatten()
+
+    hist, _, _ = np.histogram2d(el, az, shape)
+
+    if normalize_area:
+        # Normalize by the area of the bins
+        areas = calculate_spherical_areas(shape)
+        hist = np.divide(hist, areas)
+        # Drop the bins at the poles where normalization is not possible
+        hist = hist[1:-1, :]
+
+    values = hist.T
+    axis_phi, axis_theta = values.shape
+
+    phi = np.linspace(0, 360, axis_phi)
+    theta = np.linspace(0, 180, axis_theta)
+
+    theta, phi = np.meshgrid(theta, phi)
+
+    data = np.stack([theta, phi, values], axis=-1)
+
+    return OrientationalVolume(
+        data=data,
+        shape=shape,
+        radii=r,
+        lattice=orientations.trajectory.get_lattice(),
+        label='orientations',
+        units=None,
+    )

tests/integration/orientations_test.py

@@ -5,7 +5,7 @@
 import numpy as np
 import pytest
 
-from gemdat.orientations import calculate_spherical_areas
+from gemdat.utils import calculate_spherical_areas
 
 matrix = np.array(
     [
@@ -95,3 +95,21 @@ def test_calculate_spherical_areas():
 
     assert isclose(area_grid[0, 1], 3.7304874810631827e-20)
     assert isclose(area_grid.mean(), 0.00018727792392184294)
+
+
+@pytest.vasporicache_available  # type: ignore
+def test_orientations_to_volume(vasp_orientations):
+    ov = vasp_orientations.to_volume()
+    assert ov.label == 'orientations'
+    assert ov.data.shape == (360, 90, 3)
+    assert isclose(ov.data[-1].mean(), 153.43703703703704)
+
+
+@pytest.vasporicache_available  # type: ignore
+def test_orientational_peaks(vasp_orientations):
+    ov = vasp_orientations.to_volume()
+    peaks = ov.orientational_peaks()
+
+    assert len(peaks) == 17
+    assert isclose(peaks[0][1], 88.98876404494382)
+    assert isclose(peaks[-1][0], 197.54874651810584)

tests/integration/plot_mpl_test.py

@@ -92,7 +92,7 @@ def test_energy_along_path(vasp_path):
     plots.energy_along_path(path=vasp_path, structure=structure)
 
 
-@image_comparison2(baseline_images=['rectilinear'])
+@image_comparison2(baseline_images=['rectilinear', 'rectilinear_wpeaks'])
 def test_rectilinear(vasp_orientations):
     matrix = np.array(
         [
@@ -104,6 +104,7 @@ def test_rectilinear(vasp_orientations):
 
     orientations = vasp_orientations.normalize().transform(matrix=matrix)
     plots.rectilinear(orientations=orientations, normalize_histo=False)
+    plots.rectilinear(orientations=orientations, normalize_histo=False, add_peaks=True)
 
 
 @image_comparison2(baseline_images=['bond_length_distribution'])

tests/orientations_test.py

@@ -6,11 +6,8 @@
 from numpy.testing import assert_allclose
 from pymatgen.core import Species
 
-from gemdat.orientations import (
-    Orientations,
-    calculate_spherical_areas,
-)
-from gemdat.utils import fft_autocorrelation
+from gemdat.orientations import Orientations
+from gemdat.utils import fft_autocorrelation, calculate_spherical_areas
 
 
 def test_orientations_init(trajectory):