# Copyright (c) 2010-2025 The Regents of the University of Michigan
# This file is from the freud project, released under the BSD 3-Clause License.
r"""
The :class:`freud.diffraction` module provides functions for computing the
diffraction pattern of particles in systems with long range order.
.. rubric:: Stability
:mod:`freud.diffraction` is **unstable**. When upgrading from version 2.x to
2.y (y > x), existing freud scripts may need to be updated. The API will be
finalized in a future release.
"""
import logging
from importlib.util import find_spec
import numpy as np
import rowan
import scipy.ndimage
import freud._diffraction
import freud.locality
import freud.util
from freud.util import _Compute
_HAS_MPL = find_spec("matplotlib") is not None
if _HAS_MPL:
import matplotlib.cm
import matplotlib.colors
import freud.plot
else:
msg_mpl = "Plotting requires matplotlib."
logger = logging.getLogger(__name__)
class _StaticStructureFactor(_Compute):
def __init__(self):
# abstract class
pass
@_Compute._computed_property
def S_k(self):
"""(:math:`N_{bins}`,) :class:`numpy.ndarray`: Static
structure factor :math:`S(k)` values."""
return self._cpp_obj.getStructureFactor().toNumpyArray()
@property
def k_max(self):
"""float: Maximum value of k at which to calculate the structure
factor."""
return self.bounds[1]
@property
def k_min(self):
"""float: Minimum value of k at which to calculate the structure
factor."""
return self.bounds[0]
@_Compute._computed_property
def min_valid_k(self):
return self._cpp_obj.getMinValidK()
def _repr_png_(self):
try:
return freud.plot._ax_to_bytes(self.plot())
except (AttributeError, ImportError):
return None
[docs]
class StaticStructureFactorDebye(_StaticStructureFactor):
r"""Computes a 1D static structure factor using the
Debye scattering equation.
This computes the static `structure factor
<https://en.wikipedia.org/wiki/Structure_factor>`__ :math:`S(k)` at given
:math:`k` values by averaging over all :math:`\vec{k}` vectors of the same
magnitude. Note that freud employs the physics convention in which
:math:`k` is used, as opposed to the crystallographic one where :math:`q`
is used. The relation is :math:`k=2 \pi q`. The static structure factor
calculation is implemented using the Debye scattering equation:
.. math::
S(k) = \frac{1}{N} \sum_{i=0}^{N} \sum_{j=0}^{N} \text{sinc}(k r_{ij})
where :math:`N` is the number of particles, :math:`\text{sinc}` function is
defined as :math:`\sin x / x` (no factor of :math:`\pi` as in some
conventions). For more information see `this Wikipedia article
<https://en.wikipedia.org/wiki/Structure_factor>`__. For a full derivation
see :cite:`Farrow2009`. Note that the definition requires :math:`S(0) = N`.
.. note::
For 2D systems freud uses the Bessel function :math:`J_0` instead of the
:math:`\text{sinc}` function in the equation above. See
:cite:`Wieder2012` for more information. For users wishing to calculate
the structure factor of quasi 2D systems (i.e. a 2D simulation is used
to model a real system such as particles on a 2D interface or similar)
the 3D formula should be used. In these cases users should use a 3D box
with its longest dimension being in the z-direction and particle
positions of the form :math:`(x, y, 0)`.
This implementation uses an evenly spaced number of :math:`k` points between
`k_min`` and ``k_max``. If ``k_min`` is set to 0 (the default behavior), the
computed structure factor will show :math:`S(0) = N`.
The Debye implementation provides a much faster algorithm, but gives worse
results than :py:attr:`freud.diffraction.StaticStructureFactorDirect`
at low :math:`k` values.
.. note::
This code assumes all particles have a form factor :math:`f` of 1.
Partial structure factors can be computed by providing a set of
``query_points`` and the total number of points in the system ``N_total`` to
the :py:meth:`compute` method. The normalization criterion is based on the
Faber-Ziman formalism. For particle types :math:`\alpha` and :math:`\beta`,
we compute the total scattering function as a sum of the partial scattering
functions as:
.. math::
S(k) - 1 = \sum_{\alpha}\sum_{\beta} \frac{N_{\alpha}
N_{\beta}}{N_{total}^2} \left(S_{\alpha \beta}(k) - 1\right)
Args:
num_k_values (unsigned int):
Number of values to use in :math:`k` space.
k_max (float):
Maximum :math:`k` value to include in the calculation.
k_min (float, optional):
Minimum :math:`k` value included in the calculation. Note that there
are practical restrictions on the validity of the calculation in the
long wavelength regime, see :py:attr:`min_valid_k` (Default value =
0).
"""
def __init__(self, num_k_values, k_max, k_min=0):
self._cpp_obj = freud._diffraction.StaticStructureFactorDebye(
num_k_values, k_max, k_min
)
@property
def num_k_values(self):
"""int: The number of k values used."""
return len(self.k_values)
@property
def k_values(self):
""":class:`numpy.ndarray`: The :math:`k` values for the calculation."""
return self._cpp_obj.getBinCenters().toNumpyArray()
@property
def bounds(self):
"""tuple: A tuple indicating the smallest and largest :math:`k` values
used."""
k_values = self.k_values
return (k_values[0], k_values[len(k_values) - 1])
[docs]
def compute(self, system, query_points=None, N_total=None, reset=True):
r"""Computes static structure factor.
Example for a single component system::
>>> box, points = freud.data.make_random_system(10, 100, seed=0)
>>> sf = freud.diffraction.StaticStructureFactorDebye(
... num_k_values=100, k_max=10, k_min=0
... )
>>> sf.compute((box, points))
freud.diffraction.StaticStructureFactorDebye(...)
Example for partial mixed structure factor for a multiple component
system with types A and B::
>>> N_particles = 100
>>> box, points = freud.data.make_random_system(
... 10, N_particles, seed=0
... )
>>> A_points = points[:N_particles//2]
>>> B_points = points[N_particles//2:]
>>> sf = freud.diffraction.StaticStructureFactorDebye(
... num_k_values=100, k_max=10, k_min=0
... )
>>> sf.compute(
... system=(box, A_points),
... query_points=B_points,
... N_total=N_particles
... )
freud.diffraction.StaticStructureFactorDebye(...)
Args:
system:
Any object that is a valid argument to
:class:`freud.locality.NeighborQuery.from_system`.
Note that box is allowed to change when accumulating
average static structure factor.
query_points ((:math:`N_{query\_points}`, 3) :class:`numpy.ndarray`, optional):
Query points used to calculate the partial structure factor.
Uses the system's points if :code:`None`. See class
documentation for information about the normalization of partial
structure factors. If :code:`None`, the full scattering is
computed. (Default value = :code:`None`).
N_total (int, optional):
Total number of points in the system. This is required if
``query_points`` are provided. See class documentation for
information about the normalization of partial structure
factors.
reset (bool, optional):
Whether to erase the previously computed values before adding
the new computation; if False, will accumulate data (Default
value: True).
""" # noqa E501
if (query_points is None) != (N_total is None):
msg = (
"If query_points are provided, N_total must also be provided "
"in order to correctly compute the normalization of the "
"partial structure factor."
)
raise ValueError(msg)
if reset:
self._reset()
nq = freud.locality.NeighborQuery.from_system(system)
if query_points is None:
query_points = nq.points
else:
query_points = freud.util._convert_array(query_points, shape=(None, 3))
if N_total is None:
N_total = query_points.shape[0]
self._cpp_obj.accumulate(nq._cpp_obj, query_points, N_total)
return self
def _reset(self):
self._cpp_obj.reset()
@_Compute._computed_property
def min_valid_k(self):
"""float: Minimum valid value of k for the computed system box, equal
to :math:`2\\pi/(L/2)=4\\pi/L` where :math:`L` is the minimum side length.
For more information see :cite:`Liu2016`."""
return self._cpp_obj.getMinValidK()
def __repr__(self):
return (
f"freud.diffraction.{type(self).__name__}(num_k_values={self.num_k_values},"
f" k_max={self.k_max}, k_min={self.k_min})"
)
[docs]
def plot(self, ax=None, **kwargs):
r"""Plot static structure factor.
.. note::
This function plots :math:`S(k)` for values above
:py:attr:`min_valid_k`.
Args:
ax (:class:`matplotlib.axes.Axes`, optional): Axis to plot on. If
:code:`None`, make a new figure and axis.
(Default value = :code:`None`)
Returns:
(:class:`matplotlib.axes.Axes`): Axis with the plot.
"""
if not _HAS_MPL:
raise ImportError(msg_mpl)
return freud.plot.line_plot(
self.k_values[self.k_values > self.min_valid_k],
self.S_k[self.k_values > self.min_valid_k],
title="Static Structure Factor",
xlabel=r"$k$",
ylabel=r"$S(k)$",
ax=ax,
)
[docs]
class StaticStructureFactorDirect(_StaticStructureFactor):
r"""Computes a 1D static structure factor by operating on a
:math:`k` space grid.
This computes the static `structure factor
<https://en.wikipedia.org/wiki/Structure_factor>`__ :math:`S(k)` at given
:math:`k` values by averaging over all :math:`\vec{k}` vectors directions of
the same magnitude. Note that freud employs the physics convention in which
:math:`k` is used, as opposed to the crystallographic one where :math:`q` is
used. The relation is :math:`k=2 \pi q`. This is implemented using the
following formula:
.. math::
S(\vec{k}) = \frac{1}{N} \sum_{i=0}^{N} \sum_{j=0}^N e^{i\vec{k} \cdot
\vec{r}_{ij}}
where :math:`N` is the number of particles. Note that the definition
requires :math:`S(0) = N`.
This implementation provides a much slower algorithm, but gives better
results than the :py:attr:`freud.diffraction.StaticStructureFactorDebye`
method at low k values.
The :math:`\vec{k}` vectors are sampled isotropically from a grid defined by
the box's reciprocal lattice vectors. This sampling of reciprocal space is
based on the MIT licensed `Dynasor library
<https://gitlab.com/materials-modeling/dynasor/>`__, modified to use
parallelized C++ and to support larger ranges of :math:`k` values.
For more information see :cite:`Fransson2021`.
.. note::
Currently 2D boxes are not supported for this method. Use Debye instead.
.. note::
This code assumes all particles have a form factor :math:`f` of 1.
Partial structure factors can be computed by providing ``query_points`` and
total number of points in the system ``N_total`` to the :py:meth:`compute`
method. The normalization criterion is based on the Faber-Ziman formalism.
The implemented normalization scheme produces partial structure factors that
converge to the corresponding particle fraction in the infinite-:math:k limit for
pure (diagonal) terms, and to zero for the mixed partial structure factor.
For particle types :math:`\alpha` and :math:`\beta`, we compute the total
scattering function as a sum of the partial scattering functions as:
.. math::
S(k) - 1 = \sum_{\alpha}\sum_{\beta} \frac{N_{\alpha}
N_{\beta}}{N_{total}^2} \left(S_{\alpha \beta}(k) - 1\right)
To convert the partial structure factor to the normalization scheme where all the
partial structure factors converge to unity in the infinite-:math:k limit, the
users can use the following formula (see :cite:`Liu2016`):
.. math::
S_{\alpha \beta, \text{norm2}} = \frac{S_{\alpha \beta}}{x_{\alpha} x_{\beta}}
+ 1 - \frac{\delta_{\alpha \beta}}{x_{\beta}}
Args:
bins (unsigned int):
Number of bins in :math:`k` space.
k_max (float):
Maximum :math:`k` value to include in the calculation.
k_min (float, optional):
Minimum :math:`k` value included in the calculation. Note that
there are practical restrictions on the validity of the
calculation in the long wavelength regime, see :py:attr:`min_valid_k`
(Default value = 0).
num_sampled_k_points (unsigned int, optional):
The desired number of :math:`\vec{k}` vectors to sample from the
reciprocal lattice grid. If set to 0, all :math:`\vec{k}` vectors
are used. If greater than 0, the :math:`\vec{k}` vectors are sampled
from the full grid with uniform radial density, resulting in a
sample of ``num_sampled_k_points`` vectors on average (Default
value = 0).
"""
def __init__(self, bins, k_max, k_min=0, num_sampled_k_points=0):
self._cpp_obj = freud._diffraction.StaticStructureFactorDirect(
int(bins), float(k_max), float(k_min), int(num_sampled_k_points)
)
@property
def nbins(self):
"""float: Number of bins in the histogram."""
return len(self.bin_centers)
@property
def bin_edges(self):
""":class:`numpy.ndarray`: The edges of each bin of :math:`k`."""
return self._cpp_obj.getBinEdges().toNumpyArray()
@property
def bin_centers(self):
""":class:`numpy.ndarray`: The centers of each bin of :math:`k`."""
return self._cpp_obj.getBinCenters().toNumpyArray()
@property
def bounds(self):
"""tuple: A tuple indicating upper and lower bounds of the
histogram."""
bin_edges = self.bin_edges
return (bin_edges[0], bin_edges[len(bin_edges) - 1])
[docs]
def compute(self, system, query_points=None, N_total=None, reset=True):
r"""Computes static structure factor.
Example for a single component system::
>>> box, points = freud.data.make_random_system(10, 100, seed=0)
>>> sf = freud.diffraction.StaticStructureFactorDirect(
... bins=100, k_max=10, k_min=0
... )
>>> sf.compute((box, points))
freud.diffraction.StaticStructureFactorDirect(...)
Example for partial mixed structure factor for multiple component
system with types A and B::
>>> N_particles = 100
>>> box, points = freud.data.make_random_system(
... 10, N_particles, seed=0
... )
>>> A_points = points[:N_particles//2]
>>> B_points = points[N_particles//2:]
>>> sf = freud.diffraction.StaticStructureFactorDirect(
... bins=100, k_max=10, k_min=0
... )
>>> sf.compute(
... (box, A_points),
... query_points=B_points,
... N_total=N_particles
... )
freud.diffraction.StaticStructureFactorDirect(...)
Args:
system:
Any object that is a valid argument to
:class:`freud.locality.NeighborQuery.from_system`. Note that box is
allowed to change when accumulating average static structure factor.
For non-orthorhombic boxes the points are wrapped into a orthorhombic
box.
query_points ((:math:`N_{query\_points}`, 3) :class:`numpy.ndarray`, optional):
Query points used to calculate the partial structure factor.
Uses the system's points if :code:`None`. See class
documentation for information about the normalization of partial
structure factors. If :code:`None`, the full scattering is
computed. (Default value = :code:`None`).
N_total (int, optional):
Total number of points in the system. This is required if
``query_points`` are provided. See class documentation for
information about the normalization of partial structure
factors.
reset (bool, optional):
Whether to erase the previously computed values before adding
the new computation; if False, will accumulate data (Default
value = True).
""" # noqa E501
if (query_points is None) != (N_total is None):
msg = (
"If query_points are provided, N_total must also be provided "
"in order to correctly compute the normalization of the "
"partial structure factor."
)
raise ValueError(msg)
if reset:
self._reset()
# Convert points to float32 to avoid errors when float64 is passed
nq = freud.locality.NeighborQuery.from_system(system)
if query_points is not None:
query_points = freud.util._convert_array(query_points)
else:
query_points = nq.points
if N_total is None:
N_total = nq.points.shape[0]
self._cpp_obj.accumulate(nq._cpp_obj, query_points, N_total)
return self
def _reset(self):
self._cpp_obj.reset()
@_Compute._computed_property
def min_valid_k(self):
"""float: Minimum valid value of k for the computed system box, equal
to :math:`2\\pi/L` where :math:`L` is the minimum side length.
For more information see :cite:`Liu2016`."""
return self._cpp_obj.getMinValidK()
@property
def num_sampled_k_points(self):
r"""int: The target number of :math:`\vec{k}` points to use when
constructing :math:`k` space grid."""
return self._cpp_obj.getNumSampledKPoints()
@_Compute._computed_property
def k_points(self):
r""":class:`numpy.ndarray`: The :math:`\vec{k}` points used in the
calculation."""
k_points = self._cpp_obj.getKPoints()
return np.asarray([[k.x, k.y, k.z] for k in k_points])
def __repr__(self):
return (
f"freud.diffraction.{type(self).__name__}(bins={self.nbins}, "
f"k_max={self.k_max}, k_min={self.k_min}, "
f"num_sampled_k_points={self.num_sampled_k_points})"
)
[docs]
def plot(self, ax=None, **kwargs):
r"""Plot static structure factor.
.. note::
This function plots :math:`S(k)` for values above
:py:attr:`min_valid_k`.
Args:
ax (:class:`matplotlib.axes.Axes`, optional): Axis to plot on. If
:code:`None`, make a new figure and axis.
(Default value = :code:`None`)
Returns:
(:class:`matplotlib.axes.Axes`): Axis with the plot.
"""
if not _HAS_MPL:
raise ImportError(msg_mpl)
return freud.plot.line_plot(
self.bin_centers[self.bin_centers > self.min_valid_k],
self.S_k[self.bin_centers > self.min_valid_k],
title="Static Structure Factor",
xlabel=r"$k$",
ylabel=r"$S(k)$",
ax=ax,
)
[docs]
class DiffractionPattern(_Compute):
r"""Computes a 2D diffraction pattern.
The diffraction image represents the scattering of incident radiation,
and is useful for identifying translational and/or rotational symmetry
present in the system. This class computes the static `structure factor
<https://en.wikipedia.org/wiki/Structure_factor>`__ :math:`S(\vec{k})` for
a plane of wavevectors :math:`\vec{k}` orthogonal to a view axis. The
view orientation :math:`(1, 0, 0, 0)` defaults to looking down the
:math:`z` axis (at the :math:`xy` plane). The points in the system are
converted to fractional coordinates, then binned into a grid whose
resolution is given by ``grid_size``. A higher ``grid_size`` will lead to
a higher resolution. The points are convolved with a Gaussian of width
:math:`\sigma`, given by ``peak_width``. This convolution is performed
as a multiplication in Fourier space. The computed diffraction pattern
can be accessed as a square array of shape ``(output_size, output_size)``.
The :math:`\vec{k}=0` peak is always located at index
``(output_size // 2, output_size // 2)`` and is normalized to have a value
of :math:`S(\vec{k}=0) = N`, per common convention. The
remaining :math:`\vec{k}` vectors are computed such that each peak in the
diffraction pattern satisfies the relationship :math:`\vec{k} \cdot
\vec{R} = 2 \pi N` for some integer :math:`N` and lattice vector of
the system :math:`\vec{R}`. See the `reciprocal lattice Wikipedia page
<https://en.wikipedia.org/wiki/Reciprocal_lattice>`__ for more information.
This method is based on the implementations in the open-source
`GIXStapose application <https://github.com/cmelab/GIXStapose>`_ and its
predecessor, diffractometer :cite:`Jankowski2017`.
Note:
freud only supports diffraction patterns for cubic boxes.
Args:
grid_size (unsigned int):
Resolution of the diffraction grid (Default value = 512).
output_size (unsigned int):
Resolution of the output diffraction image, uses ``grid_size`` if
not provided or ``None`` (Default value = :code:`None`).
"""
def __init__(self, grid_size=512, output_size=None):
super().__init__()
self._grid_size = int(grid_size)
self._output_size = int(grid_size) if output_size is None else int(output_size)
# Cache these because they are system-independent.
self._k_values_orig = np.empty(self.output_size)
self._k_vectors_orig = np.empty((self.output_size, self.output_size, 3))
# Store these computed arrays which are exposed as properties.
self._k_values = np.empty_like(self._k_values_orig)
self._k_vectors = np.empty_like(self._k_vectors_orig)
self._diffraction = np.zeros((self.output_size, self.output_size))
self._frame_counter = 0
def _calc_proj(self, view_orientation, box):
"""Calculate the inverse shear matrix from finding the projected box
vectors whose area of parallogram is the largest.
Args:
view_orientation ((:math:`4`) :class:`numpy.ndarray`):
View orientation as a quaternion.
box (:class:`~.box.Box`):
Simulation box.
Returns:
(2, 2) :class:`numpy.ndarray`:
Inverse shear matrix.
"""
# Rotate the box matrix by the view orientation.
box_matrix = rowan.rotate(view_orientation, box.to_matrix())
# Compute normals for each box face.
# The area of the face is the length of the vector.
box_face_normals = np.cross(
np.roll(box_matrix, 1, axis=-1), np.roll(box_matrix, -1, axis=-1), axis=0
)
# Compute view axis projections.
projections = np.abs(box_face_normals.T @ np.array([0.0, 0.0, 1.0]))
# Determine the largest projection area along the view axis and use
# that face for the projection into 2D.
best_projection_axis = np.argmax(projections)
secondary_axes = (
np.array([best_projection_axis + 1, best_projection_axis + 2]) % 3
)
# Figure out appropriate shear matrix
shear = box_matrix[np.ix_([0, 1], secondary_axes)]
# Return the inverse shear matrix
return np.linalg.inv(shear)
def _transform(self, img, box, inv_shear, zoom):
"""Zoom, shear, and scale diffraction intensities.
Args:
img ((``grid_size, grid_size``) :class:`numpy.ndarray`):
Array of diffraction intensities.
box (:class:`~.box.Box`):
Simulation box.
inv_shear ((2, 2) :class:`numpy.ndarray`):
Inverse shear matrix.
zoom (float):
Scaling factor for incident wavevectors.
Returns:
(``output_size, output_size``) :class:`numpy.ndarray`:
Transformed array of diffraction intensities.
"""
# The adjustments to roll and roll_shift ensure that the peak
# corresponding to k=0 is located at exactly
# (output_size//2, output_size//2), regardless of whether the grid_size
# and output_size are odd or even. This keeps the peak aligned at the
# center of a single pixel, which should always have the maximum value.
roll = img.shape[0] / 2
if img.shape[0] % 2 == 1:
roll -= 0.5
roll_shift = self.output_size / zoom / 2
if self.output_size % 2 == 1:
roll_shift -= 0.5 / zoom
box_matrix = box.to_matrix()
ss = np.max(box_matrix) * inv_shear
shift_matrix = np.array([[1, 0, -roll], [0, 1, -roll], [0, 0, 1]])
# Translation for [roll_shift, roll_shift]
# Then shift using ss
shear_matrix = np.array(
[
[ss[1, 0], ss[0, 0], roll_shift],
[ss[1, 1], ss[0, 1], roll_shift],
[0, 0, 1],
]
)
zoom_matrix = np.diag((zoom, zoom, 1))
# This matrix uses homogeneous coordinates. It is a 3x3 matrix that
# transforms 2D points and adds an offset.
inverse_transform = np.linalg.inv(zoom_matrix @ shear_matrix @ shift_matrix)
return scipy.ndimage.affine_transform(
input=img,
matrix=inverse_transform,
output_shape=(self.output_size, self.output_size),
order=1,
mode="constant",
)
[docs]
def compute(
self,
system,
view_orientation=None,
zoom=4,
peak_width=1,
weights=None,
reset=True,
):
r"""Computes diffraction pattern.
Args:
system:
Any object that is a valid argument to
:class:`freud.locality.NeighborQuery.from_system`.
view_orientation ((:math:`4`) :class:`numpy.ndarray`, optional):
View orientation. Uses :math:`(1, 0, 0, 0)` if not provided
or :code:`None` (Default value = :code:`None`).
zoom (float, optional):
Scaling factor for incident wavevectors (Default value = 4).
peak_width (float, optional):
Width of Gaussian convolved with points, in system length units
(Default value = 1).
weights ((:math:`N`,) :class:`numpy.ndarray`, optional):
Scattering length density (SLD) weights for each point in the system.
If provided, the diffraction pattern will be normalized such
that the sum of the weights is equal to the number of points in
the system, resulting in :math:`S(0) = N`. If not provided,
all points are assumed to have a weight of 1.
(Default value = :code:`None`).
reset (bool, optional):
Whether to erase the previously computed values before adding
the new computations; if False, will accumulate data (Default
value = True).
"""
if reset:
self._diffraction = np.zeros((self.output_size, self.output_size))
self._frame_counter = 0
system = freud.locality.NeighborQuery.from_system(system)
if not system.box.cubic:
msg = "freud.diffraction.DiffractionPattern only supports cubic boxes"
raise ValueError(msg)
if view_orientation is None:
view_orientation = np.array([1.0, 0.0, 0.0, 0.0])
view_orientation = freud.util._convert_array(view_orientation, (4,), np.double)
# Compute the box projection matrix
inv_shear = self._calc_proj(view_orientation, system.box)
# Rotate points by the view quaternion and shear by the box projection
xy = rowan.rotate(view_orientation, system.points)[:, 0:2]
xy = xy @ inv_shear.T
# Normalize weights s.t. S(0) = N
if weights is not None:
weights = weights * len(weights) / np.sum(weights)
# Map positions to [0, 1] and compute a histogram "image"
# Use grid_size+1 bin edges so that there are grid_size bins
xy += 0.5
xy %= 1
im, _, _ = np.histogram2d(
xy[:, 0],
xy[:, 1],
bins=np.linspace(0, 1, self.grid_size + 1),
weights=weights,
)
# Compute FFT and convolve with Gaussian
diffraction_fft = np.fft.fft2(im)
diffraction_fft = scipy.ndimage.fourier_gaussian(
diffraction_fft, peak_width / zoom
)
diffraction_fft = np.fft.fftshift(diffraction_fft)
# Compute the squared modulus of the FFT, which is S(k)
diffraction_frame = np.real(diffraction_fft * np.conjugate(diffraction_fft))
# Transform the image (scale, shear, zoom) and normalize S(k) by the
# number of points
self._N_points = len(system.points)
diffraction_frame = (
self._transform(diffraction_frame, system.box, inv_shear, zoom)
/ self._N_points
)
# Add to the diffraction pattern and increment the frame counter
self._diffraction += np.asarray(diffraction_frame)
self._frame_counter += 1
# Compute a cached array of k-vectors that can be rotated and scaled
if not self._called_compute:
# Create a 1D axis of k-vector magnitudes
self._k_values_orig = np.fft.fftshift(np.fft.fftfreq(n=self.output_size))
# Create a 3D meshgrid of k-vectors with shape
# (output_size, output_size, 3)
self._k_vectors_orig = np.asarray(
np.meshgrid(self._k_values_orig, self._k_values_orig, [0])
).T[0]
# Cache the view orientation and box matrix scale factor for
# lazy evaluation of k-values and k-vectors
self._box_matrix_scale_factor = np.max(system.box.to_matrix())
self._view_orientation = view_orientation
self._k_scale_factor = (
2 * np.pi * self.output_size / (self._box_matrix_scale_factor * zoom)
)
self._k_values_cached = False
self._k_vectors_cached = False
return self
@property
def grid_size(self):
"""int: Resolution of the diffraction grid."""
return self._grid_size
@property
def output_size(self):
"""int: Resolution of the output diffraction image."""
return self._output_size
@_Compute._computed_property
def diffraction(self):
"""
(``output_size``, ``output_size``) :class:`numpy.ndarray`:
Diffraction pattern.
"""
return np.asarray(self._diffraction) / self._frame_counter
@_Compute._computed_property
def N_points(self):
"""int: Number of points used in the last computation."""
return self._N_points
@_Compute._computed_property
def k_values(self):
"""(``output_size``,) :class:`numpy.ndarray`: k-values."""
if not self._k_values_cached:
self._k_values = np.asarray(self._k_values_orig) * self._k_scale_factor
self._k_values_cached = True
return np.asarray(self._k_values)
@_Compute._computed_property
def k_vectors(self):
"""(``output_size``, ``output_size``, 3) :class:`numpy.ndarray`: \
k-vectors."""
if not self._k_vectors_cached:
self._k_vectors = (
rowan.rotate(self._view_orientation, self._k_vectors_orig)
* self._k_scale_factor
)
self._k_vectors_cached = True
return np.asarray(self._k_vectors)
def __repr__(self):
return (
f"freud.diffraction.{type(self).__name__}(grid_size={self.grid_size}, "
f"output_size={self.output_size})"
)
[docs]
def to_image(self, cmap="afmhot", vmin=None, vmax=None):
"""Generates image of diffraction pattern.
Args:
cmap (str, optional):
Colormap name to use (Default value = :code:`'afmhot'`).
vmin (float):
Minimum of the color scale. Uses :code:`4e-6 * N_points` if
not provided or :code:`None` (Default value = :code:`None`).
vmax (float):
Maximum of the color scale. Uses :code:`0.7 * N_points` if
not provided or :code:`None` (Default value = :code:`None`).
Returns:
((output_size, output_size, 4) :class:`numpy.ndarray`):
RGBA array of pixels.
"""
if not _HAS_MPL:
raise ImportError(msg_mpl)
if vmin is None:
vmin = 4e-6 * self.N_points
if vmax is None:
vmax = 0.7 * self.N_points
norm = matplotlib.colors.LogNorm(vmin=vmin, vmax=vmax)
cmap = matplotlib.colormaps[cmap]
image = cmap(norm(np.clip(self.diffraction, vmin, vmax)))
return (image * 255).astype(np.uint8)
[docs]
def plot(self, ax=None, cmap="afmhot", vmin=None, vmax=None):
"""Plot Diffraction Pattern.
Args:
ax (:class:`matplotlib.axes.Axes`, optional): Axis to plot on. If
:code:`None`, make a new figure and axis.
(Default value = :code:`None`)
cmap (str, optional):
Colormap name to use (Default value = :code:`'afmhot'`).
vmin (float):
Minimum of the color scale. Uses :code:`4e-6 * N_points` if
not provided or :code:`None` (Default value = :code:`None`).
vmax (float):
Maximum of the color scale. Uses :code:`0.7 * N_points` if
not provided or :code:`None` (Default value = :code:`None`).
Returns:
(:class:`matplotlib.axes.Axes`): Axis with the plot.
"""
if not _HAS_MPL:
raise ImportError(msg_mpl)
if vmin is None:
vmin = 4e-6 * self.N_points
if vmax is None:
vmax = 0.7 * self.N_points
return freud.plot.diffraction_plot(
self.diffraction, self.k_values, self.N_points, ax, cmap, vmin, vmax
)
def _repr_png_(self):
try:
return freud.plot._ax_to_bytes(self.plot())
except (AttributeError, ImportError):
return None