From 9bccc5eb2897f3c1e6b4d46fd892d6343fda79b8 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Thu, 9 May 2024 11:07:41 -0500 Subject: [PATCH 01/70] Update: Update versions --- doc/_static/switcher.json | 6 +++++- pyxem/conftest.py | 2 +- pyxem/release_info.py | 2 +- 3 files changed, 7 insertions(+), 3 deletions(-) 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/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/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 From 5cbcc2245348c7e17728eb5390564d22468d0204 Mon Sep 17 00:00:00 2001 From: Sivert Dagenborg Date: Mon, 27 May 2024 15:59:25 +0200 Subject: [PATCH 02/70] Fixed indexing error --- pyxem/signals/diffraction2d.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pyxem/signals/diffraction2d.py b/pyxem/signals/diffraction2d.py index 54e793bad..a6a8c649a 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, From 35c2d326648c5d0f25c7776b5c7bd5962b299f36 Mon Sep 17 00:00:00 2001 From: Sivert Dagenborg Date: Tue, 28 May 2024 10:40:32 +0200 Subject: [PATCH 03/70] Updated changelog --- CHANGELOG.rst | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 7adde2ce0..d55b9a028 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -8,6 +8,14 @@ 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) + + + 2024-05-08 - version 0.18.0 ========== Fixed From a0e857c97a00fa6fc51979856c3862038e7d60d1 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 07:57:23 -0500 Subject: [PATCH 04/70] Testing: Add Test for mask+center_of_mass --- pyxem/tests/signals/test_diffraction2d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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): From c9ccddd9650ee98f322aa52f857ac86d64820b1d Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 12:57:03 -0500 Subject: [PATCH 05/70] BugFix: Fix "Broken" link --- doc/bibliography.bib | 2 +- doc/conf.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/doc/bibliography.bib b/doc/bibliography.bib index 3e247826b..4800f1a58 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} } 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", ] From 011e4d27b8fd468a43b57934add19b29bafdc8d9 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 29 Mar 2024 12:38:02 -0500 Subject: [PATCH 06/70] New Feature: Add Template Matching with dilated windowed template --- pyxem/signals/diffraction2d.py | 8 ++++ pyxem/utils/expt_utils.py | 87 ++++++++++++++++++++++++++++++++++ 2 files changed, 95 insertions(+) diff --git a/pyxem/signals/diffraction2d.py b/pyxem/signals/diffraction2d.py index a6a8c649a..3297758e1 100644 --- a/pyxem/signals/diffraction2d.py +++ b/pyxem/signals/diffraction2d.py @@ -1008,6 +1008,8 @@ def template_match_disk(self, disk_r=4, inplace=False, **kwargs): ---------- disk_r : scalar, optional Radius of the disk. Default 4. + circular_background : bool, optional + If True, the windowed normalization will use a circular window lazy_output : bool, default True If True, will return a LazyDiffraction2D object. If False, will compute the result and return a Diffraction2D object. @@ -1050,6 +1052,12 @@ 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. + circular_background : bool, optional + If True, the windowed normalization will use a circular window + 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. diff --git a/pyxem/utils/expt_utils.py b/pyxem/utils/expt_utils.py index f9ae27f53..711d0a146 100644 --- a/pyxem/utils/expt_utils.py +++ b/pyxem/utils/expt_utils.py @@ -24,3 +24,90 @@ ) from pyxem.utils.diffraction import * + + +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 preforms very well + for image with different background intensities. This is a slower version of the skimage `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' + """ + + if image.ndim < template.ndim: + 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)): + 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.mean() + template_volume = math.prod(template.shape) + 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)] From 3dad15d160b230ab5d46babc91c3d8bd507d5d59 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 29 Mar 2024 17:19:56 -0500 Subject: [PATCH 07/70] Documentation: Add example and docstrings --- examples/processing/template_matching.py | 84 ++++++++++++++++++++++++ 1 file changed, 84 insertions(+) create mode 100644 examples/processing/template_matching.py diff --git a/examples/processing/template_matching.py b/examples/processing/template_matching.py new file mode 100644 index 000000000..cea589db5 --- /dev/null +++ b/examples/processing/template_matching.py @@ -0,0 +1,84 @@ +""" +Template Matching +================= + +This example shows how the template matching is done in pyxem to find peaks in a diffraction pattern. +""" + +from skimage.morphology import disk +import numpy as np +import hyperspy.api as hs + +import pyxem as pxm +import pyxem.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. +data = mdtd.generate_4d_data( + image_size_x=256, + image_size_y=256, + disk_x=128, + disk_y=128, + ring_r=45, + ring_x=128, + ring_y=128, + disk_I=10, + ring_lw=6, + ring_I=2, +) +amorphous_data = data + s + + +template_normal = amorphous_data.template_match_disk(disk_r=6) +template_circular = s.template_match_disk(disk_r=6, dilated_template_window=True) + + +ind = (5, 5) +hs.plot.plot_images( + [amorphous_data.inav[ind], template_normal.inav[ind], template_circular.inav[ind]], + label=["Signal", "Square Window", "Large Window"], + tight_layout=True, + per_row=3, +) From 5240c54eeb98a3109d1b8fc43cefc323699c2044 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 10 Apr 2024 07:31:57 -0500 Subject: [PATCH 08/70] BugFix: Update how volume is computed. --- pyxem/utils/expt_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyxem/utils/expt_utils.py b/pyxem/utils/expt_utils.py index 711d0a146..bef4c6547 100644 --- a/pyxem/utils/expt_utils.py +++ b/pyxem/utils/expt_utils.py @@ -82,7 +82,7 @@ def match_template_dilate( )[1:-1, 1:-1] template_mean = template.mean() - template_volume = math.prod(template.shape) + 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] From f0e3259f767935628b88290fa9ec6b7422fae4f1 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 10 Apr 2024 07:46:36 -0500 Subject: [PATCH 09/70] BugFix: Update how volume is computed. --- pyxem/utils/expt_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyxem/utils/expt_utils.py b/pyxem/utils/expt_utils.py index bef4c6547..b98640305 100644 --- a/pyxem/utils/expt_utils.py +++ b/pyxem/utils/expt_utils.py @@ -89,8 +89,8 @@ def match_template_dilate( 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) + np.multiply(image_window_sum, image_window_sum, out=image_window_sum) denominator -= image_window_sum denominator *= template_ssd np.maximum(denominator, 0, out=denominator) # sqrt of negative number not allowed From d39ae2a66b7e697db4ba609c5606c28e4a48c484 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 19 Apr 2024 15:38:47 -0500 Subject: [PATCH 10/70] BugFix: Fix normalization bug --- pyxem/utils/expt_utils.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pyxem/utils/expt_utils.py b/pyxem/utils/expt_utils.py index b98640305..f60f35459 100644 --- a/pyxem/utils/expt_utils.py +++ b/pyxem/utils/expt_utils.py @@ -81,7 +81,9 @@ def match_template_dilate( image**2, dilated_template[::-1, ::-1], mode="valid" )[1:-1, 1:-1] - template_mean = template.mean() + 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) @@ -89,8 +91,8 @@ def match_template_dilate( numerator = xcorr - image_window_sum * template_mean denominator = image_window_sum2 - np.divide(image_window_sum, template_volume, out=image_window_sum) 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 From 0ec6839b62e5a13a1d5ef3059ef4f7ecf85358bf Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Sat, 20 Apr 2024 10:47:06 -0500 Subject: [PATCH 11/70] BugFix: Float 64 dtype --- pyxem/utils/expt_utils.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyxem/utils/expt_utils.py b/pyxem/utils/expt_utils.py index f60f35459..e34db82cd 100644 --- a/pyxem/utils/expt_utils.py +++ b/pyxem/utils/expt_utils.py @@ -49,6 +49,7 @@ def match_template_dilate( Value to pad the image with if mode is 'constant' """ + image = image.astype(float, copy=False) if image.ndim < template.ndim: raise ValueError( "Dimensionality of template must be less than or " From 4495a3f9c2a54d26718e5980730910248ae2abb3 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 28 May 2024 09:55:56 -0500 Subject: [PATCH 12/70] Refactor: Move match_template_dilate --- pyxem/utils/diffraction.py | 90 ++++++++++++++++++++++++++++++++++++++ pyxem/utils/expt_utils.py | 90 -------------------------------------- 2 files changed, 90 insertions(+), 90 deletions(-) diff --git a/pyxem/utils/diffraction.py b/pyxem/utils/diffraction.py index e1c38771a..f7e63036f 100644 --- a/pyxem/utils/diffraction.py +++ b/pyxem/utils/diffraction.py @@ -859,3 +859,93 @@ def normalize_template_match(z, template, subtract_min=True, pad_input=True, **k if subtract_min: template_match = template_match - np.min(template_match) return template_match + + +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 preforms very well + for image with different background intensities. This is a slower version of the skimage `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' + """ + + image = image.astype(float, copy=False) + if image.ndim < template.ndim: + 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)): + 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/expt_utils.py b/pyxem/utils/expt_utils.py index e34db82cd..f9ae27f53 100644 --- a/pyxem/utils/expt_utils.py +++ b/pyxem/utils/expt_utils.py @@ -24,93 +24,3 @@ ) from pyxem.utils.diffraction import * - - -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 preforms very well - for image with different background intensities. This is a slower version of the skimage `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' - """ - - image = image.astype(float, copy=False) - if image.ndim < template.ndim: - 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)): - 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)] From d4d1518ffe1b1cff39df1619ee40521599c6adfb Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 28 May 2024 13:42:46 -0500 Subject: [PATCH 13/70] BugFix: Fix labeled vectors to markers. --- pyxem/signals/labeled_diffraction_vectors2d.py | 2 +- pyxem/utils/vectors.py | 18 +++++++++--------- 2 files changed, 10 insertions(+), 10 deletions(-) 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/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 From 6667df6dabf54ebefc0f4cb10bcff916187b1463 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 13:51:27 -0500 Subject: [PATCH 14/70] BugFix: Fix Template Matching with circular window. --- examples/processing/template_matching.py | 28 +++++++++++++++++++----- pyxem/utils/diffraction.py | 7 +++++- 2 files changed, 28 insertions(+), 7 deletions(-) diff --git a/examples/processing/template_matching.py b/examples/processing/template_matching.py index cea589db5..1348a8b48 100644 --- a/examples/processing/template_matching.py +++ b/examples/processing/template_matching.py @@ -10,7 +10,7 @@ import hyperspy.api as hs import pyxem as pxm -import pyxem.dummy_data.make_diffraction_test_data as mdtd +import pyxem.data.dummy_data.make_diffraction_test_data as mdtd s = pxm.data.tilt_boundary_data() @@ -55,7 +55,8 @@ # 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. +# 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=256, image_size_y=256, @@ -65,20 +66,35 @@ ring_x=128, ring_y=128, disk_I=10, - ring_lw=6, + 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 = amorphous_data.template_match_disk(disk_r=6) -template_circular = s.template_match_disk(disk_r=6, dilated_template_window=True) +template_normal.data[:, :, mask] = 0 +template_circular.data[:, :, mask] = 0 +template_normal.data[:, :, mask2] = 0 +template_circular.data[:, :, mask2] = 0 +import matplotlib.pyplot as plt +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", "Large Window"], + label=["Signal", "Square Window", "Circular Window"], tight_layout=True, per_row=3, + vmin=[0, 0.7, 0.7], + fig=f, ) + +# %% diff --git a/pyxem/utils/diffraction.py b/pyxem/utils/diffraction.py index f7e63036f..ab07331f2 100644 --- a/pyxem/utils/diffraction.py +++ b/pyxem/utils/diffraction.py @@ -26,10 +26,12 @@ 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 from skimage.draw import ellipse_perimeter +from skimage._shared.utils import _supported_float_type from skimage.registration import phase_cross_correlation from tqdm import tqdm from packaging.version import Version @@ -855,7 +857,10 @@ 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 From 6e357a5620fad1b3f39fdce4d59d9ccd8eb60b5f Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 14:23:35 -0500 Subject: [PATCH 15/70] Documentation: Add to CHANGELOG.rst --- CHANGELOG.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index d55b9a028..f64133782 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -14,6 +14,10 @@ Fixed ----- - Fixed indexing error in :meth:`~pyxem.signals.Diffraction2D.get_direct_beam_position` (#1080) +Added +----- +- Added `circular_background` to :meth:`~pyxem.signals.Diffraction2D.template_match_disk` to account for + an amorphous circular background when template matching (#1084) 2024-05-08 - version 0.18.0 From b67045e9c09c53558f7ce93d820afb94d98a1f0d Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Sun, 26 May 2024 11:35:02 -0500 Subject: [PATCH 16/70] Addition: Add Extra files for tutorials. --- pyxem/data/__init__.py | 4 +++ pyxem/data/_data.py | 57 ++++++++++++++++++++++++++++++++++++++++- pyxem/data/_registry.py | 2 ++ 3 files changed, 62 insertions(+), 1 deletion(-) diff --git a/pyxem/data/__init__.py b/pyxem/data/__init__.py index 85bbba5e4..adc336f80 100644 --- a/pyxem/data/__init__.py +++ b/pyxem/data/__init__.py @@ -37,6 +37,8 @@ twinned_nanowire, sample_with_g, mgo_nanocrystals, + organic_semiconductor, + cuag_orientations, ) __all__ = [ @@ -47,4 +49,6 @@ "sample_with_g", "mgo_nanocrystals", "tilt_boundary_data", + "cuag_orientations", + "organic_semiconductor", ] diff --git a/pyxem/data/_data.py b/pyxem/data/_data.py index 9774ecb16..251357b69 100644 --- a/pyxem/data/_data.py +++ b/pyxem/data/_data.py @@ -137,7 +137,7 @@ 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 @@ -162,6 +162,61 @@ 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 can be acessed from https://zenodo.org/records/8429302 + 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("cuag_orientations.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 can be acessed from https://zenodo.org/records/8429302 + 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 a larger d + """ + + def mgo_nanocrystals(allow_download=False, **kwargs): # pragma: no cover """A small 4-D STEM dataset of overlapping MgO nanocrystals diff --git a/pyxem/data/_registry.py b/pyxem/data/_registry.py index c8e108fd3..9d500206a 100644 --- a/pyxem/data/_registry.py +++ b/pyxem/data/_registry.py @@ -35,6 +35,8 @@ "sample_with_g.zspy": "md5:65c0a17387a10ca21ebf1812bab67568", "twinned_nanowire.hdf5": "md5:2765fa855db8bc011252dbb425facc81", "ZrNbPercipitate.zspy": "md5:7012a2bdf564f4fb25299af05412723f", + "cuzipProcessed.zspy": "md5:829ecb8f765acb6e9a22092339c6a268", + "colorwheel.txt": "md5:1555136e42cae858be0716f007bda4e4", } # add to _urls and to _file_names_hash # if you want to download the data from a different location From 03cb67878cba01e10587dc4868018e3990554947 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 28 May 2024 14:09:37 -0500 Subject: [PATCH 17/70] New Feature: Add new datasets to data repo from zenodo. --- pyxem/data/__init__.py | 6 +++ pyxem/data/_data.py | 109 ++++++++++++++++++++++++++++++++++++---- pyxem/data/_registry.py | 5 +- 3 files changed, 110 insertions(+), 10 deletions(-) diff --git a/pyxem/data/__init__.py b/pyxem/data/__init__.py index adc336f80..342859fa6 100644 --- a/pyxem/data/__init__.py +++ b/pyxem/data/__init__.py @@ -39,6 +39,9 @@ mgo_nanocrystals, organic_semiconductor, cuag_orientations, + feal_stripes, + sped_ag, + pdcusi_insitu, ) __all__ = [ @@ -51,4 +54,7 @@ "tilt_boundary_data", "cuag_orientations", "organic_semiconductor", + "feal_stripes", + "sped_ag", + "pdcusi_insitu", ] diff --git a/pyxem/data/_data.py b/pyxem/data/_data.py index 251357b69..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 ---------- @@ -139,7 +142,9 @@ 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 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 ---------- @@ -166,7 +171,10 @@ 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 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 ---------- allow_download: bool @@ -181,7 +189,7 @@ def cuag_orientations(allow_download=False, **kwargs): # pragma: no cover """ import zarr - sample = Dataset("cuag_orientations.zspy") + 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) @@ -190,7 +198,10 @@ def cuag_orientations(allow_download=False, **kwargs): # pragma: no cover def organic_semiconductor(allow_download=False, **kwargs): # pragma: no cover """A small 4-D STEM dataset of an organic semiconductor for orientation mapping. - 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 ---------- allow_download: bool @@ -213,14 +224,94 @@ def organic_semiconductor(allow_download=False, **kwargs): # pragma: no cover 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 a larger d + 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 9d500206a..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", @@ -37,6 +37,9 @@ "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 From 81725758e3867f045c918ad3b09778c1c6e2af28 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 28 May 2024 14:45:40 -0500 Subject: [PATCH 18/70] Documentation: Add to CHANGELOG.rst --- CHANGELOG.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index d55b9a028..2fd35a8c6 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -7,6 +7,12 @@ 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 +========== +Added +----- +- Added new datasets of in situ crystalization, Ag SPED, + Organic Semiconductor Orientation mapping, Orientation Mapping, and DPC (#1081) Unreleased ========== From 9f010332d7ab5dec63e926454d133d809abce561 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 31 May 2024 13:48:49 -0500 Subject: [PATCH 19/70] Testing: Add testing for circular Background --- pyxem/tests/signals/test_pixelated_stem_class.py | 7 +++++++ pyxem/utils/diffraction.py | 4 ++-- 2 files changed, 9 insertions(+), 2 deletions(-) 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 ab07331f2..dd7afb460 100644 --- a/pyxem/utils/diffraction.py +++ b/pyxem/utils/diffraction.py @@ -890,12 +890,12 @@ def match_template_dilate( """ image = image.astype(float, copy=False) - if image.ndim < template.ndim: + 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)): + if np.any(np.less(image.shape, template.shape)): # pragma: no cover raise ValueError("Image must be larger than template.") image_shape = image.shape From 1a25638f2a13715a3e893f6a4684d17c991dc8dc Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 31 May 2024 14:42:24 -0500 Subject: [PATCH 20/70] Refactor: Changes per @hakonanes --- examples/processing/template_matching.py | 28 ++++++----- pyxem/signals/diffraction2d.py | 10 ++-- pyxem/utils/diffraction.py | 60 +++++++++++++++++++++--- 3 files changed, 73 insertions(+), 25 deletions(-) diff --git a/examples/processing/template_matching.py b/examples/processing/template_matching.py index 1348a8b48..d3272fef9 100644 --- a/examples/processing/template_matching.py +++ b/examples/processing/template_matching.py @@ -5,18 +5,21 @@ This example shows how the template matching is done in pyxem to find peaks in a diffraction pattern. """ -from skimage.morphology import disk 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. @@ -45,26 +48,26 @@ ) # %% - # 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=256, - image_size_y=256, - disk_x=128, - disk_y=128, - ring_r=45, - ring_x=128, - ring_y=128, + 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, @@ -84,7 +87,6 @@ template_normal.data[:, :, mask2] = 0 template_circular.data[:, :, mask2] = 0 -import matplotlib.pyplot as plt f = plt.figure(figsize=(15, 5)) ind = (5, 5) diff --git a/pyxem/signals/diffraction2d.py b/pyxem/signals/diffraction2d.py index 3297758e1..b9f4ba4ae 100644 --- a/pyxem/signals/diffraction2d.py +++ b/pyxem/signals/diffraction2d.py @@ -1008,13 +1008,12 @@ def template_match_disk(self, disk_r=4, inplace=False, **kwargs): ---------- disk_r : scalar, optional Radius of the disk. Default 4. - circular_background : bool, optional - If True, the windowed normalization will use a circular window lazy_output : bool, default True 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 ------- @@ -1052,8 +1051,6 @@ 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. - circular_background : bool, optional - If True, the windowed normalization will use a circular window lazy_output : bool, default True If True, will return a LazyDiffraction2D object. If False, will compute the result and return a Diffraction2D object. @@ -1062,7 +1059,8 @@ def template_match_ring(self, r_inner=5, r_outer=7, inplace=False, **kwargs): 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 ------- diff --git a/pyxem/utils/diffraction.py b/pyxem/utils/diffraction.py index dd7afb460..1a5042f8e 100644 --- a/pyxem/utils/diffraction.py +++ b/pyxem/utils/diffraction.py @@ -31,7 +31,6 @@ from skimage.feature import match_template from skimage import morphology, filters from skimage.draw import ellipse_perimeter -from skimage._shared.utils import _supported_float_type from skimage.registration import phase_cross_correlation from tqdm import tqdm from packaging.version import Version @@ -50,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. @@ -866,14 +877,46 @@ def normalize_template_match(z, template, subtract_min=True, pad_input=True, **k 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 preforms very well - for image with different background intensities. This is a slower version of the skimage `match_template` - but performs better for images with circular variations in background intensity, specifically accounting - for an amorphous halo around the diffraction pattern. + """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 ---------- @@ -887,6 +930,11 @@ def match_template_dilate( 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) From c826c02c8de956754e8f04fa1ba7162eb82d8375 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 09:02:27 -0500 Subject: [PATCH 21/70] Documentation: Add additional information about averaging neighboring patterns. --- examples/processing/filtering_data.py | 16 ++++++++++++++-- pyxem/signals/diffraction2d.py | 4 ++++ 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/examples/processing/filtering_data.py b/examples/processing/filtering_data.py index e54953d62..93ad10702 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 diff --git a/pyxem/signals/diffraction2d.py b/pyxem/signals/diffraction2d.py index b9f4ba4ae..9d158ba59 100644 --- a/pyxem/signals/diffraction2d.py +++ b/pyxem/signals/diffraction2d.py @@ -1101,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 From 812418072e35cd337ae57f7273cdf1220b5cc71a Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 09:45:02 -0500 Subject: [PATCH 22/70] Documentation: Add Hough Transform example. --- .../processing/circular_hough_transform.py | 72 +++++++++++++++++++ examples/processing/filtering_data.py | 15 ++-- 2 files changed, 77 insertions(+), 10 deletions(-) create mode 100644 examples/processing/circular_hough_transform.py diff --git a/examples/processing/circular_hough_transform.py b/examples/processing/circular_hough_transform.py new file mode 100644 index 000000000..0377982ba --- /dev/null +++ b/examples/processing/circular_hough_transform.py @@ -0,0 +1,72 @@ +""" +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) diff --git a/examples/processing/filtering_data.py b/examples/processing/filtering_data.py index 93ad10702..76c770d88 100644 --- a/examples/processing/filtering_data.py +++ b/examples/processing/filtering_data.py @@ -43,10 +43,8 @@ ) # %% -""" -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): @@ -63,12 +61,9 @@ def custom_filter(array): 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( From 7576cc8bc132f75cad7c90498694f137476133c5 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 10:32:55 -0500 Subject: [PATCH 23/70] Documentation: Add CHANGELOG.rst --- CHANGELOG.rst | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 285b47020..354deaeb1 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -7,24 +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 -========== -Added ------ -- Added new datasets of in situ crystalization, Ag SPED, - Organic Semiconductor Orientation mapping, Orientation Mapping, and DPC (#1081) - 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 ========== From f5d7bca407d6c137cb40dd0cc120ea6fadd88ff1 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 12:01:42 -0500 Subject: [PATCH 24/70] Refactor: Clean up Examples --- .../processing/circular_hough_transform.py | 3 ++ .../processing/determining_ellipticity.py | 38 +++++++++---------- examples/standards/pixel_coodinates.py | 1 + examples/vectors/clustering_vectors.py | 9 +++++ examples/vectors/slicing_vectors.py | 1 + 5 files changed, 31 insertions(+), 21 deletions(-) diff --git a/examples/processing/circular_hough_transform.py b/examples/processing/circular_hough_transform.py index 0377982ba..f10385919 100644 --- a/examples/processing/circular_hough_transform.py +++ b/examples/processing/circular_hough_transform.py @@ -70,3 +70,6 @@ def hough_circle_single_rad(img, radius, **kwargs): 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..d1341bbfc 100644 --- a/examples/processing/determining_ellipticity.py +++ b/examples/processing/determining_ellipticity.py @@ -34,14 +34,13 @@ # %% -""" -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 @@ -57,12 +56,11 @@ 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) @@ -99,14 +97,12 @@ # %% -""" -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) diff --git a/examples/standards/pixel_coodinates.py b/examples/standards/pixel_coodinates.py index b665f442b..65a9d21dc 100644 --- a/examples/standards/pixel_coodinates.py +++ b/examples/standards/pixel_coodinates.py @@ -63,3 +63,4 @@ az.plot() # %% +# sphinx_gallery_thumbnail_number = 4 diff --git a/examples/vectors/clustering_vectors.py b/examples/vectors/clustering_vectors.py index 4e0230a8d..1670314b7 100644 --- a/examples/vectors/clustering_vectors.py +++ b/examples/vectors/clustering_vectors.py @@ -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/slicing_vectors.py b/examples/vectors/slicing_vectors.py index a525511b6..5697ad466 100644 --- a/examples/vectors/slicing_vectors.py +++ b/examples/vectors/slicing_vectors.py @@ -59,3 +59,4 @@ s.add_marker([all_vectors, slic_vectors]) s.add_marker([all_vectors, slic_vectors]) # %% +# sphinx_gallery_thumbnail_number = 8 From 62ab7c900bb2f6cc2cdc0e862af7732780cd18df Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 12:34:27 -0500 Subject: [PATCH 25/70] Refactor: Remove Long Running Examples --- examples/processing/azimuthal_integration.py | 82 +++++++++----------- examples/processing/filtering_data.py | 9 ++- examples/vectors/clustering_vectors.py | 2 +- 3 files changed, 43 insertions(+), 50 deletions(-) diff --git a/examples/processing/azimuthal_integration.py b/examples/processing/azimuthal_integration.py index 7877c9abd..6a1bfa4b5 100644 --- a/examples/processing/azimuthal_integration.py +++ b/examples/processing/azimuthal_integration.py @@ -13,95 +13,87 @@ 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() +s1d.sum().plot() # %% -""" -Similarly, the `get_azimuthal_integral2d` method will return a `Polar2D` signal. -""" +# 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() +s_polar = s.get_azimuthal_integral2d(npt=100, npt_azim=360, inplace=False) +s_polar.sum().plot() # %% -""" -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. +# 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. +# 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! +# 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! -""" +# We only show the 1D case here, but the same applies for the 2D case as well! -nano_crystals.calibration.detector( +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/filtering_data.py b/examples/processing/filtering_data.py index 76c770d88..3ce0f8322 100644 --- a/examples/processing/filtering_data.py +++ b/examples/processing/filtering_data.py @@ -25,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 @@ -36,7 +36,7 @@ ) # 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", @@ -55,7 +55,7 @@ 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", @@ -71,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/vectors/clustering_vectors.py b/examples/vectors/clustering_vectors.py index 1670314b7..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) From dd05e44e10fa1e397af767722db6632c83d3d943 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 14:32:42 -0500 Subject: [PATCH 26/70] Documentation: Fix Formatting --- examples/processing/azimuthal_integration.py | 15 ++++++--------- examples/processing/centering_the_zero_beam.py | 5 ++++- examples/processing/circular_hough_transform.py | 4 ---- examples/processing/determining_ellipticity.py | 9 ++++----- examples/processing/vector_finding.py | 4 +--- examples/standards/pixel_coodinates.py | 3 --- examples/vectors/glass_symmetry.py | 9 ++++++++- examples/vectors/slicing_vectors.py | 2 ++ .../creating_a_segmented_detector.py | 3 ++- .../virtual_imaging/interactive_virtual_images.py | 3 ++- examples/virtual_imaging/other_virtual_imaging.py | 2 +- 11 files changed, 30 insertions(+), 29 deletions(-) diff --git a/examples/processing/azimuthal_integration.py b/examples/processing/azimuthal_integration.py index 6a1bfa4b5..02772307a 100644 --- a/examples/processing/azimuthal_integration.py +++ b/examples/processing/azimuthal_integration.py @@ -22,27 +22,26 @@ s1d = s.get_azimuthal_integral1d(npt=100, inplace=False) s1d.sum().plot() -# %% +# %% # Similarly, the `get_azimuthal_integral2d` method will return a `Polar2D` signal. s_polar = s.get_azimuthal_integral2d(npt=100, npt_azim=360, inplace=False) s_polar.sum().plot() # %% - # 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( @@ -67,12 +66,10 @@ 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. @@ -82,8 +79,8 @@ ) # Just a random affine transformation for illustration 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. 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 index f10385919..a3fbd1d2a 100644 --- a/examples/processing/circular_hough_transform.py +++ b/examples/processing/circular_hough_transform.py @@ -24,7 +24,6 @@ # %% - # Canny Filter # --------------------- # First we have to apply a Canny filter to the dataset to get a binary image of the @@ -42,8 +41,6 @@ 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 @@ -60,7 +57,6 @@ def hough_circle_single_rad(img, radius, **kwargs): circular_hough.plot() # %% - # Finding Peaks # ------------- # Finding peaks is fairly easy from this point and uses the ``get_diffraction_vectors`` function. diff --git a/examples/processing/determining_ellipticity.py b/examples/processing/determining_ellipticity.py index d1341bbfc..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,7 +33,6 @@ s.plot() # %% - # Using RANSAC Ellipse Determination # ---------------------------------- # The RANSAC algorithm is a robust method for fitting an ellipse to points with outliers. Here we @@ -54,8 +53,8 @@ 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 @@ -95,8 +94,6 @@ colorbar=None, ) # %% - - # Ellipse from Points # ------------------- # This workflow will likely change slightly in the future to add better parllelism, @@ -114,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/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 65a9d21dc..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 # --------------- # 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 5697ad466..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,5 +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) From a7cedeb4053731a78438c18d9a5e03993c1f71a4 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Thu, 8 Feb 2024 05:55:50 -0600 Subject: [PATCH 27/70] Refactor: Moved orientation mapping to polar_diffraction2d.py --- pyxem/signals/polar_diffraction2d.py | 49 ++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/pyxem/signals/polar_diffraction2d.py b/pyxem/signals/polar_diffraction2d.py index 502581de0..ad62903ca 100644 --- a/pyxem/signals/polar_diffraction2d.py +++ b/pyxem/signals/polar_diffraction2d.py @@ -23,6 +23,10 @@ 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, +) from pyxem.utils._background_subtraction import ( _polar_subtract_radial_median, @@ -304,6 +308,51 @@ def subtract_diffraction_background( **kwargs, ) + def get_orientation( + self, simulation, n_keep=1, frac_keep=0.1, n_best=1, normalize_templates=True + ): + """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. + + Returns + ------- + orientation : BaseSignal + A signal with the orientation at each navigation position. + """ + ( + r_templates, + theta_templates, + intensities_templates, + ) = simulation.get_polar_templates() + radius = self.axes_manager.signal_axes[0].size # number radial pixels + integrated_templates = _get_integrated_polar_templates( + radius, r_templates, intensities_templates, normalize_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, + ) + + return orientation + class LazyPolarDiffraction2D(LazySignal, PolarDiffraction2D): pass From dc77b799a72d1d697e85ca2cf4bf70e22b5e1bdc Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 14 Feb 2024 09:08:57 -0600 Subject: [PATCH 28/70] New Feature: Added OrientationMap class --- pyxem/signals/indexation_results.py | 53 ++++++++++++++++++++++++++++ pyxem/signals/polar_diffraction2d.py | 13 ++++++- 2 files changed, 65 insertions(+), 1 deletion(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index cc6d78a49..9b7b07286 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -26,6 +26,7 @@ from transforms3d.euler import mat2euler 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 @@ -147,6 +148,58 @@ def _get_second_best_phase(z): return -1 +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" + + def __init__(self): + super().__init__() + self._signal_type = "orientation_map" + + @property + def simulation(self): + return self.metadata.get_item("simulation") + + @simulation.setter + def simulation(self, value): + self.metadata.set_item("simulation", value) + + def to_crystal_map(self): + pass + + def to_markers(self): + pass + + def to_navigator(self): + pass + + def plot_over_signal(self, annotate=False, **kwargs): + """ + Parameters + ---------- + annotate + kwargs + + Returns + ------- + + """ + pass + + class GenericMatchingResults: def __init__(self, data): self.data = hs.signals.Signal2D(data) diff --git a/pyxem/signals/polar_diffraction2d.py b/pyxem/signals/polar_diffraction2d.py index ad62903ca..0ae5f956f 100644 --- a/pyxem/signals/polar_diffraction2d.py +++ b/pyxem/signals/polar_diffraction2d.py @@ -324,6 +324,14 @@ def get_orientation( ---------- 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. Returns ------- @@ -334,7 +342,10 @@ def get_orientation( r_templates, theta_templates, intensities_templates, - ) = simulation.get_polar_templates() + ) = simulation.polar_flatten_simulations( + radial_axes=self.axes_manager.signal_axes[0].axis, + azimuthal_axes=self.axes_manager.signal_axes[1].axis, + ) radius = self.axes_manager.signal_axes[0].size # number radial pixels integrated_templates = _get_integrated_polar_templates( radius, r_templates, intensities_templates, normalize_templates From ad9eb9be2e732e7d45ba9e8eb68834a4836d7f70 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 14 Feb 2024 09:21:45 -0600 Subject: [PATCH 29/70] Refactor: Clean up `OrientationMap` class --- pyxem/hyperspy_extension.yaml | 12 +++++++++ pyxem/signals/indexation_results.py | 39 ++++++++++++++++++++++------ pyxem/signals/polar_diffraction2d.py | 4 +++ 3 files changed, 47 insertions(+), 8 deletions(-) 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/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 9b7b07286..97b502739 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -19,6 +19,7 @@ from warnings import warn import hyperspy.api as hs +from hyperspy._signals.lazy import LazySignal from hyperspy.signal import BaseSignal import numpy as np from orix.crystal_map import CrystalMap @@ -178,27 +179,45 @@ def simulation(self, value): self.metadata.set_item("simulation", value) def to_crystal_map(self): + """Convert the orientation map to an `orix.CrystalMap` object""" pass - def to_markers(self): + def to_markers(self, annotate=False, **kwargs): + """Convert the orientation map to a set of markers for plotting. + + Parameters + ---------- + annotate : bool + If True, the euler rotation and the correlation will be annotated on the plot using + the `Texts` class from hyperspy. + """ + pass def to_navigator(self): + """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. + """ pass - def plot_over_signal(self, annotate=False, **kwargs): - """ + def plot_over_signal(self, signal, annotate=False, **kwargs): + """Convenience method to plot the orientation map and the n-best matches over the signal. + Parameters ---------- - annotate - kwargs - - Returns - ------- + signal : BaseSignal + The signal to plot the orientation map over. + annotate: bool + If True, the euler rotation and the correlation will be annotated on the plot using + the `Texts` class from hyperspy. """ pass + def plot_inplane_rotation(self, **kwargs): + """Plot the in-plane rotation of the orientation map as a 2D map.""" + pass + class GenericMatchingResults: def __init__(self, data): @@ -252,6 +271,10 @@ def to_crystal_map(self): ) +class LazyOrientationMap(LazySignal, OrientationMap): + pass + + class VectorMatchingResults(BaseSignal): """Vector matching results containing the top n best matching crystal phase and orientation at each navigation position with associated metrics. diff --git a/pyxem/signals/polar_diffraction2d.py b/pyxem/signals/polar_diffraction2d.py index 0ae5f956f..153c5869a 100644 --- a/pyxem/signals/polar_diffraction2d.py +++ b/pyxem/signals/polar_diffraction2d.py @@ -362,6 +362,10 @@ def get_orientation( inplace=False, ) + 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 From d6dc3a0b05bbd3807df6e4b1daa6eb4009b976b1 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 14 Feb 2024 09:25:13 -0600 Subject: [PATCH 30/70] Bugfix: Fix `OrientationMap` error --- pyxem/signals/indexation_results.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 97b502739..92f7c5288 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -166,10 +166,6 @@ class OrientationMap(DiffractionVectors2D): _signal_type = "orientation_map" - def __init__(self): - super().__init__() - self._signal_type = "orientation_map" - @property def simulation(self): return self.metadata.get_item("simulation") From 7f4afc9b43d99cdfb37a16cd4d8a08a046d29cf0 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 14 Feb 2024 15:10:29 -0600 Subject: [PATCH 31/70] Bugfix: transpose polar image to match --- pyxem/signals/polar_diffraction2d.py | 12 ++++++++---- pyxem/utils/indexation_utils.py | 21 +++++++++++++++++++++ 2 files changed, 29 insertions(+), 4 deletions(-) diff --git a/pyxem/signals/polar_diffraction2d.py b/pyxem/signals/polar_diffraction2d.py index 153c5869a..669b6f985 100644 --- a/pyxem/signals/polar_diffraction2d.py +++ b/pyxem/signals/polar_diffraction2d.py @@ -26,6 +26,7 @@ from pyxem.utils.indexation_utils import ( _mixed_matching_lib_to_polar, _get_integrated_polar_templates, + _norm_rows, ) from pyxem.utils._background_subtraction import ( @@ -309,7 +310,7 @@ def subtract_diffraction_background( ) def get_orientation( - self, simulation, n_keep=1, frac_keep=0.1, n_best=1, normalize_templates=True + self, simulation, n_keep=None, frac_keep=0.1, n_best=1, normalize_templates=True ): """Match the orientation with some simulated diffraction patterns using an accelerated orientation mapping algorithm. @@ -343,13 +344,15 @@ def get_orientation( theta_templates, intensities_templates, ) = simulation.polar_flatten_simulations( - radial_axes=self.axes_manager.signal_axes[0].axis, - azimuthal_axes=self.axes_manager.signal_axes[1].axis, + radial_axes=self.axes_manager.signal_axes[1].axis, + azimuthal_axes=self.axes_manager.signal_axes[0].axis, ) - radius = self.axes_manager.signal_axes[0].size # number radial pixels + 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, @@ -360,6 +363,7 @@ def get_orientation( frac_keep=frac_keep, n_best=n_best, inplace=False, + transpose=True, ) orientation.set_signal_type("orientation_map") diff --git a/pyxem/utils/indexation_utils.py b/pyxem/utils/indexation_utils.py index a77e117c7..ef77a3ec5 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,8 @@ 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 dispatcher = get_array_module(polar_image) # remove templates we don't care about with a fast match ( @@ -1099,6 +1102,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) From 9ad2e38f3f4271af1d6ed84c7c0baa7698538223 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Sat, 17 Feb 2024 17:20:15 -0600 Subject: [PATCH 32/70] Testing: Add tests for orientation mapping --- .../tests/signals/test_polar_diffraction2d.py | 23 +++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/pyxem/tests/signals/test_polar_diffraction2d.py b/pyxem/tests/signals/test_polar_diffraction2d.py index 639b7ec1c..50d76e868 100644 --- a/pyxem/tests/signals/test_polar_diffraction2d.py +++ b/pyxem/tests/signals/test_polar_diffraction2d.py @@ -251,6 +251,29 @@ def test_decomposition_class_assignment(self, diffraction_pattern): assert isinstance(s, PolarDiffraction2D) +class TestOrientationMap: + @pytest.fixture + def polar_image(self, diffraction_pattern): + from pyxem.signals import PolarDiffraction2D + + pol_image = np.zeros((2, 2, 20, 60)) + pol_image[:, :, 4:6, 9:11] = 1 + pol_image[:, :, 4:6, 29:31] = 1 + pol_image[:, :, 4:6, 49:51] = 1 + + pol_image[:, :, 9:11, 9:11] = 1 + pol_image[:, :, 9:11, 29:31] = 1 + pol_image[:, :, 9:11, 49:51] = 1 + pol = PolarDiffraction2D(pol_image) + pol.axes_manager[3].scale = 0.5 + return pol + + def test_orientation_map_inplace(self, diffraction_pattern): + s = PolarDiffraction2D(diffraction_pattern) + s.orientation_map(inplace=True) + assert s.learning_results is not None + + class TestSubtractingDiffractionBackground: @pytest.fixture def noisy_data(self): From bf23e6886237581717b5af1dfc2e35de548bfe6a Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Mon, 26 Feb 2024 08:35:04 -0600 Subject: [PATCH 33/70] Refactor: Avoid testing for signal size --- pyxem/signals/polar_diffraction2d.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/pyxem/signals/polar_diffraction2d.py b/pyxem/signals/polar_diffraction2d.py index 669b6f985..4f23369ef 100644 --- a/pyxem/signals/polar_diffraction2d.py +++ b/pyxem/signals/polar_diffraction2d.py @@ -310,7 +310,13 @@ def subtract_diffraction_background( ) def get_orientation( - self, simulation, n_keep=None, frac_keep=0.1, n_best=1, normalize_templates=True + self, + simulation, + n_keep=None, + frac_keep=0.1, + n_best=1, + normalize_templates=True, + **kwargs ): """Match the orientation with some simulated diffraction patterns using an accelerated orientation mapping algorithm. @@ -333,6 +339,8 @@ def get_orientation( The number of best matching orientations to keep. normalize_templates : bool Normalize the templates to the same intensity. + kwargs : dict + Any additional options for the :meth:`~hyperspy.signal.BaseSignal.map` function. Returns ------- @@ -364,6 +372,8 @@ def get_orientation( n_best=n_best, inplace=False, transpose=True, + output_signal_size=(n_best, 4), + output_dtype=float, ) orientation.set_signal_type("orientation_map") From 15260f8ab8bd44789735002b1734fc1f42220fa6 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 7 May 2024 17:09:52 -0500 Subject: [PATCH 34/70] BugFix: Fix regression with rebase --- pyxem/signals/polar_diffraction2d.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/pyxem/signals/polar_diffraction2d.py b/pyxem/signals/polar_diffraction2d.py index 4f23369ef..bff6aa80a 100644 --- a/pyxem/signals/polar_diffraction2d.py +++ b/pyxem/signals/polar_diffraction2d.py @@ -316,17 +316,14 @@ def get_orientation( frac_keep=0.1, n_best=1, normalize_templates=True, - **kwargs + **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 @@ -341,7 +338,6 @@ def get_orientation( Normalize the templates to the same intensity. kwargs : dict Any additional options for the :meth:`~hyperspy.signal.BaseSignal.map` function. - Returns ------- orientation : BaseSignal From 22d7fb3243ac868a5ae9cb292e9faeb4f4a36858 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 7 May 2024 17:22:58 -0500 Subject: [PATCH 35/70] BugFix: Remove test --- pyxem/tests/signals/test_polar_diffraction2d.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/pyxem/tests/signals/test_polar_diffraction2d.py b/pyxem/tests/signals/test_polar_diffraction2d.py index 50d76e868..e708356ca 100644 --- a/pyxem/tests/signals/test_polar_diffraction2d.py +++ b/pyxem/tests/signals/test_polar_diffraction2d.py @@ -268,11 +268,6 @@ def polar_image(self, diffraction_pattern): pol.axes_manager[3].scale = 0.5 return pol - def test_orientation_map_inplace(self, diffraction_pattern): - s = PolarDiffraction2D(diffraction_pattern) - s.orientation_map(inplace=True) - assert s.learning_results is not None - class TestSubtractingDiffractionBackground: @pytest.fixture From 178a815558b00c626871ff318c3eec1856daf652 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 14 Feb 2024 09:08:57 -0600 Subject: [PATCH 36/70] New Feature: Added OrientationMap class --- pyxem/signals/indexation_results.py | 38 +++++++++-------------------- 1 file changed, 12 insertions(+), 26 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 92f7c5288..d01622ecb 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -166,6 +166,10 @@ class OrientationMap(DiffractionVectors2D): _signal_type = "orientation_map" + def __init__(self): + super().__init__() + self._signal_type = "orientation_map" + @property def simulation(self): return self.metadata.get_item("simulation") @@ -175,43 +179,25 @@ def simulation(self, value): self.metadata.set_item("simulation", value) def to_crystal_map(self): - """Convert the orientation map to an `orix.CrystalMap` object""" pass - def to_markers(self, annotate=False, **kwargs): - """Convert the orientation map to a set of markers for plotting. - - Parameters - ---------- - annotate : bool - If True, the euler rotation and the correlation will be annotated on the plot using - the `Texts` class from hyperspy. - """ - + def to_markers(self): pass def to_navigator(self): - """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. - """ pass - def plot_over_signal(self, signal, annotate=False, **kwargs): - """Convenience method to plot the orientation map and the n-best matches over the signal. - + def plot_over_signal(self, annotate=False, **kwargs): + """ Parameters ---------- - signal : BaseSignal - The signal to plot the orientation map over. - annotate: bool - If True, the euler rotation and the correlation will be annotated on the plot using - the `Texts` class from hyperspy. + annotate + kwargs - """ - pass + Returns + ------- - def plot_inplane_rotation(self, **kwargs): - """Plot the in-plane rotation of the orientation map as a 2D map.""" + """ pass From 44b0faa03324725175d35fe356a2f7975bfa79c7 Mon Sep 17 00:00:00 2001 From: Viljar Femoen Date: Sat, 20 Apr 2024 11:21:03 +0200 Subject: [PATCH 37/70] Implement `to_markers` --- pyxem/signals/indexation_results.py | 116 +++++++++++++++++++++++++--- 1 file changed, 107 insertions(+), 9 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index d01622ecb..cfff2f01b 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -179,27 +179,125 @@ def simulation(self, value): self.metadata.set_item("simulation", value) def to_crystal_map(self): + """Convert the orientation map to an `orix.CrystalMap` object""" pass - def to_markers(self): - pass + def to_markers(self, n_best: int = 1, annotate=False, **kwargs): + """Convert the orientation map to a set of markers for plotting. + + Parameters + ---------- + annotate : bool + If True, the euler rotation and the correlation will be annotated on the plot using + the `Texts` class from hyperspy. + """ + def marker_generator_factory(n_best_entry: int): + def marker_generator(entry): + # Get data + index, correlation, rotation, factor = entry[n_best_entry] + # Get coordinates of reflections + _, _, coords = self.simulation.get_simulation(int(index)) + # Mirror data if necessary + coords.data[:, 1] *= factor + # Rotation matrix for the in-plane rotation + T = Rotation.from_euler((rotation, 0, 0), degrees=True).to_matrix().squeeze() + coords = coords.data @ T + # x and y needs to swap, and we don't want z. Therefore, use slice(1, 0, -1) + return coords[:, 1::-1] + return marker_generator + + def reciprocal_lattice_vector_to_text(vec): + def add_bar(i: int) -> str: + if i < 0: + return f"$\\bar{{{abs(i)}}}$" + else: + return f"{i}" + out = [] + for hkl in vec.hkl: + hkl = np.round(hkl).astype(np.int16) + out.append(f"({add_bar(hkl[0])} {add_bar(hkl[1])} {add_bar(hkl[2])})") + return out + + def text_generator_factory(n_best_entry: int): + def text_generator(entry): + # Get data + index, correlation, rotation, factor = entry[n_best_entry] + _, _, vecs = self.simulation.get_simulation(int(index)) + return reciprocal_lattice_vector_to_text(vecs) + return text_generator + + for n in range(n_best): + markers_signal = self.map( + marker_generator_factory(n), + inplace=False, + ragged=True, + lazy_output=True, + ) + markers = hs.plot.markers.Points.from_signal(markers_signal, color="red") + yield markers + if annotate: + text_signal = self.map( + text_generator_factory(n), + inplace=False, + ragged=True, + lazy_output=False, + ) + texts = np.empty(self.axes_manager.navigation_shape, dtype=object) + for i in range(texts.shape[0]): + for j in range(texts.shape[1]): + texts[i, j] = ["a", str(i) + str(j), "a b c"] + text_markers = hs.plot.markers.Texts.from_signal(markers_signal, texts=np.swapaxes(text_signal.data, 0, 1), color="black") + yield text_markers + + + def to_polar_markers(self, n_best: int = 1): + r_templates, theta_templates, intensities_templates = self.simulation.polar_flatten_simulations() + + def marker_generator_factory(n_best_entry: int): + def marker_generator(entry): + index, correlation, rotation, factor = entry[n_best_entry] + r = r_templates[int(index)] + theta = theta_templates[int(index)] + theta += 2*np.pi + np.deg2rad(rotation) + theta %= 2*np.pi + theta -= np.pi + return np.array((theta, r)).T + return marker_generator + + for n in range(n_best): + markers_signal = self.map( + marker_generator_factory(n), + inplace=False, + ragged=True, + lazy_output=True, + ) + markers = hs.plot.markers.Points.from_signal(markers_signal) + yield markers def to_navigator(self): + """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. + """ pass - def plot_over_signal(self, annotate=False, **kwargs): - """ + def plot_over_signal(self, signal, annotate=False, **kwargs): + """Convenience method to plot the orientation map and the n-best matches over the signal. + Parameters ---------- - annotate - kwargs - - Returns - ------- + signal : BaseSignal + The signal to plot the orientation map over. + annotate: bool + If True, the euler rotation and the correlation will be annotated on the plot using + the `Texts` class from hyperspy. """ pass + def plot_inplane_rotation(self, **kwargs): + """Plot the in-plane rotation of the orientation map as a 2D map.""" + pass + class GenericMatchingResults: def __init__(self, data): From adb33a1a155e0e67017317112d8b97af1af3cd60 Mon Sep 17 00:00:00 2001 From: Viljar Femoen Date: Sat, 20 Apr 2024 11:21:03 +0200 Subject: [PATCH 38/70] Type hints and colors --- pyxem/signals/indexation_results.py | 30 ++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index cfff2f01b..31cc1eaba 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -17,6 +17,7 @@ # along with pyXem. If not, see . from warnings import warn +from typing import Iterator, Union, Literal import hyperspy.api as hs from hyperspy._signals.lazy import LazySignal @@ -24,6 +25,9 @@ import numpy as np from orix.crystal_map import CrystalMap from orix.quaternion import Rotation +from orix.vector import Vector3d +from orix.plot import IPFColorKeyTSL +from diffsims.simulations import Simulation2D from transforms3d.euler import mat2euler from pyxem.utils.indexation_utils import get_nth_best_solution @@ -171,7 +175,7 @@ def __init__(self): self._signal_type = "orientation_map" @property - def simulation(self): + def simulation(self) -> Simulation2D: return self.metadata.get_item("simulation") @simulation.setter @@ -182,7 +186,7 @@ def to_crystal_map(self): """Convert the orientation map to an `orix.CrystalMap` object""" pass - def to_markers(self, n_best: int = 1, annotate=False, **kwargs): + def to_markers(self, n_best: int = 1, annotate=False, marker_color: str = "red", text_color: str = "black"): """Convert the orientation map to a set of markers for plotting. Parameters @@ -190,6 +194,10 @@ def to_markers(self, n_best: int = 1, annotate=False, **kwargs): 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`. """ def marker_generator_factory(n_best_entry: int): def marker_generator(entry): @@ -233,7 +241,7 @@ def text_generator(entry): ragged=True, lazy_output=True, ) - markers = hs.plot.markers.Points.from_signal(markers_signal, color="red") + markers = hs.plot.markers.Points.from_signal(markers_signal, color=marker_color) yield markers if annotate: text_signal = self.map( @@ -242,15 +250,11 @@ def text_generator(entry): ragged=True, lazy_output=False, ) - texts = np.empty(self.axes_manager.navigation_shape, dtype=object) - for i in range(texts.shape[0]): - for j in range(texts.shape[1]): - texts[i, j] = ["a", str(i) + str(j), "a b c"] - text_markers = hs.plot.markers.Texts.from_signal(markers_signal, texts=np.swapaxes(text_signal.data, 0, 1), color="black") + text_markers = hs.plot.markers.Texts.from_signal(markers_signal, texts=text_signal.data.T, color=text_color) yield text_markers - def to_polar_markers(self, n_best: int = 1): + def to_polar_markers(self, n_best: int = 1) -> Iterator[hs.plot.markers.Markers]: r_templates, theta_templates, intensities_templates = self.simulation.polar_flatten_simulations() def marker_generator_factory(n_best_entry: int): @@ -290,9 +294,13 @@ def plot_over_signal(self, signal, annotate=False, **kwargs): annotate: bool If True, the euler rotation and the correlation will be annotated on the plot using the `Texts` class from hyperspy. - + + Notes + ----- + The kwargs are passed to the `signal.plot` function call """ - pass + signal.plot(**kwargs) + signal.add_marker(self.to_markers(1, annotate=annotate)) def plot_inplane_rotation(self, **kwargs): """Plot the in-plane rotation of the orientation map as a 2D map.""" From 7e538f7ccbeebeefd9fea8f563c48b522505194d Mon Sep 17 00:00:00 2001 From: Viljar Femoen Date: Sat, 20 Apr 2024 11:42:37 +0200 Subject: [PATCH 39/70] Implement to_navigator and to_single_phase_orientations --- pyxem/signals/indexation_results.py | 49 +++++++++++++++++++++++++---- 1 file changed, 43 insertions(+), 6 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 31cc1eaba..ecd752314 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -23,11 +23,10 @@ from hyperspy._signals.lazy import LazySignal from hyperspy.signal import BaseSignal import numpy as np -from orix.crystal_map import CrystalMap -from orix.quaternion import Rotation +from orix.crystal_map import CrystalMap, Phase +from orix.quaternion import Rotation, Orientation from orix.vector import Vector3d from orix.plot import IPFColorKeyTSL -from diffsims.simulations import Simulation2D from transforms3d.euler import mat2euler from pyxem.utils.indexation_utils import get_nth_best_solution @@ -182,7 +181,36 @@ def simulation(self) -> Simulation2D: def simulation(self, value): self.metadata.set_item("simulation", value) - def to_crystal_map(self): + def to_single_phase_orientations(self) -> Orientation: + """Convert the orientation map to an `Orientation`-object, + given a single-phase simulation. + """ + if not isinstance(self.simulation.phases, Phase): + raise ValueError("Multiple phases found in simulation") + + # 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) + + def rotation_from_orientation_map(result, rotations): + index, _, rotation, mirror = result.T + index = index.astype(int) + ori = rotations[index] + euler = Orientation(ori).to_euler(degrees=True) * mirror[..., np.newaxis] + euler[:, 0] = rotation + ori = Orientation.from_euler(euler, degrees=True).data + return ori + + return Orientation( + self.map( + rotation_from_orientation_map, + rotations=rotations, + inplace=False, + ), + symmetry=self.simulation.phases.point_group, + ) + + def to_crystal_map(self) -> CrystalMap: """Convert the orientation map to an `orix.CrystalMap` object""" pass @@ -278,11 +306,20 @@ def marker_generator(entry): markers = hs.plot.markers.Points.from_signal(markers_signal) yield markers - def to_navigator(self): + def to_navigator(self, direction: Vector3d = Vector3d.zvector()): """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. """ - pass + 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") + + return s def plot_over_signal(self, signal, annotate=False, **kwargs): """Convenience method to plot the orientation map and the n-best matches over the signal. From b2031ca40ea5a16912a28df9cbcb60e396928ab2 Mon Sep 17 00:00:00 2001 From: Viljar Femoen Date: Sat, 20 Apr 2024 11:44:00 +0200 Subject: [PATCH 40/70] Change units of in-plane angle to degrees --- pyxem/signals/polar_diffraction2d.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/pyxem/signals/polar_diffraction2d.py b/pyxem/signals/polar_diffraction2d.py index bff6aa80a..7d4fc1758 100644 --- a/pyxem/signals/polar_diffraction2d.py +++ b/pyxem/signals/polar_diffraction2d.py @@ -19,6 +19,7 @@ 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 @@ -372,6 +373,20 @@ def get_orientation( output_dtype=float, ) + # Translate in-plane rotation from index to degrees + # by using the calibration of the axis + def rotation_index_to_degrees(data): + data = data.copy() + axis = self.axes_manager.signal_axes[0].axis + ind = data[:, 2].astype(int) + data[:, 2] = rad2deg(axis[ind]) + # Account for the 90 degree discrepancy between the X_L axis + # and the origin of the azimuthal axis of the polar signal + data[:, 2] -= 90 + return data + + orientation.map(rotation_index_to_degrees) + orientation.set_signal_type("orientation_map") orientation.simulation = simulation orientation.column_names = ["index", "correlation", "rotation", "factor"] From 5e21eb20caaabda13287c97c4bca85baa359bda4 Mon Sep 17 00:00:00 2001 From: Viljar Femoen Date: Sat, 20 Apr 2024 15:25:27 +0200 Subject: [PATCH 41/70] Refactor to_marker --- pyxem/signals/indexation_results.py | 134 ++++++++++++++++++---------- 1 file changed, 85 insertions(+), 49 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index ecd752314..a4937d181 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -28,6 +28,7 @@ from orix.vector import Vector3d from orix.plot import IPFColorKeyTSL from transforms3d.euler import mat2euler +from diffsims.crystallography import ReciprocalLatticeVector from pyxem.utils.indexation_utils import get_nth_best_solution from pyxem.signals.diffraction_vectors2d import DiffractionVectors2D @@ -174,7 +175,7 @@ def __init__(self): self._signal_type = "orientation_map" @property - def simulation(self) -> Simulation2D: + def simulation(self): return self.metadata.get_item("simulation") @simulation.setter @@ -185,7 +186,7 @@ def to_single_phase_orientations(self) -> Orientation: """Convert the orientation map to an `Orientation`-object, given a single-phase simulation. """ - if not isinstance(self.simulation.phases, Phase): + if self.simulation.has_multiple_phases: raise ValueError("Multiple phases found in simulation") # Use the quaternion data from rotations to support 2D rotations, @@ -210,15 +211,65 @@ def rotation_from_orientation_map(result, rotations): symmetry=self.simulation.phases.point_group, ) + def to_single_phase_vectors(self, n_best_index: int = 0) -> 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 + """ + + 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) + + def extract_vectors_from_orientation_map(result, all_vectors): + index, _, rotation, mirror = result[n_best_index, :].T + index = index.astype(int) + vectors = all_vectors[index] + # Copy manually, as deepcopy adds a lot of overhead with the phase + vectors = ReciprocalLatticeVector(vectors.phase, xyz=vectors.data.copy()) + # Flip y, as discussed in https://github.com/pyxem/pyxem/issues/925 + vectors.y = -vectors.y + # Mirror if necessary + vectors.y = mirror * vectors.y + rotation_matrix = ( + Rotation.from_euler( + (rotation, 0, 0), degrees=True, direction="crystal2lab" + ) + ).to_matrix() + vectors = vectors.rotate_from_matrix(rotation_matrix.squeeze()) + + return vectors + + return self.map( + extract_vectors_from_orientation_map, + all_vectors=vectors_signal, + inplace=False, + ) + def to_crystal_map(self) -> CrystalMap: """Convert the orientation map to an `orix.CrystalMap` object""" pass - def to_markers(self, n_best: int = 1, annotate=False, marker_color: str = "red", text_color: str = "black"): + def to_single_phase_markers( + self, + n_best: int = 1, + annotate: bool = False, + marker_color: str = "red", + text_color: str = "black", + lazy: bool = True, + ): """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. @@ -227,73 +278,57 @@ def to_markers(self, n_best: int = 1, annotate=False, marker_color: str = "red", text_color: str, optional The color used for the text annotations for reflections. Does nothing if `annotate` is `False`. """ - def marker_generator_factory(n_best_entry: int): - def marker_generator(entry): - # Get data - index, correlation, rotation, factor = entry[n_best_entry] - # Get coordinates of reflections - _, _, coords = self.simulation.get_simulation(int(index)) - # Mirror data if necessary - coords.data[:, 1] *= factor - # Rotation matrix for the in-plane rotation - T = Rotation.from_euler((rotation, 0, 0), degrees=True).to_matrix().squeeze() - coords = coords.data @ T - # x and y needs to swap, and we don't want z. Therefore, use slice(1, 0, -1) - return coords[:, 1::-1] - return marker_generator - def reciprocal_lattice_vector_to_text(vec): + def vectors_to_coordinates(vectors): + return np.vstack((vectors.x, vectors.y)).T + + def vectors_to_text(vectors): def add_bar(i: int) -> str: if i < 0: return f"$\\bar{{{abs(i)}}}$" else: return f"{i}" + out = [] - for hkl in vec.hkl: - hkl = np.round(hkl).astype(np.int16) - out.append(f"({add_bar(hkl[0])} {add_bar(hkl[1])} {add_bar(hkl[2])})") + 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 text_generator_factory(n_best_entry: int): - def text_generator(entry): - # Get data - index, correlation, rotation, factor = entry[n_best_entry] - _, _, vecs = self.simulation.get_simulation(int(index)) - return reciprocal_lattice_vector_to_text(vecs) - return text_generator - for n in range(n_best): - markers_signal = self.map( - marker_generator_factory(n), - inplace=False, - ragged=True, - lazy_output=True, + vectors = self.to_single_phase_vectors(n) + coords = vectors.map( + vectors_to_coordinates, inplace=False, lazy_output=lazy, ragged=True ) - markers = hs.plot.markers.Points.from_signal(markers_signal, color=marker_color) + markers = hs.plot.markers.Points.from_signal(coords, color=marker_color) yield markers + if annotate: - text_signal = self.map( - text_generator_factory(n), - inplace=False, - ragged=True, - lazy_output=False, + texts = vectors.map( + vectors_to_text, inplace=False, lazy_output=lazy, ragged=True + ) + text_markers = hs.plot.markers.Texts.from_signal( + coords, texts=texts.data.T, color=text_color ) - text_markers = hs.plot.markers.Texts.from_signal(markers_signal, texts=text_signal.data.T, color=text_color) yield text_markers - def to_polar_markers(self, n_best: int = 1) -> Iterator[hs.plot.markers.Markers]: - r_templates, theta_templates, intensities_templates = self.simulation.polar_flatten_simulations() + ( + r_templates, + theta_templates, + intensities_templates, + ) = self.simulation.polar_flatten_simulations() def marker_generator_factory(n_best_entry: int): def marker_generator(entry): index, correlation, rotation, factor = entry[n_best_entry] r = r_templates[int(index)] theta = theta_templates[int(index)] - theta += 2*np.pi + np.deg2rad(rotation) - theta %= 2*np.pi + theta += 2 * np.pi + np.deg2rad(rotation) + theta %= 2 * np.pi theta -= np.pi return np.array((theta, r)).T + return marker_generator for n in range(n_best): @@ -321,7 +356,7 @@ def to_navigator(self, direction: Vector3d = Vector3d.zvector()): return s - def plot_over_signal(self, signal, annotate=False, **kwargs): + def plot_over_signal(self, signal, annotate=False, **plot_kwargs): """Convenience method to plot the orientation map and the n-best matches over the signal. Parameters @@ -331,13 +366,14 @@ def plot_over_signal(self, signal, annotate=False, **kwargs): annotate: bool If True, the euler rotation and the correlation will be annotated on the plot using the `Texts` class from hyperspy. - + Notes ----- The kwargs are passed to the `signal.plot` function call """ - signal.plot(**kwargs) - signal.add_marker(self.to_markers(1, annotate=annotate)) + nav = self.to_navigator() + signal.plot(navigator=nav, **plot_kwargs) + signal.add_marker(self.to_single_phase_markers(1, annotate=annotate)) def plot_inplane_rotation(self, **kwargs): """Plot the in-plane rotation of the orientation map as a 2D map.""" From 3690c0993ce43b660db22d548fd251ce7c201afb Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 7 May 2024 17:52:33 -0500 Subject: [PATCH 42/70] BugFix: Fix regression from rebase --- pyxem/signals/indexation_results.py | 24 +++++++++++++----------- pyxem/signals/polar_diffraction2d.py | 14 -------------- 2 files changed, 13 insertions(+), 25 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index a4937d181..501a4529f 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -170,10 +170,6 @@ class OrientationMap(DiffractionVectors2D): _signal_type = "orientation_map" - def __init__(self): - super().__init__() - self._signal_type = "orientation_map" - @property def simulation(self): return self.metadata.get_item("simulation") @@ -182,6 +178,12 @@ def simulation(self): def simulation(self, value): self.metadata.set_item("simulation", value) + def deepcopy(self): + """Deepcopy the signal""" + self.simulation._phase_slider = None + self.simulation._rotation_slider = None + return super().deepcopy() + def to_single_phase_orientations(self) -> Orientation: """Convert the orientation map to an `Orientation`-object, given a single-phase simulation. @@ -237,12 +239,10 @@ def extract_vectors_from_orientation_map(result, all_vectors): vectors.y = -vectors.y # Mirror if necessary vectors.y = mirror * vectors.y - rotation_matrix = ( - Rotation.from_euler( - (rotation, 0, 0), degrees=True, direction="crystal2lab" - ) - ).to_matrix() - vectors = vectors.rotate_from_matrix(rotation_matrix.squeeze()) + rotation = Rotation.from_euler( + (rotation, 0, 0), degrees=True, direction="crystal2lab" + ) + vectors = ~rotation * vectors return vectors @@ -300,7 +300,9 @@ def add_bar(i: int) -> str: coords = vectors.map( vectors_to_coordinates, inplace=False, lazy_output=lazy, ragged=True ) - markers = hs.plot.markers.Points.from_signal(coords, color=marker_color) + markers = hs.plot.markers.Points.from_signal( + coords, facecolor="none", edgecolor=marker_color, sizes=(20,) + ) yield markers if annotate: diff --git a/pyxem/signals/polar_diffraction2d.py b/pyxem/signals/polar_diffraction2d.py index 7d4fc1758..40571828b 100644 --- a/pyxem/signals/polar_diffraction2d.py +++ b/pyxem/signals/polar_diffraction2d.py @@ -373,20 +373,6 @@ def get_orientation( output_dtype=float, ) - # Translate in-plane rotation from index to degrees - # by using the calibration of the axis - def rotation_index_to_degrees(data): - data = data.copy() - axis = self.axes_manager.signal_axes[0].axis - ind = data[:, 2].astype(int) - data[:, 2] = rad2deg(axis[ind]) - # Account for the 90 degree discrepancy between the X_L axis - # and the origin of the azimuthal axis of the polar signal - data[:, 2] -= 90 - return data - - orientation.map(rotation_index_to_degrees) - orientation.set_signal_type("orientation_map") orientation.simulation = simulation orientation.column_names = ["index", "correlation", "rotation", "factor"] From 2edfc3af4b105e015d6cbfd87c5f7843dce21114 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 7 May 2024 17:53:06 -0500 Subject: [PATCH 43/70] Documentation: Add example --- examples/orientation_mapping/README.rst | 3 + .../on_zone_orientation.py | 43 +++++++++++ pyxem/data/__init__.py | 3 + pyxem/data/simulated_si.py | 77 +++++++++++++++++++ 4 files changed, 126 insertions(+) create mode 100644 examples/orientation_mapping/README.rst create mode 100644 examples/orientation_mapping/on_zone_orientation.py create mode 100644 pyxem/data/simulated_si.py 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/on_zone_orientation.py b/examples/orientation_mapping/on_zone_orientation.py new file mode 100644 index 000000000..725b2019f --- /dev/null +++ b/examples/orientation_mapping/on_zone_orientation.py @@ -0,0 +1,43 @@ +""" +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. +""" + +from pyxem.data import si_tilt, si_phase +from diffsims.generators.simulation_generator import SimulationGenerator +from orix.quaternion import Rotation + +simulated_si_tilt = si_tilt() + +# %% +# 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.calibrate.center = None +polar_si_tilt = simulated_si_tilt.get_azimuthal_integral2d( + npt=100, npt_azim=360, inplace=False, mean=True +) +polar_si_tilt.plot() + +# %% + +# Now we can get make the orientation map +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, +) +orientation_map = polar_si_tilt.get_orientation(sim) +orientation_map.plot_over_signal(simulated_si_tilt, annotate=True) diff --git a/pyxem/data/__init__.py b/pyxem/data/__init__.py index 342859fa6..5d5d84e5e 100644 --- a/pyxem/data/__init__.py +++ b/pyxem/data/__init__.py @@ -30,6 +30,7 @@ """ from pyxem.data.simulated_tilt import tilt_boundary_data +from pyxem.data.simulated_si import si_phase, si_tilt from pyxem.data._data import ( au_grating, pdnip_glass, @@ -52,6 +53,8 @@ "sample_with_g", "mgo_nanocrystals", "tilt_boundary_data", + "si_phase", + "si_tilt", "cuag_orientations", "organic_semiconductor", "feal_stripes", diff --git a/pyxem/data/simulated_si.py b/pyxem/data/simulated_si.py new file mode 100644 index 000000000..6f93034dd --- /dev/null +++ b/pyxem/data/simulated_si.py @@ -0,0 +1,77 @@ +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 scipy.ndimage import gaussian_filter1d +from pyxem.signals import Diffraction1D, Diffraction2D + + +def si_phase(): + 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 = sim.get_diffraction_pattern( + sigma=5, + ) + dp2 = sim.irot[1].get_diffraction_pattern(sigma=5) + + 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 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 From dc7772b8210002287d1079b7d1bfb63076dfd964 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 7 May 2024 19:27:33 -0500 Subject: [PATCH 44/70] Documentation: Clean up Example for on_zone_orientation.py --- examples/orientation_mapping/on_zone_orientation.py | 7 +++++-- pyxem/signals/indexation_results.py | 12 ++++++++++-- pyxem/utils/indexation_utils.py | 1 + 3 files changed, 16 insertions(+), 4 deletions(-) diff --git a/examples/orientation_mapping/on_zone_orientation.py b/examples/orientation_mapping/on_zone_orientation.py index 725b2019f..bf6661d6e 100644 --- a/examples/orientation_mapping/on_zone_orientation.py +++ b/examples/orientation_mapping/on_zone_orientation.py @@ -37,7 +37,10 @@ degrees=True, ), max_excitation_error=0.1, - reciprocal_radius=1, + reciprocal_radius=1.5, ) orientation_map = polar_si_tilt.get_orientation(sim) -orientation_map.plot_over_signal(simulated_si_tilt, annotate=True) +simulated_si_tilt.plot(vmax="99th") +simulated_si_tilt.add_marker( + orientation_map.to_single_phase_markers(annotate=True, text_color="w") +) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 501a4529f..3a9131b7a 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -17,7 +17,7 @@ # along with pyXem. If not, see . from warnings import warn -from typing import Iterator, Union, Literal +from typing import Union, Literal, Sequence, Iterator import hyperspy.api as hs from hyperspy._signals.lazy import LazySignal @@ -263,6 +263,7 @@ def to_single_phase_markers( marker_color: str = "red", text_color: str = "black", lazy: bool = True, + annotation_shift: Sequence[float] = None, ): """Convert the orientation map to a set of markers for plotting. @@ -277,7 +278,11 @@ def to_single_phase_markers( 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] """ + if annotation_shift is None: + annotation_shift = [0, -0.15] def vectors_to_coordinates(vectors): return np.vstack((vectors.x, vectors.y)).T @@ -309,8 +314,11 @@ def add_bar(i: int) -> str: texts = vectors.map( vectors_to_text, inplace=False, lazy_output=lazy, ragged=True ) + coords.map(lambda x: x + annotation_shift, inplace=True) text_markers = hs.plot.markers.Texts.from_signal( - coords, texts=texts.data.T, color=text_color + coords, + texts=texts.data.T, + color=text_color, ) yield text_markers diff --git a/pyxem/utils/indexation_utils.py b/pyxem/utils/indexation_utils.py index ef77a3ec5..ebbefcbbe 100644 --- a/pyxem/utils/indexation_utils.py +++ b/pyxem/utils/indexation_utils.py @@ -779,6 +779,7 @@ def _mixed_matching_lib_to_polar( """ 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 ( From 00d463318c6c01cfe65f5b928b22713ea094a303 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 7 May 2024 21:32:06 -0500 Subject: [PATCH 45/70] Documentation: Add examples for simulating multiple phases. --- .../mulit_phase_orientation.py | 57 ++++++++++++ .../single_phase_orientation.py | 47 ++++++++++ pyxem/data/__init__.py | 7 +- pyxem/data/simulated_si.py | 40 +++++++- pyxem/data/simulation_fe.py | 91 +++++++++++++++++++ pyxem/signals/indexation_results.py | 5 +- 6 files changed, 244 insertions(+), 3 deletions(-) create mode 100644 examples/orientation_mapping/mulit_phase_orientation.py create mode 100644 examples/orientation_mapping/single_phase_orientation.py create mode 100644 pyxem/data/simulation_fe.py diff --git a/examples/orientation_mapping/mulit_phase_orientation.py b/examples/orientation_mapping/mulit_phase_orientation.py new file mode 100644 index 000000000..a48235f43 --- /dev/null +++ b/examples/orientation_mapping/mulit_phase_orientation.py @@ -0,0 +1,57 @@ +""" +Mulit 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. +""" + +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.calibrate.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() + +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, +) +orientation_map = polar_multi.get_orientation(sim) + +# mulit_phase.add_marker(orientation_map.to_single_phase_markers(annotate=True, text_color="w")) + +# %% +# sphinx_gallery_thumbnail_number = 3 diff --git a/examples/orientation_mapping/single_phase_orientation.py b/examples/orientation_mapping/single_phase_orientation.py new file mode 100644 index 000000000..cd2503e07 --- /dev/null +++ b/examples/orientation_mapping/single_phase_orientation.py @@ -0,0 +1,47 @@ +""" +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. +""" + +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() + +# %% +# 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.calibrate.center = None +polar_si = simulated_si.get_azimuthal_integral2d( + npt=100, npt_azim=360, inplace=False, mean=True +) +polar_si.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` +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 +) +orientation_map = polar_si.get_orientation(sim) +navigator = orientation_map.to_navigator() + +simulated_si.plot(navigator=navigator, vmax="99th") +simulated_si.add_marker( + orientation_map.to_single_phase_markers(annotate=True, text_color="w") +) + +# %% +# sphinx_gallery_thumbnail_number = 3 diff --git a/pyxem/data/__init__.py b/pyxem/data/__init__.py index 5d5d84e5e..961693902 100644 --- a/pyxem/data/__init__.py +++ b/pyxem/data/__init__.py @@ -30,7 +30,8 @@ """ from pyxem.data.simulated_tilt import tilt_boundary_data -from pyxem.data.simulated_si import si_phase, si_tilt +from pyxem.data.simulated_si import si_phase, si_tilt, si_grains +from pyxem.data.simulation_fe import fe_bcc_phase, fe_fcc_phase, fe_multi_phase_grains from pyxem.data._data import ( au_grating, pdnip_glass, @@ -55,6 +56,10 @@ "tilt_boundary_data", "si_phase", "si_tilt", + "si_grains", + "fe_multi_phase_grains", + "fe_fcc_phase", + "fe_bcc_phase", "cuag_orientations", "organic_semiconductor", "feal_stripes", diff --git a/pyxem/data/simulated_si.py b/pyxem/data/simulated_si.py index 6f93034dd..4e4a29b3f 100644 --- a/pyxem/data/simulated_si.py +++ b/pyxem/data/simulated_si.py @@ -5,9 +5,14 @@ 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 = [] @@ -25,7 +30,9 @@ def si_phase(): def si_tilt(): p = si_phase() gen = SimulationGenerator() - rotations = Rotation.from_euler([[0, 0, 0], [10, 0, 0]], degrees=True) + rotations = Rotation.from_euler( + [[0, 0, 0], [10, 0, 0]], degrees=True, direction="crystal2lab" + ) sim = gen.calculate_diffraction2d( phase=p, rotation=rotations, reciprocal_radius=1.5, max_excitation_error=0.1 ) @@ -53,6 +60,37 @@ def si_tilt(): return tilt +def si_grains(num_grains=4, seed=2, size=20, recip_pixels=128): + p = si_phase() + gen = SimulationGenerator() + rotations = Rotation.random(num_grains) + sim = gen.calculate_diffraction2d( + phase=p, rotation=rotations, reciprocal_radius=1.5, max_excitation_error=0.1 + ) + dps = [ + 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 simulated1dsi( num_points=200, accelerating_voltage=200, diff --git a/pyxem/data/simulation_fe.py b/pyxem/data/simulation_fe.py new file mode 100644 index 000000000..670b34dfc --- /dev/null +++ b/pyxem/data/simulation_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/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 3a9131b7a..50adb08af 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -232,7 +232,10 @@ def to_single_phase_vectors(self, n_best_index: int = 0) -> hs.signals.Signal1D: def extract_vectors_from_orientation_map(result, all_vectors): index, _, rotation, mirror = result[n_best_index, :].T index = index.astype(int) - vectors = all_vectors[index] + 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 vectors = ReciprocalLatticeVector(vectors.phase, xyz=vectors.data.copy()) # Flip y, as discussed in https://github.com/pyxem/pyxem/issues/925 From d7483f67cec7f38f66cc0675592bc0a941601905 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Thu, 9 May 2024 13:53:09 -0500 Subject: [PATCH 46/70] Update: Update diffsims version --- examples/orientation_mapping/mulit_phase_orientation.py | 2 +- examples/orientation_mapping/on_zone_orientation.py | 2 +- examples/orientation_mapping/single_phase_orientation.py | 2 +- setup.py | 4 ++-- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/orientation_mapping/mulit_phase_orientation.py b/examples/orientation_mapping/mulit_phase_orientation.py index a48235f43..5daf6ee13 100644 --- a/examples/orientation_mapping/mulit_phase_orientation.py +++ b/examples/orientation_mapping/mulit_phase_orientation.py @@ -21,7 +21,7 @@ # 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.calibrate.center = None +mulit_phase.calibration.center = None polar_multi = mulit_phase.get_azimuthal_integral2d( npt=100, npt_azim=360, inplace=False, mean=True ) diff --git a/examples/orientation_mapping/on_zone_orientation.py b/examples/orientation_mapping/on_zone_orientation.py index bf6661d6e..84159d73a 100644 --- a/examples/orientation_mapping/on_zone_orientation.py +++ b/examples/orientation_mapping/on_zone_orientation.py @@ -19,7 +19,7 @@ # 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.calibrate.center = None +simulated_si_tilt.calibration.center = None polar_si_tilt = simulated_si_tilt.get_azimuthal_integral2d( npt=100, npt_azim=360, inplace=False, mean=True ) diff --git a/examples/orientation_mapping/single_phase_orientation.py b/examples/orientation_mapping/single_phase_orientation.py index cd2503e07..a89dece22 100644 --- a/examples/orientation_mapping/single_phase_orientation.py +++ b/examples/orientation_mapping/single_phase_orientation.py @@ -19,7 +19,7 @@ # 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.calibrate.center = None +simulated_si.calibration.center = None polar_si = simulated_si.get_azimuthal_integral2d( npt=100, npt_azim=360, inplace=False, mean=True ) diff --git a/setup.py b/setup.py index 2d3ab83ca..d74797e6c 100644 --- a/setup.py +++ b/setup.py @@ -87,7 +87,7 @@ extras_require=extra_feature_requirements, install_requires=[ "dask", - "diffsims >= 0.5", + "diffsims @ git+https://github.com/pyxem/diffsims.git@main", "hyperspy >= 2.0", "h5py", "lmfit >= 0.9.12", @@ -95,7 +95,7 @@ "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 From 691ed5bc315c6fbd17d412f91740bfe4e7cf2b57 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Thu, 9 May 2024 16:41:05 -0500 Subject: [PATCH 47/70] Update: Update tests for grain orientation. --- pyxem/data/simulated_si.py | 28 +++--- pyxem/signals/__init__.py | 3 +- pyxem/signals/indexation_results.py | 17 ++-- .../tests/signals/test_indexation_results.py | 86 ++++++++++++++++++- 4 files changed, 114 insertions(+), 20 deletions(-) diff --git a/pyxem/data/simulated_si.py b/pyxem/data/simulated_si.py index 4e4a29b3f..177eb4b4e 100644 --- a/pyxem/data/simulated_si.py +++ b/pyxem/data/simulated_si.py @@ -1,6 +1,6 @@ import numpy as np from orix.crystal_map import Phase -from orix.quaternion import Rotation +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 @@ -31,15 +31,16 @@ def si_tilt(): p = si_phase() gen = SimulationGenerator() rotations = Rotation.from_euler( - [[0, 0, 0], [10, 0, 0]], degrees=True, direction="crystal2lab" + [[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 = sim.get_diffraction_pattern( - sigma=5, - ) - dp2 = sim.irot[1].get_diffraction_pattern(sigma=5) + 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, @@ -60,15 +61,19 @@ def si_tilt(): return tilt -def si_grains(num_grains=4, seed=2, size=20, recip_pixels=128): +def si_grains(num_grains=4, seed=2, size=20, recip_pixels=128, return_rotations=False): p = si_phase() gen = SimulationGenerator() - rotations = Rotation.random(num_grains) + rotations = Orientation.random(num_grains, symmetry=p.point_group) sim = gen.calculate_diffraction2d( phase=p, rotation=rotations, reciprocal_radius=1.5, max_excitation_error=0.1 ) dps = [ - sim.irot[i].get_diffraction_pattern(sigma=5, shape=(recip_pixels, recip_pixels)) + 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) @@ -88,7 +93,10 @@ def si_grains(num_grains=4, seed=2, size=20, recip_pixels=128): 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 + if return_rotations: + return grains, rotations + else: + return grains def simulated1dsi( 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/indexation_results.py b/pyxem/signals/indexation_results.py index 50adb08af..c66e944c9 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -195,11 +195,18 @@ def to_single_phase_orientations(self) -> Orientation: # i.e. unique rotations for each navigation position rotations = hs.signals.Signal2D(self.simulation.rotations.data) - def rotation_from_orientation_map(result, rotations): + 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 = rotations[index] - euler = Orientation(ori).to_euler(degrees=True) * mirror[..., np.newaxis] + 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 @@ -207,7 +214,7 @@ def rotation_from_orientation_map(result, rotations): return Orientation( self.map( rotation_from_orientation_map, - rotations=rotations, + rots=rotations, inplace=False, ), symmetry=self.simulation.phases.point_group, @@ -245,7 +252,7 @@ def extract_vectors_from_orientation_map(result, all_vectors): rotation = Rotation.from_euler( (rotation, 0, 0), degrees=True, direction="crystal2lab" ) - vectors = ~rotation * vectors + vectors = ~rotation * vectors.to_miller() return vectors diff --git a/pyxem/tests/signals/test_indexation_results.py b/pyxem/tests/signals/test_indexation_results.py index 5af586e6f..de7621d46 100644 --- a/pyxem/tests/signals/test_indexation_results.py +++ b/pyxem/tests/signals/test_indexation_results.py @@ -25,12 +25,14 @@ 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 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 @pytest.fixture @@ -143,3 +145,79 @@ 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 + 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, + ) + orientation_map = polar_si_tilt.get_orientation(sim) + return orientation_map + + @pytest.fixture + 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 + ) + orientations = polar.get_orientation(sims) + return orientations, r + + 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, multi_rot_orientation_result): + orientations, rotations = multi_rot_orientation_result + assert isinstance(rotations, Orientation) + assert isinstance(orientations, OrientationMap) + orients = orientations.to_single_phase_orientations() + # Check that the orientations are within 1 degree of the expected value + + degrees_between = orients.angle_with(rotations, degrees=True) + assert np.all(np.min(degrees_between, axis=2) <= 1) + + def test_orientation_result(self, orientation_result): + assert isinstance(orientation_result[0], OrientationMap) + + def test_to_orientation(self, orientation_result): + orientation_result[0].to_single_phase_orientations() From 5a82d192cb5b81006a0f11cfc118835320623ea5 Mon Sep 17 00:00:00 2001 From: Viljar Femoen Date: Fri, 10 May 2024 18:19:58 +0200 Subject: [PATCH 48/70] Add low-index simulations --- pyxem/data/simulated_si.py | 35 +++++++++++++++++++++++++++++++---- 1 file changed, 31 insertions(+), 4 deletions(-) diff --git a/pyxem/data/simulated_si.py b/pyxem/data/simulated_si.py index 177eb4b4e..8f4be1e57 100644 --- a/pyxem/data/simulated_si.py +++ b/pyxem/data/simulated_si.py @@ -60,13 +60,13 @@ def si_tilt(): tilt.axes_manager.signal_axes[1].scale = 0.01 return tilt - -def si_grains(num_grains=4, seed=2, size=20, recip_pixels=128, return_rotations=False): +def si_grains_from_orientations(oris: Orientation, seed: int = 2, size: int = 20, recip_pixels: int = 128): p = si_phase() gen = SimulationGenerator() - rotations = Orientation.random(num_grains, symmetry=p.point_group) + num_grains = oris.size + sim = gen.calculate_diffraction2d( - phase=p, rotation=rotations, reciprocal_radius=1.5, max_excitation_error=0.1 + phase=p, rotation=oris, reciprocal_radius=1.5, max_excitation_error=0.1 ) dps = [ np.flipud( @@ -93,11 +93,38 @@ def si_grains(num_grains=4, seed=2, size=20, recip_pixels=128, return_rotations= 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_random(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 simulated1dsi( num_points=200, From cb4aba4b7d43398c03d4554403b3030cce15d8ca Mon Sep 17 00:00:00 2001 From: Viljar Femoen Date: Fri, 10 May 2024 18:20:15 +0200 Subject: [PATCH 49/70] Translate angle index to degrees --- pyxem/signals/polar_diffraction2d.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/pyxem/signals/polar_diffraction2d.py b/pyxem/signals/polar_diffraction2d.py index 40571828b..cb19f1c96 100644 --- a/pyxem/signals/polar_diffraction2d.py +++ b/pyxem/signals/polar_diffraction2d.py @@ -373,6 +373,18 @@ def get_orientation( 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.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"] From 438c2ed91c56a4a98348af9771469e0282b2939a Mon Sep 17 00:00:00 2001 From: Viljar Femoen Date: Fri, 10 May 2024 18:25:42 +0200 Subject: [PATCH 50/70] Update test to use `si_phase_simple` --- pyxem/data/__init__.py | 3 +- pyxem/data/simulated_si.py | 28 ++++++++++++------ .../tests/signals/test_indexation_results.py | 29 +++++++++++++++---- 3 files changed, 45 insertions(+), 15 deletions(-) diff --git a/pyxem/data/__init__.py b/pyxem/data/__init__.py index 961693902..5f7524ccb 100644 --- a/pyxem/data/__init__.py +++ b/pyxem/data/__init__.py @@ -30,7 +30,7 @@ """ from pyxem.data.simulated_tilt import tilt_boundary_data -from pyxem.data.simulated_si import si_phase, si_tilt, si_grains +from pyxem.data.simulated_si import si_phase, si_tilt, si_grains, si_grains_simple from pyxem.data.simulation_fe import fe_bcc_phase, fe_fcc_phase, fe_multi_phase_grains from pyxem.data._data import ( au_grating, @@ -57,6 +57,7 @@ "si_phase", "si_tilt", "si_grains", + "si_grains_simple", "fe_multi_phase_grains", "fe_fcc_phase", "fe_bcc_phase", diff --git a/pyxem/data/simulated_si.py b/pyxem/data/simulated_si.py index 8f4be1e57..b4f09e1e9 100644 --- a/pyxem/data/simulated_si.py +++ b/pyxem/data/simulated_si.py @@ -60,7 +60,10 @@ def si_tilt(): 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): + +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 @@ -95,37 +98,44 @@ def si_grains_from_orientations(oris: Orientation, seed: int = 2, size: int = 20 grains.axes_manager.signal_axes[1].scale = 0.01 return grains -def si_grains_random(num_grains=4, seed=2, size=20, recip_pixels=128, return_rotations=False): + +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) - + 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, 0, 0], # [0 0 1] + [0, 45, 0], # [0 1 1] [0, 54.7, 45], # [1 1 1] - [0, 35, 45], # [1 1 2] + [0, 35, 45], # [1 1 2] ], degrees=True, - symmetry=p.point_group + symmetry=p.point_group, + ) + grains = si_grains_from_orientations( + rotations, seed=seed, size=size, recip_pixels=recip_pixels ) - grains = si_grains_from_orientations(rotations, seed=seed, size=size, recip_pixels=recip_pixels) if return_rotations: return grains, rotations else: return grains + def simulated1dsi( num_points=200, accelerating_voltage=200, diff --git a/pyxem/tests/signals/test_indexation_results.py b/pyxem/tests/signals/test_indexation_results.py index de7621d46..3fa7a1079 100644 --- a/pyxem/tests/signals/test_indexation_results.py +++ b/pyxem/tests/signals/test_indexation_results.py @@ -32,7 +32,7 @@ from pyxem.generators import TemplateIndexationGenerator from pyxem.signals import VectorMatchingResults, DiffractionVectors, OrientationMap from pyxem.utils.indexation_utils import OrientationResult -from pyxem.data import si_grains, si_phase, si_tilt +from pyxem.data import si_grains, si_phase, si_tilt, si_grains_simple @pytest.fixture @@ -191,6 +191,24 @@ def multi_rot_orientation_result(self): orientations = polar.get_orientation(sims) return orientations, r + @pytest.fixture + 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 + ) + orientations = polar.get_orientation(sims) + return orientations, r + 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() @@ -206,15 +224,16 @@ def test_tilt_orientation_result(self, single_rot_orientation_result): ) assert np.all(degrees_between[:, 5:] <= 1) - def test_grain_orientation_result(self, multi_rot_orientation_result): - orientations, rotations = multi_rot_orientation_result + def test_grain_orientation_result(self, simple_multi_rot_orientation_result): + orientations, rotations = 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 1 degree of the expected value + # 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) <= 1) + assert np.all(np.min(degrees_between, axis=2) <= 2) def test_orientation_result(self, orientation_result): assert isinstance(orientation_result[0], OrientationMap) From 9f8e68668a17ee3eff7cea68906cda03f7887c61 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Mon, 13 May 2024 18:24:10 -0500 Subject: [PATCH 51/70] Testing: Remove unused tests --- pyxem/tests/signals/test_indexation_results.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/pyxem/tests/signals/test_indexation_results.py b/pyxem/tests/signals/test_indexation_results.py index 3fa7a1079..8d955b4be 100644 --- a/pyxem/tests/signals/test_indexation_results.py +++ b/pyxem/tests/signals/test_indexation_results.py @@ -234,9 +234,3 @@ def test_grain_orientation_result(self, simple_multi_rot_orientation_result): # 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_orientation_result(self, orientation_result): - assert isinstance(orientation_result[0], OrientationMap) - - def test_to_orientation(self, orientation_result): - orientation_result[0].to_single_phase_orientations() From 0c735a024a95ad4a42b185262239248115b03f60 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 14 May 2024 13:07:07 -0500 Subject: [PATCH 52/70] Plotting: Add plotting for diffracting intensity --- pyxem/signals/indexation_results.py | 151 +++++++++++++++++++--------- 1 file changed, 101 insertions(+), 50 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index c66e944c9..369b65d49 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -28,7 +28,7 @@ from orix.vector import Vector3d from orix.plot import IPFColorKeyTSL from transforms3d.euler import mat2euler -from diffsims.crystallography import ReciprocalLatticeVector +from diffsims.crystallography._diffracting_vector import DiffractingVector from pyxem.utils.indexation_utils import get_nth_best_solution from pyxem.signals.diffraction_vectors2d import DiffractionVectors2D @@ -184,7 +184,7 @@ def deepcopy(self): self.simulation._rotation_slider = None return super().deepcopy() - def to_single_phase_orientations(self) -> Orientation: + def to_single_phase_orientations(self, **kwargs) -> Orientation: """Convert the orientation map to an `Orientation`-object, given a single-phase simulation. """ @@ -216,6 +216,7 @@ def rotation_from_orientation_map(result, rots): rotation_from_orientation_map, rots=rotations, inplace=False, + **kwargs, ), symmetry=self.simulation.phases.point_group, ) @@ -244,7 +245,8 @@ def extract_vectors_from_orientation_map(result, all_vectors): else: vectors = all_vectors[index] # Copy manually, as deepcopy adds a lot of overhead with the phase - vectors = ReciprocalLatticeVector(vectors.phase, xyz=vectors.data.copy()) + 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 # Mirror if necessary @@ -253,6 +255,9 @@ def extract_vectors_from_orientation_map(result, all_vectors): (rotation, 0, 0), degrees=True, direction="crystal2lab" ) vectors = ~rotation * vectors.to_miller() + vectors = DiffractingVector( + vectors.phase, xyz=vectors.data.copy(), intensity=intensity + ) return vectors @@ -266,14 +271,25 @@ def to_crystal_map(self) -> CrystalMap: """Convert the orientation map to an `orix.CrystalMap` object""" pass + def to_ipf_markers(self): + """Convert the orientation map to a set of inverse pole figure + markers which visualizes the best matching orientations in the + reduced S2 space.""" + + pass + def to_single_phase_markers( self, n_best: int = 1, annotate: bool = False, - marker_color: str = "red", + marker_colors: str = ("red", "blue", "green", "orange", "purple"), text_color: str = "black", - lazy: bool = True, + lazy_output: bool = None, annotation_shift: Sequence[float] = None, + text_kwargs: dict = None, + include_intensity: bool = False, + intesity_scale: float = 1, + **kwargs, ): """Convert the orientation map to a set of markers for plotting. @@ -290,45 +306,46 @@ def to_single_phase_markers( 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. """ + if text_kwargs is None: + text_kwargs = dict() if annotation_shift is None: annotation_shift = [0, -0.15] - def vectors_to_coordinates(vectors): - return np.vstack((vectors.x, vectors.y)).T - - def vectors_to_text(vectors): - 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 - for n in range(n_best): vectors = self.to_single_phase_vectors(n) + color = marker_colors[n % len(marker_colors)] + if include_intensity: + intensity = vectors.map( + vectors_to_intensity, + scale=intesity_scale, + inplace=False, + lazy_output=lazy_output, + ragged=True, + ).data + kwargs["sizes"] = intensity + coords = vectors.map( - vectors_to_coordinates, inplace=False, lazy_output=lazy, ragged=True + vectors_to_coordinates, + inplace=False, + lazy_output=lazy_output, + ragged=True, ) markers = hs.plot.markers.Points.from_signal( - coords, facecolor="none", edgecolor=marker_color, sizes=(20,) + coords, facecolor="none", edgecolor=color, **kwargs ) yield markers if annotate: texts = vectors.map( - vectors_to_text, inplace=False, lazy_output=lazy, ragged=True + vectors_to_text, inplace=False, lazy_output=lazy_output, ragged=True ) coords.map(lambda x: x + annotation_shift, inplace=True) text_markers = hs.plot.markers.Texts.from_signal( - coords, - texts=texts.data.T, - color=text_color, + coords, texts=texts.data.T, color=text_color, **text_kwargs ) yield text_markers @@ -503,32 +520,66 @@ def get_crystallographic_map(self, *args, **kwargs): crystal_map = _transfer_navigation_axes(crystal_map, self) return crystal_map - def get_indexed_diffraction_vectors( - self, vectors, overwrite=False, *args, **kwargs - ): - """Obtain an indexed diffraction vectors object. - Parameters - ---------- - vectors : pyxem.signals.DiffractionVectors - A diffraction vectors object to be indexed. +def get_indexed_diffraction_vectors(self, vectors, overwrite=False, *args, **kwargs): + """Obtain an indexed diffraction vectors object. - Returns - ------- - indexed_vectors : pyxem.signals.DiffractionVectors - An indexed diffraction vectors object. + Parameters + ---------- + vectors : pyxem.signals.DiffractionVectors + A diffraction vectors object to be indexed. - """ - if overwrite is False: - if vectors.hkls is not None: - warn( - "The vectors supplied are already associated with hkls set " - "overwrite=True to replace these hkls." - ) - else: - vectors.hkls = self.hkls + Returns + ------- + indexed_vectors : pyxem.signals.DiffractionVectors + An indexed diffraction vectors object. - elif overwrite is True: + """ + if overwrite is False: + if vectors.hkls is not None: + warn( + "The vectors supplied are already associated with hkls set " + "overwrite=True to replace these hkls." + ) + else: vectors.hkls = self.hkls - return vectors + elif overwrite is True: + vectors.hkls = self.hkls + return vectors + + +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 From 82c0e0113968d645b45132545402028324786d0f Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 14 May 2024 15:12:30 -0500 Subject: [PATCH 53/70] Plotting: Add plotting of ipf --- pyxem/data/__init__.py | 9 +- pyxem/data/simulated_si.py | 23 +++++ pyxem/signals/indexation_results.py | 132 ++++++++++++++++++++++------ 3 files changed, 138 insertions(+), 26 deletions(-) diff --git a/pyxem/data/__init__.py b/pyxem/data/__init__.py index 5f7524ccb..0a73806c9 100644 --- a/pyxem/data/__init__.py +++ b/pyxem/data/__init__.py @@ -30,7 +30,13 @@ """ from pyxem.data.simulated_tilt import tilt_boundary_data -from pyxem.data.simulated_si import si_phase, si_tilt, si_grains, si_grains_simple +from pyxem.data.simulated_si import ( + si_phase, + si_tilt, + si_grains, + si_grains_simple, + si_rotations_line, +) from pyxem.data.simulation_fe import fe_bcc_phase, fe_fcc_phase, fe_multi_phase_grains from pyxem.data._data import ( au_grating, @@ -58,6 +64,7 @@ "si_tilt", "si_grains", "si_grains_simple", + "si_rotations_line", "fe_multi_phase_grains", "fe_fcc_phase", "fe_bcc_phase", diff --git a/pyxem/data/simulated_si.py b/pyxem/data/simulated_si.py index b4f09e1e9..7f6697af0 100644 --- a/pyxem/data/simulated_si.py +++ b/pyxem/data/simulated_si.py @@ -136,6 +136,29 @@ def si_grains_simple(seed=2, size=20, recip_pixels=128, return_rotations=False): 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, diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 369b65d49..3e5e3fc73 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -29,6 +29,12 @@ 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 +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 @@ -274,9 +280,84 @@ def to_crystal_map(self) -> CrystalMap: def to_ipf_markers(self): """Convert the orientation map to a set of inverse pole figure markers which visualizes the best matching orientations in the - reduced S2 space.""" + reduced S2 space. + """ + 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") - pass + orients = self.to_single_phase_orientations() + sector = self.simulation.phases.point_group.fundamental_sector + labels = _get_ipf_axes_labels( + sector.vertices, symmetry=self.simulation.phases.point_group + ) + s = StereographicProjection() + vectors = orients * Vector3d.zvector() + edges = _closed_edges_in_hemisphere(sector.edges, sector) + vectors = vectors.in_fundamental_sector(self.simulation.phases.point_group) + x, y = s.vector2xy(vectors) + ex, ey = s.vector2xy(edges) + tx, ty = s.vector2xy(sector.vertices) + + 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) + original_offset = np.vstack((ex, ey)).T + texts_offset = np.vstack((tx, ty)).T + + mins, maxes = original_offset.min(axis=0), original_offset.max(axis=0) + + original_offset = ( + (original_offset - ((maxes + mins) / 2)) / (maxes - mins) * 0.2 + ) + original_offset = original_offset + 0.85 + + texts_offset = (texts_offset - ((maxes + mins) / 2)) / (maxes - mins) * 0.2 + texts_offset = texts_offset + 0.85 + for i in np.ndindex(offsets.shape): + off = np.vstack((x[i], y[i])).T + norm_points = (off - ((maxes + mins) / 2)) / (maxes - mins) * 0.2 + norm_points = norm_points + 0.85 + offsets[i] = norm_points + correlation[i] = cor[i] / np.max(cor[i]) * 0.5 + + square = hs.plot.markers.Squares( + offsets=[[0.85, 0.85]], + widths=(0.3,), + units="width", + offset_transform="axes", + facecolor="white", + edgecolor="black", + ) + polygon_sector = hs.plot.markers.Polygons( + verts=original_offset[np.newaxis], + transform="axes", + alpha=1, + facecolor="none", + ) + + best_points = hs.plot.markers.Points( + offsets=offsets.T, + sizes=(4,), + offset_transform="axes", + alpha=correlation.T, + facecolor="green", + ) + + texts = hs.plot.markers.Texts( + texts=labels, + offsets=texts_offset, + sizes=(1,), + offset_transform="axes", + facecolor="k", + ) + return square, polygon_sector, best_points, texts def to_single_phase_markers( self, @@ -325,7 +406,7 @@ def to_single_phase_markers( inplace=False, lazy_output=lazy_output, ragged=True, - ).data + ).data.T kwargs["sizes"] = intensity coords = vectors.map( @@ -520,33 +601,34 @@ def get_crystallographic_map(self, *args, **kwargs): crystal_map = _transfer_navigation_axes(crystal_map, self) return crystal_map + def get_indexed_diffraction_vectors( + self, vectors, overwrite=False, *args, **kwargs + ): + """Obtain an indexed diffraction vectors object. -def get_indexed_diffraction_vectors(self, vectors, overwrite=False, *args, **kwargs): - """Obtain an indexed diffraction vectors object. + Parameters + ---------- + vectors : pyxem.signals.DiffractionVectors + A diffraction vectors object to be indexed. - Parameters - ---------- - vectors : pyxem.signals.DiffractionVectors - A diffraction vectors object to be indexed. + Returns + ------- + indexed_vectors : pyxem.signals.DiffractionVectors + An indexed diffraction vectors object. - Returns - ------- - indexed_vectors : pyxem.signals.DiffractionVectors - An indexed diffraction vectors object. + """ + if overwrite is False: + if vectors.hkls is not None: + warn( + "The vectors supplied are already associated with hkls set " + "overwrite=True to replace these hkls." + ) + else: + vectors.hkls = self.hkls - """ - if overwrite is False: - if vectors.hkls is not None: - warn( - "The vectors supplied are already associated with hkls set " - "overwrite=True to replace these hkls." - ) - else: + elif overwrite is True: vectors.hkls = self.hkls - - elif overwrite is True: - vectors.hkls = self.hkls - return vectors + return vectors def vectors_to_coordinates(vectors): From 73491332bbca9f9187bb910952eb4d3cc09010b0 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 24 May 2024 11:28:01 -0500 Subject: [PATCH 54/70] Refactor:Add kwargs --- pyxem/signals/indexation_results.py | 30 ++++++++++++++++++++++++----- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 3e5e3fc73..7e00f3fce 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -227,7 +227,9 @@ def rotation_from_orientation_map(result, rots): symmetry=self.simulation.phases.point_group, ) - def to_single_phase_vectors(self, n_best_index: int = 0) -> hs.signals.Signal1D: + 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. @@ -271,6 +273,9 @@ def extract_vectors_from_orientation_map(result, all_vectors): extract_vectors_from_orientation_map, all_vectors=vectors_signal, inplace=False, + output_signal_size=(), + output_dtype=object, + **kwargs, ) def to_crystal_map(self) -> CrystalMap: @@ -395,25 +400,35 @@ def to_single_phase_markers( 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 for n in range(n_best): - vectors = self.to_single_phase_vectors(n) + vectors = self.to_single_phase_vectors( + lazy_output=lazy_output, 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, - lazy_output=lazy_output, ragged=True, + output_dtype=object, + output_signal_size=(), + navigation_chunks=navigation_chunks, ).data.T kwargs["sizes"] = intensity coords = vectors.map( vectors_to_coordinates, inplace=False, - lazy_output=lazy_output, ragged=True, + output_dtype=object, + output_signal_size=(), + navigation_chunks=navigation_chunks, ) markers = hs.plot.markers.Points.from_signal( coords, facecolor="none", edgecolor=color, **kwargs @@ -422,7 +437,12 @@ def to_single_phase_markers( if annotate: texts = vectors.map( - vectors_to_text, inplace=False, lazy_output=lazy_output, ragged=True + vectors_to_text, + inplace=False, + lazy_output=lazy_output, + ragged=True, + output_dtype=object, + output_signal_size=(), ) coords.map(lambda x: x + annotation_shift, inplace=True) text_markers = hs.plot.markers.Texts.from_signal( From fabe79f784939f5dcf2e110e25525c7bcf592351 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Sun, 26 May 2024 21:38:06 -0500 Subject: [PATCH 55/70] New Feature: Add to `CrystalMap` --- pyxem/signals/indexation_results.py | 163 +++++++++++++++++++++------- 1 file changed, 126 insertions(+), 37 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 7e00f3fce..9e60bb65c 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -23,7 +23,7 @@ from hyperspy._signals.lazy import LazySignal from hyperspy.signal import BaseSignal import numpy as np -from orix.crystal_map import CrystalMap, Phase +from orix.crystal_map import CrystalMap, Phase, PhaseList from orix.quaternion import Rotation, Orientation from orix.vector import Vector3d from orix.plot import IPFColorKeyTSL @@ -159,6 +159,64 @@ 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 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 @@ -190,6 +248,38 @@ def deepcopy(self): self.simulation._rotation_slider = None return super().deepcopy() + def to_rotation(self, flatten=False): + """ + Convert the orientation map to a set of `orix.Quaternion.Rotation` objects. + Returns + ------- + """ + if self._lazy: + warn("Computing the s ") + all_rotations = Rotation.stack(self.simulation.rotations).flatten() + rotations = self.map( + rotation_from_orientation_map, + rots=all_rotations, + inplace=False, + lazy_output=False, + ) + + 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 + """ + 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. @@ -280,6 +370,41 @@ def extract_vectors_from_orientation_map(result, all_vectors): def to_crystal_map(self) -> CrystalMap: """Convert the orientation map to an `orix.CrystalMap` object""" + 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)) + else: + phases = PhaseList(self.simulation.phases) + + return CrystalMap( + rotations=rotations, + x=xx.flatten(), + phase_id=phase_index, + y=yy.flatten(), + scan_unit=scan_unit, + phase_list=phases, + ) + + def to_polar_mask(self, reciporical_radius: float = 0.1, **kwargs): + """Convert the orientation map to a mask in polar coordinate for + finding multiple overlapping phases. This is useful for indexing + multiple phases by finding the best matching orientation and then + masking the data to find the next best matching orientation for + the second phase. + """ pass def to_ipf_markers(self): @@ -649,39 +774,3 @@ def get_indexed_diffraction_vectors( elif overwrite is True: vectors.hkls = self.hkls return vectors - - -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 From de23d145cd9d2344b9ed6bb1f21e76e3f60f499d Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Sun, 26 May 2024 22:10:18 -0500 Subject: [PATCH 56/70] New Feature: Add to `CrystalMap` --- examples/orientation_mapping/mulit_phase_orientation.py | 3 ++- pyxem/signals/indexation_results.py | 3 +++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/examples/orientation_mapping/mulit_phase_orientation.py b/examples/orientation_mapping/mulit_phase_orientation.py index 5daf6ee13..d25bf6d3f 100644 --- a/examples/orientation_mapping/mulit_phase_orientation.py +++ b/examples/orientation_mapping/mulit_phase_orientation.py @@ -51,7 +51,8 @@ ) orientation_map = polar_multi.get_orientation(sim) -# mulit_phase.add_marker(orientation_map.to_single_phase_markers(annotate=True, text_color="w")) +cmap = orientation_map.to_crystal_map() +cmap.plot() # %% # sphinx_gallery_thumbnail_number = 3 diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 9e60bb65c..7b7f56b63 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -386,6 +386,9 @@ def to_crystal_map(self) -> CrystalMap: 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) From 3a245d18f7bbbc570690981df6f19e7d7375be63 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Mon, 27 May 2024 08:41:47 -0500 Subject: [PATCH 57/70] Refactor: Refactor indexation_results.py --- pyxem/signals/indexation_results.py | 85 ++++++++++++----------------- 1 file changed, 35 insertions(+), 50 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 7b7f56b63..7f62c4f84 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -212,6 +212,31 @@ def rotation_from_orientation_map(result, rots): return ori +def extract_vectors_from_orientation_map(result, all_vectors): + 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 + # Mirror if necessary + vectors.y = mirror * vectors.y + rotation = Rotation.from_euler( + (rotation, 0, 0), degrees=True, direction="crystal2lab" + ) + vectors = ~rotation * vectors.to_miller() + vectors = DiffractingVector( + vectors.phase, xyz=vectors.data.copy(), intensity=intensity + ) + + return vectors + + def orientation2phase(orientations, sizes): o = orientations[:, 0] return np.searchsorted(sizes, o) @@ -273,6 +298,12 @@ def to_rotation(self, flatten=False): 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]) @@ -285,28 +316,14 @@ def to_single_phase_orientations(self, **kwargs) -> Orientation: given a single-phase simulation. """ if self.simulation.has_multiple_phases: - raise ValueError("Multiple phases found in simulation") + 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) - 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 - return Orientation( self.map( rotation_from_orientation_map, @@ -335,36 +352,13 @@ def to_single_phase_vectors( # Use vector data as signal in case of different vectors per navigation position vectors_signal = hs.signals.Signal1D(self.simulation.coordinates) - def extract_vectors_from_orientation_map(result, all_vectors): - 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 - # Mirror if necessary - vectors.y = mirror * vectors.y - rotation = Rotation.from_euler( - (rotation, 0, 0), degrees=True, direction="crystal2lab" - ) - vectors = ~rotation * vectors.to_miller() - vectors = DiffractingVector( - vectors.phase, xyz=vectors.data.copy(), intensity=intensity - ) - - return vectors - return self.map( extract_vectors_from_orientation_map, all_vectors=vectors_signal, inplace=False, output_signal_size=(), output_dtype=object, + show_progressbar=False, **kwargs, ) @@ -401,15 +395,6 @@ def to_crystal_map(self) -> CrystalMap: phase_list=phases, ) - def to_polar_mask(self, reciporical_radius: float = 0.1, **kwargs): - """Convert the orientation map to a mask in polar coordinate for - finding multiple overlapping phases. This is useful for indexing - multiple phases by finding the best matching orientation and then - masking the data to find the next best matching orientation for - the second phase. - """ - pass - def to_ipf_markers(self): """Convert the orientation map to a set of inverse pole figure markers which visualizes the best matching orientations in the From 994eba806256da09e317f7be2547f06019ad66f3 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Mon, 27 May 2024 09:23:06 -0500 Subject: [PATCH 58/70] Testing: Test `to_crystal_map` function --- pyxem/signals/indexation_results.py | 17 +++++++++++++---- pyxem/tests/signals/test_indexation_results.py | 8 ++++++++ 2 files changed, 21 insertions(+), 4 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 7f62c4f84..4f3bdf674 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -196,8 +196,8 @@ def add_bar(i: int) -> str: def rotation_from_orientation_map(result, rots): - # if rots.ndim == 1: - # rots = rots[np.newaxis, :] + if rots.ndim == 1: + rots = rots[np.newaxis, :] index, _, rotation, mirror = result.T index = index.astype(int) ori = rots[index] @@ -273,6 +273,10 @@ def deepcopy(self): self.simulation._rotation_slider = None return super().deepcopy() + @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. @@ -280,13 +284,16 @@ def to_rotation(self, flatten=False): ------- """ if self._lazy: - warn("Computing the s ") + warn("Computing the signal") + self.compute() all_rotations = Rotation.stack(self.simulation.rotations).flatten() rotations = self.map( rotation_from_orientation_map, - rots=all_rotations, + rots=all_rotations.data, inplace=False, lazy_output=False, + output_signal_size=(self.num_rows, 4), + output_dtype=float, ) rots = Rotation(rotations) @@ -329,6 +336,8 @@ def to_single_phase_orientations(self, **kwargs) -> Orientation: 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, diff --git a/pyxem/tests/signals/test_indexation_results.py b/pyxem/tests/signals/test_indexation_results.py index 8d955b4be..2c5014ab9 100644 --- a/pyxem/tests/signals/test_indexation_results.py +++ b/pyxem/tests/signals/test_indexation_results.py @@ -28,11 +28,13 @@ 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, OrientationMap from pyxem.utils.indexation_utils import OrientationResult from pyxem.data import si_grains, si_phase, si_tilt, si_grains_simple +import hyperspy.api as hs @pytest.fixture @@ -234,3 +236,9 @@ def test_grain_orientation_result(self, simple_multi_rot_orientation_result): # 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 = simple_multi_rot_orientation_result + crystal_map = orientations.to_crystal_map() + assert isinstance(crystal_map, CrystalMap) + assert np.all(crystal_map.phase_id == 0) From 9e6afda90b2d53520ff649c706ea96a7144fa291 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Mon, 27 May 2024 16:46:32 -0500 Subject: [PATCH 59/70] Testing: Test `to_crystal_map` function for multiple phases --- pyxem/signals/polar_diffraction2d.py | 3 ++ .../tests/signals/test_indexation_results.py | 41 ++++++++++++++++++- 2 files changed, 43 insertions(+), 1 deletion(-) diff --git a/pyxem/signals/polar_diffraction2d.py b/pyxem/signals/polar_diffraction2d.py index cb19f1c96..476cd7f4e 100644 --- a/pyxem/signals/polar_diffraction2d.py +++ b/pyxem/signals/polar_diffraction2d.py @@ -381,6 +381,9 @@ def rotation_index_to_degrees(data, axis): 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 ) diff --git a/pyxem/tests/signals/test_indexation_results.py b/pyxem/tests/signals/test_indexation_results.py index 2c5014ab9..6688ea2db 100644 --- a/pyxem/tests/signals/test_indexation_results.py +++ b/pyxem/tests/signals/test_indexation_results.py @@ -33,7 +33,15 @@ from pyxem.generators import TemplateIndexationGenerator 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 +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 @@ -211,6 +219,32 @@ def simple_multi_rot_orientation_result(self): orientations = polar.get_orientation(sims) return orientations, r + @pytest.fixture + 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) + 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() @@ -242,3 +276,8 @@ def test_to_crystal_map(self, 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) From 2d03b93b7f5dcad0cbb830ae03d3f23aaca160b4 Mon Sep 17 00:00:00 2001 From: Viljar Femoen Date: Mon, 27 May 2024 14:35:51 +0200 Subject: [PATCH 60/70] Try to ensure correct conventions is followed --- pyxem/signals/indexation_results.py | 83 +++++++++++++++++++++++------ 1 file changed, 67 insertions(+), 16 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 4f3bdf674..f6002609f 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -22,6 +22,7 @@ import hyperspy.api as hs from hyperspy._signals.lazy import LazySignal from hyperspy.signal import BaseSignal +from hyperspy.axes import AxesManager import numpy as np from orix.crystal_map import CrystalMap, Phase, PhaseList from orix.quaternion import Rotation, Orientation @@ -196,8 +197,8 @@ def add_bar(i: int) -> str: def rotation_from_orientation_map(result, rots): - if rots.ndim == 1: - rots = rots[np.newaxis, :] + # if rots.ndim == 1: + # rots = rots[np.newaxis, :] index, _, rotation, mirror = result.T index = index.astype(int) ori = rots[index] @@ -361,6 +362,34 @@ def to_single_phase_vectors( # Use vector data as signal in case of different vectors per navigation position vectors_signal = hs.signals.Signal1D(self.simulation.coordinates) + def extract_vectors_from_orientation_map(result, all_vectors): + 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 + return self.map( extract_vectors_from_orientation_map, all_vectors=vectors_signal, @@ -566,39 +595,61 @@ def to_single_phase_markers( output_dtype=object, output_signal_size=(), ) - coords.map(lambda x: x + annotation_shift, inplace=True) + # 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) text_markers = hs.plot.markers.Texts.from_signal( - coords, texts=texts.data.T, color=text_color, **text_kwargs + text_coords, texts=texts.data.T, color=text_color, **text_kwargs ) yield text_markers - def to_polar_markers(self, n_best: int = 1) -> Iterator[hs.plot.markers.Markers]: + def to_single_phase_polar_markers( + self, signal_axes: AxesManager, n_best: int = 1, marker_color: str = "red" + ) -> Iterator[hs.plot.markers.Markers]: ( r_templates, theta_templates, intensities_templates, - ) = self.simulation.polar_flatten_simulations() + ) = self.simulation.polar_flatten_simulations( + signal_axes[1].axis, signal_axes[0].axis + ) + + 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_factory(n_best_entry: int): def marker_generator(entry): - index, correlation, rotation, factor = entry[n_best_entry] - r = r_templates[int(index)] - theta = theta_templates[int(index)] - theta += 2 * np.pi + np.deg2rad(rotation) - theta %= 2 * np.pi - theta -= np.pi - return np.array((theta, r)).T + 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 for n in range(n_best): markers_signal = self.map( - marker_generator_factory(n), + marker_generator_factory(n, signal_axes[1].axis, signal_axes[0].axis), inplace=False, ragged=True, lazy_output=True, ) - markers = hs.plot.markers.Points.from_signal(markers_signal) + markers = hs.plot.markers.Points.from_signal( + markers_signal, facecolor="none", edgecolor=marker_color, sizes=(10,) + ) yield markers def to_navigator(self, direction: Vector3d = Vector3d.zvector()): From 20bbfe5772d7e02c78002d637ca39c551631a725 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 11:52:16 -0500 Subject: [PATCH 61/70] Rebase: Follow up for Rebase --- pyxem/signals/indexation_results.py | 38 ++++++----------------------- 1 file changed, 7 insertions(+), 31 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index f6002609f..b4cb4e07a 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -225,16 +225,20 @@ def extract_vectors_from_orientation_map(result, all_vectors): vectors = DiffractingVector(vectors.phase, xyz=vectors.data.copy()) # Flip y, as discussed in https://github.com/pyxem/pyxem/issues/925 vectors.y = -vectors.y - # Mirror if necessary - vectors.y = mirror * vectors.y + rotation = Rotation.from_euler( - (rotation, 0, 0), degrees=True, direction="crystal2lab" + (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 @@ -362,34 +366,6 @@ def to_single_phase_vectors( # Use vector data as signal in case of different vectors per navigation position vectors_signal = hs.signals.Signal1D(self.simulation.coordinates) - def extract_vectors_from_orientation_map(result, all_vectors): - 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 - return self.map( extract_vectors_from_orientation_map, all_vectors=vectors_signal, From eb86caa087961f50bd4883b2dae1905646cb83c0 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 12:49:12 -0500 Subject: [PATCH 62/70] Refactor: Fix Tests --- pyxem/signals/indexation_results.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index b4cb4e07a..4d87495e7 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -197,8 +197,8 @@ def add_bar(i: int) -> str: def rotation_from_orientation_map(result, rots): - # if rots.ndim == 1: - # rots = rots[np.newaxis, :] + if rots.ndim == 1: + rots = rots[np.newaxis, :] index, _, rotation, mirror = result.T index = index.astype(int) ori = rots[index] From 5642625071f922bb4bd2050fba166bc65ac7c080 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 16:26:43 -0500 Subject: [PATCH 63/70] BugFix: Fix Examples --- .../mulit_phase_orientation.py | 7 +++- .../on_zone_orientation.py | 27 ++++++++++++-- .../single_phase_orientation.py | 37 +++++++++++++++---- pyxem/signals/indexation_results.py | 6 ++- pyxem/signals/polar_diffraction2d.py | 8 ++++ 5 files changed, 72 insertions(+), 13 deletions(-) diff --git a/examples/orientation_mapping/mulit_phase_orientation.py b/examples/orientation_mapping/mulit_phase_orientation.py index d25bf6d3f..296f51df0 100644 --- a/examples/orientation_mapping/mulit_phase_orientation.py +++ b/examples/orientation_mapping/mulit_phase_orientation.py @@ -28,12 +28,16 @@ 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( @@ -48,6 +52,7 @@ rotation=[rotations_bcc, rotations_fcc], max_excitation_error=0.1, reciprocal_radius=2, + with_direct_beam=False, ) orientation_map = polar_multi.get_orientation(sim) diff --git a/examples/orientation_mapping/on_zone_orientation.py b/examples/orientation_mapping/on_zone_orientation.py index 84159d73a..762f2e55d 100644 --- a/examples/orientation_mapping/on_zone_orientation.py +++ b/examples/orientation_mapping/on_zone_orientation.py @@ -9,10 +9,13 @@ 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 @@ -27,7 +30,14 @@ # %% -# Now we can get make the orientation map +# 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( @@ -38,9 +48,20 @@ ), 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) -simulated_si_tilt.plot(vmax="99th") +simulated_si_tilt.plot(vmax="99th", navigator=orientation_map.isig[2, 0]) simulated_si_tilt.add_marker( - orientation_map.to_single_phase_markers(annotate=True, text_color="w") + orientation_map.to_single_phase_markers( + annotate=True, text_color="w", lw=4, include_intensity=True, intesity_scale=40 + ) ) +# %% +# sphinx_gallery_thumbnail_number = 4 diff --git a/examples/orientation_mapping/single_phase_orientation.py b/examples/orientation_mapping/single_phase_orientation.py index a89dece22..189faeac3 100644 --- a/examples/orientation_mapping/single_phase_orientation.py +++ b/examples/orientation_mapping/single_phase_orientation.py @@ -13,6 +13,8 @@ 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 @@ -26,22 +28,41 @@ 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` +# 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 -) + 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) +ipf_markers = orientation_map.to_ipf_markers() navigator = orientation_map.to_navigator() - -simulated_si.plot(navigator=navigator, vmax="99th") +simulated_si.plot(navigator=navigator, vmax="95th") simulated_si.add_marker( - orientation_map.to_single_phase_markers(annotate=True, text_color="w") + orientation_map.to_single_phase_markers( + annotate=True, text_color="w", include_intensity=True, intesity_scale=40, lw=4 + ) ) +simulated_si.add_marker(ipf_markers) # %% -# sphinx_gallery_thumbnail_number = 3 +# sphinx_gallery_thumbnail_number = 4 diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 4d87495e7..176898070 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -18,6 +18,7 @@ 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 @@ -213,7 +214,7 @@ def rotation_from_orientation_map(result, rots): return ori -def extract_vectors_from_orientation_map(result, all_vectors): +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: @@ -373,6 +374,7 @@ def to_single_phase_vectors( output_signal_size=(), output_dtype=object, show_progressbar=False, + n_best_index=n_best_index, **kwargs, ) @@ -399,6 +401,8 @@ def to_crystal_map(self) -> CrystalMap: phase_index = phase_index.flatten() else: phases = PhaseList(self.simulation.phases) + if scan_unit is Undefined: + scan_unit = "px" return CrystalMap( rotations=rotations, diff --git a/pyxem/signals/polar_diffraction2d.py b/pyxem/signals/polar_diffraction2d.py index 476cd7f4e..a3c215b82 100644 --- a/pyxem/signals/polar_diffraction2d.py +++ b/pyxem/signals/polar_diffraction2d.py @@ -316,6 +316,7 @@ def get_orientation( n_keep=None, frac_keep=0.1, n_best=1, + gamma=0.5, normalize_templates=True, **kwargs, ): @@ -337,6 +338,11 @@ def get_orientation( 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. kwargs : dict Any additional options for the :meth:`~hyperspy.signal.BaseSignal.map` function. Returns @@ -344,6 +350,8 @@ def get_orientation( orientation : BaseSignal A signal with the orientation at each navigation position. """ + if gamma != 1: + self.data = self.data**gamma ( r_templates, theta_templates, From 050c332a4bb53b7805c93572e5fdc782d05bbda5 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 29 May 2024 16:39:08 -0500 Subject: [PATCH 64/70] Version: diffsims == 0.6.rc1 --- .github/workflows/build.yml | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index f7de899de..b6c3bd56c 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.6 orix==0.9.0 scikit-image==0.19.0 scikit-learn==1.0.0 LABEL: -oldest steps: - uses: actions/checkout@v3 diff --git a/setup.py b/setup.py index d74797e6c..a5b077dba 100644 --- a/setup.py +++ b/setup.py @@ -87,7 +87,7 @@ extras_require=extra_feature_requirements, install_requires=[ "dask", - "diffsims @ git+https://github.com/pyxem/diffsims.git@main", + "diffsims == 0.6.rc1", "hyperspy >= 2.0", "h5py", "lmfit >= 0.9.12", From 2d6789041bf2be53aaaa1e9b1ae15f85ed9113f1 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 31 May 2024 08:51:37 -0500 Subject: [PATCH 65/70] Visualization: Add marker for IPF visualization --- pyxem/signals/indexation_results.py | 203 +++++++++++++----- .../tests/signals/test_indexation_results.py | 24 +++ pyxem/utils/signal.py | 29 +++ 3 files changed, 202 insertions(+), 54 deletions(-) diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 176898070..326f955c5 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -35,12 +35,22 @@ 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): @@ -273,11 +283,16 @@ def simulation(self): def simulation(self, value): self.metadata.set_item("simulation", value) - def deepcopy(self): + def deepcopy(self, deepcopy_simulation=False): """Deepcopy the signal""" - self.simulation._phase_slider = None - self.simulation._rotation_slider = None - return super().deepcopy() + simulation = self.simulation + simulation._phase_slider = None + simulation._rotation_slider = None + self.simulation = None + new = super().deepcopy() + new.simulation = simulation + self.simulation = simulation + return new @property def num_rows(self): @@ -366,8 +381,7 @@ def to_single_phase_vectors( # Use vector data as signal in case of different vectors per navigation position vectors_signal = hs.signals.Signal1D(self.simulation.coordinates) - - return self.map( + v = self.map( extract_vectors_from_orientation_map, all_vectors=vectors_signal, inplace=False, @@ -377,6 +391,8 @@ def to_single_phase_vectors( 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 `orix.CrystalMap` object""" @@ -413,7 +429,7 @@ def to_crystal_map(self) -> CrystalMap: phase_list=phases, ) - def to_ipf_markers(self): + def to_ipf_markers(self, offset=0.85, scale=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. @@ -424,59 +440,36 @@ def to_ipf_markers(self): ) 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() - sector = self.simulation.phases.point_group.fundamental_sector - labels = _get_ipf_axes_labels( - sector.vertices, symmetry=self.simulation.phases.point_group - ) - s = StereographicProjection() vectors = orients * Vector3d.zvector() - edges = _closed_edges_in_hemisphere(sector.edges, sector) vectors = vectors.in_fundamental_sector(self.simulation.phases.point_group) + s = StereographicProjection() x, y = s.vector2xy(vectors) - ex, ey = s.vector2xy(edges) - tx, ty = s.vector2xy(sector.vertices) - 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) - original_offset = np.vstack((ex, ey)).T - texts_offset = np.vstack((tx, ty)).T - mins, maxes = original_offset.min(axis=0), original_offset.max(axis=0) - - original_offset = ( - (original_offset - ((maxes + mins) / 2)) / (maxes - mins) * 0.2 - ) - original_offset = original_offset + 0.85 - - texts_offset = (texts_offset - ((maxes + mins) / 2)) / (maxes - mins) * 0.2 - texts_offset = texts_offset + 0.85 for i in np.ndindex(offsets.shape): off = np.vstack((x[i], y[i])).T - norm_points = (off - ((maxes + mins) / 2)) / (maxes - mins) * 0.2 - norm_points = norm_points + 0.85 + 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=[[0.85, 0.85]], - widths=(0.3,), + offsets=[[offset, offset]], + widths=(scale + scale / 2,), units="width", offset_transform="axes", facecolor="white", edgecolor="black", ) - polygon_sector = hs.plot.markers.Polygons( - verts=original_offset[np.newaxis], - transform="axes", - alpha=1, - facecolor="none", - ) best_points = hs.plot.markers.Points( offsets=offsets.T, @@ -486,13 +479,6 @@ def to_ipf_markers(self): facecolor="green", ) - texts = hs.plot.markers.Texts( - texts=labels, - offsets=texts_offset, - sizes=(1,), - offset_transform="axes", - facecolor="k", - ) return square, polygon_sector, best_points, texts def to_single_phase_markers( @@ -526,7 +512,11 @@ def to_single_phase_markers( 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. """ + if lazy_output is None: + lazy_output = self._lazy if text_kwargs is None: text_kwargs = dict() if annotation_shift is None: @@ -535,10 +525,10 @@ def to_single_phase_markers( 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=lazy_output, navigation_chunks=navigation_chunks + lazy_output=True, navigation_chunks=navigation_chunks ) color = marker_colors[n % len(marker_colors)] if include_intensity: @@ -550,6 +540,7 @@ def to_single_phase_markers( output_dtype=object, output_signal_size=(), navigation_chunks=navigation_chunks, + lazy_output=True, ).data.T kwargs["sizes"] = intensity @@ -560,27 +551,40 @@ def to_single_phase_markers( 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 ) - yield markers + all_markers.append(markers) if annotate: texts = vectors.map( vectors_to_text, inplace=False, - lazy_output=lazy_output, + 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) + 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 ) - yield text_markers + 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: AxesManager, n_best: int = 1, marker_color: str = "red" @@ -632,7 +636,88 @@ def marker_generator(entry): ) yield markers - def to_navigator(self, direction: Vector3d = Vector3d.zvector()): + def _get_ipf_outline(self, include_labels=True, offset=0.85, scale=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.""" + # 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=0.85, scale=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.""" + + 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. """ @@ -644,7 +729,17 @@ def to_navigator(self, direction: Vector3d = Vector3d.zvector()): 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, annotate=False, **plot_kwargs): @@ -723,7 +818,7 @@ def to_crystal_map(self): ) -class LazyOrientationMap(LazySignal, OrientationMap): +class LazyOrientationMap(OrientationMap, LazySignal): pass diff --git a/pyxem/tests/signals/test_indexation_results.py b/pyxem/tests/signals/test_indexation_results.py index 6688ea2db..4f3c9d032 100644 --- a/pyxem/tests/signals/test_indexation_results.py +++ b/pyxem/tests/signals/test_indexation_results.py @@ -281,3 +281,27 @@ 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) + + def test_to_markers(self, simple_multi_rot_orientation_result): + orientations, rotations = simple_multi_rot_orientation_result + markers = orientations.to_single_phase_markers(lazy_output=True) + assert isinstance(markers[0], hs.plot.markers.Markers) + + def test_to_ipf_markers(self, simple_multi_rot_orientation_result): + orientations, rotations = 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 = 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 = 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 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 From 772dafa6a52893ceac15f2379a5978980b452a69 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 31 May 2024 10:37:42 -0500 Subject: [PATCH 66/70] Testing: Added new tests to cover new class and fixed up documentation. --- .../mulit_phase_orientation.py | 2 +- .../on_zone_orientation.py | 9 +- .../single_phase_orientation.py | 11 +- pyxem/signals/indexation_results.py | 223 +++++++++++++++--- .../tests/signals/test_indexation_results.py | 96 ++++++-- .../tests/signals/test_polar_diffraction2d.py | 18 -- 6 files changed, 268 insertions(+), 91 deletions(-) diff --git a/examples/orientation_mapping/mulit_phase_orientation.py b/examples/orientation_mapping/mulit_phase_orientation.py index 296f51df0..c40599302 100644 --- a/examples/orientation_mapping/mulit_phase_orientation.py +++ b/examples/orientation_mapping/mulit_phase_orientation.py @@ -1,5 +1,5 @@ """ -Mulit Phase Orientation Mapping +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 diff --git a/examples/orientation_mapping/on_zone_orientation.py b/examples/orientation_mapping/on_zone_orientation.py index 762f2e55d..e703d7add 100644 --- a/examples/orientation_mapping/on_zone_orientation.py +++ b/examples/orientation_mapping/on_zone_orientation.py @@ -57,11 +57,4 @@ # rotation as a navigator or plot it directly. orientation_map = polar_si_tilt.get_orientation(sim) -simulated_si_tilt.plot(vmax="99th", navigator=orientation_map.isig[2, 0]) -simulated_si_tilt.add_marker( - orientation_map.to_single_phase_markers( - annotate=True, text_color="w", lw=4, include_intensity=True, intesity_scale=40 - ) -) -# %% -# sphinx_gallery_thumbnail_number = 4 +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 index 189faeac3..6687ac366 100644 --- a/examples/orientation_mapping/single_phase_orientation.py +++ b/examples/orientation_mapping/single_phase_orientation.py @@ -53,16 +53,9 @@ # 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) -ipf_markers = orientation_map.to_ipf_markers() -navigator = orientation_map.to_navigator() -simulated_si.plot(navigator=navigator, vmax="95th") -simulated_si.add_marker( - orientation_map.to_single_phase_markers( - annotate=True, text_color="w", include_intensity=True, intesity_scale=40, lw=4 - ) -) -simulated_si.add_marker(ipf_markers) +orientation_map.plot_over_signal(simulated_si, vmax="96th") # %% # sphinx_gallery_thumbnail_number = 4 diff --git a/pyxem/signals/indexation_results.py b/pyxem/signals/indexation_results.py index 326f955c5..e35c428f3 100644 --- a/pyxem/signals/indexation_results.py +++ b/pyxem/signals/indexation_results.py @@ -23,7 +23,7 @@ import hyperspy.api as hs from hyperspy._signals.lazy import LazySignal from hyperspy.signal import BaseSignal -from hyperspy.axes import AxesManager +from hyperspy.axes import AxesManager, BaseDataAxis import numpy as np from orix.crystal_map import CrystalMap, Phase, PhaseList from orix.quaternion import Rotation, Orientation @@ -277,6 +277,9 @@ class OrientationMap(DiffractionVectors2D): @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 @@ -284,14 +287,19 @@ def simulation(self, value): self.metadata.set_item("simulation", value) def deepcopy(self, deepcopy_simulation=False): - """Deepcopy the signal""" - simulation = self.simulation - simulation._phase_slider = None - simulation._rotation_slider = None - self.simulation = None + """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() - new.simulation = simulation - self.simulation = simulation + if not deepcopy_simulation: + new.simulation = simulation + self.simulation = simulation return new @property @@ -301,12 +309,20 @@ def num_rows(self): 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: - warn("Computing the signal") - self.compute() + 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, @@ -342,6 +358,11 @@ def to_phase_index(self): 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( @@ -367,13 +388,14 @@ def to_single_phase_orientations(self, **kwargs) -> Orientation: 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. + """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: @@ -395,7 +417,13 @@ def to_single_phase_vectors( return v def to_crystal_map(self) -> CrystalMap: - """Convert the orientation map to an `orix.CrystalMap` object""" + """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 " @@ -429,10 +457,17 @@ def to_crystal_map(self) -> CrystalMap: phase_list=phases, ) - def to_ipf_markers(self, offset=0.85, scale=0.2): + 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( @@ -481,7 +516,7 @@ def to_ipf_markers(self, offset=0.85, scale=0.2): return square, polygon_sector, best_points, texts - def to_single_phase_markers( + def to_markers( self, n_best: int = 1, annotate: bool = False, @@ -493,7 +528,7 @@ def to_single_phase_markers( 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 @@ -514,6 +549,12 @@ def to_single_phase_markers( 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 @@ -587,8 +628,36 @@ def to_single_phase_markers( return all_markers def to_single_phase_polar_markers( - self, signal_axes: AxesManager, n_best: int = 1, marker_color: str = "red" + 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, @@ -597,6 +666,9 @@ def to_single_phase_polar_markers( 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() @@ -624,21 +696,52 @@ def marker_generator(entry): 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=marker_color, sizes=(10,) + markers_signal, facecolor="none", edgecolor=color, **kwargs ) - yield markers - def _get_ipf_outline(self, include_labels=True, offset=0.85, scale=0.2): + 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.""" + 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() @@ -676,10 +779,27 @@ def _get_ipf_outline(self, include_labels=True, offset=0.85, scale=0.2): ) return polygon_sector, texts, maxes, mins - def get_ipf_annotation_markers(self, offset=0.85, scale=0.2): + 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.""" + 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) @@ -720,6 +840,17 @@ def to_ipf_colormap( ): """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) @@ -742,28 +873,44 @@ def to_ipf_colormap( ) return s - def plot_over_signal(self, signal, annotate=False, **plot_kwargs): + 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. - annotate: bool - If True, the euler rotation and the correlation will be annotated on the plot using - the `Texts` class from hyperspy. - - Notes - ----- - The kwargs are passed to the `signal.plot` function call + 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_navigator() - signal.plot(navigator=nav, **plot_kwargs) - signal.add_marker(self.to_single_phase_markers(1, annotate=annotate)) - - def plot_inplane_rotation(self, **kwargs): - """Plot the in-plane rotation of the orientation map as a 2D map.""" - pass + 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: diff --git a/pyxem/tests/signals/test_indexation_results.py b/pyxem/tests/signals/test_indexation_results.py index 4f3c9d032..3f3096afa 100644 --- a/pyxem/tests/signals/test_indexation_results.py +++ b/pyxem/tests/signals/test_indexation_results.py @@ -20,6 +20,7 @@ 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 @@ -162,7 +163,7 @@ class TestOrientationResult: examples provided in the documentation. """ - @pytest.fixture + @pytest.fixture(scope="class") def single_rot_orientation_result(self): s = si_tilt() s.calibration.center = None @@ -179,11 +180,12 @@ def single_rot_orientation_result(self): ), 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 + @pytest.fixture(scope="class") def multi_rot_orientation_result(self): s, r = si_grains(return_rotations=True) s.calibration.center = None @@ -196,12 +198,17 @@ def multi_rot_orientation_result(self): resolution=1, point_group=phase.point_group ) sims = generator.calculate_diffraction2d( - phase, rotation=rotations, max_excitation_error=0.1, reciprocal_radius=2 + 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 + @pytest.fixture(scope="class") def simple_multi_rot_orientation_result(self): s, r = si_grains_simple(return_rotations=True) s.calibration.center = None @@ -214,12 +221,16 @@ def simple_multi_rot_orientation_result(self): resolution=1, point_group=phase.point_group ) sims = generator.calculate_diffraction2d( - phase, rotation=rotations, max_excitation_error=0.1, reciprocal_radius=2 + phase, + rotation=rotations, + max_excitation_error=0.1, + reciprocal_radius=2, + with_direct_beam=True, ) - orientations = polar.get_orientation(sims) - return orientations, r + orientations = polar.get_orientation(sims, alpha=1) # in this case we + return orientations, r, s - @pytest.fixture + @pytest.fixture(scope="class") def multi_phase_orientation_result(self): s = fe_multi_phase_grains() s.calibration.center = None @@ -242,7 +253,7 @@ def multi_phase_orientation_result(self): max_excitation_error=0.1, reciprocal_radius=2, ) - orientations = polar.get_orientation(sims) + orientations = polar.get_orientation(sims, n_best=3) return orientations def test_tilt_orientation_result(self, single_rot_orientation_result): @@ -261,7 +272,7 @@ def test_tilt_orientation_result(self, single_rot_orientation_result): assert np.all(degrees_between[:, 5:] <= 1) def test_grain_orientation_result(self, simple_multi_rot_orientation_result): - orientations, rotations = 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() @@ -272,7 +283,7 @@ def test_grain_orientation_result(self, simple_multi_rot_orientation_result): assert np.all(np.min(degrees_between, axis=2) <= 2) def test_to_crystal_map(self, simple_multi_rot_orientation_result): - orientations, rotations = 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) @@ -282,26 +293,77 @@ def test_to_crystal_map_multi_phase(self, multi_phase_orientation_result): assert isinstance(crystal_map, CrystalMap) assert np.all(crystal_map.phase_id < 2) - def test_to_markers(self, simple_multi_rot_orientation_result): - orientations, rotations = simple_multi_rot_orientation_result - markers = orientations.to_single_phase_markers(lazy_output=True) + @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 = 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 = 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 = simple_multi_rot_orientation_result + 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_polar_diffraction2d.py b/pyxem/tests/signals/test_polar_diffraction2d.py index e708356ca..639b7ec1c 100644 --- a/pyxem/tests/signals/test_polar_diffraction2d.py +++ b/pyxem/tests/signals/test_polar_diffraction2d.py @@ -251,24 +251,6 @@ def test_decomposition_class_assignment(self, diffraction_pattern): assert isinstance(s, PolarDiffraction2D) -class TestOrientationMap: - @pytest.fixture - def polar_image(self, diffraction_pattern): - from pyxem.signals import PolarDiffraction2D - - pol_image = np.zeros((2, 2, 20, 60)) - pol_image[:, :, 4:6, 9:11] = 1 - pol_image[:, :, 4:6, 29:31] = 1 - pol_image[:, :, 4:6, 49:51] = 1 - - pol_image[:, :, 9:11, 9:11] = 1 - pol_image[:, :, 9:11, 29:31] = 1 - pol_image[:, :, 9:11, 49:51] = 1 - pol = PolarDiffraction2D(pol_image) - pol.axes_manager[3].scale = 0.5 - return pol - - class TestSubtractingDiffractionBackground: @pytest.fixture def noisy_data(self): From c479e07c36742520838eee7ec307a2e751950099 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 31 May 2024 10:50:52 -0500 Subject: [PATCH 67/70] Testing: Update Oldest --- .github/workflows/build.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index b6c3bd56c..206089871 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.6.rc1 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.6 orix==0.12.1 scikit-image==0.19.0 scikit-learn==1.0.0 LABEL: -oldest steps: - uses: actions/checkout@v3 From 80f542880f317583f3d661a003ad4ae4e519f6b1 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 31 May 2024 11:01:59 -0500 Subject: [PATCH 68/70] Refactor: rename simulated_fe.py --- pyxem/data/__init__.py | 2 +- pyxem/data/{simulation_fe.py => simulated_fe.py} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename pyxem/data/{simulation_fe.py => simulated_fe.py} (100%) diff --git a/pyxem/data/__init__.py b/pyxem/data/__init__.py index 0a73806c9..6f392494a 100644 --- a/pyxem/data/__init__.py +++ b/pyxem/data/__init__.py @@ -37,7 +37,7 @@ si_grains_simple, si_rotations_line, ) -from pyxem.data.simulation_fe import fe_bcc_phase, fe_fcc_phase, fe_multi_phase_grains +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, diff --git a/pyxem/data/simulation_fe.py b/pyxem/data/simulated_fe.py similarity index 100% rename from pyxem/data/simulation_fe.py rename to pyxem/data/simulated_fe.py From 2cf1c95d965deb1aff3a354118727c5bef5f3a1b Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 31 May 2024 11:07:25 -0500 Subject: [PATCH 69/70] Testing: ++ min matplotlib version --- .github/workflows/build.yml | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 206089871..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.6.rc1 hyperspy~=2.0.rc0 lmfit==0.9.12 matplotlib==3.6 orix==0.12.1 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/setup.py b/setup.py index a5b077dba..bcfd0159b 100644 --- a/setup.py +++ b/setup.py @@ -91,7 +91,7 @@ "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 From 2d9383769e1375d6dd2dca811f7c872fe4ffc1ba Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 31 May 2024 15:36:04 -0500 Subject: [PATCH 70/70] Refactor: References --- doc/bibliography.bib | 9 +++++++++ examples/orientation_mapping/mulit_phase_orientation.py | 2 ++ examples/orientation_mapping/on_zone_orientation.py | 2 ++ examples/orientation_mapping/single_phase_orientation.py | 2 ++ pyxem/signals/polar_diffraction2d.py | 7 ++++++- 5 files changed, 21 insertions(+), 1 deletion(-) diff --git a/doc/bibliography.bib b/doc/bibliography.bib index 4800f1a58..1b41a0dc7 100644 --- a/doc/bibliography.bib +++ b/doc/bibliography.bib @@ -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/examples/orientation_mapping/mulit_phase_orientation.py b/examples/orientation_mapping/mulit_phase_orientation.py index c40599302..ef3463b3c 100644 --- a/examples/orientation_mapping/mulit_phase_orientation.py +++ b/examples/orientation_mapping/mulit_phase_orientation.py @@ -4,6 +4,8 @@ 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 diff --git a/examples/orientation_mapping/on_zone_orientation.py b/examples/orientation_mapping/on_zone_orientation.py index e703d7add..383b1a555 100644 --- a/examples/orientation_mapping/on_zone_orientation.py +++ b/examples/orientation_mapping/on_zone_orientation.py @@ -4,6 +4,8 @@ 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 diff --git a/examples/orientation_mapping/single_phase_orientation.py b/examples/orientation_mapping/single_phase_orientation.py index 6687ac366..f111f265d 100644 --- a/examples/orientation_mapping/single_phase_orientation.py +++ b/examples/orientation_mapping/single_phase_orientation.py @@ -4,6 +4,8 @@ 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 diff --git a/pyxem/signals/polar_diffraction2d.py b/pyxem/signals/polar_diffraction2d.py index a3c215b82..9d620f498 100644 --- a/pyxem/signals/polar_diffraction2d.py +++ b/pyxem/signals/polar_diffraction2d.py @@ -342,13 +342,18 @@ def get_orientation( 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. + 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