diff --git a/tests/test_signals/test_ebsd_master_pattern.py b/tests/test_signals/test_ebsd_master_pattern.py index 114fa48b..8448b447 100644 --- a/tests/test_signals/test_ebsd_master_pattern.py +++ b/tests/test_signals/test_ebsd_master_pattern.py @@ -14,8 +14,10 @@ # # You should have received a copy of the GNU General Public License # along with kikuchipy. If not, see . +from copy import deepcopy import dask.array as da +from diffsims.crystallography import ReciprocalLatticeVector from hyperspy._signals.signal2d import Signal2D import hyperspy.api as hs import numpy as np @@ -24,6 +26,7 @@ import pytest import kikuchipy as kp +from kikuchipy._rotation import _rotate_vector from kikuchipy.signals.util._master_pattern import ( _get_direction_cosines_for_fixed_pc, _get_direction_cosines_for_varying_pc, @@ -771,3 +774,487 @@ def test_vector2xy(self): ] assert np.allclose(_vector2lambert.py_func(xyz), lambert_xy) assert np.allclose(_vector2lambert(xyz), lambert_xy) + + +class TestFitPatternDetectorOrientation: + """ + Test the fit between an EBSD pattern generated + from a master pattern and an associated + GeometricalKikuchiPatternSimulation for different + detector orientations. + """ + + detector = kp.detectors.EBSDDetector( + shape=(480, 640), px_size=50, pc=(20, 20, 15000), convention="emsoft4" + ) + + phase = kp.data.nickel_ebsd_master_pattern_small().phase + + pqr_p = [(1, 1, 1), (2, 0, 0), (2, 2, 0), (3, 1, 1)] + + ref = ReciprocalLatticeVector(phase=phase, hkl=pqr_p) + ref = ref.symmetrise() + + simulator = kp.simulations.KikuchiPatternSimulator(ref) + + rotations = Rotation.from_euler([23, 14, 5], degrees=True) + + # Transformation from CSs to cartesian crystal reference frame CSc + u_o = rotations.to_matrix().squeeze() + + # Transformation from CSc to direct crystal reference frame CSk + u_a = phase.structure.lattice.base + + def setup_detector_sim_and_u_os(self, tilt, azimuthal, euler_2): + """Setup the detector with the correct tilt, azimuthal and + euler_2 angles, get the GeometricalKikuchiPatternSimulation from + the simulator and calculate the orientation matrix + transforming from the detector reference frame CSd to + the crystal cartesian reference frame CSc. Return these + three items. + + Parameters + ---------- + tilt : Float + The detector tilt angle in degrees (i.e. detector.tilt). + Detector Euler angle PHI (EBSDDetector.euler[1]) == 90 + tilt + azimuthal : Float + The detector azimuthal angle in degrees (i.e. detector.azimuthal). + Detector Euler angle phi1 (EBSDDetector.euler[0]) == azimuthal + euler_2 : Float + The detector Euler angle phi2 (EBSDDetector.euler[2]) in deg. + """ + det = self.detector.deepcopy() + det.euler = np.array([azimuthal, 90 + tilt, euler_2]) + + sim_lines = self.simulator.on_detector(det, self.rotations) + + u_os = np.matmul(self.u_o, det.u_s) + + return det, sim_lines, u_os + + def setup_za_and_cds(self, zone_axes, sim_lines, det): + """Find the indices of the zone_axes in the + GeometricalKikuchiPatternSimulation (sim_lines). + Find the gnomonic coordinates of these zone axes. + Then obtain convert these into pixel coordinates + on the detector and round to the nearest pixel. + Return the indices of the zone axes, the + conversion dict from pix to gn and back, and + the coordinates of the nearest detector pixels + to the three zone axes. + """ + za_idx = [ + index_row_in_array(sim_lines._zone_axes.vector.uvw, i) for i in zone_axes + ] + + x_gn = sim_lines._zone_axes.x_gnomonic[0, za_idx] + y_gn = sim_lines._zone_axes.y_gnomonic[0, za_idx] + cds_gn = np.stack((x_gn, y_gn), axis=1) + + conversions = get_coordinate_conversions(det) + cds_px = convert_coordinates(cds_gn, "gn_to_pix", conversions) + + cds_px_int = np.around(cds_px, decimals=0).astype(int) + + return za_idx, conversions, cds_px_int + + def get_a_angles(self, sim_lines, za_idx, u_os): + """Check that the GeometricalKikuchiPatternSimulation + is self-consistent. We do this by taking the vectors + in the detector reference frame of 3 zone axes in + the simulation, transforming these vectors from CSd + to the cartesian crystal reference frame, CSc, and + measuring the angle between these and vectors + transformed from the uvw indices of the zone axes + into the cartesian crystal reference frame + i.e. CSk to CSc. These angles (a_ang) should all be + zero if the GeometricalKikuchiPatternSimulation is + self-consistent. This is tested later. + The a_ang are returned. + """ + # CSk to CSc: + uvw = sim_lines._zone_axes.vector.uvw[za_idx] + d_vec = np.matmul(uvw, self.u_a) + N1 = np.sqrt(np.sum(np.square(d_vec), axis=1)) + d_vec_n = d_vec / np.expand_dims(N1, 1) + + # CSd to CSc: + za_vec = sim_lines._zone_axes.vector_detector[0, za_idx].data + za_vec_trans = np.matmul(za_vec, np.linalg.inv(u_os)).squeeze() + N2 = np.sqrt(np.sum(np.square(za_vec_trans), axis=1)) + za_vec_trans_n = za_vec_trans / np.expand_dims(N2, 1) + + # angles between the two sets of vectors: + a_ang = np.array( + [ + np.rad2deg( + np.arccos( + np.around( + np.dot(za_vec_trans_n[i, :], d_vec_n[i, :]), decimals=8 + ) + ) + ) + for i in range(3) + ] + ) + + return a_ang, d_vec_n + + def get_d_ang(self, cds_px_int, conversions, u_os, d_vec_n, det): + """Find the gnomonic coordinates of the nearest + pixel centre to each of the 3 zone axes. Turn them + into vectors and transform them from CSd to CSc using + the same transformation as the + GeometricalKikuchiPatternSimulation. Then + calculate the angles (d_ang) between these vectors + and the vectors representing the zone axes in CSc but + calculated from CSk to CSc (d_vec_n). + """ + # Here we add 0.5 because pixel centres are used by the direction + # cosines method and we need to take the same coords for this + # alternative approach. + cds_gn_int = convert_coordinates(cds_px_int + 0.5, "pix_to_gn", conversions) + vecs = np.hstack( + ( + cds_gn_int * det.pcz, + np.repeat(np.atleast_2d(det.pcz), cds_gn_int.shape[0], axis=0), + ) + ) + N3 = np.sqrt(np.sum(np.square(vecs), axis=1)) + vecs_n = vecs / np.expand_dims(N3, 1) + + # CSd to CSc: + dddd = np.matmul(vecs_n, np.linalg.inv(u_os)).squeeze() + + d_ang = np.array( + [ + np.rad2deg( + np.arccos(np.around(np.dot(dddd[i, :], d_vec_n[i, :]), decimals=8)) + ) + for i in range(3) + ] + ) + + return d_ang, dddd + + def get_r_ang_and_n_ang(self, cds_px_int, det, d_vec_n, dddd): + """Calculate the direction cosines of all the + detector pixels then rotate them to account + for the crystal orientation. The resulting + vectors are in CSc. Select the vectors + corresponding to the pixel centres representing the 3 + zone axes. Then calculate the angles (r_ang) between + these vectors and the vectors representing the + zone axes in CSc but calculated from CSk to CSc + (d_vec_n). Finally calculate the angles (n_ang) between + the two sets of vectors representing the centres of the + nearest pixels to the zone axes (one set is from the direction + cosines of the detecotr, the other is from the + GeometricalKikuchiPatternSimulation). The angles are + zero if the transformations used for the direction + cosines and for the GeometricalKikuchiPatternSimulation + are the same (this is tested later). + """ + # swap columns for i,j array indexing: + cds_px_int_ij = cds_px_int[:, ::-1] + + # all detector pixels: + r_g_array = _get_direction_cosines_from_detector(det) + r_g_array_rot = _rotate_vector(self.rotations.data[0], r_g_array) + rgarrrot_reshaped = r_g_array_rot.reshape((*self.detector.shape, 3)) + + # select vectors corresponding to the nearest pixels to the chosen zone axes + rgrar_vec = rgarrrot_reshaped[cds_px_int_ij[:, 0], cds_px_int_ij[:, 1]] + + r_ang = np.array( + [ + np.rad2deg( + np.arccos( + np.around(np.dot(d_vec_n[i, :], rgrar_vec[i, :]), decimals=8) + ) + ) + for i in range(3) + ] + ) + + n_ang = np.array( + [ + np.rad2deg( + np.arccos( + np.around(np.dot(dddd[i, :], rgrar_vec[i, :]), decimals=8) + ) + ) + for i in range(3) + ] + ) + + return r_ang, n_ang + + def calculate_fit(self, tilt, azimuthal, euler_2, zone_axes): + """ + Calculates four sets of angles with which the fit + between the EBSD pattern simulated from a master + pattern, and the GeometricalKikuchiPatternSimulation + generated using the same parameters, can be evaluated. + The function can be tested with different values of + the detector tilt, azimuthal and euler_2 angles and + appropriate zone axes (indices uvw) which appear on + the pattern under those conditions. + + The approach has several steps: + + 1. Check that the GeometricalKikuchiPatternSimulation + is self-consistent. We do this by taking the vectors + in the detector reference frame of 3 zone axes in + the simulation, transforming these vectors from CSd + to the cartesian crystal reference frame, CSc, and + measuring the angle between these and vectors + transformed from the uvw indices of the zone axes + into the cartesian crystal reference frame + i.e. CSk to CSc. These angles (a_ang) should all be + zero if the GeometricalKikuchiPatternSimulation is + self-consistent. + + 2. Find the gnomonic coordinates of the nearest + pixel centre to the 3 zone axes. This enables us + to check the vectors corresponding to the detector + pixels. We do two things with these coordinates. + a) turn them into vectors and transform them + from CSd to CSc using the same transformation as + the GeometricalKikuchiPatternSimulation. Then + calculate the angles (d_ang) between these vectors + and the vectors representing the zone axes in CSc. + b) calculate the direction cosines of all the + detector pixels then rotate them to account + for the crystal orientation. The resulting + vectors are in CSc but calculated using + different functions. Select the vectors + corresponding to the pixels representing the 3 + zone axes. Then calculate the angles (r_ang) between + these vectors and the vectors representing the + zone axes in CSc. + These angles should be the same for a) and b), + meaning that the angle between the zone axes + and the centre of the nearest pixel is the + same for both transformation routes. + + 3. Finally calculate the angles (n_ang) between the two + sets of vectors representing the centres of the + nearest pixels to the zone axes. The angles are + zero if the transformations used for the direction + cosines and for the GeometricalKikuchiPatternSimulation + are the same. + + + Parameters + ---------- + tilt : Float + The detector tilt angle in degrees (i.e. detector.tilt). + Detector Euler angle PHI (EBSDDetector.euler[1]) == 90 + tilt + azimuthal : Float + The detector azimuthal angle in degrees (i.e. detector.azimuthal). + Detector Euler angle phi1 (EBSDDetector.euler[0]) == azimuthal + euler_2 : Float + The detector Euler angle phi2 (EBSDDetector.euler[2]) in deg. + zone_axes : List or np.ndarray + List/array containing three lists, each containing a set + of uvw indices describing a zone axis on the pattern, + e.g. [[0,0,1], [1,1,0], [1,2,3]]. + + Returns + ------- + a_ang : np.ndarray + The angles in degrees, between vectors in CSc calculated + 1) from uvw indeces in CSk and 2) from vectors in CSd. + These should all be zero for self-consistency of + the GeometricalKikuchiPatternSimulation. + d_ang : np.ndarray + The angles in degrees, between vectors representing + 3 zone axes on the pattern and vectors representing + the nearest pixel centres to the same zone axes. + Both sets of vectors are transformed into CSc. + The transformation for both sets was the one used + by the GeometricalKikuchiPatternSimulation. + r_ang : np.ndarray + The angles in degrees, between vectors representing + 3 zone axes on the pattern and vectors representing + the nearest pixel centres to the same zone axes but + this time, the pixel centre vectors use the + transformation of the direction cosines of the + detector pixels. + n_ang : np.ndarray + The angles in degrees between the two sets of vectors + representing the centres of the nearest pixels to 3 + zone axes on the pattern. Both sets of vectors are + in CSc. The transformation for one set was the one + used by the GeometricalKikuchiPatternSimulation, + and the one used by the other set was for the + direction cosines of the detector. These angles + should all be zero if the pattern and simulation + match. + """ + det, sim_lines, u_os = self.setup_detector_sim_and_u_os( + tilt, azimuthal, euler_2 + ) + za_idx, conversions, cds_px_int = self.setup_za_and_cds( + zone_axes, sim_lines, det + ) + a_ang, d_vec_n = self.get_a_angles(sim_lines, za_idx, u_os) + d_ang, dddd = self.get_d_ang(cds_px_int, conversions, u_os, d_vec_n, det) + r_ang, n_ang = self.get_r_ang_and_n_ang(cds_px_int, det, d_vec_n, dddd) + + return a_ang, d_ang, r_ang, n_ang + + @pytest.mark.parametrize( + "tilt, azimuthal, euler_2, zone_axes", + [ + (0.0, 0.0, 0.0, [[1, 0, 1], [0, 0, 1], [1, 1, 2]]), + (0.0, 0.0, 1.2, [[1, 0, 1], [0, 0, 1], [1, 1, 2]]), + (40.0, 0.0, 0.0, [[1, 0, 1], [1, 0, 0], [1, -2, 1]]), + (40.0, 0.0, 1.2, [[1, 0, 1], [1, 0, 0], [1, -2, 1]]), + (0.0, 40.0, 0.0, [[1, 0, 1], [1, 1, 0], [1, 2, 1]]), + (0.0, 40.0, 1.2, [[1, 0, 1], [1, 1, 0], [1, 2, 1]]), + (40.0, 40.0, 0.0, [[1, 0, 1], [1, 0, 0], [3, 1, 0]]), + (40.0, 40.0, 1.2, [[1, 0, 1], [1, 0, 0], [3, 1, 0]]), + ], + ) + def test_fit_detector_orientation(self, tilt, azimuthal, euler_2, zone_axes): + """ + Check that the EBSD pattern simulated from a master + pattern and the associated + GeometricalKikuchiPatternSimulation match perfectly, + for various detector orientations. + + 4 sets of angles are returned by self.calculate_fit(). + See the doctstring of that function for details. + + Here we assert that the first set of angles are all + zero, that the second and third sets are equal, and + that the fourth set are all zero. If these conditions + are all met, the GeometricalKikuchiPatternSimulation + should match the EBSD pattern simulated from a + master pattern perfectly for the given detector + orientations. + + + Parameters + ---------- + tilt : Float + The detector tilt angle in degrees (i.e. detector.tilt). + Detector Euler angle PHI (EBSDDetector.euler[1]) == 90 + tilt + azimuthal : Float + The detector azimuthal angle in degrees (i.e. detector.azimuthal). + Detector Euler angle phi1 (EBSDDetector.euler[0]) == azimuthal + euler_2 : Float + The detector Euler angle phi2 (EBSDDetector.euler[2]) in deg. + zone_axes : List or np.ndarray + List/array containing three lists, each containing a set + of uvw indices describing a zone axis on the pattern, + e.g. [[0,0,1], [1,1,0], [1,2,3]]. + + Returns + ------- + None. + + """ + angles = self.calculate_fit(tilt, azimuthal, euler_2, zone_axes) + + assert np.allclose(angles[0], 0.0) + assert np.allclose(angles[1], angles[2]) + assert np.allclose(angles[3], 0.0) + + +def index_row_in_array(myarray, myrow): + """Check if the row "myrow" is present in the array "myarray". + If it is, return an int containing the row index of the first + occurrence. If the row is not present, return None. + """ + loc = np.where((myarray == myrow).all(-1))[0] + if len(loc) > 0: + return loc[0] + return None + + +def get_coordinate_conversions(detector): + """Return a dict 'conversions' containing the keys + "pix_to_gn" and "gn_to_pix". Under each of these keys + is a further dict with the keys: "m_x", "c_x", "m_y" + and "c_y". These are the slope (m) and y-intercept (c) + corresponding to y = mx + c, which describes the + conversion of the coordinates from pixels to + gnomonic or vice versa. + + """ + m_pix_to_gn_x = ( + detector.gnomonic_bounds[..., 1] - detector.gnomonic_bounds[..., 0] + ) / (detector.bounds[1] + 1) + c_pix_to_gn_x = deepcopy(detector.gnomonic_bounds[..., 0]) + + m_pix_to_gn_y = ( + detector.gnomonic_bounds[..., 2] - detector.gnomonic_bounds[..., 3] + ) / (detector.bounds[3] + 1) + c_pix_to_gn_y = deepcopy(detector.gnomonic_bounds[..., 3]) + + m_gn_to_pix_x = 1 / m_pix_to_gn_x + c_gn_to_pix_x = -c_pix_to_gn_x / m_pix_to_gn_x + + m_gn_to_pix_y = 1 / m_pix_to_gn_y + c_gn_to_pix_y = -c_pix_to_gn_y / m_pix_to_gn_y + + conversions = { + "pix_to_gn": { + "m_x": m_pix_to_gn_x, + "c_x": c_pix_to_gn_x, + "m_y": m_pix_to_gn_y, + "c_y": c_pix_to_gn_y, + }, + "gn_to_pix": { + "m_x": m_gn_to_pix_x, + "c_x": c_gn_to_pix_x, + "m_y": m_gn_to_pix_y, + "c_y": c_gn_to_pix_y, + }, + } + + return conversions + + +def convert_coordinates( + coords: np.array, + direction: str, + conversions: dict, +) -> np.ndarray: + """ + Convert coordinates from gnomonic to pixels or from pixels + to gnomonic. + + Parameters + ---------- + coords + An array of coordinates whereby the x and y coordinates + to be converted are stored in the last axis + direction + Either "pix_to_gn" or "gn_to_pix", depending on the + direction of conversion needed. + conversions + dict + + + Returns + ------- + coords_out + Array of coords but with values converted as specified + by direction. The shape is either the same as the input + or is the naviagation shape then the shape of the input. + + """ + coords_out = np.zeros_like(coords, dtype=np.float32) + + coords_out[..., 0] = ( + conversions[direction]["m_x"] * coords[..., 0] + conversions[direction]["c_x"] + ) + coords_out[..., 1] = ( + conversions[direction]["m_y"] * coords[..., 1] + conversions[direction]["c_y"] + ) + return coords_out