diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index f7de899de..833f74083 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -22,7 +22,7 @@ jobs: - os: ubuntu-latest python-version: 3.8 - DEPENDENCIES: diffsims==0.5.0 hyperspy~=2.0.rc0 lmfit==0.9.12 matplotlib==3.6 orix==0.9.0 scikit-image==0.19.0 scikit-learn==1.0.0 + DEPENDENCIES: diffsims==0.6.rc1 hyperspy~=2.0.rc0 lmfit==0.9.12 matplotlib==3.7.5 orix==0.12.1 scikit-image==0.19.0 scikit-learn==1.0.0 LABEL: -oldest steps: - uses: actions/checkout@v3 diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 7adde2ce0..354deaeb1 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -7,6 +7,23 @@ All notable changes to this project will be documented in this file. The format is based on `Keep a Changelog `_, and this project adheres to `Semantic Versioning `_. +Unreleased +========== +Fixed +----- +- Fixed indexing error in :meth:`~pyxem.signals.Diffraction2D.get_direct_beam_position` (#1080) + +Added +----- +- Added Examples for doing a Circular Hough Transform and Increased Documentation for Filtering Data (#1082) +Added +----- +- Added `circular_background` to :meth:`~pyxem.signals.Diffraction2D.template_match_disk` to account for + an amorphous circular background when template matching (#1084) + +- Added new datasets of in situ crystalization, Ag SPED, + Organic Semiconductor Orientation mapping, Orientation Mapping, and DPC (#1081) + 2024-05-08 - version 0.18.0 ========== diff --git a/doc/_static/switcher.json b/doc/_static/switcher.json index 817f117af..f9848cfba 100644 --- a/doc/_static/switcher.json +++ b/doc/_static/switcher.json @@ -6,9 +6,13 @@ }, { "name": "stable", - "version": "0.17.0", + "version": "0.18.0", "url": "https://pyxem.readthedocs.io/en/stable/" }, + { + "version": "0.17.0", + "url": "https://pyxem.readthedocs.io/en/v0.17.0/" + }, { "version": "0.16.0", "url": "https://pyxem.readthedocs.io/en/v0.16.0/" diff --git a/doc/bibliography.bib b/doc/bibliography.bib index 3e247826b..1b41a0dc7 100644 --- a/doc/bibliography.bib +++ b/doc/bibliography.bib @@ -42,7 +42,7 @@ @article{Pekin:17 keywords = { Strain measurement,Nanobeam electron diffraction}, pages = {170--176}, title = {{Optimizing disk registration algorithms for nanobeam electron diffraction strain mapping}}, -url = {https://www.sciencedirect.com/science/article/pii/S0304399116304065}, +url = {https://doi.org/10.1016/j.ultramic.2016.12.021}, volume = {176}, year = {2017} } @@ -90,4 +90,13 @@ @article{Viladot:13 volume = {252}, journal = {Journal of microscopy}, doi = {10.1111/jmi.12065} +} +@article{pyxemorientationmapping2022, + title={Free, flexible and fast: Orientation mapping using the multi-core and GPU-accelerated template matching capabilities in the python-based open source 4D-STEM analysis toolbox Pyxem}, + author={Cautaerts, Niels and Crout, Phillip and {\AA}nes, H{\aa}kon Wiik and Prestat, Eric and Jeong, Jiwon and Dehm, Gerhard and Liebscher, Christian H}, + journal={Ultramicroscopy}, + pages={113517}, + year={2022}, + publisher={Elsevier}, + doi={10.1016/j.ultramic.2022.113517} } \ No newline at end of file diff --git a/doc/conf.py b/doc/conf.py index 39fb48dd4..86ffd5938 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -81,6 +81,7 @@ "http://dx.doi.org/10.1088/0965-0393/23/8/083501", # 403 Client Error: Forbidden for url "https://doi.org/10.1111/jmi.12065", # 403 Client Error: Forbidden for url "https://doi.org/10.1111/j.0022-2720.2004.01293.x", # 403 Client Error: Forbidden for url + "https://doi.org/10.1016/j.ultramic.2016.12.021", ] diff --git a/examples/orientation_mapping/README.rst b/examples/orientation_mapping/README.rst new file mode 100644 index 000000000..7260f5880 --- /dev/null +++ b/examples/orientation_mapping/README.rst @@ -0,0 +1,3 @@ +Orientation Mapping +=================== +Below is a gallery of examples on how to do orientation mapping in pyxem. \ No newline at end of file diff --git a/examples/orientation_mapping/mulit_phase_orientation.py b/examples/orientation_mapping/mulit_phase_orientation.py new file mode 100644 index 000000000..ef3463b3c --- /dev/null +++ b/examples/orientation_mapping/mulit_phase_orientation.py @@ -0,0 +1,65 @@ +""" +Multi Phase Orientation Mapping +=============================== +You can also calculate the orientation of the grains for multiple phases using the +:meth:`pyxem.signals.PolarSignal2D.get_orientation` method. This requires that you +simulate the entire S2 space for the phase and then compare to the simulated diffraction. + +For more information on the orientation mapping process see :cite:`pyxemorientationmapping2022` +""" + +import pyxem as pxm +from pyxem.data import fe_multi_phase_grains, fe_bcc_phase, fe_fcc_phase +from diffsims.generators.simulation_generator import SimulationGenerator +from orix.quaternion import Rotation +from orix.sampling import get_sample_reduced_fundamental + +mulit_phase = fe_multi_phase_grains() + +# %% +# First we center the diffraction patterns and get a polar signal +# Increasing the number of npt_azim with give better polar sampling +# but will take longer to compute the orientation map +# The mean=True argument will return the mean pixel value in each bin rather than the sum +# this makes the high k values more visible + +mulit_phase.calibration.center = None +polar_multi = mulit_phase.get_azimuthal_integral2d( + npt=100, npt_azim=360, inplace=False, mean=True +) +polar_multi.plot() + +# %% +# Now we can get make a simulation. In this case we want to set a minimum_intensity which removes the low intensity reflections. +# we also sample the S2 space using the :func`orix.sampling.get_sample_reduced_fundamental` +# We have two phases here so we can make a simulation object with both of the phases. + +bcc = fe_bcc_phase() +fcc = fe_fcc_phase() +bcc.name = "BCC Phase" +fcc.name = "FCC Phase" +fcc.color = "red" +bcc.color = "blue" + +generator = SimulationGenerator(200, minimum_intensity=0.05) +rotations_bcc = get_sample_reduced_fundamental( + resolution=1, point_group=bcc.point_group +) +rotations_fcc = get_sample_reduced_fundamental( + resolution=1, point_group=fcc.point_group +) + +sim = generator.calculate_diffraction2d( + [bcc, fcc], + rotation=[rotations_bcc, rotations_fcc], + max_excitation_error=0.1, + reciprocal_radius=2, + with_direct_beam=False, +) +orientation_map = polar_multi.get_orientation(sim) + + +cmap = orientation_map.to_crystal_map() +cmap.plot() +# %% +# sphinx_gallery_thumbnail_number = 3 diff --git a/examples/orientation_mapping/on_zone_orientation.py b/examples/orientation_mapping/on_zone_orientation.py new file mode 100644 index 000000000..383b1a555 --- /dev/null +++ b/examples/orientation_mapping/on_zone_orientation.py @@ -0,0 +1,62 @@ +""" +On Zone Orientation +=================== +Sometimes you have a tilt boundary and you might want to know the orientation of the +grains on each side of the boundary. This can be done using the +:meth:`pyxem.signals.PolarSignal2D.get_orientation` method. + +For more information on the orientation mapping process see :cite:`pyxemorientationmapping2022` +""" + +from pyxem.data import si_tilt, si_phase +from diffsims.generators.simulation_generator import SimulationGenerator +from orix.quaternion import Rotation +from orix.vector import Vector3d + +simulated_si_tilt = si_tilt() + +# %% +# Pre-Processing +# -------------- +# First we center the diffraction patterns and get a polar signal +# Increasing the number of npt_azim with give better polar sampling +# but will take longer to compute the orientation map +# The mean=True argument will return the mean pixel value in each bin rather than the sum +# this makes the high k values more visible + +simulated_si_tilt.calibration.center = None +polar_si_tilt = simulated_si_tilt.get_azimuthal_integral2d( + npt=100, npt_azim=360, inplace=False, mean=True +) +polar_si_tilt.plot() + +# %% + +# Building a Simulation +# --------------------- +# Now we can get make the orientation map. In this case we have aligned the tilt axis with the z-axis +# so we can use the :func:`orix.vector.Vector3d.from_euler` method to get the rotation axis. +# As always ``with_direct_beam=False`` is important to make sure that the center +# beam does not affect the orientation mapping. + + +phase = si_phase() +generator = SimulationGenerator(200) +sim = generator.calculate_diffraction2d( + phase, + rotation=Rotation.from_euler( + [[0, 0, 0], [0, 0, 0]], + degrees=True, + ), + max_excitation_error=0.1, + reciprocal_radius=1.5, + with_direct_beam=False, +) +# Getting the Orientation +# ----------------------- +# This should be fairly good at finding the orientation of the grains on each side of the tilt boundary. +# The rotation is stored in the rotation column of the orientation map or .isg[2,0] if you want to use the +# rotation as a navigator or plot it directly. + +orientation_map = polar_si_tilt.get_orientation(sim) +orientation_map.plot_over_signal(simulated_si_tilt) diff --git a/examples/orientation_mapping/single_phase_orientation.py b/examples/orientation_mapping/single_phase_orientation.py new file mode 100644 index 000000000..f111f265d --- /dev/null +++ b/examples/orientation_mapping/single_phase_orientation.py @@ -0,0 +1,63 @@ +""" +Single Phase Orientation Mapping +================================ +You can also calculate the orientation of the grains in a single phase sample using the +:meth:`pyxem.signals.PolarSignal2D.get_orientation` method. This requires that you +simulate the entire S2 space for the phase and then compare to the simulated diffraction. + +For more information on the orientation mapping process see :cite:`pyxemorientationmapping2022` +""" + +from pyxem.data import si_phase, si_grains +from diffsims.generators.simulation_generator import SimulationGenerator +from orix.sampling import get_sample_reduced_fundamental + +simulated_si = si_grains() + +# %% +# Pre-Processing +# -------------- +# First we center the diffraction patterns and get a polar signal +# Increasing the number of npt_azim with give better polar sampling +# but will take longer to compute the orientation map +# The mean=True argument will return the mean pixel value in each bin rather than the sum +# this makes the high k values more visible + +simulated_si.calibration.center = None +polar_si = simulated_si.get_azimuthal_integral2d( + npt=100, npt_azim=360, inplace=False, mean=True +) +polar_si.plot() + +# %% +# Building a Simulation +# --------------------- +# Now we can get make a simulation. In this case we want to set a minimum_intensity which removes the low intensity reflections. +# we also sample the S2 space using the :func`orix.sampling.get_sample_reduced_fundamental`. Make sure that you set +# ``with_direct_beam=False`` or the orientation mapping will be unduely affected by the center beam. +phase = si_phase() +generator = SimulationGenerator(200, minimum_intensity=0.05) +rotations = get_sample_reduced_fundamental(resolution=1, point_group=phase.point_group) +sim = generator.calculate_diffraction2d( + phase, + rotation=rotations, + max_excitation_error=0.1, + reciprocal_radius=2, + with_direct_beam=False, +) # Make sure that with_direct_beam ==False + +# %% +# Getting the Orientation +# ----------------------- +# By default the `get_orientation` function uses a gamma correction equilivent to polar_si**0.5. For noisy datasets +# it might be a good idea to reduce the noise (Maybe by averaging neighboring patterns?) or simple background subtraction, +# otherwise the gamma correction will increase the effects of noise on the data. This tries to focus on "Is the Bragg vector +# there?" rather than "Is the Bragg vector the right intensity?" patially because the intensity of the Bragg vector might have +# many different effects. + + +orientation_map = polar_si.get_orientation(sim) +orientation_map.plot_over_signal(simulated_si, vmax="96th") + +# %% +# sphinx_gallery_thumbnail_number = 4 diff --git a/examples/processing/azimuthal_integration.py b/examples/processing/azimuthal_integration.py index 7877c9abd..02772307a 100644 --- a/examples/processing/azimuthal_integration.py +++ b/examples/processing/azimuthal_integration.py @@ -13,95 +13,84 @@ import hyperspy.api as hs import numpy as np -nano_crystals = pxm.data.mgo_nanocrystals(lazy=True) -nano_crystals.calibration( +s = pxm.data.tilt_boundary_data() +s.calibration( center=None ) # set the center to None to use center of the diffraction patterns -nano_crystals1d = nano_crystals.get_azimuthal_integral1d(npt=100, inplace=False) +s.calibration.units = "/AA^{-1}" +s.calibration.scale = 0.03 # To angstroms +s1d = s.get_azimuthal_integral1d(npt=100, inplace=False) -nano_crystals1d.sum().plot() -# %% - -""" -Similarly, the `get_azimuthal_integral2d` method will return a `Polar2D` signal. -""" - -nano_crystals_polar = nano_crystals.get_azimuthal_integral2d( - npt=100, npt_azim=360, inplace=False -) -nano_crystals_polar.sum().plot() +s1d.sum().plot() # %% +# Similarly, the `get_azimuthal_integral2d` method will return a `Polar2D` signal. -""" -There are also other things you can account for with azimuthal integration, such as the -effects of the Ewald sphere. This can be done by calibrating with a known detector distance, -and beam energy. - -Here we just show the effect of just calibrating with the first peak vs. calibrating -with the known beam energy and detector distance. For things like accurate template matching good -calibration can be important when matching to high diffraction vectors. The calibration example gives -more information on how to get the correct values for your microscope/setup. - -If you are doing x-ray diffraction please raise an issue on the pyxem github to let us know! The same -assumptions should apply for each case, but it would be good to test! +s_polar = s.get_azimuthal_integral2d(npt=100, npt_azim=360, inplace=False) +s_polar.sum().plot() -We only show the 1D case here, but the same applies for the 2D case as well! -""" - -nano_crystals.calibration.detector( +# %% +# There are also other things you can account for with azimuthal integration, such as the +# effects of the Ewald sphere. This can be done by calibrating with a known detector distance, +# and beam energy. +# +# Here we just show the effect of just calibrating with the first peak vs. calibrating +# with the known beam energy and detector distance. For things like accurate template matching good +# calibration can be important when matching to high diffraction vectors. The calibration example gives +# more information on how to get the correct values for your microscope/setup. +# +# If you are doing x-ray diffraction please raise an issue on the pyxem github to let us know! The same +# assumptions should apply for each case, but it would be good to test! +# +# We only show the 1D case here, but the same applies for the 2D case as well! + +s.calibration.detector( pixel_size=0.001, detector_distance=0.125, beam_energy=200, center=None, units="k_A^-1", ) # set the center= None to use the center of the diffraction patterns -nano_crystals1d_200 = nano_crystals.get_azimuthal_integral1d(npt=100, inplace=False) -nano_crystals.calibration.detector( +s1d_200 = s.get_azimuthal_integral1d(npt=100, inplace=False) +s.calibration.detector( pixel_size=0.001, detector_distance=0.075, beam_energy=80, center=None, units="k_A^-1", ) # These are just made up pixel sizes and detector distances for illustration -nano_crystals1d_80 = nano_crystals.get_azimuthal_integral1d(npt=100, inplace=False) +s1d_80 = s.get_azimuthal_integral1d(npt=100, inplace=False) hs.plot.plot_spectra( - [nano_crystals1d.sum(), nano_crystals1d_200.sum(), nano_crystals1d_80.sum()], + [s1d.sum(), s1d_200.sum(), s1d_80.sum()], legend=["Flat Ewald Sphere Assumption", "200keV Corrected", "80keV Corrected"], ) # %% +# At times you may want to use a mask to exclude certain pixels from the azimuthal integration or apply an affine +# transformation to the diffraction patterns before azimuthal integration. This can be done using the `mask` and +# `affine` parameters of the `Calibration` object. +# +# Here we just show a random affine transformation for illustration. -""" -At times you may want to use a mask to exclude certain pixels from the azimuthal integration or apply an affine -transformation to the diffraction patterns before azimuthal integration. This can be done using the `mask` and -`affine` parameters of the `Calibration` object. -Here we just show a random affine transformation for illustration. -""" - -mask = nano_crystals.get_direct_beam_mask(radius=20) # Mask the direct beam +mask = s.get_direct_beam_mask(radius=20) # Mask the direct beam affine = np.array( [[0.9, 0.1, 0], [0.1, 0.9, 0], [0, 0, 1]] ) # Just a random affine transformation for illustration -nano_crystals.calibration(mask=mask, affine=affine) -nano_crystals.get_azimuthal_integral2d( - npt=100, npt_azim=360, inplace=False -).sum().plot() -# %% +s.calibration(mask=mask, affine=affine) +s.get_azimuthal_integral2d(npt=100, npt_azim=360, inplace=False).sum().plot() +# %% # The `azimuth_range`-argument lets you choose what angular range to calculate the azimuthal integral for. # The range can be increasing, decreasing, and does not need to be a multiple of pi. -pol1 = nano_crystals.get_azimuthal_integral2d(npt=100, azimuth_range=(-np.pi, np.pi)) +pol1 = s.get_azimuthal_integral2d(npt=100, azimuth_range=(-np.pi, np.pi)) -pol2 = nano_crystals.get_azimuthal_integral2d(npt=100, azimuth_range=(0, 1)) +pol2 = s.get_azimuthal_integral2d(npt=100, azimuth_range=(0, 1)) -pol3 = nano_crystals.get_azimuthal_integral2d( - npt=100, npt_azim=720, azimuth_range=(0, 4 * np.pi) -) +pol3 = s.get_azimuthal_integral2d(npt=100, npt_azim=720, azimuth_range=(0, 4 * np.pi)) -pol4 = nano_crystals.get_azimuthal_integral2d(npt=100, azimuth_range=(np.pi, 0)) +pol4 = s.get_azimuthal_integral2d(npt=100, azimuth_range=(np.pi, 0)) hs.plot.plot_images( [pol1.sum(), pol2.sum(), pol3.sum(), pol4.sum()], diff --git a/examples/processing/centering_the_zero_beam.py b/examples/processing/centering_the_zero_beam.py index 3e70b39d9..48631f7bd 100644 --- a/examples/processing/centering_the_zero_beam.py +++ b/examples/processing/centering_the_zero_beam.py @@ -10,22 +10,25 @@ s = pxm.data.tilt_boundary_data(correct_pivot_point=False) +# %% # Getting the Position of the Zero beam # ------------------------------------- # The zero beam position can be obtained using the :meth:`get_direct_beam_position` method. shifts = s.get_direct_beam_position(method="blur", sigma=5, half_square_width=20) hs.plot.plot_images(shifts.T) # Plotting the shifts -# %% +# %% # Making a Linear Plane # --------------------- # In many instances the zero beam position will vary systematically with the scan position. # This can be corrected by fitting a linear plane to the zero beam position using the # :meth:`make_linear_plane` method. + shifts.make_linear_plane() # Making a linear plane to remove the systematic shift hs.plot.plot_images(shifts.T) # Plotting the shifts after making a linear plane +# %% # Centering the Zero Beam # ----------------------- # The zero beam can be centered using the :meth:`center_direct_beam` method. diff --git a/examples/processing/circular_hough_transform.py b/examples/processing/circular_hough_transform.py new file mode 100644 index 000000000..a3fbd1d2a --- /dev/null +++ b/examples/processing/circular_hough_transform.py @@ -0,0 +1,71 @@ +""" +Circular Hough Transform Peak Finding +===================================== + +The Circular Hough Transform is a method to detect circular features in an image. It can be +used to detect rings (or disks) in the diffraction pattern. For strain mapping there is +some evidence that the Circular Hough transform is more accurate than typical template matching +and can provide subpixel accuracy. + +That being said, the Windowed Template Matching that pyxem uses is quite robust, (both to noise +and in homogeniety in the diffraction disks) although a complete comparison hasn't been fully +studied. +""" + +# %% +# Making a Dummy Dataset +import hyperspy.api as hs +import pyxem as pxm +from skimage.transform import hough_circle +from skimage.feature import canny +import numpy as np + +s = pxm.data.tilt_boundary_data(correct_pivot_point=True) + + +# %% +# Canny Filter +# --------------------- +# First we have to apply a Canny filter to the dataset to get a binary image of the +# outlines of the disks. This is basically a 1st diriviative in reciporical space + +# Filter the image with a Canny filter. +canny_img = s.map( + canny, + sigma=2, + low_threshold=0.6, + high_threshold=0.8, + inplace=False, + use_quantiles=True, +) + +canny_img.plot() # Plotting canny filtered image with outlines +# %% +# Computing the Hough Transform +# ----------------------------- +# We can then compute the hough Transform using the radius of the disk. It is possible to +# have multiple radii but that will return multiple signals for each radii + + +def hough_circle_single_rad(img, radius, **kwargs): + return hough_circle(img, radius, **kwargs)[ + 0 + ] # Otherwise multiple radii are returned + + +circular_hough = canny_img.map(hough_circle_single_rad, radius=6, inplace=False) +circular_hough.plot() + +# %% +# Finding Peaks +# ------------- +# Finding peaks is fairly easy from this point and uses the ``get_diffraction_vectors`` function. + +dv = circular_hough.get_diffraction_vectors(threshold_abs=0.4, min_distance=4) + +m = dv.to_markers(facecolor="none", edgecolor="w") +circular_hough.plot() +circular_hough.add_marker(m) + +# %% +# sphinx_gallery_thumbnail_number = 6 diff --git a/examples/processing/determining_ellipticity.py b/examples/processing/determining_ellipticity.py index 31f600a12..573d294ef 100644 --- a/examples/processing/determining_ellipticity.py +++ b/examples/processing/determining_ellipticity.py @@ -13,7 +13,7 @@ 1.0.0 is developed/released. """ -######################################################################################## +# %% # First, we import the necessary packages and make a single test diffraction pattern. import pyxem.data.dummy_data.make_diffraction_test_data as mdtd import pyxem as pxm @@ -33,15 +33,13 @@ s.plot() # %% +# Using RANSAC Ellipse Determination +# ---------------------------------- +# The RANSAC algorithm is a robust method for fitting an ellipse to points with outliers. Here we +# use it to determine the elliptic distortion of the diffraction pattern and exclude the zero +# beam as "outliers". We can also manually exclude other points from the fit by using the +# mask parameter. -""" -Using RANSAC Ellipse Determination ----------------------------------- -The RANSAC algorithm is a robust method for fitting an ellipse to points with outliers. Here we -use it to determine the elliptic distortion of the diffraction pattern and exclude the zero -beam as "outliers". We can also manually exclude other points from the fit by using the -mask parameter. -""" center, affine, params, pos, inlier = pxm.utils.ransac_ellipse_tools.determine_ellipse( s, use_ransac=True, return_params=True @@ -55,14 +53,13 @@ s.add_marker(in_points, plot_marker=True) s.add_marker(el, plot_marker=True) s.add_marker(out_points, plot_marker=True) + # %% +# Using Manual Ellipse Determination +# ---------------------------------- +# Sometimes it is useful to force the ellipse to fit certain points. For example, here we +# can force the ellipse to fit the first ring by masking the zero beam. -""" -Using Manual Ellipse Determination ----------------------------------- -Sometimes it is useful to force the ellipse to fit certain points. For example, here we -can force the ellipse to fit the first ring by masking the zero beam. -""" mask = s.get_direct_beam_mask(radius=20) @@ -97,16 +94,12 @@ colorbar=None, ) # %% - - -""" -Ellipse from Points -------------------- -This workflow will likely change slightly in the future to add better parllelism, -but here is how to determine the ellipticity from a set of points. -""" - +# Ellipse from Points +# ------------------- +# This workflow will likely change slightly in the future to add better parllelism, +# but here is how to determine the ellipticity from a set of points. # making a test elliptical diffraction pattern + xf, yf, a, b, r, nt = 100, 115, 45, 35, 0, 15 data_points = ret.make_ellipse_data_points(xf, yf, a, b, r, nt) image = np.zeros(shape=(200, 210), dtype=np.float32) @@ -118,7 +111,9 @@ data[:, :] = image.T s = Diffraction2D(data) +# %% # Finding the peaks + s_t = s.template_match_disk(disk_r=5) peak_array = s_t.find_peaks( method="difference_of_gaussian", diff --git a/examples/processing/filtering_data.py b/examples/processing/filtering_data.py index e54953d62..3ce0f8322 100644 --- a/examples/processing/filtering_data.py +++ b/examples/processing/filtering_data.py @@ -1,10 +1,22 @@ """ -Filtering Data -============== +Filtering Data (Averaging Neighboring Patterns etc.) +==================================================== If you have a low number of counts in your data, you may want to filter the data to remove noise. This can be done using the `filter` function which applies some function to the entire dataset and returns a filtered dataset of the same shape. + +In this example, we will use the `filter` function to apply a Gaussian filter in +only the real space dimensions of the dataset. This is useful for removing noise +as it acts like a spatial bandpass filter for high frequencies (noise) in the dataset. + +Because the STEM probe is Gaussian-like, and the convolution of two Gaussian functions +is another Gaussian function, the Gaussian filter is a good choice for filtering and is +equivalent to aquiring the data with a larger probe size equal to +:math:`\sigma_{Filtererd} = \sqrt{\sigma_{Beam}^2 + \sigma_{Filter}^2}`. The benefit is +that the total aquisition time is much shorter than if the probe size was increased to +reach the same number of counts. For detectors with low saturation counts, this can be +a significant advantage. """ from scipy.ndimage import gaussian_filter @@ -13,7 +25,7 @@ import hyperspy.api as hs import numpy as np -s = pxm.data.mgo_nanocrystals(allow_download=True) # MgO nanocrystals dataset +s = pxm.data.tilt_boundary_data() s_filtered = s.filter( gaussian_filter, sigma=1.0, inplace=False @@ -24,17 +36,15 @@ ) # Only filter in real space hs.plot.plot_images( - [s.inav[10, 10], s_filtered.inav[10, 10], s_filtered2.inav[10, 10]], + [s.inav[3, 3], s_filtered.inav[3, 3], s_filtered2.inav[3, 3]], label=["Original", "GaussFilt(all)", "GaussFilt(real space)"], tight_layout=True, vmax="99th", ) # %% -""" -The `filter` function can also be used with a custom function as long as the function -takes a numpy array as input and returns a numpy array of the same shape. -""" +# The `filter` function can also be used with a custom function as long as the function +# takes a numpy array as input and returns a numpy array of the same shape. def custom_filter(array): @@ -45,18 +55,15 @@ def custom_filter(array): s_filtered3 = s.filter(custom_filter, inplace=False) # Custom filter hs.plot.plot_images( - [s.inav[10, 10], s_filtered3.inav[10, 10]], + [s.inav[3, 3], s_filtered3.inav[3, 3]], label=["Original", "GaussFilt(Custom)"], tight_layout=True, vmax="99th", ) # %% - -""" -For lazy datasets, functions which operate on dask arrays can be used. For example, -the `gaussian_filter` function from `scipy.ndimage` is replaced with the `dask_image` -version which operates on dask arrays. -""" +# For lazy datasets, functions which operate on dask arrays can be used. For example, +# the `gaussian_filter` function from `scipy.ndimage` is replaced with the `dask_image` +# version which operates on dask arrays. s = s.as_lazy() # Convert to lazy dataset s_filtered4 = s.filter( @@ -64,9 +71,10 @@ def custom_filter(array): ) # Gaussian filter with sigma=1.0 hs.plot.plot_images( - [s_filtered.inav[10, 10], s_filtered4.inav[10, 10]], + [s_filtered.inav[3, 3], s_filtered4.inav[3, 3]], label=["GaussFilt", "GaussFilt(Lazy)"], tight_layout=True, vmax="99th", ) + # %% diff --git a/examples/processing/template_matching.py b/examples/processing/template_matching.py new file mode 100644 index 000000000..d3272fef9 --- /dev/null +++ b/examples/processing/template_matching.py @@ -0,0 +1,102 @@ +""" +Template Matching +================= + +This example shows how the template matching is done in pyxem to find peaks in a diffraction pattern. +""" + +import numpy as np +import hyperspy.api as hs +import matplotlib.pyplot as plt +from skimage.morphology import disk + + +import pyxem as pxm +import pyxem.data.dummy_data.make_diffraction_test_data as mdtd + +s = pxm.data.tilt_boundary_data() + +# %% +# How Template Matching Works +# =========================== +# +# Pyxem uses a window-normalized cross-correlation to find the peaks. This is much better for finding +# both strongly and weakly scattering peaks but sometimes if the window is too small, too large, or the +# wrong shape, the behavior can be unexpected. + +template = disk(5) +padded_template = np.pad(template, 3) # padding a template increased the window size +padded_template_large = np.pad( + template, 20 +) # padding a template increased the window size + +template_normal = s.template_match(template) +template_padded = s.template_match(padded_template) +template_padded_large = s.template_match(padded_template_large) + +ind = (5, 5) +hs.plot.plot_images( + [ + s.inav[ind], + template_normal.inav[ind], + template_padded.inav[ind], + template_padded_large.inav[ind], + ], + label=["Signal", "Normal Window", "Large Window", "Very Large Window"], + tight_layout=True, + per_row=2, +) + +# %% +# In the very large window case you can see that in the middle the strong zero beam is __included__ in the window +# and the intensity for those pixels is suppressed in the template matching. We also see this occurs in a +# square pattern because even though our template is a disk, the window is a square! +# +# Template Matching with an Amorphous Halo +# ========================================= +# +# Sometimes the template matching can be thrown off by an amorphous halo around the diffraction spots. This +# can be seen in the following example. In this case we can use a dilated (circular) window to reduce the +# effect of the imposed square window. This can be seen in the intensity of the diffraction vectors at +# +y,+x and -y,-x. + +data = mdtd.generate_4d_data( + image_size_x=s.axes_manager.signal_axes[0].size, + image_size_y=s.axes_manager.signal_axes[1].size, + disk_x=s.axes_manager.signal_axes[0].size // 2, + disk_y=s.axes_manager.signal_axes[0].size // 2, + ring_r=45 / 128 * s.axes_manager.signal_axes[0].size // 2, + ring_x=s.axes_manager.signal_axes[0].size // 2, + ring_y=s.axes_manager.signal_axes[0].size // 2, + disk_I=10, + ring_lw=8, + ring_I=2, +) +amorphous_data = data + s + +template_normal = amorphous_data.template_match_disk(disk_r=6, subtract_min=False) +template_circular = amorphous_data.template_match_disk( + disk_r=6, circular_background=True, template_dilation=5, subtract_min=False +) +mask = template_normal.get_direct_beam_mask(35) + +mask2 = ~template_normal.get_direct_beam_mask(55) + +template_normal.data[:, :, mask] = 0 +template_circular.data[:, :, mask] = 0 + +template_normal.data[:, :, mask2] = 0 +template_circular.data[:, :, mask2] = 0 + +f = plt.figure(figsize=(15, 5)) +ind = (5, 5) +hs.plot.plot_images( + [amorphous_data.inav[ind], template_normal.inav[ind], template_circular.inav[ind]], + label=["Signal", "Square Window", "Circular Window"], + tight_layout=True, + per_row=3, + vmin=[0, 0.7, 0.7], + fig=f, +) + +# %% diff --git a/examples/processing/vector_finding.py b/examples/processing/vector_finding.py index fee5ca65b..e794d5b81 100644 --- a/examples/processing/vector_finding.py +++ b/examples/processing/vector_finding.py @@ -43,7 +43,6 @@ vectors = temp.get_diffraction_vectors(threshold_abs=0.4, min_distance=5) # %% - # Plotting Peaks # ============== # We can plot the peaks using hyperSpy's markers and DiffractionVectors. @@ -52,10 +51,9 @@ s.add_marker(vectors.to_markers(color="red", sizes=10, alpha=0.5)) # %% - # Subpixel Peak Fitting # ===================== - +# # The template matching is done on the pixel grid. To find the peak position more accurately the correlation # can be up-sampled using the :func:`pyxem.signals.DiffractionVectors.subpixel_refine` method. This method takes a # `DiffractionSignal2D` object and uses that to refine the peak positions. diff --git a/examples/standards/pixel_coodinates.py b/examples/standards/pixel_coodinates.py index b665f442b..5c7bfca41 100644 --- a/examples/standards/pixel_coodinates.py +++ b/examples/standards/pixel_coodinates.py @@ -24,7 +24,6 @@ s.plot(axes_ticks=True) # %% - # From the plot above you can see that hyperspy automatically sets the axes ticks to be centered # on each pixel. This means that for a 21x21 pixel image, the center is at (-10, -10) in pixel coordinates. # if we change the scale using the calibration function it will automatically adjust the center. Here it is @@ -37,7 +36,6 @@ # %% - # Azimuthal Integration # --------------------- # @@ -47,7 +45,6 @@ az.plot() # %% - # Non-Linear Axes # --------------- # @@ -63,3 +60,4 @@ az.plot() # %% +# sphinx_gallery_thumbnail_number = 4 diff --git a/examples/vectors/clustering_vectors.py b/examples/vectors/clustering_vectors.py index 4e0230a8d..7a3d3c024 100644 --- a/examples/vectors/clustering_vectors.py +++ b/examples/vectors/clustering_vectors.py @@ -13,7 +13,7 @@ from sklearn.cluster import DBSCAN # Getting the vectors for some dataset -s = pxm.data.mgo_nanocrystals() +s = pxm.data.mgo_nanocrystals(allow_download=True) s.data[s.data < 120] = 1 s.filter(gaussian_filter, sigma=(0.5, 0.5, 0, 0), inplace=True) # only in real space s.template_match_disk(disk_r=3, subtract_min=False, inplace=True) @@ -24,6 +24,9 @@ vectors.flatten_diffraction_vectors() ) # flatten the vectors into a 2D array scan = DBSCAN(eps=1.0, min_samples=2) +# %% +# Clustering the Vectors +# ====================== # It is very important that we first normalize the real and reciprocal space distances # The column scale factors map the real space and reciprocal space distances to the same scale # Here this means that the clustering algorithm operates on 10 nm in real space and .1 nm^-1 in @@ -58,6 +61,10 @@ clustered2 = clustered.cluster_labeled_vectors(method=clusterer) m, p = clustered2.to_markers(s, alpha=0.8, get_polygons=True) +# %% +# Visualizing the Clustering +# ========================== +# # This clustering is decent. It shows that there might be some small tilt boundaries in the data # which segment some of the nano-crystals into different clusters. It also shows the effect of using # a phosphor screen which has some pretty severe after glow. This results in a smearing of the @@ -66,3 +73,5 @@ s.add_marker(m) s.add_marker(p, plot_on_signal=False) # %% + +# sphinx_gallery_thumbnail_number = 3 diff --git a/examples/vectors/glass_symmetry.py b/examples/vectors/glass_symmetry.py index 612558464..6d545dafb 100644 --- a/examples/vectors/glass_symmetry.py +++ b/examples/vectors/glass_symmetry.py @@ -15,6 +15,7 @@ import matplotlib.pyplot as plt import numpy as np +# %% # First we load the data and do some basic processing s = pxm.data.pdnip_glass(allow_download=True) s.axes_manager.signal_axes[0].offset = -23.7 @@ -25,11 +26,15 @@ vectors = s.get_diffraction_vectors(threshold_abs=0.5, min_distance=3) +# %% # Now we can convert to polar vectors + pol = vectors.to_polar() +# %% # This function gets the inscribed angle # accept_threshold is the maximum difference between the two angles subtended by the 3 vectors + ins = pol.get_angles(min_angle=0.05, min_k=0.3, accept_threshold=0.1) flat_vect = ins.flatten_diffraction_vectors() @@ -53,8 +58,8 @@ ) # %% - # cycle through colors in groups of 3 for each symmetry cluster + points = ins.to_markers( color=["b", "b", "b", "g", "g", "g", "y", "y", "y", "r", "r", "r"] ) @@ -64,3 +69,5 @@ s.plot(vmin=0.0) s.add_marker(points) s.add_marker(original_points) + +# %% diff --git a/examples/vectors/slicing_vectors.py b/examples/vectors/slicing_vectors.py index a525511b6..1eab5198d 100644 --- a/examples/vectors/slicing_vectors.py +++ b/examples/vectors/slicing_vectors.py @@ -20,6 +20,7 @@ vectors = s.get_diffraction_vectors(threshold_abs=0.4, min_distance=5) +# %% # Plotting all the vectors s.plot() @@ -58,4 +59,6 @@ s.plot() s.add_marker([all_vectors, slic_vectors]) s.add_marker([all_vectors, slic_vectors]) + # %% +# sphinx_gallery_thumbnail_number = 8 diff --git a/examples/virtual_imaging/creating_a_segmented_detector.py b/examples/virtual_imaging/creating_a_segmented_detector.py index ab40cf97b..c33d33a55 100644 --- a/examples/virtual_imaging/creating_a_segmented_detector.py +++ b/examples/virtual_imaging/creating_a_segmented_detector.py @@ -19,8 +19,9 @@ dp.calibration.center = None # Center the diffraction patterns dp.calibration.scale = 0.1 # Scale the diffraction patterns in reciprocal space - +# %% # Visualizing the virtual detector + cp = _get_control_points( 1, npt_azim=8, diff --git a/examples/virtual_imaging/interactive_virtual_images.py b/examples/virtual_imaging/interactive_virtual_images.py index 032fef8bf..ad610a3d7 100644 --- a/examples/virtual_imaging/interactive_virtual_images.py +++ b/examples/virtual_imaging/interactive_virtual_images.py @@ -16,9 +16,10 @@ s = pxm.data.dummy_data.get_cbed_signal() circle = hs.roi.CircleROI(cx=26, cy=74, r=5.5, r_inner=0) s.plot_integrated_intensity(circle) -# %% +# %% # Also we can do the same with a 1D signal + s = pxm.data.dummy_data.get_cbed_signal() s.calibration(center=None) s1d = s.get_azimuthal_integral1d(npt=100, mean=True) diff --git a/examples/virtual_imaging/other_virtual_imaging.py b/examples/virtual_imaging/other_virtual_imaging.py index 7aeda5149..9595b8f29 100644 --- a/examples/virtual_imaging/other_virtual_imaging.py +++ b/examples/virtual_imaging/other_virtual_imaging.py @@ -24,8 +24,8 @@ vdfs = c.get_virtual_image(rois) vdfs.plot() -# %% +# %% # We can also use multiple different ROIs and combine them into a virtual image. r1 = hs.roi.RectangularROI(-35, -35, 35, 35) diff --git a/pyxem/conftest.py b/pyxem/conftest.py index fcc74c537..34c739924 100644 --- a/pyxem/conftest.py +++ b/pyxem/conftest.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -# Copyright 2016-2022 The pyXem developers +# Copyright 2016-2024 The pyXem developers # # This file is part of pyXem. # diff --git a/pyxem/data/__init__.py b/pyxem/data/__init__.py index 85bbba5e4..6f392494a 100644 --- a/pyxem/data/__init__.py +++ b/pyxem/data/__init__.py @@ -30,6 +30,14 @@ """ from pyxem.data.simulated_tilt import tilt_boundary_data +from pyxem.data.simulated_si import ( + si_phase, + si_tilt, + si_grains, + si_grains_simple, + si_rotations_line, +) +from pyxem.data.simulated_fe import fe_bcc_phase, fe_fcc_phase, fe_multi_phase_grains from pyxem.data._data import ( au_grating, pdnip_glass, @@ -37,6 +45,11 @@ twinned_nanowire, sample_with_g, mgo_nanocrystals, + organic_semiconductor, + cuag_orientations, + feal_stripes, + sped_ag, + pdcusi_insitu, ) __all__ = [ @@ -47,4 +60,17 @@ "sample_with_g", "mgo_nanocrystals", "tilt_boundary_data", + "si_phase", + "si_tilt", + "si_grains", + "si_grains_simple", + "si_rotations_line", + "fe_multi_phase_grains", + "fe_fcc_phase", + "fe_bcc_phase", + "cuag_orientations", + "organic_semiconductor", + "feal_stripes", + "sped_ag", + "pdcusi_insitu", ] diff --git a/pyxem/data/_data.py b/pyxem/data/_data.py index 9774ecb16..cb261c1c4 100644 --- a/pyxem/data/_data.py +++ b/pyxem/data/_data.py @@ -39,6 +39,9 @@ def au_grating(allow_download=False, **kwargs): """An au_grating 4-D STEM dataset used to show calibration. + Data can be acessed from https://zenodo.org/records/11284654 + + Parameters ---------- **kwargs @@ -60,7 +63,7 @@ def au_grating(allow_download=False, **kwargs): def pdnip_glass(allow_download=False, **kwargs): # pragma: no cover """A small PdNiP glass 4-D STEM dataset. - Data can be acessed from https://zenodo.org/records/8429302 + Data can be acessed from https://zenodo.org/records/11284654 Data liscenced under CC BY 4.0 @@ -90,7 +93,7 @@ def zrnb_precipitate(allow_download=False, **kwargs): # pragma: no cover Data liscenced under CC BY 4.0 - Data can be acessed from https://zenodo.org/records/8429302 + Data can be acessed from https://zenodo.org/records/11284654 Parameters ---------- @@ -116,7 +119,7 @@ def zrnb_precipitate(allow_download=False, **kwargs): # pragma: no cover def twinned_nanowire(allow_download=False, **kwargs): # pragma: no cover """A small 4-D STEM dataset of a twinned nanowire for orientation mapping. - Data can be acessed from https://zenodo.org/records/8429302 + Data can be acessed from https://zenodo.org/records/11284654 Parameters ---------- @@ -137,9 +140,11 @@ def twinned_nanowire(allow_download=False, **kwargs): # pragma: no cover def sample_with_g(allow_download=False, **kwargs): # pragma: no cover - """A small 4-D STEM dataset of a twinned nanowire for orientation mapping. + """A small 4-D STEM dataset for orientation mapping with mulitple overlapping grains. - Data can be acessed from https://zenodo.org/records/8429302 + Data liscenced under CC BY 4.0 + + Data can be acessed from https://zenodo.org/records/11284654 Parameters ---------- @@ -162,10 +167,151 @@ def sample_with_g(allow_download=False, **kwargs): # pragma: no cover return hs.load(store, **kwargs) +def cuag_orientations(allow_download=False, **kwargs): # pragma: no cover + """A small 4-D STEM dataset of CuAg with multiple orientations of the Cu FCC phase + for orientation mapping. + + Data liscenced under CC BY 4.0 + + Data can be acessed from https://zenodo.org/records/11284654 + + Parameters + ---------- + allow_download: bool + If True, the file will be downloaded from the repository to the local cache. + **kwargs + Keyword arguments passed to :func:`~hyperspy.api.load`. + Examples + -------- + >>> import pyxem as pxm + >>> s = pxm.data.cuag_orientations(allow_download=True) + >>> print(s) + """ + import zarr + + sample = Dataset("cuzipProcessed.zspy") + file_path = sample.fetch_file_path(allow_download=allow_download) + store = zarr.ZipStore(file_path) # for loading from zip + return hs.load(store, **kwargs) + + +def organic_semiconductor(allow_download=False, **kwargs): # pragma: no cover + """A small 4-D STEM dataset of an organic semiconductor for orientation mapping. + + Data liscenced under CC BY 4.0 + + Data can be acessed from https://zenodo.org/records/11284654 + + Parameters + ---------- + allow_download: bool + If True, the file will be downloaded from the repository to the local cache. + **kwargs + Keyword arguments passed to :func:`~hyperspy.api.load`. + Examples + -------- + >>> import pyxem as pxm + >>> s = pxm.data.organic_semiconductor(allow_download=True) + >>> print(s) + """ + import zarr + + sample = Dataset("data_processed.zspy") + file_path = sample.fetch_file_path(allow_download=allow_download) + store = zarr.ZipStore(file_path) # for loading from zip + return hs.load(store, **kwargs) + + +def pdcusi_insitu(allow_download=False, **kwargs): # pragma: no cover + """A decently sized 5D STEM dataset of a PdCuSi alloy held at 390C. This + dataset is sliced from a larger dataset and is used to demonstrate how to operate + on a in situ 4D STEM dataset. Note that this dataset is ~6GB compressed and about + 30GB uncompressed so it might take a while to download. + + Data liscenced under CC BY 4.0 + + Data can be acessed from https://zenodo.org/records/11284654 + + Parameters + ---------- + allow_download: bool + If True, the file will be downloaded from the repository to the local cache. + **kwargs + Keyword arguments passed to :func:`~hyperspy.api.load`. + + Examples + -------- + >>> import pyxem as pxm + >>> s = pxm.data.pdcusi_insitu(allow_download=True) + >>> print(s) + + """ + import zarr + + sample = Dataset("PdCuSiCrystalization-zip.zspy") + file_path = sample.fetch_file_path(allow_download=allow_download) + store = zarr.ZipStore(file_path) # for loading from zip + return hs.load(store, **kwargs) + + +def feal_stripes(allow_download=False, **kwargs): # pragma: no cover + """A small 4-D STEM dataset for doing DPC on a set of magnetic FeAl stripes + + Data can be acessed from https://zenodo.org/records/11284654 + + Parameters + ---------- + allow_download: bool + If True, the file will be downloaded from the repository to the local cache. + **kwargs + Keyword arguments passed to :func:`~hyperspy.api.load`. + + Examples + -------- + >>> import pyxem as pxm + >>> s = pxm.data.feal_stripes() + >>> print(s) + """ + import zarr + + sample = Dataset("FeAl_stripes.zspy") + file_path = sample.fetch_file_path(allow_download=allow_download) + store = zarr.ZipStore(file_path) # for loading from zip + return hs.load(store, **kwargs) + + +def sped_ag(allow_download=False, **kwargs): # pragma: no cover + """A small 4-D STEM dataset of a Ag sample which is good for demonstrating the + ML capabilities of pyxem. + + Data can be acessed from https://zenodo.org/records/11284654 + + Parameters + ---------- + allow_download: bool + If True, the file will be downloaded from the repository to the local cache. + **kwargs + Keyword arguments passed to :func:`~hyperspy.api.load`. + + Examples + -------- + >>> import pyxem as pxm + >>> s = pxm.data.sped_ag() + >>> print(s) + + """ + import zarr + + sample = Dataset("SPED-Ag.zspy") + file_path = sample.fetch_file_path(allow_download=allow_download) + store = zarr.ZipStore(file_path) # for loading from zip + return hs.load(store, **kwargs) + + def mgo_nanocrystals(allow_download=False, **kwargs): # pragma: no cover """A small 4-D STEM dataset of overlapping MgO nanocrystals - Data can be acessed from https://zenodo.org/records/8429302 + Data can be acessed from https://zenodo.org/records/11284654 Parameters ---------- diff --git a/pyxem/data/_registry.py b/pyxem/data/_registry.py index c8e108fd3..1a3cc4204 100644 --- a/pyxem/data/_registry.py +++ b/pyxem/data/_registry.py @@ -22,7 +22,7 @@ """ -_zenodo_url = "https://zenodo.org/record/10459209/files" +_zenodo_url = "https://zenodo.org/records/11284654/files" # file name : hash _file_names_hash = { "au_xgrating_100kX.hspy": "md5:b0733af9d0a272fc1a47e2f56a324fe5", @@ -35,6 +35,11 @@ "sample_with_g.zspy": "md5:65c0a17387a10ca21ebf1812bab67568", "twinned_nanowire.hdf5": "md5:2765fa855db8bc011252dbb425facc81", "ZrNbPercipitate.zspy": "md5:7012a2bdf564f4fb25299af05412723f", + "cuzipProcessed.zspy": "md5:829ecb8f765acb6e9a22092339c6a268", + "colorwheel.txt": "md5:1555136e42cae858be0716f007bda4e4", + "SPED-Ag.zspy": "md5:8556346543fc19f0ef9bdc0f4a6619b5", + "PdCuSiCrystalization-zip.zspy": "md5:80ec7f95ec250106c586debf5d814325", + "FeAl_stripes.zspy": "md5:702cb0c8ff75062c0cb23b3722c2f859", } # add to _urls and to _file_names_hash # if you want to download the data from a different location diff --git a/pyxem/data/simulated_fe.py b/pyxem/data/simulated_fe.py new file mode 100644 index 000000000..670b34dfc --- /dev/null +++ b/pyxem/data/simulated_fe.py @@ -0,0 +1,91 @@ +import numpy as np +from orix.crystal_map import Phase +from orix.quaternion import Rotation +from diffpy.structure import Atom, Lattice, Structure +from diffsims.generators.simulation_generator import SimulationGenerator +from pyxem.signals import Diffraction2D +import skimage + + +def fe_fcc_phase(): + """Create a fe fcc phase with space group 225. This is a fcc structure with + a lattice parameter of 3.571 Å. + + """ + a = 3.571 + latt = Lattice(a, a, a, 90, 90, 90) + atom_list = [] + for coords in [[0, 0, 0], [0.5, 0, 0.5], [0, 0.5, 0.5], [0.5, 0.5, 0]]: + x, y, z = coords[0], coords[1], coords[2] + atom_list.append(Atom(atype="Fe", xyz=[x, y, z], lattice=latt)) # Motif part A + struct = Structure(atoms=atom_list, lattice=latt) + p = Phase(structure=struct, space_group=225) + return p + + +def fe_bcc_phase(): + """Create a fe bcc phase with space group 229. This is a bcc structure with + a lattice parameter of 2.866 Å. + + """ + a = 2.866 + latt = Lattice(a, a, a, 90, 90, 90) + atom_list = [] + for coords in [[0, 0, 0], [0.5, 0.5, 0.5]]: + x, y, z = coords[0], coords[1], coords[2] + atom_list.append(Atom(atype="Fe", xyz=[x, y, z], lattice=latt)) # Motif part A + struct = Structure(atoms=atom_list, lattice=latt) + p = Phase(structure=struct, space_group=229) + return p + + +def fe_multi_phase_grains(num_grains=2, seed=2, size=20, recip_pixels=128): + bcc = fe_bcc_phase() + fcc = fe_fcc_phase() + gen = SimulationGenerator() + rotations1 = Rotation.random(num_grains) + rotations2 = Rotation.random(num_grains) + + sim = gen.calculate_diffraction2d( + phase=[bcc, fcc], + rotation=[rotations1, rotations2], + reciprocal_radius=2.5, + max_excitation_error=0.1, + ) + dps_bcc = [ + sim.iphase[0] + .irot[i] + .get_diffraction_pattern( + sigma=3, shape=(recip_pixels, recip_pixels), calibration=0.02 + ) + for i in range(num_grains) + ] + dps_fcc = [ + sim.iphase[1] + .irot[i] + .get_diffraction_pattern( + sigma=3, shape=(recip_pixels, recip_pixels), calibration=0.02 + ) + for i in range(num_grains) + ] + rng = np.random.default_rng(seed) + x = rng.integers(0, size, size=num_grains * 2) + y = rng.integers(0, size, size=num_grains * 2) + navigator = np.zeros((size, size)) + for i in range(num_grains * 2): + navigator[x[i], y[i]] = i + 1 + navigator = skimage.segmentation.expand_labels(navigator, distance=size) + grain_data = np.empty((size, size, recip_pixels, recip_pixels)) + for i in range(num_grains * 2): + if i < num_grains: + grain_data[navigator == i + 1] = dps_bcc[i][np.newaxis] + else: + grain_data[navigator == i + 1] = dps_fcc[i % num_grains][np.newaxis] + grains = Diffraction2D(grain_data) + grains.axes_manager.signal_axes[0].name = "kx" + grains.axes_manager.signal_axes[1].name = "kx" + grains.axes_manager.signal_axes[0].units = r"$\AA^{-1}$" + grains.axes_manager.signal_axes[1].units = r"$\AA^{-1}$" + grains.axes_manager.signal_axes[0].scale = 0.02 + grains.axes_manager.signal_axes[1].scale = 0.02 + return grains diff --git a/pyxem/data/simulated_si.py b/pyxem/data/simulated_si.py new file mode 100644 index 000000000..7f6697af0 --- /dev/null +++ b/pyxem/data/simulated_si.py @@ -0,0 +1,183 @@ +import numpy as np +from orix.crystal_map import Phase +from orix.quaternion import Rotation, Orientation +from diffpy.structure import Atom, Lattice, Structure +from diffsims.generators.simulation_generator import SimulationGenerator +from scipy.ndimage import gaussian_filter1d +from pyxem.signals import Diffraction1D, Diffraction2D +import skimage + + +def si_phase(): + """Create a silicon phase with space group 227. This is a diamond cubic structure with + a lattice parameter of 5.431 Å. + + """ + a = 5.431 + latt = Lattice(a, a, a, 90, 90, 90) + atom_list = [] + for coords in [[0, 0, 0], [0.5, 0, 0.5], [0, 0.5, 0.5], [0.5, 0.5, 0]]: + x, y, z = coords[0], coords[1], coords[2] + atom_list.append(Atom(atype="Si", xyz=[x, y, z], lattice=latt)) # Motif part A + atom_list.append( + Atom(atype="Si", xyz=[x + 0.25, y + 0.25, z + 0.25], lattice=latt) + ) # Motif part B + struct = Structure(atoms=atom_list, lattice=latt) + p = Phase(structure=struct, space_group=227) + return p + + +def si_tilt(): + p = si_phase() + gen = SimulationGenerator() + rotations = Rotation.from_euler( + [[0, 0, 0], [10, 0, 0]], + degrees=True, + ) + sim = gen.calculate_diffraction2d( + phase=p, rotation=rotations, reciprocal_radius=1.5, max_excitation_error=0.1 + ) + dp1 = np.flipud( + sim.get_diffraction_pattern(sigma=5, shape=(128, 128)) + ) # flip up/down to go from scatter to diffraction + dp2 = np.flipud(sim.irot[1].get_diffraction_pattern(sigma=5, shape=(128, 128))) + + top = np.tile( + dp1, + (10, 5, 1, 1), + ) + bottom = np.tile( + dp2, + (10, 5, 1, 1), + ) + total = np.hstack((top, bottom)) + tilt = Diffraction2D(total) + tilt.axes_manager.signal_axes[0].name = "kx" + tilt.axes_manager.signal_axes[1].name = "kx" + tilt.axes_manager.signal_axes[0].units = r"$\AA^{-1}$" + tilt.axes_manager.signal_axes[1].units = r"$\AA^{-1}$" + tilt.axes_manager.signal_axes[0].scale = 0.01 + tilt.axes_manager.signal_axes[1].scale = 0.01 + return tilt + + +def si_grains_from_orientations( + oris: Orientation, seed: int = 2, size: int = 20, recip_pixels: int = 128 +): + p = si_phase() + gen = SimulationGenerator() + num_grains = oris.size + + sim = gen.calculate_diffraction2d( + phase=p, rotation=oris, reciprocal_radius=1.5, max_excitation_error=0.1 + ) + dps = [ + np.flipud( + sim.irot[i].get_diffraction_pattern( + sigma=5, shape=(recip_pixels, recip_pixels) + ) + ) + for i in range(num_grains) + ] + rng = np.random.default_rng(seed) + x = rng.integers(0, size, size=num_grains) + y = rng.integers(0, size, size=num_grains) + navigator = np.zeros((size, size)) + for i in range(num_grains): + navigator[x[i], y[i]] = i + 1 + navigator = skimage.segmentation.expand_labels(navigator, distance=size) + grain_data = np.empty((size, size, recip_pixels, recip_pixels)) + for i in range(num_grains): + grain_data[navigator == i + 1] = dps[i][np.newaxis] + grains = Diffraction2D(grain_data) + grains.axes_manager.signal_axes[0].name = "kx" + grains.axes_manager.signal_axes[1].name = "kx" + grains.axes_manager.signal_axes[0].units = r"$\AA^{-1}$" + grains.axes_manager.signal_axes[1].units = r"$\AA^{-1}$" + grains.axes_manager.signal_axes[0].scale = 0.01 + grains.axes_manager.signal_axes[1].scale = 0.01 + return grains + + +def si_grains(num_grains=4, seed=2, size=20, recip_pixels=128, return_rotations=False): + """Generate a simulated dataset with grains in random orientations""" + p = si_phase() + rotations = Orientation.random(num_grains, symmetry=p.point_group) + grains = si_grains_from_orientations( + rotations, seed=seed, size=size, recip_pixels=recip_pixels + ) + + if return_rotations: + return grains, rotations + else: + return grains + + +def si_grains_simple(seed=2, size=20, recip_pixels=128, return_rotations=False): + """Generate a simulated dataset with low-index zone axes""" + p = si_phase() + rotations = Orientation.from_euler( + [ + [0, 0, 0], # [0 0 1] + [0, 45, 0], # [0 1 1] + [0, 54.7, 45], # [1 1 1] + [0, 35, 45], # [1 1 2] + ], + degrees=True, + symmetry=p.point_group, + ) + grains = si_grains_from_orientations( + rotations, seed=seed, size=size, recip_pixels=recip_pixels + ) + + if return_rotations: + return grains, rotations + else: + return grains + + +def si_rotations_line(): + from orix.sampling import get_sample_reduced_fundamental + + p = si_phase() + gen = SimulationGenerator() + rotations = get_sample_reduced_fundamental(resolution=3, point_group=p.point_group) + sim = gen.calculate_diffraction2d( + phase=p, rotation=rotations, max_excitation_error=0.1, reciprocal_radius=2 + ) + dps = [] + for i in range(rotations.size): + dp = np.flipud(sim.irot[i].get_diffraction_pattern(sigma=5, shape=(256, 256))) + dps.append(dp) + line = Diffraction2D(np.array(dps)) + line.axes_manager.signal_axes[0].name = "kx" + line.axes_manager.signal_axes[1].name = "kx" + line.axes_manager.signal_axes[0].units = r"$\AA^{-1}$" + line.axes_manager.signal_axes[1].units = r"$\AA^{-1}$" + line.axes_manager.signal_axes[0].scale = 0.01 + line.axes_manager.signal_axes[1].scale = 0.01 + return line + + +def simulated1dsi( + num_points=200, + accelerating_voltage=200, + reciporical_radius=1, + sigma=2, +): + p = si_phase() + gen = SimulationGenerator( + accelerating_voltage=accelerating_voltage, + ) + x = np.zeros(num_points) + sim = gen.calculate_diffraction1d(phase=p, reciprocal_radius=reciporical_radius) + int_points = np.round( + np.array(sim.reciprocal_spacing) * num_points / reciporical_radius + ).astype(int) + x[int_points] = sim.intensities + x = gaussian_filter1d(x, sigma) + s = Diffraction1D(x) + s.axes_manager[0].name = "k" + s.axes_manager[0].units = r"$\AA^{-1}$" + s.axes_manager[0].scale = reciporical_radius / num_points + return s diff --git a/pyxem/hyperspy_extension.yaml b/pyxem/hyperspy_extension.yaml index ec15642fc..cd7c63a93 100644 --- a/pyxem/hyperspy_extension.yaml +++ b/pyxem/hyperspy_extension.yaml @@ -125,6 +125,18 @@ signals: dtype: real lazy: False module: pyxem.signals.indexation_results + OrientationMap: + signal_type: orientation_map + signal_dimension: 2 + dtype: real + lazy: False + module: pyxem.signals.indexation_results + LazyOrientationMap: + signal_type: orientation_map + signal_dimension: 2 + dtype: real + lazy: True + module: pyxem.signals.indexation_results DisplacementGradientMap: signal_type: tensor_field signal_dimension: 2 diff --git a/pyxem/release_info.py b/pyxem/release_info.py index 0e328bbcb..3c79a81db 100644 --- a/pyxem/release_info.py +++ b/pyxem/release_info.py @@ -1,5 +1,5 @@ name = "pyxem" -version = "0.18.0" +version = "0.19.dev0" author = "Duncan Johnstone, Phillip Crout, Carter Francis, Magnus Nord" copyright = "Copyright 2016-2024, The pyxem developers" # Contributors listed by original committer, maintainers, then other diff --git a/pyxem/signals/__init__.py b/pyxem/signals/__init__.py index 167da8cc4..1df774012 100644 --- a/pyxem/signals/__init__.py +++ b/pyxem/signals/__init__.py @@ -41,7 +41,7 @@ from .polar_vectors import PolarVectors, LazyPolarVectors from .electron_diffraction1d import ElectronDiffraction1D, LazyElectronDiffraction1D from .electron_diffraction2d import ElectronDiffraction2D, LazyElectronDiffraction2D -from .indexation_results import VectorMatchingResults +from .indexation_results import VectorMatchingResults, OrientationMap from .pair_distribution_function1d import PairDistributionFunction1D from .polar_diffraction2d import PolarDiffraction2D, LazyPolarDiffraction2D from .power2d import Power2D, LazyPower2D @@ -99,4 +99,5 @@ "VirtualDarkFieldImage", "InSituDiffraction2D", "LabeledDiffractionVectors2D", + "OrientationMap", ] diff --git a/pyxem/signals/diffraction2d.py b/pyxem/signals/diffraction2d.py index 54e793bad..9d158ba59 100644 --- a/pyxem/signals/diffraction2d.py +++ b/pyxem/signals/diffraction2d.py @@ -654,7 +654,7 @@ def get_direct_beam_position( Half the side length of square that captures the direct beam in all scans. Means that the centering algorithm is stable against diffracted spots brighter than the direct beam. Crops the diffraction - pattern to `half_square_width` pixels around th center of the diffraction + pattern to `half_square_width` pixels around the center of the diffraction pattern. Only one of `half_square_width` or signal_slice can be defined. **kwargs: Additional arguments accepted by :func:`pyxem.utils.diffraction.find_beam_center_blur`, @@ -751,9 +751,10 @@ def get_direct_beam_position( shifts = -centers + origin_coordinates elif method == "center_of_mass": if "mask" in kwargs and signal_slice is not None: + # Shifts mask into coordinate space of sliced signal x, y, r = kwargs["mask"] x = x - signal_slice[0] - y = y - signal_slice[1] + y = y - signal_slice[2] kwargs["mask"] = (x, y, r) centers = signal.center_of_mass( lazy_result=lazy_output, @@ -1011,7 +1012,8 @@ def template_match_disk(self, disk_r=4, inplace=False, **kwargs): If True, will return a LazyDiffraction2D object. If False, will compute the result and return a Diffraction2D object. kwargs : - Passed to :func:`pyxem.utils.diffraction.normalize_template_match` + Passed to :func:`pyxem.utils.diffraction.normalize_template_match` or + the :func:`hyperspy.api.BaseSignal.map` method. Returns ------- @@ -1049,11 +1051,16 @@ def template_match_ring(self, r_inner=5, r_outer=7, inplace=False, **kwargs): ---------- r_inner, r_outer : scalar, optional Inner and outer radius of the rings. + lazy_output : bool, default True + If True, will return a LazyDiffraction2D object. If False, + will compute the result and return a Diffraction2D object. + show_progressbar : bool, default True inplace : bool, optional If True, the data is replaced by the filtered data. If False, a new signal is returned. Default False. kwargs : - Passed to :func:`pyxem.utils.diffraction.normalize_template_match` + Passed to :func:`pyxem.utils.diffraction.normalize_template_match` or + the :func:`hyperspy.api.BaseSignal.map` method. Returns ------- @@ -1094,6 +1101,10 @@ def filter(self, func, inplace=False, **kwargs): numpy or dask array as output which has the same shape, and axes as the input. + This can be useful for adjusting the data inplace or for filtering the + data using something like a Gaussian Kernel to average neighboring + diffraction patterns. + Parameters ---------- func : function diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index cc6d78a49..e35c428f3 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -17,16 +17,40 @@ # along with pyXem. If not, see . from warnings import warn +from typing import Union, Literal, Sequence, Iterator +from traits.api import Undefined import hyperspy.api as hs +from hyperspy._signals.lazy import LazySignal from hyperspy.signal import BaseSignal +from hyperspy.axes import AxesManager, BaseDataAxis import numpy as np -from orix.crystal_map import CrystalMap -from orix.quaternion import Rotation +from orix.crystal_map import CrystalMap, Phase, PhaseList +from orix.quaternion import Rotation, Orientation +from orix.vector import Vector3d +from orix.plot import IPFColorKeyTSL from transforms3d.euler import mat2euler +from diffsims.crystallography._diffracting_vector import DiffractingVector +from orix.vector import Vector3d +from orix.projections import StereographicProjection +from orix.plot.inverse_pole_figure_plot import _get_ipf_axes_labels +from orix.vector.fundamental_sector import _closed_edges_in_hemisphere +from orix.vector import Vector3d +from orix.projections import StereographicProjection +from orix.plot.inverse_pole_figure_plot import _get_ipf_axes_labels +from orix.vector.fundamental_sector import _closed_edges_in_hemisphere +from orix.plot import IPFColorKeyTSL, DirectionColorKeyTSL +import hyperspy.api as hs +import numpy as np +from matplotlib.collections import QuadMesh + +import numpy as np +import hyperspy.api as hs from pyxem.utils.indexation_utils import get_nth_best_solution +from pyxem.signals.diffraction_vectors2d import DiffractionVectors2D from pyxem.utils._signals import _transfer_navigation_axes +from pyxem.utils.signal import compute_markers def crystal_from_vector_matching(z_matches): @@ -147,6 +171,748 @@ def _get_second_best_phase(z): return -1 +def vectors_to_coordinates(vectors): + """ + Convert a set of diffraction vectors to coordinates. For use with the map + function and making markers. + """ + return np.vstack((vectors.x, vectors.y)).T + + +def vectors_to_intensity(vectors, scale=1): + """ + Convert a set of diffraction vectors to coordinates. For use with the map + function and making markers. + """ + + return (vectors.intensity / np.max(vectors.intensity)) * scale + + +def vectors_to_text(vectors): + """ + Convert a set of diffraction vectors to text. For use with the map function + and making text markers. + """ + + def add_bar(i: int) -> str: + if i < 0: + return f"$\\bar{{{abs(i)}}}$" + else: + return f"{i}" + + out = [] + for hkl in vectors.hkl: + h, k, l = np.round(hkl).astype(np.int16) + out.append(f"({add_bar(h)} {add_bar(k)} {add_bar(l)})") + return out + + +def rotation_from_orientation_map(result, rots): + if rots.ndim == 1: + rots = rots[np.newaxis, :] + index, _, rotation, mirror = result.T + index = index.astype(int) + ori = rots[index] + euler = ( + Orientation(ori).to_euler( + degrees=True, + ) + * mirror[..., np.newaxis] + ) + euler[:, 0] = rotation + ori = Orientation.from_euler(euler, degrees=True).data + return ori + + +def extract_vectors_from_orientation_map(result, all_vectors, n_best_index=0): + index, _, rotation, mirror = result[n_best_index, :].T + index = index.astype(int) + if all_vectors.ndim == 0: + vectors = all_vectors + else: + vectors = all_vectors[index] + # Copy manually, as deepcopy adds a lot of overhead with the phase + intensity = vectors.intensity + vectors = DiffractingVector(vectors.phase, xyz=vectors.data.copy()) + # Flip y, as discussed in https://github.com/pyxem/pyxem/issues/925 + vectors.y = -vectors.y + + rotation = Rotation.from_euler( + (mirror * rotation, 0, 0), degrees=True, direction="crystal2lab" + ) + vectors = ~rotation * vectors.to_miller() + vectors = DiffractingVector( + vectors.phase, xyz=vectors.data.copy(), intensity=intensity + ) + + # Mirror if necessary. + # Mirroring, in this case, means casting (r, theta) to (r, -theta). + # This is equivalent to (x, y) -> (x, -y) + vectors.y = mirror * vectors.y + + return vectors + + +def orientation2phase(orientations, sizes): + o = orientations[:, 0] + return np.searchsorted(sizes, o) + + +class OrientationMap(DiffractionVectors2D): + """Signal class for orientation maps. Note that this is a special case where + for each navigation position, the signal contains the top n best matches in the form + of a nx4 array with columns [index,correlation, in-plane rotation, mirror(factor)] + + The Simulation is saved in the metadata but can be accessed using the .simulation attribute. + + Parameters + ---------- + *args + See :class:`~hyperspy._signals.signal2d.Signal2D`. + **kwargs + See :class:`~hyperspy._signals.signal2d.Signal2D` + """ + + _signal_type = "orientation_map" + + @property + def simulation(self): + """Simulation object used to generate the orientation map. This is stored in the metadata + but can be accessed using the ``.simulation`` attribute. The simulation object is a + :class:`diffsims.simulation.Simulation2D` object.""" + return self.metadata.get_item("simulation") + + @simulation.setter + def simulation(self, value): + self.metadata.set_item("simulation", value) + + def deepcopy(self, deepcopy_simulation=False): + """Deepcopy the signal. Deepcopying the simulation is optional and is computationally + expensive. If the simulation is not deepcopied, the simulation attribute of the + copied signal is passed along""" + + if not deepcopy_simulation: + simulation = self.simulation + simulation._phase_slider = None + simulation._rotation_slider = None + self.simulation = None + new = super().deepcopy() + if not deepcopy_simulation: + new.simulation = simulation + self.simulation = simulation + return new + + @property + def num_rows(self): + return self.axes_manager.signal_axes[1].size + + def to_rotation(self, flatten=False): + """ + Convert the orientation map to a set of `orix.Quaternion.Rotation` objects. + + Parameters + ---------- + flatten : bool + If True, the rotations will be flattened to a 2D array. This is useful for passing + to a `CrystalMap` object. + Returns + ------- + orix.Quaternion.Rotation + """ + if self._lazy: + raise ValueError( + "Cannot create rotation from lazy signal. Please compute the signal first." + ) + all_rotations = Rotation.stack(self.simulation.rotations).flatten() + rotations = self.map( + rotation_from_orientation_map, + rots=all_rotations.data, + inplace=False, + lazy_output=False, + output_signal_size=(self.num_rows, 4), + output_dtype=float, + ) + + rots = Rotation(rotations) + if flatten: + shape1 = np.prod(rots.shape[:-1]) + rots = rots.reshape((shape1, rots.shape[-1])) + return rots + + def to_phase_index(self): + """ + Convert the orientation map to a set of phase ids + + Returns + ------- + np.ndarray or None + The phase ids for each pixel in the orientation map or None if the simulation + does not have multiple phases. + """ + if self.simulation.has_multiple_phases: + sizes = np.cumsum([i.size for i in self.simulation.rotations]) + return self.map(orientation2phase, sizes=sizes, inplace=False).data + else: + return None + + def to_single_phase_orientations(self, **kwargs) -> Orientation: + """Convert the orientation map to an `Orientation`-object, + given a single-phase simulation. + + Parameters + ---------- + **kwargs + Additional keyword arguments to pass to the :meth:`hyperspy.api.signals.BaseSignal.map` function. + """ + if self.simulation.has_multiple_phases: + raise ValueError( + "Multiple phases found in simulation (use to_crystal_map instead)" + ) + + # Use the quaternion data from rotations to support 2D rotations, + # i.e. unique rotations for each navigation position + rotations = hs.signals.Signal2D(self.simulation.rotations.data) + + return Orientation( + self.map( + rotation_from_orientation_map, + rots=rotations, + inplace=False, + output_signal_size=(self.num_rows, 4), + output_dtype=float, + **kwargs, + ), + symmetry=self.simulation.phases.point_group, + ) + + def to_single_phase_vectors( + self, n_best_index: int = 0, **kwargs + ) -> hs.signals.Signal1D: + """Get the reciprocal lattice vectors for a single-phase simulation. + + Parameters + ---------- + n_best_index: int + The index into the `n_best` matchese + **kwargs + Additional keyword arguments to pass to the :meth:`hyperspy.api.signals.BaseSignal.map` function. + """ + + if self.simulation.has_multiple_phases: + raise ValueError("Multiple phases found in simulation") + + # Use vector data as signal in case of different vectors per navigation position + vectors_signal = hs.signals.Signal1D(self.simulation.coordinates) + v = self.map( + extract_vectors_from_orientation_map, + all_vectors=vectors_signal, + inplace=False, + output_signal_size=(), + output_dtype=object, + show_progressbar=False, + n_best_index=n_best_index, + **kwargs, + ) + v.metadata.simulation = None + return v + + def to_crystal_map(self) -> CrystalMap: + """Convert the orientation map to an :class:`orix.crystal_map.CrystalMap` object + + Returns + ------- + orix.crystal_map.CrystalMap + A crystal map object with the phase id as the highest matching phase. + """ + if self.axes_manager.navigation_dimension != 2: + raise ValueError( + "Only 2D navigation supported. Please raise an issue if you are interested in " + "support for 3+ navigation dimensions." + ) + + x, y = [ax.axis for ax in self.axes_manager.navigation_axes] + xx, yy = np.meshgrid(x, y) + xx = xx - np.min(xx) + yy = yy - np.min(yy) + scan_unit = self.axes_manager.navigation_axes[0].units + rotations = self.to_rotation(flatten=True) + phase_index = self.to_phase_index() + + if self.simulation.has_multiple_phases: + phases = PhaseList(list(self.simulation.phases)) + if phase_index.ndim == 3: + phase_index = phase_index[..., 0] + phase_index = phase_index.flatten() + else: + phases = PhaseList(self.simulation.phases) + if scan_unit is Undefined: + scan_unit = "px" + + return CrystalMap( + rotations=rotations, + x=xx.flatten(), + phase_id=phase_index, + y=yy.flatten(), + scan_unit=scan_unit, + phase_list=phases, + ) + + def to_ipf_markers(self, offset: float = 0.85, scale: float = 0.2): + """Convert the orientation map to a set of inverse pole figure + markers which visualizes the best matching orientations in the + reduced S2 space. + + Parameters + ---------- + offset : float + The offset of the markers from the center of the plot + scale : float + The scale (as a fraction of the axis) for the markers. + """ + if self._lazy: + raise ValueError( + "Cannot create markers from lazy signal. Please compute the signal first." + ) + if self.simulation.has_multiple_phases: + raise ValueError("Multiple phases found in simulation") + polygon_sector, texts, maxes, mins = self._get_ipf_outline( + offset=offset, scale=scale + ) + + orients = self.to_single_phase_orientations() + vectors = orients * Vector3d.zvector() + vectors = vectors.in_fundamental_sector(self.simulation.phases.point_group) + s = StereographicProjection() + x, y = s.vector2xy(vectors) + x = x.reshape(vectors.shape) + y = y.reshape(vectors.shape) + cor = self.data[..., 1] + offsets = np.empty(shape=vectors.shape[:-1], dtype=object) + correlation = np.empty(shape=vectors.shape[:-1], dtype=object) + + for i in np.ndindex(offsets.shape): + off = np.vstack((x[i], y[i])).T + norm_points = (off - ((maxes + mins) / 2)) / (maxes - mins) * scale + norm_points = norm_points + offset + offsets[i] = norm_points + correlation[i] = cor[i] / np.max(cor[i]) * 0.5 + + square = hs.plot.markers.Squares( + offsets=[[offset, offset]], + widths=(scale + scale / 2,), + units="width", + offset_transform="axes", + facecolor="white", + edgecolor="black", + ) + + best_points = hs.plot.markers.Points( + offsets=offsets.T, + sizes=(4,), + offset_transform="axes", + alpha=correlation.T, + facecolor="green", + ) + + return square, polygon_sector, best_points, texts + + def to_markers( + self, + n_best: int = 1, + annotate: bool = False, + marker_colors: str = ("red", "blue", "green", "orange", "purple"), + text_color: str = "black", + lazy_output: bool = None, + annotation_shift: Sequence[float] = None, + text_kwargs: dict = None, + include_intensity: bool = False, + intesity_scale: float = 1, + **kwargs, + ) -> Sequence[hs.plot.markers.Markers]: + """Convert the orientation map to a set of markers for plotting. + + Parameters + ---------- + n_best: int + The amount of solutions to plot + annotate : bool + If True, the euler rotation and the correlation will be annotated on the plot using + the `Texts` class from hyperspy. + marker_color: str, optional + The color of the point markers used for simulated reflections + text_color: str, optional + The color used for the text annotations for reflections. Does nothing if `annotate` is `False`. + annotation_shift: List[float,float], optional + The shift to apply to the annotation text. Default is [0,-0.1] + include_intensity: bool + If True, the intensity of the diffraction spot will be displayed with more intense peaks + having a larger marker size. + lazy_output: bool + If True, the output will be a lazy signal. If None, the output will be lazy if the input is lazy. + + Returns + ------- + all_markers : Sequence[hs.plot.markers.Markers] + A list of markers for each of the n_best solutions + + """ + if lazy_output is None: + lazy_output = self._lazy + if text_kwargs is None: + text_kwargs = dict() + if annotation_shift is None: + annotation_shift = [0, -0.15] + if not self._lazy: + navigation_chunks = (5,) * self.axes_manager.navigation_dimension + else: + navigation_chunks = None + all_markers = [] + for n in range(n_best): + vectors = self.to_single_phase_vectors( + lazy_output=True, navigation_chunks=navigation_chunks + ) + color = marker_colors[n % len(marker_colors)] + if include_intensity: + intensity = vectors.map( + vectors_to_intensity, + scale=intesity_scale, + inplace=False, + ragged=True, + output_dtype=object, + output_signal_size=(), + navigation_chunks=navigation_chunks, + lazy_output=True, + ).data.T + kwargs["sizes"] = intensity + + coords = vectors.map( + vectors_to_coordinates, + inplace=False, + ragged=True, + output_dtype=object, + output_signal_size=(), + navigation_chunks=navigation_chunks, + lazy_output=True, + ) + markers = hs.plot.markers.Points.from_signal( + coords, facecolor="none", edgecolor=color, **kwargs + ) + all_markers.append(markers) + + if annotate: + texts = vectors.map( + vectors_to_text, + inplace=False, + lazy_output=True, + ragged=True, + output_dtype=object, + output_signal_size=(), + ) + # New signal for offset coordinates, as using inplace=True shifts the point markers too + text_coords = coords.map( + lambda x: x + annotation_shift, + inplace=False, + lazy_output=True, + ) + text_coords = coords.map( + lambda x: x + annotation_shift, + inplace=False, + lazy_output=True, + ) + text_markers = hs.plot.markers.Texts.from_signal( + text_coords, texts=texts.data.T, color=text_color, **text_kwargs + ) + all_markers.append(text_markers) + if lazy_output is False: # Compute all at once (as it is faster) + all_markers = compute_markers(all_markers) + return all_markers + + def to_single_phase_polar_markers( + self, + signal_axes: Sequence[BaseDataAxis], + n_best: int = 1, + marker_colors: str = ("red", "blue", "green", "orange", "purple"), + lazy_output: bool = None, + **kwargs, + ) -> Iterator[hs.plot.markers.Markers]: + """ + Convert the orientation map to a set of markers for plotting in polar coordinates. + + Parameters + ---------- + signal_axes: Sequence[BaseDataAxis] + The signal axes for the orientation map. The first axis should be the azimuthal axis + and the second axis should be the radial axis. + n_best: int + The number of best fit solutions to return as markers + marker_colors: + The colors to use for the markers. If there are more than 5 solutions, the colors will repeat. + lazy_output: + If True, the output will be a set of lazy markers. If False the output will be a set of computed markers. + If None, the output will be lazy if the input is lazy or not lazy if the input is not lazy. + kwargs: + Additional keyword arguments to pass to the hyperspy.plot.markers.Points.from_signal function. + Returns + ------- + all_markers : Sequence[hs.plot.markers.Markers] + An list of markers for each of the n_best solutions + + """ + ( + r_templates, + theta_templates, + intensities_templates, + ) = self.simulation.polar_flatten_simulations( + signal_axes[1].axis, signal_axes[0].axis + ) + + if lazy_output is None: + lazy_output = self._lazy + + def marker_generator_factory(n_best_entry: int, r_axis, theta_axis): + theta_min, theta_max = theta_axis.min(), theta_axis.max() + + def marker_generator(entry): + index, _, rotation, mirror = entry[n_best_entry] + index = index.astype(int) + mirror = mirror.astype(int) + + r_ind = r_templates[index] + r = r_axis[r_ind] + + theta_ind = theta_templates[index] + theta = theta_axis[::mirror][theta_ind] + np.deg2rad(rotation) + + # Rotate as per https://github.com/pyxem/pyxem/issues/925 + theta += np.pi + + # handle wrap-around theta + theta -= theta_min + theta %= theta_max - theta_min + theta += theta_min + + mask = r != 0 + return np.vstack((theta[mask], r[mask])).T + + return marker_generator + + all_markers = [] + for n in range(n_best): + color = marker_colors[n % len(marker_colors)] + markers_signal = self.map( + marker_generator_factory(n, signal_axes[1].axis, signal_axes[0].axis), + inplace=False, + ragged=True, + lazy_output=True, + ) + if "sizes" not in kwargs: + kwargs["sizes"] = 15 + markers = hs.plot.markers.Points.from_signal( + markers_signal, facecolor="none", edgecolor=color, **kwargs + ) + + all_markers.append(markers) + if lazy_output is False: # Compute all at once (as it is faster) + all_markers = compute_markers(all_markers) + return all_markers + + def _get_ipf_outline( + self, include_labels: bool = True, offset: float = 0.85, scale: float = 0.2 + ): + """Get the outline of the IPF for the orientation map as a marker in the + upper right hand corner including labels if desired. + + Parameters + ---------- + include_labels : bool + If True, the labels for the axes will be included. + offset : float + The offset of the markers from the lower left of the plot (as a fraction of the axis). + scale : float + The scale (as a fraction of the axis) for the markers. + + Returns + ------- + polygon_sector : hs.plot.markers.Polygons + The outline of the IPF as a marker + texts : hs.plot.markers.Texts + The text labels for the IPF axes + maxes : np.ndarray + The maximum values for the axes + mins : np.ndarray + The minimum values for the axes + """ + # Creating Lines around QuadMesh + sector = self.simulation.phases.point_group.fundamental_sector + s = StereographicProjection() + + edges = _closed_edges_in_hemisphere(sector.edges, sector) + ex, ey = s.vector2xy(edges) + original_offset = np.vstack((ex, ey)).T + mins, maxes = original_offset.min(axis=0), original_offset.max(axis=0) + original_offset = ( + (original_offset - ((maxes + mins) / 2)) / (maxes - mins) * scale + ) + original_offset = original_offset + offset + polygon_sector = hs.plot.markers.Polygons( + verts=original_offset[np.newaxis], + transform="axes", + alpha=1, + facecolor="none", + ) + if include_labels: + labels = _get_ipf_axes_labels( + sector.vertices, symmetry=self.simulation.phases.point_group + ) + tx, ty = s.vector2xy(sector.vertices) + texts_offset = np.vstack((tx, ty)).T + texts_offset = ( + (texts_offset - ((maxes + mins) / 2)) / (maxes - mins) * scale + ) + texts_offset = texts_offset + offset + texts = hs.plot.markers.Texts( + texts=labels, + offsets=texts_offset, + sizes=(1,), + offset_transform="axes", + facecolor="k", + ) + return polygon_sector, texts, maxes, mins + + def get_ipf_annotation_markers(self, offset: float = 0.85, scale: float = 0.2): + """Get the outline of the IPF for the orientation map as a marker in the + upper right hand corner including labels if desired. As well as the color + mesh for the IPF. + + Parameters + ---------- + offset : float + The offset of the markers from the lower left of the plot (as a fraction of the axis). + scale : float + The scale (as a fraction of the axis) for the markers. + + Returns + ------- + polygon_sector : hs.plot.markers.Polygons + The outline of the IPF as a marker + texts : hs.plot.markers.Texts + The text labels for the IPF axes + mesh : hs.plot.markers.Markers + The color mesh for the IPF (using :class:`matplotlib.collections.QuadMesh`) + """ + + polygon_sector, texts, _, _ = self._get_ipf_outline(offset=offset, scale=scale) + + # Create Color Mesh + color_key = DirectionColorKeyTSL(symmetry=self.simulation.phases.point_group) + g, ext = color_key._create_rgba_grid(return_extent=True) + + max_x = np.max(ext[1]) + min_x = np.min(ext[1]) + + max_y = np.max(ext[0]) + min_y = np.min(ext[0]) + + # center extent: + y = np.linspace(ext[0][0], ext[0][1], g.shape[1] + 1) - ((max_y + min_y) / 2) + + y = y / (max_y - min_y) * scale + offset + + x = np.linspace(ext[1][1], ext[1][0], g.shape[0] + 1) - ((max_x + min_x) / 2) + + x = x / (max_x - min_x) * scale + offset + xx, yy = np.meshgrid(y, x) + + mesh = hs.plot.markers.Markers( + collection=QuadMesh, + coordinates=np.stack((xx, yy), axis=-1), + array=g, + transform="axes", + offset_transform="display", + offsets=[[0, 0]], + ) + return polygon_sector, texts, mesh + + def to_ipf_colormap( + self, + direction: Vector3d = Vector3d.zvector(), + add_markers: bool = True, + ): + """Create a colored navigator and a legend (in the form of a marker) which can be passed as the + navigator argument to the `plot` method of some signal. + + Parameters + ---------- + direction : Vector3d + The direction to plot the IPF in + add_markers : bool + If True, the markers for the IPF will be added to the navigator as permanent markers. + + Returns + ------- + hs.signals.BaseSignal + """ + oris = self.to_single_phase_orientations()[:, :, 0] + ipfcolorkey = IPFColorKeyTSL(oris.symmetry, direction) + + float_rgb = ipfcolorkey.orientation2color(oris) + int_rgb = (float_rgb * 255).astype(np.uint8) + + s = hs.signals.Signal1D(int_rgb) + s.change_dtype("rgb8") + s = s.T + + if add_markers: + annotations = self.get_ipf_annotation_markers() + s.add_marker( + annotations, + permanent=True, + plot_signal=False, + plot_marker=False, + plot_on_signal=False, + ) + return s + + def plot_over_signal( + self, + signal, + add_vector_markers=True, + add_ipf_markers=True, + add_ipf_colorkey=True, + vector_kwargs=None, + **kwargs, + ): + """Convenience method to plot the orientation map and the n-best matches over the signal. + + Parameters + ---------- + signal : BaseSignal + The signal to plot the orientation map over. + add_vector_markers : bool + If True, the vector markers will be added to the signal. + add_ipf_markers : bool + If True, the IPF best fit will be added to the signal in an overlay + add_ipf_colorkey : bool + If True, the IPF colorkey will be added to the signal + vector_kwargs : dict + Additional keyword arguments to pass to the `to_single_phase_markers` method + kwargs + Additional keyword arguments to pass to the + :meth:`hyperspy.api.signals.Signal2D.plot` method + """ + nav = self.to_ipf_colormap(add_markers=False) + if vector_kwargs is None: + vector_kwargs = dict() + signal.plot(navigator=nav, **kwargs) + if add_vector_markers: + signal.add_marker(self.to_markers(1, **vector_kwargs)) + if add_ipf_markers: + ipf_markers = self.to_ipf_markers() + signal.add_marker(ipf_markers) + if add_ipf_colorkey: + signal.add_marker(self.get_ipf_annotation_markers(), plot_on_signal=False) + + class GenericMatchingResults: def __init__(self, data): self.data = hs.signals.Signal2D(data) @@ -199,6 +965,10 @@ def to_crystal_map(self): ) +class LazyOrientationMap(OrientationMap, LazySignal): + pass + + class VectorMatchingResults(BaseSignal): """Vector matching results containing the top n best matching crystal phase and orientation at each navigation position with associated metrics. @@ -273,5 +1043,4 @@ def get_indexed_diffraction_vectors( elif overwrite is True: vectors.hkls = self.hkls - return vectors diff --git a/pyxem/signals/labeled_diffraction_vectors2d.py b/pyxem/signals/labeled_diffraction_vectors2d.py index 775ee79b8..edc1b7700 100644 --- a/pyxem/signals/labeled_diffraction_vectors2d.py +++ b/pyxem/signals/labeled_diffraction_vectors2d.py @@ -255,7 +255,7 @@ def to_markers(self, signal, get_polygons=False, num_points=10, **kwargs): """ offsets, colors, colors_by_index = convert_to_markers(self, signal) - points = hs.plot.markers.Points(offsets=offsets.T, color=colors.T, **kwargs) + points = hs.plot.markers.Points(offsets=offsets.T, edgecolor=colors.T, **kwargs) if get_polygons: verts = self.map_vectors( points_to_polygon, num_points=num_points, dtype=object diff --git a/pyxem/signals/polar_diffraction2d.py b/pyxem/signals/polar_diffraction2d.py index 502581de0..9d620f498 100644 --- a/pyxem/signals/polar_diffraction2d.py +++ b/pyxem/signals/polar_diffraction2d.py @@ -19,10 +19,16 @@ from hyperspy.signals import Signal2D from hyperspy._signals.lazy import LazySignal +from numpy import rad2deg from pyxem.signals.common_diffraction import CommonDiffraction from pyxem.utils._correlations import _correlation, _power, _pearson_correlation from pyxem.utils._deprecated import deprecated +from pyxem.utils.indexation_utils import ( + _mixed_matching_lib_to_polar, + _get_integrated_polar_templates, + _norm_rows, +) from pyxem.utils._background_subtraction import ( _polar_subtract_radial_median, @@ -304,6 +310,103 @@ def subtract_diffraction_background( **kwargs, ) + def get_orientation( + self, + simulation, + n_keep=None, + frac_keep=0.1, + n_best=1, + gamma=0.5, + normalize_templates=True, + **kwargs, + ): + """Match the orientation with some simulated diffraction patterns using + an accelerated orientation mapping algorithm. + The details of the algorithm are described in the paper: + "Free, flexible and fast: Orientation mapping using the multi-core and + GPU-accelerated template matching capabilities in the python-based open + source 4D-STEM analysis toolbox Pyxem" + Parameters + ---------- + simulation : DiffractionSimulation + The diffraction simulation object to use for indexing. + n_keep : int + The number of orientations to keep for each diffraction pattern. + frac_keep : float + The fraction of the best matching orientations to keep. + n_best : int + The number of best matching orientations to keep. + normalize_templates : bool + Normalize the templates to the same intensity. + gamma : float + The gamma correction applied to the diffraction patterns. The default + value is 0.5 which takes the square root of the diffraction patterns to + increase the intensity of the low intensity reflections and decrease the + intensity of the high intensity reflections. In most cases gamma<1 is a + good starting point. See :cite:`pyxemorientationmapping2022` for more information. + kwargs : dict + Any additional options for the :meth:`~hyperspy.signal.BaseSignal.map` function. + Returns + ------- + orientation : BaseSignal + A signal with the orientation at each navigation position. + + References + ---------- + .. bibliography:: + """ + if gamma != 1: + self.data = self.data**gamma + ( + r_templates, + theta_templates, + intensities_templates, + ) = simulation.polar_flatten_simulations( + radial_axes=self.axes_manager.signal_axes[1].axis, + azimuthal_axes=self.axes_manager.signal_axes[0].axis, + ) + radius = self.axes_manager.signal_axes[1].size # number radial pixels + integrated_templates = _get_integrated_polar_templates( + radius, r_templates, intensities_templates, normalize_templates + ) + if normalize_templates: + intensities_templates = _norm_rows(intensities_templates) + orientation = self.map( + _mixed_matching_lib_to_polar, + integrated_templates=integrated_templates, + r_templates=r_templates, + theta_templates=theta_templates, + intensities_templates=intensities_templates, + n_keep=n_keep, + frac_keep=frac_keep, + n_best=n_best, + inplace=False, + transpose=True, + output_signal_size=(n_best, 4), + output_dtype=float, + ) + + # Translate in-plane rotation from index to degrees + # by using the calibration of the axis + def rotation_index_to_degrees(data, axis): + data = data.copy() + ind = data[:, 2].astype(int) + data[:, 2] = rad2deg(axis[ind]) + return data + + orientation.axes_manager.signal_axes[0].name = "n-best" + orientation.axes_manager.signal_axes[1].name = "columns" + + orientation.map( + rotation_index_to_degrees, axis=self.axes_manager.signal_axes[0].axis + ) + + orientation.set_signal_type("orientation_map") + orientation.simulation = simulation + orientation.column_names = ["index", "correlation", "rotation", "factor"] + orientation.units = ["a.u.", "a.u.", "deg", "a.u."] + return orientation + class LazyPolarDiffraction2D(LazySignal, PolarDiffraction2D): pass diff --git a/pyxem/tests/signals/test_diffraction2d.py b/pyxem/tests/signals/test_diffraction2d.py index ed5d5dd60..93fd6c33d 100644 --- a/pyxem/tests/signals/test_diffraction2d.py +++ b/pyxem/tests/signals/test_diffraction2d.py @@ -857,7 +857,7 @@ def setup_method(self): "kind": "nearest", }, ), - ("center_of_mass", {}), + ("center_of_mass", {"mask": (10, 13, 10)}), ], ) def test_get_direct_beam(self, method, sig_slice, kwargs): diff --git a/pyxem/tests/signals/test_indexation_results.py b/pyxem/tests/signals/test_indexation_results.py index 5af586e6f..3f3096afa 100644 --- a/pyxem/tests/signals/test_indexation_results.py +++ b/pyxem/tests/signals/test_indexation_results.py @@ -20,17 +20,30 @@ import pytest from matplotlib import pyplot as plt from transforms3d.euler import euler2mat +import dask.array as da from diffsims.libraries.structure_library import StructureLibrary from diffsims.generators.diffraction_generator import DiffractionGenerator from diffsims.generators.library_generator import DiffractionLibraryGenerator +from diffsims.generators.simulation_generator import SimulationGenerator +from orix.sampling import get_sample_reduced_fundamental +from orix.quaternion import Rotation, Orientation +from orix.crystal_map import CrystalMap + from pyxem.generators import TemplateIndexationGenerator -from pyxem.signals import ( - VectorMatchingResults, - DiffractionVectors, -) +from pyxem.signals import VectorMatchingResults, DiffractionVectors, OrientationMap from pyxem.utils.indexation_utils import OrientationResult +from pyxem.data import ( + si_grains, + si_phase, + si_tilt, + si_grains_simple, + fe_multi_phase_grains, + fe_bcc_phase, + fe_fcc_phase, +) +import hyperspy.api as hs @pytest.fixture @@ -143,3 +156,214 @@ def test_vector_get_indexed_diffraction_vectors_warn(): with pytest.warns(Warning): match_results.get_indexed_diffraction_vectors(vectors) np.testing.assert_allclose(vectors.hkls, [0, 0, 0]) + + +class TestOrientationResult: + """Testing the OrientationMap class for valid outputs. These tests are based on the + examples provided in the documentation. + """ + + @pytest.fixture(scope="class") + def single_rot_orientation_result(self): + s = si_tilt() + s.calibration.center = None + polar_si_tilt = s.get_azimuthal_integral2d( + npt=100, npt_azim=360, inplace=False, mean=True + ) + phase = si_phase() + generator = SimulationGenerator(200) + sim = generator.calculate_diffraction2d( + phase, + rotation=Rotation.from_euler( + [0, 0, 0], + degrees=True, + ), + max_excitation_error=0.1, + reciprocal_radius=1.5, + with_direct_beam=False, + ) + orientation_map = polar_si_tilt.get_orientation(sim) + return orientation_map + + @pytest.fixture(scope="class") + def multi_rot_orientation_result(self): + s, r = si_grains(return_rotations=True) + s.calibration.center = None + polar = s.get_azimuthal_integral2d( + npt=100, npt_azim=180, inplace=False, mean=True + ) + phase = si_phase() + generator = SimulationGenerator(200, minimum_intensity=0.05) + rotations = get_sample_reduced_fundamental( + resolution=1, point_group=phase.point_group + ) + sims = generator.calculate_diffraction2d( + phase, + rotation=rotations, + max_excitation_error=0.1, + reciprocal_radius=2, + with_direct_beam=False, + ) + orientations = polar.get_orientation(sims) + + return orientations, r + + @pytest.fixture(scope="class") + def simple_multi_rot_orientation_result(self): + s, r = si_grains_simple(return_rotations=True) + s.calibration.center = None + polar = s.get_azimuthal_integral2d( + npt=100, npt_azim=180, inplace=False, mean=True + ) + phase = si_phase() + generator = SimulationGenerator(200, minimum_intensity=0.05) + rotations = get_sample_reduced_fundamental( + resolution=1, point_group=phase.point_group + ) + sims = generator.calculate_diffraction2d( + phase, + rotation=rotations, + max_excitation_error=0.1, + reciprocal_radius=2, + with_direct_beam=True, + ) + orientations = polar.get_orientation(sims, alpha=1) # in this case we + return orientations, r, s + + @pytest.fixture(scope="class") + def multi_phase_orientation_result(self): + s = fe_multi_phase_grains() + s.calibration.center = None + polar = s.get_azimuthal_integral2d( + npt=100, npt_azim=180, inplace=False, mean=True + ) + phase = fe_fcc_phase() + phase2 = fe_bcc_phase() + + generator = SimulationGenerator(200, minimum_intensity=0.05) + rotations = get_sample_reduced_fundamental( + resolution=1, point_group=phase.point_group + ) + rotations2 = get_sample_reduced_fundamental( + resolution=1, point_group=phase2.point_group + ) + sims = generator.calculate_diffraction2d( + [phase, phase2], + rotation=[rotations, rotations2], + max_excitation_error=0.1, + reciprocal_radius=2, + ) + orientations = polar.get_orientation(sims, n_best=3) + return orientations + + def test_tilt_orientation_result(self, single_rot_orientation_result): + assert isinstance(single_rot_orientation_result, OrientationMap) + orients = single_rot_orientation_result.to_single_phase_orientations() + # Check that the orientations are within 1 degree of the expected value + degrees_between = orients.angle_with( + Orientation.from_euler([0, 0, 0]), degrees=True + ) + assert np.all( + degrees_between[:, :5] <= 1 + ) # off by 1 degree (due to pixelation?) + degrees_between = orients.angle_with( + Orientation.from_euler([10, 0, 0], degrees=True), degrees=True + ) + assert np.all(degrees_between[:, 5:] <= 1) + + def test_grain_orientation_result(self, simple_multi_rot_orientation_result): + orientations, rotations, s = simple_multi_rot_orientation_result + assert isinstance(rotations, Orientation) + assert isinstance(orientations, OrientationMap) + orients = orientations.to_single_phase_orientations() + + # Check that the orientations are within 2 degrees of the expected value. + # Use 2 degrees since that is the angular resolution of the polar dataset + degrees_between = orients.angle_with(rotations, degrees=True) + assert np.all(np.min(degrees_between, axis=2) <= 2) + + def test_to_crystal_map(self, simple_multi_rot_orientation_result): + orientations, rotations, s = simple_multi_rot_orientation_result + crystal_map = orientations.to_crystal_map() + assert isinstance(crystal_map, CrystalMap) + assert np.all(crystal_map.phase_id == 0) + + def test_to_crystal_map_multi_phase(self, multi_phase_orientation_result): + crystal_map = multi_phase_orientation_result.to_crystal_map() + assert isinstance(crystal_map, CrystalMap) + assert np.all(crystal_map.phase_id < 2) + + @pytest.mark.parametrize("annotate", [True, False]) + @pytest.mark.parametrize("lazy_output", [True, False]) + @pytest.mark.parametrize("add_intensity", [True, False]) + def test_to_markers( + self, simple_multi_rot_orientation_result, annotate, lazy_output, add_intensity + ): + orientations, rotations, s = simple_multi_rot_orientation_result + markers = orientations.to_markers( + lazy_output=lazy_output, annotate=annotate, include_intensity=add_intensity + ) + assert isinstance(markers[0], hs.plot.markers.Markers) + + def test_to_markers_lazy(self, simple_multi_rot_orientation_result): + orientations, rotations, s = simple_multi_rot_orientation_result + orientations = orientations.as_lazy() + markers = orientations.to_markers() + assert isinstance(markers[0].kwargs["offsets"], da.Array) + + def test_to_markers_polar(self, simple_multi_rot_orientation_result): + orientations, rotations, s = simple_multi_rot_orientation_result + polar = s.get_azimuthal_integral2d( + npt=100, npt_azim=180, inplace=False, mean=True + ) + markers = orientations.to_single_phase_polar_markers( + polar.axes_manager.signal_axes + ) + assert isinstance(markers[0], hs.plot.markers.Markers) + + def test_to_ipf_markers(self, simple_multi_rot_orientation_result): + orientations, rotations, s = simple_multi_rot_orientation_result + markers = orientations.to_ipf_markers() + assert isinstance(markers[0], hs.plot.markers.Markers) + + def test_to_ipf_annotation(self, simple_multi_rot_orientation_result): + orientations, rotations, s = simple_multi_rot_orientation_result + annotations = orientations.get_ipf_annotation_markers() + for a in annotations: + assert isinstance(a, hs.plot.markers.Markers) + + @pytest.mark.parametrize("add_markers", [True, False]) + def test_to_ipf_map(self, simple_multi_rot_orientation_result, add_markers): + orientations, rotations, s = simple_multi_rot_orientation_result + navigator = orientations.to_ipf_colormap(add_markers=add_markers) + assert isinstance(navigator, hs.signals.BaseSignal) + if add_markers: + assert len(navigator.metadata.Markers) == 3 + + def test_multi_phase_errors(self, multi_phase_orientation_result): + with pytest.raises(ValueError): + multi_phase_orientation_result.to_ipf_colormap() + with pytest.raises(ValueError): + multi_phase_orientation_result.to_single_phase_vectors() + with pytest.raises(ValueError): + multi_phase_orientation_result.to_single_phase_orientations() + with pytest.raises(ValueError): + multi_phase_orientation_result.to_ipf_markers() + + def test_lazy_error(self, simple_multi_rot_orientation_result): + orientations, rotations, s = simple_multi_rot_orientation_result + orientations = orientations.as_lazy() + with pytest.raises(ValueError): + rotations = orientations.to_rotation() + with pytest.raises(ValueError): + rotations = orientations.to_ipf_markers() + + def test_to_crystal_map_error(self, simple_multi_rot_orientation_result): + orientations, rotations, s = simple_multi_rot_orientation_result + orientations = hs.stack((orientations, orientations)) + with pytest.raises(ValueError): + rotations = orientations.to_crystal_map() + + def test_plot_over_signal(self, simple_multi_rot_orientation_result): + orientations, rotations, s = simple_multi_rot_orientation_result + orientations.plot_over_signal(s) diff --git a/pyxem/tests/signals/test_pixelated_stem_class.py b/pyxem/tests/signals/test_pixelated_stem_class.py index a56d89c5e..a7e061b3f 100644 --- a/pyxem/tests/signals/test_pixelated_stem_class.py +++ b/pyxem/tests/signals/test_pixelated_stem_class.py @@ -683,6 +683,13 @@ def test_disk_radius(self): assert s.data.shape == s_template1.data.shape assert not (s_template0.data == s_template1.data).all() + @pytest.mark.parametrize("mode", ["constant", "reflect"]) + def test_template_circular(self, mode): + s = Diffraction2D(np.random.randint(100, size=(5, 5, 20, 20))) + s_template = s.template_match_disk(circular_background=True, mode=mode) + assert s.data.shape == s_template.data.shape + assert not s_template._lazy + class TestDiffraction2DTemplateMatchRing: def test_simple(self): diff --git a/pyxem/utils/diffraction.py b/pyxem/utils/diffraction.py index e1c38771a..1a5042f8e 100644 --- a/pyxem/utils/diffraction.py +++ b/pyxem/utils/diffraction.py @@ -26,6 +26,7 @@ import pyxem as pxm # for ElectronDiffraction2D from scipy.interpolate import interp1d +from scipy.signal import fftconvolve from skimage import transform as tf from skimage.feature import match_template from skimage import morphology, filters @@ -48,6 +49,18 @@ cp = None ndigpu = None +new_float_type = { + # preserved types + np.float32().dtype.char: np.float32, + np.float64().dtype.char: np.float64, + np.complex64().dtype.char: np.complex64, + np.complex128().dtype.char: np.complex128, + # altered types + np.float16().dtype.char: np.float32, + "g": np.float64, # np.float128 ; doesn't exist on windows + "G": np.complex128, # np.complex256 ; doesn't exist on windows +} + def _index_coords(z, origin=None): """Creates x & y coords for the indices in a numpy array. @@ -855,7 +868,137 @@ def normalize_template_match(z, template, subtract_min=True, pad_input=True, **k **kwargs : Keyword arguments to be passed to :func:`skimage.feature.match_template` """ - template_match = match_template(z, template, pad_input=pad_input, **kwargs) + if kwargs.pop("circular_background", False): + template_match = match_template_dilate(z, template, **kwargs) + else: + template_match = match_template(z, template, pad_input=pad_input, **kwargs) if subtract_min: template_match = template_match - np.min(template_match) return template_match + + +def _supported_float_type(input_dtype, allow_complex=False): + """Return an appropriate floating-point dtype for a given dtype. + + float32, float64, complex64, complex128 are preserved. + float16 is promoted to float32. + complex256 is demoted to complex128. + Other types are cast to float64. + + Parameters + ---------- + input_dtype : np.dtype or tuple of np.dtype + The input dtype. If a tuple of multiple dtypes is provided, each + dtype is first converted to a supported floating point type and the + final dtype is then determined by applying `np.result_type` on the + sequence of supported floating point types. + allow_complex : bool, optional + If False, raise a ValueError on complex-valued inputs. + + Returns + ------- + float_type : dtype + Floating-point dtype for the image. + """ + if isinstance(input_dtype, tuple): + return np.result_type(*(_supported_float_type(d) for d in input_dtype)) + input_dtype = np.dtype(input_dtype) + if not allow_complex and input_dtype.kind == "c": + raise ValueError("complex valued input is not supported") + return new_float_type.get(input_dtype.char, np.float64) + + +def match_template_dilate( + image, template, template_dilation=2, mode="constant", constant_values=0 +): + """Matches a template with an image using a window normalized cross-correlation. + + This performs very well for image with different background intensities. This is a slower version + of the skimage :func:`skimage.feature.match_template` but performs better for images with circular + variations in background intensity, specifically accounting for an amorphous halo around the + diffraction pattern. + + Parameters + ---------- + image : np.array + Image to be matched + template : np.array + Template to preform the normalized cross-correlation with + template_dilation : int + The number of pixels to dilate the template by for the windowed cross-correlation + mode : str + Padding mode for the image. Options are 'constant', 'edge', 'wrap', 'reflect' + constant_values : int + Value to pad the image with if mode is 'constant' + + Returns + ------- + response : np.array + The windowed cross-correlation of the image and template + """ + + image = image.astype(float, copy=False) + if image.ndim < template.ndim: # pragma: no cover + raise ValueError( + "Dimensionality of template must be less than or " + "equal to the dimensionality of image." + ) + if np.any(np.less(image.shape, template.shape)): # pragma: no cover + raise ValueError("Image must be larger than template.") + + image_shape = image.shape + float_dtype = _supported_float_type(image.dtype) + + template = np.pad(template, template_dilation) + pad_width = tuple((width, width) for width in template.shape) + if mode == "constant": + image = np.pad( + image, pad_width=pad_width, mode=mode, constant_values=constant_values + ) + else: + image = np.pad(image, pad_width=pad_width, mode=mode) + + dilated_template = morphology.dilation( + template, footprint=morphology.disk(template_dilation) + ) + # Use special case for 2-D images for much better performance in + # computation of integral images + image_window_sum = fftconvolve(image, dilated_template[::-1, ::-1], mode="valid")[ + 1:-1, 1:-1 + ] + image_window_sum2 = fftconvolve( + image**2, dilated_template[::-1, ::-1], mode="valid" + )[1:-1, 1:-1] + + template_mean = template[ + dilated_template.astype(bool) + ].mean() # only consider the pixels in the dilated template + template_volume = np.sum(dilated_template) + template_ssd = np.sum((template - template_mean) ** 2) + + xcorr = fftconvolve(image, template[::-1, ::-1], mode="valid")[1:-1, 1:-1] + numerator = xcorr - image_window_sum * template_mean + + denominator = image_window_sum2 + np.multiply(image_window_sum, image_window_sum, out=image_window_sum) + np.divide(image_window_sum, template_volume, out=image_window_sum) + denominator -= image_window_sum + denominator *= template_ssd + np.maximum(denominator, 0, out=denominator) # sqrt of negative number not allowed + np.sqrt(denominator, out=denominator) + + response = np.zeros_like(xcorr, dtype=float_dtype) + + # avoid zero-division + mask = denominator > np.finfo(float_dtype).eps + + response[mask] = numerator[mask] / denominator[mask] + + slices = [] + for i in range(template.ndim): + d0 = (template.shape[i] - 1) // 2 + d1 = d0 + image_shape[i] + + slices.append(slice(d0, d1)) + + return response[tuple(slices)] diff --git a/pyxem/utils/indexation_utils.py b/pyxem/utils/indexation_utils.py index a77e117c7..ebbefcbbe 100644 --- a/pyxem/utils/indexation_utils.py +++ b/pyxem/utils/indexation_utils.py @@ -742,6 +742,7 @@ def _mixed_matching_lib_to_polar( n_keep, frac_keep, n_best, + transpose=False, ): """ Match a polar image to a filtered subset of polar templates @@ -776,6 +777,9 @@ def _mixed_matching_lib_to_polar( of the best fitting template, where factor is 1 if the direct template is matched and -1 if the mirror template is matched """ + if transpose: + polar_image = polar_image.T + polar_image = np.nan_to_num(polar_image) dispatcher = get_array_module(polar_image) # remove templates we don't care about with a fast match ( @@ -1099,6 +1103,24 @@ def _prefilter_templates( frac_keep, n_keep, ): + """ + Pre-filter the templates to reduce the number of templates to do a full match on + + Parameters + ---------- + polar_image: 2D numpy.ndarray + The image converted to polar coordinates in the form (theta, r) + r: + theta + intensities + integrated_templates + frac_keep + n_keep + + Returns + ------- + + """ dispatcher = get_array_module(polar_image) max_keep = _get_max_n(r.shape[0], n_keep, frac_keep) template_indexes = dispatcher.arange(r.shape[0], dtype=np.int32) diff --git a/pyxem/utils/signal.py b/pyxem/utils/signal.py index 41071c132..c59b81556 100644 --- a/pyxem/utils/signal.py +++ b/pyxem/utils/signal.py @@ -21,6 +21,7 @@ # This code is not used internally and so can safely be deleted at the end of the deprecation period from pyxem.utils._deprecated import deprecated +import dask.array as da @deprecated(since="0.18.0", removal="0.19.0") @@ -163,3 +164,31 @@ def transfer_navigation_axes_to_signal_axes(new_signal, old_signal): ax_new.units = ax_old.units return new_signal + + +def compute_markers(markers): + """Compute the markers for a list of markers quickly by merging all of the dask arrays + + Parameters + ---------- + markers : list + List of markers. + + Returns + ------- + markers : list + List of markers now computed + """ + to_compute = [] + for m in markers: + for k, i in m.kwargs.items(): + if isinstance(i, da.Array): + to_compute.append(i) + computed = da.compute(to_compute) + ind = 0 + for m in markers: + for k, i in m.kwargs.items(): + if isinstance(i, da.Array): + m.kwargs[k] = computed[0][ind] + ind += 1 + return markers diff --git a/pyxem/utils/vectors.py b/pyxem/utils/vectors.py index 291476fd5..faa3aa6d8 100644 --- a/pyxem/utils/vectors.py +++ b/pyxem/utils/vectors.py @@ -840,13 +840,13 @@ def points_to_polygon(points, num_points=50): """ if len(points) == 0: return np.array([[0, 0], [0, 0], [0, 0]]) - sorted_points = points[np.argsort(points[:, 0])] - min_x = sorted_points[0, 0] - max_x = sorted_points[-1, 0] + sorted_points = points[np.argsort(points[:, 1])] + min_x = sorted_points[0, 1] + max_x = sorted_points[-1, 1] sorted_index = np.linspace(min_x, max_x, num_points) - lo = np.searchsorted(sorted_points[:, 0], sorted_index[:-1], side="left") - hi = np.searchsorted(sorted_points[:, 0], sorted_index[1:], side="right") + lo = np.searchsorted(sorted_points[:, 1], sorted_index[:-1], side="left") + hi = np.searchsorted(sorted_points[:, 1], sorted_index[1:], side="right") min_points = [] for l, h, x in zip(lo, hi, sorted_index): @@ -855,7 +855,7 @@ def points_to_polygon(points, num_points=50): else: min_points.append( [ - np.min(sorted_points[l:h][:, 1]), + np.min(sorted_points[l:h][:, 0]), x, ] ) @@ -865,7 +865,7 @@ def points_to_polygon(points, num_points=50): if l == h: pass else: - max_points.append([np.max(sorted_points[l:h][:, 1]), x]) + max_points.append([np.max(sorted_points[l:h][:, 0]), x]) all_points = min_points + max_points[::-1] all_points = np.array(all_points) return all_points @@ -916,6 +916,6 @@ def convert_to_markers( x_values = x_inds[lo_y:hi_y, 2] y_values = x_inds[lo_y:hi_y, 3] labels = np.array(x_inds[lo_y:hi_y, -1], dtype=int) - by_ind_peaks[i, j] = np.stack((y_values, x_values), axis=1) - by_ind_colors[i, j] = colors_by_index[labels] + by_ind_peaks[j, i] = np.stack((y_values, x_values), axis=1) + by_ind_colors[j, i] = colors_by_index[labels] return by_ind_peaks, by_ind_colors, colors_by_index diff --git a/setup.py b/setup.py index 2d3ab83ca..bcfd0159b 100644 --- a/setup.py +++ b/setup.py @@ -87,15 +87,15 @@ extras_require=extra_feature_requirements, install_requires=[ "dask", - "diffsims >= 0.5", + "diffsims == 0.6.rc1", "hyperspy >= 2.0", "h5py", "lmfit >= 0.9.12", - "matplotlib >= 3.6", + "matplotlib >= 3.7.5", "numba", "numpy", "numexpr != 2.8.6", # bug in 2.8.6 for greek letters need for pyfai - "orix >= 0.9", + "orix >= 0.12.1", "pooch", "psutil", "pyfai <= 2023.9.0", # breaking changes