# Copyright 2023 Mechanics of Microstructures Group
# at The University of Manchester
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import numpy as np
from skimage import morphology as mph
import networkx as nx
import copy
from warnings import warn
from defdap.utils import Datastore
from defdap.file_readers import EBSDDataLoader
from defdap.file_writers import EBSDDataWriter
from defdap.quat import Quat
from defdap import base
from defdap._accelerated import flood_fill
from defdap import defaults
from defdap.plotting import MapPlot
from defdap.utils import report_progress
[docs]class Map(base.Map):
"""
Class to encapsulate an EBSD map and useful analysis and plotting
methods.
Attributes
----------
step_size : float
Step size in micron.
phases : list of defdap.crystal.Phase
List of phases.
mis_ori : numpy.ndarray
Map of misorientation.
mis_ori_axis : list of numpy.ndarray
Map of misorientation axis components.
origin : tuple(int)
Map origin (x, y). Used by linker class where origin is a
homologue point of the maps.
data : defdap.utils.Datastore
Must contain after loading data (maps):
phase : numpy.ndarray
1-based, 0 is non-indexed points
euler_angle : numpy.ndarray
stored as (3, y_dim, x_dim) in radians
Generated data:
orientation : numpy.ndarray of defdap.quat.Quat
Quaterion for each point of map. Shape (y_dim, x_dim).
grain_boundaries : BoundarySet
phase_boundaries : BoundarySet
grains : numpy.ndarray of int
Map of grains. Grain numbers start at 1 here but everywhere else
grainID starts at 0. Regions that are smaller than the minimum
grain size are given value -2. Remnant boundary points are -1.
KAM : numpy.ndarray
Kernal average misorientaion map.
GND : numpy.ndarray
GND scalar map.
Nye_tensor : numpy.ndarray
3x3 Nye tensor at each point.
Derived data:
Grain list data to map data from all grains
"""
MAPNAME = 'ebsd'
def __init__(self, *args, **kwargs):
"""
Initialise class and load EBSD data.
Parameters
----------
*args, **kwarg
Passed to base constructor
"""
# Initialise variables
self.step_size = None
self.phases = []
# Call base class constructor
super(Map, self).__init__(*args, **kwargs)
self.mis_ori = None
self.mis_ori_axis = None
self.origin = (0, 0)
# Phase used for the maps crystal structure and c_over_a. So old
# functions still work for the 'main' phase in the map. 0-based
self.primary_phase_id = 0
# Use euler map for defining homologous points
self.plot_default = self.plot_euler_map
self.homog_map_name = 'band_contrast'
self.highlight_alpha = 1
self.data.add_generator(
'orientation', self.calc_quat_array, unit='', type='map',
order=0, default_component='IPF_x',
)
self.data.add_generator(
('phase_boundaries', 'grain_boundaries'), self.find_boundaries,
type='boundaries',
)
self.data.add_generator(
'grains', self.find_grains, unit='', type='map', order=0
)
self.data.add_generator(
'KAM', self.calc_kam, unit='rad', type='map', order=0,
plot_params={
'plot_colour_bar': True,
'clabel': 'KAM',
}
)
self.data.add_generator(
('GND', 'Nye_tensor'), self.calc_nye,
unit='', type='map',
metadatas=({
'order': 0,
'plot_params': {
'plot_colour_bar': True,
'clabel': 'GND content',
}
}, {
'order': 2,
'save': False,
'default_component': (0, 0),
'plot_params': {
'plot_colour_bar': True,
'clabel': 'Nye tensor',
}
})
)
[docs] @report_progress("loading EBSD data")
def load_data(self, file_name, data_type=None):
"""Load in EBSD data from file.
Parameters
----------
file_name : pathlib.Path
Path to EBSD file
data_type : str, {'OxfordBinary', 'OxfordText', 'EdaxAng', 'PythonDict'}
Format of EBSD data file.
"""
data_loader = EBSDDataLoader.get_loader(data_type, file_name)
data_loader.load(file_name)
metadata_dict = data_loader.loaded_metadata
self.shape = metadata_dict['shape']
self.step_size = metadata_dict['step_size']
self.phases = metadata_dict['phases']
self.data.update(data_loader.loaded_data)
# write final status
yield (f"Loaded EBSD data (dimensions: {self.x_dim} x {self.y_dim} "
f"pixels, step size: {self.step_size} um)")
[docs] def save(self, file_name, data_type=None, file_dir=""):
"""Save EBSD map to file.
Parameters
----------
file_name : str
Name of file to save to, it must not already exist.
data_type : str, {'OxfordText'}
Format of EBSD data file to save.
file_dir : str
Directory to save the file to.
"""
data_writer = EBSDDataWriter.get_writer(data_type)
data_writer.metadata['shape'] = self.shape
data_writer.metadata['step_size'] = self.step_size
data_writer.metadata['phases'] = self.phases
data_writer.data['phase'] = self.data.phase
data_writer.data['quat'] = self.data.orientation
data_writer.data['band_contrast'] = self.data.band_contrast
data_writer.write(file_name, file_dir=file_dir)
@property
def crystal_sym(self):
"""Crystal symmetry of the primary phase.
Returns
-------
str
Crystal symmetry
"""
return self.primary_phase.crystal_structure.name
@property
def c_over_a(self):
"""C over A ratio of the primary phase
Returns
-------
float or None
C over A ratio if hexagonal crystal structure otherwise None
"""
return self.primary_phase.c_over_a
@property
def num_phases(self):
return len(self.phases) or None
@property
def primary_phase(self):
"""Primary phase of the EBSD map.
Returns
-------
defdap.crystal.Phase
Primary phase
"""
return self.phases[self.primary_phase_id]
@property
def scale(self):
return self.step_size
[docs] @report_progress("rotating EBSD data")
def rotate_data(self):
"""Rotate map by 180 degrees and transform quats accordingly.
"""
self.data.euler_angle = self.data.euler_angle[:, ::-1, ::-1]
self.data.band_contrast = self.data.band_contrast[::-1, ::-1]
self.data.band_slope = self.data.band_slope[::-1, ::-1]
self.data.phase = self.data.phase[::-1, ::-1]
self.calc_quat_array()
# Rotation from old coord system to new
transform_quat = Quat.from_axis_angle(np.array([0, 0, 1]), np.pi).conjugate
# Perform vectorised multiplication
quats = Quat.multiply_many_quats(self.data.orientation.flatten(), transform_quat)
self.data.orientation = np.array(quats).reshape(self.shape)
yield 1.
[docs] def calc_euler_colour(self, map_data, phases=None, bg_colour=None):
if phases is None:
phases = self.phases
phase_ids = range(len(phases))
else:
phase_ids = phases
phases = [self.phases[i] for i in phase_ids]
if bg_colour is None:
bg_colour = np.array([0., 0., 0.])
map_colours = np.tile(bg_colour, self.shape + (1,))
for phase, phase_id in zip(phases, phase_ids):
if phase.crystal_structure.name == 'cubic':
norm = np.array([2 * np.pi, np.pi / 2, np.pi / 2])
elif phase.crystal_structure.name == 'hexagonal':
norm = np.array([np.pi, np.pi, np.pi / 3])
else:
ValueError("Only hexagonal and cubic symGroup supported")
# Apply normalisation for each phase
phase_mask = self.data.phase == phase_id + 1
map_colours[phase_mask] = map_data[:, phase_mask].T / norm
return map_colours
[docs] def calc_ipf_colour(self, map_data, direction, phases=None,
bg_colour=None):
if phases is None:
phases = self.phases
phase_ids = range(len(phases))
else:
phase_ids = phases
phases = [self.phases[i] for i in phase_ids]
if bg_colour is None:
bg_colour = np.array([0., 0., 0.])
map_colours = np.tile(bg_colour, self.shape + (1,))
for phase, phase_id in zip(phases, phase_ids):
# calculate IPF colours for phase
phase_mask = self.data.phase == phase_id + 1
map_colours[phase_mask] = Quat.calc_ipf_colours(
map_data[phase_mask], direction, phase.crystal_structure.name
).T
return map_colours
[docs] def plot_euler_map(self, phases=None, bg_colour=None, **kwargs):
"""Plot an orientation map in Euler colouring
Parameters
----------
phases : list of int
Which phases to plot for
kwargs
All arguments are passed to :func:`defdap.plotting.MapPlot.create`.
Returns
-------
defdap.plotting.MapPlot
"""
# Set default plot parameters then update with any input
plot_params = {}
plot_params.update(kwargs)
map_colours = self.calc_euler_colour(
self.data.euler_angle, phases=phases, bg_colour=bg_colour
)
return MapPlot.create(self, map_colours, **plot_params)
[docs] def plot_ipf_map(self, direction, phases=None, bg_colour=None, **kwargs):
"""
Plot a map with points coloured in IPF colouring,
with respect to a given sample direction.
Parameters
----------
direction : np.array len 3
Sample direction.
phases : list of int
Which phases to plot IPF data for.
bg_colour : np.array len 3
Colour of background (i.e. for phases not plotted).
kwargs
Other arguments passed to :func:`defdap.plotting.MapPlot.create`.
Returns
-------
defdap.plotting.MapPlot
"""
# Set default plot parameters then update with any input
plot_params = {}
plot_params.update(kwargs)
map_colours = self.calc_ipf_colour(
self.data.orientation, direction, phases=phases,
bg_colour=bg_colour
)
return MapPlot.create(self, map_colours, **plot_params)
[docs] def plot_phase_map(self, **kwargs):
"""Plot a phase map.
Parameters
----------
kwargs
All arguments passed to :func:`defdap.plotting.MapPlot.create`.
Returns
-------
defdap.plotting.MapPlot
"""
# Set default plot parameters then update with any input
plot_params = {
'vmin': 0,
'vmax': self.num_phases
}
plot_params.update(kwargs)
plot = MapPlot.create(self, self.data.phase, **plot_params)
# add a legend to the plot
phase_ids = list(range(0, self.num_phases + 1))
phase_names = ["Non-indexed"] + [phase.name for phase in self.phases]
plot.add_legend(phase_ids, phase_names, loc=2, borderaxespad=0.)
return plot
[docs] @report_progress("calculating KAM")
def calc_kam(self):
"""
Calculates Kernel Average Misorientaion (KAM) for the EBSD map,
based on a 3x3 kernel. Crystal symmetric equivalences are not
considered. Stores result as `KAM`.
"""
quat_comps = np.empty((4, ) + self.shape)
for i, row in enumerate(self.data.orientation):
for j, quat in enumerate(row):
quat_comps[:, i, j] = quat.quat_coef
kam = np.empty(self.shape)
# Start with rows. Calculate misorientation with neighbouring rows.
# First and last row only in one direction
kam[0] = abs(np.einsum("ij,ij->j",
quat_comps[:, 0], quat_comps[:, 1]))
kam[-1] = abs(np.einsum("ij,ij->j",
quat_comps[:, -1], quat_comps[:, -2]))
for i in range(1, self.y_dim - 1):
kam[i] = (abs(np.einsum("ij,ij->j",
quat_comps[:, i], quat_comps[:, i + 1])) +
abs(np.einsum("ij,ij->j",
quat_comps[:, i], quat_comps[:, i - 1]))
) / 2
kam[kam > 1] = 1
# Do the same for columns
kam[:, 0] += abs(np.einsum("ij,ij->j",
quat_comps[:, :, 0], quat_comps[:, :, 1]))
kam[:, -1] += abs(np.einsum("ij,ij->j",
quat_comps[:, :, -1], quat_comps[:, :, -2]))
for i in range(1, self.x_dim - 1):
kam[:, i] += (abs(np.einsum("ij,ij->j",
quat_comps[:, :, i],
quat_comps[:, :, i + 1])) +
abs(np.einsum("ij,ij->j",
quat_comps[:, :, i],
quat_comps[:, :, i - 1]))
) / 2
kam /= 2
kam[kam > 1] = 1
yield 1.
return 2 * np.arccos(kam)
[docs] @report_progress("calculating Nye tensor")
def calc_nye(self):
"""
Calculates Nye tensor and related GND density for the EBSD map.
Stores result as `Nye_tensor` and `GND`. Uses the crystal
symmetry of the primary phase.
"""
syms = self.primary_phase.crystal_structure.symmetries
num_syms = len(syms)
# array to store quat components of initial and symmetric equivalents
quat_comps = np.empty((num_syms, 4) + self.shape)
# populate with initial quat components
for i, row in enumerate(self.data.orientation):
for j, quat in enumerate(row):
quat_comps[0, :, i, j] = quat.quat_coef
# loop of over symmetries and apply to initial quat components
# (excluding first symmetry as this is the identity transformation)
for i, sym in enumerate(syms[1:], start=1):
# sym[i] * quat for all points (* is quaternion product)
quat_comps[i, 0] = (quat_comps[0, 0] * sym[0] - quat_comps[0, 1] * sym[1] -
quat_comps[0, 2] * sym[2] - quat_comps[0, 3] * sym[3])
quat_comps[i, 1] = (quat_comps[0, 0] * sym[1] + quat_comps[0, 1] * sym[0] -
quat_comps[0, 2] * sym[3] + quat_comps[0, 3] * sym[2])
quat_comps[i, 2] = (quat_comps[0, 0] * sym[2] + quat_comps[0, 2] * sym[0] -
quat_comps[0, 3] * sym[1] + quat_comps[0, 1] * sym[3])
quat_comps[i, 3] = (quat_comps[0, 0] * sym[3] + quat_comps[0, 3] * sym[0] -
quat_comps[0, 1] * sym[2] + quat_comps[0, 2] * sym[1])
# swap into positive hemisphere if required
quat_comps[i, :, quat_comps[i, 0] < 0] *= -1
# Arrays to store neighbour misorientation in positive x and y direction
mis_ori_x = np.zeros((num_syms,) + self.shape)
mis_ori_y = np.zeros((num_syms, ) + self.shape)
# loop over symmetries calculating misorientation to initial
for i in range(num_syms):
for j in range(self.x_dim - 1):
mis_ori_x[i, :, j] = abs(np.einsum("ij,ij->j", quat_comps[0, :, :, j], quat_comps[i, :, :, j + 1]))
for j in range(self.y_dim - 1):
mis_ori_y[i, j, :] = abs(np.einsum("ij,ij->j", quat_comps[0, :, j, :], quat_comps[i, :, j + 1, :]))
mis_ori_x[mis_ori_x > 1] = 1
mis_ori_y[mis_ori_y > 1] = 1
# find min misorientation (max here as misorientaion is cos of this)
arg_mis_ori_x = np.argmax(mis_ori_x, axis=0)
arg_mis_ori_y = np.argmax(mis_ori_y, axis=0)
mis_ori_x = np.max(mis_ori_x, axis=0)
mis_ori_y = np.max(mis_ori_y, axis=0)
# convert to misorientation in degrees
mis_ori_x = 360 * np.arccos(mis_ori_x) / np.pi
mis_ori_y = 360 * np.arccos(mis_ori_y) / np.pi
# calculate relative elastic distortion tensors at each point in the two directions
betaderx = np.zeros((3, 3) + self.shape)
betadery = betaderx
for i in range(self.x_dim - 1):
for j in range(self.y_dim - 1):
q0x = Quat(quat_comps[0, 0, j, i], quat_comps[0, 1, j, i],
quat_comps[0, 2, j, i], quat_comps[0, 3, j, i])
qix = Quat(quat_comps[arg_mis_ori_x[j, i], 0, j, i + 1],
quat_comps[arg_mis_ori_x[j, i], 1, j, i + 1],
quat_comps[arg_mis_ori_x[j, i], 2, j, i + 1],
quat_comps[arg_mis_ori_x[j, i], 3, j, i + 1])
misoquatx = qix.conjugate * q0x
# change stepsize to meters
betaderx[:, :, j, i] = (Quat.rot_matrix(misoquatx) - np.eye(3)) / self.step_size / 1e-6
q0y = Quat(quat_comps[0, 0, j, i], quat_comps[0, 1, j, i],
quat_comps[0, 2, j, i], quat_comps[0, 3, j, i])
qiy = Quat(quat_comps[arg_mis_ori_y[j, i], 0, j + 1, i],
quat_comps[arg_mis_ori_y[j, i], 1, j + 1, i],
quat_comps[arg_mis_ori_y[j, i], 2, j + 1, i],
quat_comps[arg_mis_ori_y[j, i], 3, j + 1, i])
misoquaty = qiy.conjugate * q0y
# change stepsize to meters
betadery[:, :, j, i] = (Quat.rot_matrix(misoquaty) - np.eye(3)) / self.step_size / 1e-6
# Calculate the Nye Tensor
alpha = np.empty((3, 3) + self.shape)
bavg = 1.4e-10 # Burgers vector
alpha[0, 2] = (betadery[0, 0] - betaderx[0, 1]) / bavg
alpha[1, 2] = (betadery[1, 0] - betaderx[1, 1]) / bavg
alpha[2, 2] = (betadery[2, 0] - betaderx[2, 1]) / bavg
alpha[:, 1] = betaderx[:, 2] / bavg
alpha[:, 0] = -1 * betadery[:, 2] / bavg
# Calculate 3 possible L1 norms of Nye tensor for total
# disloction density
alpha_total3 = np.empty(self.shape)
alpha_total5 = np.empty(self.shape)
alpha_total9 = np.empty(self.shape)
alpha_total3 = 30 / 10. * (
abs(alpha[0, 2]) + abs(alpha[1, 2]) + abs(alpha[2, 2])
)
alpha_total5 = 30 / 14. * (
abs(alpha[0, 2]) + abs(alpha[1, 2]) + abs(alpha[2, 2]) +
abs(alpha[1, 0]) + abs(alpha[0, 1])
)
alpha_total9 = 30 / 20. * (
abs(alpha[0, 2]) + abs(alpha[1, 2]) + abs(alpha[2, 2]) +
abs(alpha[0, 0]) + abs(alpha[1, 0]) + abs(alpha[2, 0]) +
abs(alpha[0, 1]) + abs(alpha[1, 1]) + abs(alpha[2, 1])
)
alpha_total3[abs(alpha_total3) < 1] = 1e12
alpha_total5[abs(alpha_total3) < 1] = 1e12
alpha_total9[abs(alpha_total3) < 1] = 1e12
# choose from the different alpha_totals according to preference;
# see Ruggles GND density paper
yield 1.
return alpha_total9, alpha
[docs] @report_progress("building quaternion array")
def calc_quat_array(self):
"""Build quaternion array
"""
# create the array of quat objects
quats = Quat.create_many_quats(self.data.euler_angle)
yield 1.
return quats
[docs] def filter_data(self, misori_tol=5):
# Kuwahara filter
print("8 quadrants")
misori_tol *= np.pi / 180
misori_tol = np.cos(misori_tol / 2)
# store quat components in array
quat_comps = np.empty((4,) + self.shape)
for idx in np.ndindex(self.shape):
quat_comps[(slice(None),) + idx] = self.data.orientation[idx].quat_coef
# misorientation in each quadrant surrounding a point
mis_oris = np.zeros((8,) + self.shape)
for i in range(2, self.shape[0] - 2):
for j in range(2, self.shape[1] - 2):
ref_quat = quat_comps[:, i, j]
quadrants = [
quat_comps[:, i - 2:i + 1, j - 2:j + 1], # UL
quat_comps[:, i - 2:i + 1, j - 1:j + 2], # UC
quat_comps[:, i - 2:i + 1, j:j + 3], # UR
quat_comps[:, i - 1:i + 2, j:j + 3], # MR
quat_comps[:, i:i + 3, j:j + 3], # LR
quat_comps[:, i:i + 3, j - 1:j + 2], # LC
quat_comps[:, i:i + 3, j - 2:j + 1], # LL
quat_comps[:, i - 1:i + 2, j - 2:j + 1] # ML
]
for k, quats in enumerate(quadrants):
mis_oris_quad = np.abs(
np.einsum("ijk,i->jk", quats, ref_quat)
)
mis_oris_quad = mis_oris_quad[mis_oris_quad > misori_tol]
mis_oris[k, i, j] = mis_oris_quad.mean()
min_mis_ori_quadrant = np.argmax(mis_oris, axis=0)
# minMisOris = np.max(mis_oris, axis=0)
# minMisOris[minMisOris > 1.] = 1.
# minMisOris = 2 * np.arccos(minMisOris)
quat_comps_new = np.copy(quat_comps)
for i in range(2, self.shape[0] - 2):
for j in range(2, self.shape[1] - 2):
# if minMisOris[i, j] < misOriTol:
# continue
ref_quat = quat_comps[:, i, j]
quadrants = [
quat_comps[:, i - 2:i + 1, j - 2:j + 1], # UL
quat_comps[:, i - 2:i + 1, j - 1:j + 2], # UC
quat_comps[:, i - 2:i + 1, j:j + 3], # UR
quat_comps[:, i - 1:i + 2, j:j + 3], # MR
quat_comps[:, i:i + 3, j:j + 3], # LR
quat_comps[:, i:i + 3, j - 1:j + 2], # LC
quat_comps[:, i:i + 3, j - 2:j + 1], # LL
quat_comps[:, i - 1:i + 2, j - 2:j + 1] # ML
]
quats = quadrants[min_mis_ori_quadrant[i, j]]
mis_oris_quad = np.abs(
np.einsum("ijk,i->jk", quats, ref_quat)
)
quats = quats[:, mis_oris_quad > misori_tol]
avOri = np.einsum("ij->i", quats)
# avOri /= np.sqrt(np.dot(avOri, avOri))
quat_comps_new[:, i, j] = avOri
quat_comps_new /= np.sqrt(np.einsum("ijk,ijk->jk", quat_comps_new, quat_comps_new))
quat_array_new = np.empty(self.shape, dtype=Quat)
for idx in np.ndindex(self.shape):
quat_array_new[idx] = Quat(quat_comps_new[(slice(None),) + idx])
self.data.orientation = quat_array_new
return quats
[docs] @report_progress("finding grain boundaries")
def find_boundaries(self, misori_tol=10):
"""Find grain and phase boundaries
Parameters
----------
misori_tol : float
Critical misorientation in degrees.
"""
# TODO: what happens with non-indexed points
# TODO: grain boundaries should be calculated per crystal structure
misori_tol *= np.pi / 180
syms = self.primary_phase.crystal_structure.symmetries
num_syms = len(syms)
# array to store quat components of initial and symmetric equivalents
quat_comps = np.empty((num_syms, 4) + self.shape)
# populate with initial quat components
for i, row in enumerate(self.data.orientation):
for j, quat in enumerate(row):
quat_comps[0, :, i, j] = quat.quat_coef
# loop of over symmetries and apply to initial quat components
# (excluding first symmetry as this is the identity transformation)
for i, sym in enumerate(syms[1:], start=1):
# sym[i] * quat for all points (* is quaternion product)
quat_comps[i, 0] = (
quat_comps[0, 0]*sym[0] - quat_comps[0, 1]*sym[1] -
quat_comps[0, 2]*sym[2] - quat_comps[0, 3]*sym[3]
)
quat_comps[i, 1] = (
quat_comps[0, 0]*sym[1] + quat_comps[0, 1]*sym[0] -
quat_comps[0, 2]*sym[3] + quat_comps[0, 3]*sym[2]
)
quat_comps[i, 2] = (
quat_comps[0, 0]*sym[2] + quat_comps[0, 2]*sym[0] -
quat_comps[0, 3]*sym[1] + quat_comps[0, 1]*sym[3]
)
quat_comps[i, 3] = (
quat_comps[0, 0]*sym[3] + quat_comps[0, 3]*sym[0] -
quat_comps[0, 1]*sym[2] + quat_comps[0, 2]*sym[1]
)
# swap into positive hemisphere if required
quat_comps[i, :, quat_comps[i, 0] < 0] *= -1
# Arrays to store neighbour misorientation in positive x and y
# directions
misori_x = np.ones((num_syms, ) + self.shape)
misori_y = np.ones((num_syms, ) + self.shape)
# loop over symmetries calculating misorientation to initial
for i in range(num_syms):
for j in range(self.shape[1] - 1):
misori_x[i, :, j] = abs(np.einsum(
"ij,ij->j", quat_comps[0, :, :, j], quat_comps[i, :, :, j+1]
))
for j in range(self.shape[0] - 1):
misori_y[i, j] = abs(np.einsum(
"ij,ij->j", quat_comps[0, :, j], quat_comps[i, :, j+1]
))
misori_x[misori_x > 1] = 1
misori_y[misori_y > 1] = 1
# find max dot product and then convert to misorientation angle
misori_x = 2 * np.arccos(np.max(misori_x, axis=0))
misori_y = 2 * np.arccos(np.max(misori_y, axis=0))
# PHASE boundary POINTS
phase_im = self.data.phase
pb_im_x = np.not_equal(phase_im, np.roll(phase_im, -1, axis=1))
pb_im_x[:, -1] = False
pb_im_y = np.not_equal(phase_im, np.roll(phase_im, -1, axis=0))
pb_im_y[-1] = False
phase_boundaries = BoundarySet.from_image(self, pb_im_x, pb_im_y)
grain_boundaries = BoundarySet.from_image(
self,
(misori_x > misori_tol) | pb_im_x,
(misori_y > misori_tol) | pb_im_y
)
yield 1.
return phase_boundaries, grain_boundaries
[docs] @report_progress("constructing neighbour network")
def build_neighbour_network(self):
# create network
nn = nx.Graph()
nn.add_nodes_from(self.grains)
points_x = self.data.grain_boundaries.points_x
points_y = self.data.grain_boundaries.points_y
total_points_x = len(points_x)
total_points = total_points_x + len(points_y)
for i, points in enumerate((points_x, points_y)):
for i_point, (x, y) in enumerate(points):
# report progress
yield (i_point + i * total_points_x) / total_points
if (x == 0 or y == 0 or x == self.shape[1] - 1 or
y == self.shape[0] - 1):
# exclude boundary pixels of map
continue
grain_id = self.data.grains[y, x] - 1
nei_grain_id = self.data.grains[y + i, x - i + 1] - 1
if nei_grain_id == grain_id:
# ignore if neighbour is same as grain
continue
if nei_grain_id < 0 or grain_id < 0:
# ignore if not a grain (boundary points -1 and
# points in small grains -2)
continue
grain = self[grain_id]
nei_grain = self[nei_grain_id]
try:
# look up boundary segment if it exists
b_seg = nn[grain][nei_grain]['boundary']
except KeyError:
# neighbour relation doesn't exist so add it
b_seg = BoundarySegment(self, grain, nei_grain)
nn.add_edge(grain, nei_grain, boundary=b_seg)
# add the boundary point
b_seg.addBoundaryPoint((x, y), i, grain)
self.neighbour_network = nn
[docs] def plot_phase_boundary_map(self, dilate=False, **kwargs):
"""Plot phase boundary map.
Parameters
----------
dilate : bool
If true, dilate boundary.
kwargs
All other arguments are passed to :func:`defdap.plotting.MapPlot.create`.
Returns
-------
defdap.plotting.MapPlot
"""
# Set default plot parameters then update with any input
plot_params = {
'vmax': 1,
'plot_colour_bar': True,
'cmap': 'gray'
}
plot_params.update(kwargs)
boundaries_image = self.data.phase_boundaries.image.astype(int)
if dilate:
boundaries_image = mph.binary_dilation(boundaries_image)
plot = MapPlot.create(self, boundaries_image, **plot_params)
return plot
[docs] def plot_boundary_map(self, **kwargs):
"""Plot grain boundary map.
Parameters
----------
kwargs
All arguments are passed to :func:`defdap.plotting.MapPlot.create`.
Returns
-------
defdap.plotting.MapPlot
"""
# Set default plot parameters then update with any input
plot_params = {
'plot_gbs': True,
'boundaryColour': 'black'
}
plot_params.update(kwargs)
plot = MapPlot.create(self, None, **plot_params)
return plot
[docs] @report_progress("finding grains")
def find_grains(self, min_grain_size=10):
"""Find grains and assign IDs.
Parameters
----------
min_grain_size : int
Minimum grain area in pixels.
"""
# Initialise the grain map
# TODO: Look at grain map compared to boundary map
grains = np.zeros(self.shape, dtype=int)
grain_list = []
boundary_im_x = self.data.grain_boundaries.image_x
boundary_im_y = self.data.grain_boundaries.image_y
# List of points where no grain has be set yet
points_left = self.data.phase != 0
total_points = points_left.sum()
found_point = 0
next_point = points_left.tobytes().find(b'\x01')
# Start counter for grains
grain_index = 1
group_id = Datastore.generate_id()
# Loop until all points (except boundaries) have been assigned
# to a grain or ignored
i = 0
coords_buffer = np.zeros((boundary_im_y.size, 2), dtype=np.intp)
while found_point >= 0:
# Flood fill first unknown point and return grain object
seed = np.unravel_index(next_point, self.shape)
grain = Grain(grain_index - 1, self, group_id)
grain.data.point = flood_fill(
(seed[1], seed[0]), grain_index, points_left, grains,
boundary_im_x, boundary_im_y, coords_buffer
)
coords_buffer = coords_buffer[len(grain.data.point):]
if len(grain) < min_grain_size:
# if grain size less than minimum, ignore grain and set
# values in grain map to -2
for point in grain.data.point:
grains[point[1], point[0]] = -2
else:
# add grain to list and increment grain index
grain_list.append(grain)
grain_index += 1
# find next search point
points_left_sub = points_left.reshape(-1)[next_point + 1:]
found_point = points_left_sub.tobytes().find(b'\x01')
next_point += found_point + 1
# report progress
i += 1
if i == defaults['find_grain_report_freq']:
yield 1. - points_left_sub.sum() / total_points
i = 0
# Assign phase to each grain
for grain in grain_list:
phase_vals = grain.grain_data(self.data.phase)
if np.max(phase_vals) != np.min(phase_vals):
warn(f"Grain {grain.grain_id} could not be assigned a "
f"phase, phase vals not constant.")
continue
phase_id = phase_vals[0] - 1
if not (0 <= phase_id < self.num_phases):
warn(f"Grain {grain.grain_id} could not be assigned a "
f"phase, invalid phase {phase_id}.")
continue
grain.phase_id = phase_id
grain.phase = self.phases[phase_id]
## TODO: this will get duplicated if find grains called again
self.data.add_derivative(
grain_list[0].data, self.grain_data_to_map, pass_ref=True,
in_props={
'type': 'list'
},
out_props={
'type': 'map'
}
)
self._grains = grain_list
return grains
[docs] def plot_grain_map(self, **kwargs):
"""Plot a map with grains coloured.
Parameters
----------
kwargs
All arguments are passed to :func:`defdap.plotting.MapPlot.create`.
Returns
-------
defdap.plotting.MapPlot
"""
# Set default plot parameters then update with any input
plot_params = {
'clabel': "Grain number"
}
plot_params.update(kwargs)
plot = MapPlot.create(self, self.data.grains, **plot_params)
return plot
[docs] @report_progress("calculating grain mean orientations")
def calc_grain_av_oris(self):
"""Calculate the average orientation of grains.
"""
numGrains = len(self)
for iGrain, grain in enumerate(self):
grain.calc_average_ori()
# report progress
yield (iGrain + 1) / numGrains
[docs] @report_progress("calculating grain misorientations")
def calc_grain_mis_ori(self, calc_axis=False):
"""Calculate the misorientation within grains.
Parameters
----------
calc_axis : bool
Calculate the misorientation axis if True.
"""
num_grains = len(self)
for i_grain, grain in enumerate(self):
grain.build_mis_ori_list(calc_axis=calc_axis)
# report progress
yield (i_grain + 1) / num_grains
[docs] def plot_mis_ori_map(self, component=0, **kwargs):
"""Plot misorientation map.
Parameters
----------
component : int, {0, 1, 2, 3}
0 gives misorientation, 1, 2, 3 gives rotation about x, y, z
kwargs
All other arguments are passed to :func:`defdap.plotting.MapPlot.create`.
Returns
-------
defdap.plotting.MapPlot
"""
if component in [1, 2, 3]:
self.mis_ori = np.zeros(self.shape)
# Calculate misorientation axis if not calculated
if np.any([grain.mis_ori_axis_list is None for grain in self]):
self.calc_grain_mis_ori(calc_axis=True)
for grain in self:
for point, mis_ori_axis in zip(grain.data.point, np.array(grain.mis_ori_axis_list)):
self.mis_ori[point[1], point[0]] = mis_ori_axis[component - 1]
mis_ori = self.mis_ori * 180 / np.pi
clabel = r"Rotation around {:} axis ($^\circ$)".format(
['X', 'Y', 'Z'][component-1]
)
else:
self.mis_ori = np.ones(self.shape)
# Calculate misorientation if not calculated
if np.any([grain.mis_ori_list is None for grain in self]):
self.calc_grain_mis_ori(calc_axis=False)
for grain in self:
for point, mis_ori in zip(grain.data.point, grain.mis_ori_list):
self.mis_ori[point[1], point[0]] = mis_ori
mis_ori = np.arccos(self.mis_ori) * 360 / np.pi
clabel = r"Grain reference orienation deviation (GROD) ($^\circ$)"
# Set default plot parameters then update with any input
plot_params = {
'plot_colour_bar': True,
'clabel': clabel
}
plot_params.update(kwargs)
plot = MapPlot.create(self, mis_ori, **plot_params)
return plot
[docs] @report_progress("calculating grain average Schmid factors")
def calc_average_grain_schmid_factors(self, load_vector, slip_systems=None):
"""
Calculates Schmid factors for all slip systems, for all grains,
based on average grain orientation.
Parameters
----------
load_vector :
Loading vector, e.g. [1, 0, 0].
slip_systems : list, optional
Slip planes to calculate Schmid factor for, maximum of all
planes calculated if not given.
"""
num_grains = len(self)
for iGrain, grain in enumerate(self):
grain.calc_average_schmid_factors(load_vector, slip_systems=slip_systems)
# report progress
yield (iGrain + 1) / num_grains
[docs] @report_progress("calculating RDR values")
def calc_rdr(self):
"""Calculates Relative Displacent Ratio values for all grains"""
num_grains = len(self)
for iGrain, grain in enumerate(self):
grain.calc_rdr()
# report progress
yield (iGrain + 1) / num_grains
[docs] def plot_average_grain_schmid_factors_map(self, planes=None, directions=None,
**kwargs):
"""
Plot maximum Schmid factor map, based on average grain
orientation (for all slip systems unless specified).
Parameters
----------
planes : list, optional
Plane ID(s) to consider. All planes considered if not given.
directions : list, optional
Direction ID(s) to consider. All directions considered if not given.
kwargs
All other arguments are passed to :func:`defdap.plotting.MapPlot.create`.
Returns
-------
defdap.plotting.MapPlot
"""
# Set default plot parameters then update with any input
plot_params = {
'vmin': 0,
'vmax': 0.5,
'cmap': 'gray',
'plot_colour_bar': True,
'clabel': "Schmid factor"
}
plot_params.update(kwargs)
if self[0].average_schmid_factors is None:
raise Exception("Run 'calc_average_grain_schmid_factors' first")
grains_sf_max = []
for grain in self:
current_sf = []
if planes is not None:
for plane in planes:
if directions is not None:
for direction in directions:
current_sf.append(
grain.average_schmid_factors[plane][direction]
)
else:
current_sf += grain.average_schmid_factors[plane]
else:
for sf_group in grain.average_schmid_factors:
current_sf += sf_group
grains_sf_max.append(max(current_sf))
plot = self.plot_grain_data_map(grain_data=grains_sf_max, bg=0.5,
**plot_params)
return plot
[docs]class Grain(base.Grain):
"""
Class to encapsulate a grain in an EBSD map and useful analysis and
plotting methods.
Attributes
----------
ebsd_map : defdap.ebsd.Map
EBSD map this grain is a member of.
owner_map : defdap.ebsd.Map
EBSD map this grain is a member of.
phase_id : int
phase : defdap.crystal.Phase
data : defdap.utils.Datastore
Must contain after creating:
point : list of tuples
(x, y)
Generated data:
GROD : numpy.ndarray
GROD_axis : numpy.ndarray
Derived data:
Map data to list data from the map the grain is part of
mis_ori_list : list
MisOri at each point in grain.
mis_ori_axis_list : list
MisOri axes at each point in grain.
ref_ori : defdap.quat.Quat
Average ori of grain
average_mis_ori
Average mis_ori of grain.
average_schmid_factors : list
List of list Schmid factors (grouped by slip plane).
slip_trace_angles : list
Slip trace angles in screen plane.
slip_trace_inclinations : list
Angle between slip plane and screen plane.
"""
def __init__(self, grain_id, ebsdMap, group_id):
# Call base class constructor
super(Grain, self).__init__(grain_id, ebsdMap, group_id)
self.ebsd_map = self.owner_map # ebsd map this grain is a member of
self.mis_ori_list = None # list of mis_ori at each point in grain
self.mis_ori_axis_list = None # list of mis_ori axes at each point in grain
self.ref_ori = None # (quat) average ori of grain
self.average_mis_ori = None # average mis_ori of grain
self.average_schmid_factors = None # list of list Schmid factors (grouped by slip plane)
self.slip_trace_angles = None # list of slip trace angles
self.slip_trace_inclinations = None
self.plot_default = self.plot_unit_cell
self.data.add_generator(
('GROD', 'GROD_axis'), self.calc_grod,
type='list',
metadatas=({
'unit': 'rad',
'order': 0,
'plot_params': {
'plot_colour_bar': True,
'clabel': 'GROD',
}
}, {
'unit': '',
'order': 1,
'default_component': 0,
'plot_params': {
'plot_colour_bar': True,
'clabel': 'GROD axis',
}
})
)
@property
def crystal_sym(self):
"""Temporary"""
return self.phase.crystal_structure.name
[docs] def calc_average_ori(self):
"""Calculate the average orientation of a grain.
"""
quat_comps_sym = Quat.calc_sym_eqvs(self.data.orientation, self.crystal_sym)
self.ref_ori = Quat.calc_average_ori(quat_comps_sym)
[docs] def build_mis_ori_list(self, calc_axis=False):
"""Calculate the misorientation within given grain.
Parameters
----------
calc_axis : bool
Calculate the misorientation axis if True.
"""
quat_comps_sym = Quat.calc_sym_eqvs(self.data.orientation, self.crystal_sym)
if self.ref_ori is None:
self.ref_ori = Quat.calc_average_ori(quat_comps_sym)
mis_ori_array, min_quat_comps = Quat.calcMisOri(quat_comps_sym, self.ref_ori)
self.average_mis_ori = mis_ori_array.mean()
self.mis_ori_list = list(mis_ori_array)
if calc_axis:
# Now for axis calculation
ref_ori_inv = self.ref_ori.conjugate
mis_ori_axis = np.empty((3, min_quat_comps.shape[1]))
dq = np.empty((4, min_quat_comps.shape[1]))
# ref_ori_inv * minQuat for all points (* is quaternion product)
# change to minQuat * ref_ori_inv
dq[0, :] = (ref_ori_inv[0] * min_quat_comps[0, :] - ref_ori_inv[1] * min_quat_comps[1, :] -
ref_ori_inv[2] * min_quat_comps[2, :] - ref_ori_inv[3] * min_quat_comps[3, :])
dq[1, :] = (ref_ori_inv[1] * min_quat_comps[0, :] + ref_ori_inv[0] * min_quat_comps[1, :] +
ref_ori_inv[3] * min_quat_comps[2, :] - ref_ori_inv[2] * min_quat_comps[3, :])
dq[2, :] = (ref_ori_inv[2] * min_quat_comps[0, :] + ref_ori_inv[0] * min_quat_comps[2, :] +
ref_ori_inv[1] * min_quat_comps[3, :] - ref_ori_inv[3] * min_quat_comps[1, :])
dq[3, :] = (ref_ori_inv[3] * min_quat_comps[0, :] + ref_ori_inv[0] * min_quat_comps[3, :] +
ref_ori_inv[2] * min_quat_comps[1, :] - ref_ori_inv[1] * min_quat_comps[2, :])
dq[:, dq[0] < 0] = -dq[:, dq[0] < 0]
# numpy broadcasting taking care of different array sizes
mis_ori_axis[:, :] = (2 * dq[1:4, :] * np.arccos(dq[0, :])) / np.sqrt(1 - np.power(dq[0, :], 2))
# hack it back into a list. Need to change self.*List to be arrays, it was a bad decision to
# make them lists in the beginning
self.mis_ori_axis_list = []
for row in mis_ori_axis.transpose():
self.mis_ori_axis_list.append(row)
[docs] def calc_grod(self):
quat_comps = Quat.calc_sym_eqvs(self.data.orientation, self.crystal_sym)
if self.ref_ori is None:
self.ref_ori = Quat.calc_average_ori(quat_comps)
misori, quat_comps = Quat.calcMisOri(quat_comps, self.ref_ori)
misori = 2 * np.arccos(misori)
ref_ori_inv = self.ref_ori.conjugate
dq = np.empty((4, len(self)))
# ref_ori_inv * quat_comps for all points
# change to quat_comps * ref_ori_inv
dq[0] = (ref_ori_inv[0]*quat_comps[0] - ref_ori_inv[1]*quat_comps[1] -
ref_ori_inv[2]*quat_comps[2] - ref_ori_inv[3]*quat_comps[3])
dq[1] = (ref_ori_inv[1]*quat_comps[0] + ref_ori_inv[0]*quat_comps[1] +
ref_ori_inv[3]*quat_comps[2] - ref_ori_inv[2]*quat_comps[3])
dq[2] = (ref_ori_inv[2]*quat_comps[0] + ref_ori_inv[0]*quat_comps[2] +
ref_ori_inv[1]*quat_comps[3] - ref_ori_inv[3]*quat_comps[1])
dq[3] = (ref_ori_inv[3]*quat_comps[0] + ref_ori_inv[0]*quat_comps[3] +
ref_ori_inv[2]*quat_comps[1] - ref_ori_inv[1]*quat_comps[2])
dq[:, dq[0] < 0] *= -1
misori_axis = (2 * dq[1:4] * np.arccos(dq[0])) / np.sqrt(1 - dq[0]**2)
return misori, misori_axis
[docs] def plot_ref_ori(self, direction=np.array([0, 0, 1]), **kwargs):
"""Plot the average grain orientation on an IPF.
Parameters
----------
direction : numpy.ndarray
Sample direction for IPF.
kwargs
All other arguments are passed to :func:`defdap.quat.Quat.plot_ipf`.
Returns
-------
defdap.plotting.PolePlot
"""
plot_params = {'marker': '+'}
plot_params.update(kwargs)
return Quat.plot_ipf([self.ref_ori], direction, self.crystal_sym,
**plot_params)
[docs] def plot_ori_spread(self, direction=np.array([0, 0, 1]), **kwargs):
"""Plot all orientations within a given grain, on an IPF.
Parameters
----------
direction : numpy.ndarray
Sample direction for IPF.
kwargs
All other arguments are passed to :func:`defdap.quat.Quat.plot_ipf`.
Returns
-------
defdap.plotting.PolePlot
"""
plot_params = {'marker': '.'}
plot_params.update(kwargs)
return Quat.plot_ipf(self.data.orientation, direction, self.crystal_sym,
**plot_params)
[docs] def plot_unit_cell(self, fig=None, ax=None, plot=None, **kwargs):
"""Plot an unit cell of the average grain orientation.
Parameters
----------
fig : matplotlib.figure.Figure
Matplotlib figure to plot on
ax : matplotlib.figure.Figure
Matplotlib figure to plot on
plot : defdap.plotting.PolePlot
defdap plot to plot the figure to.
kwargs
All other arguments are passed to :func:`defdap.quat.Quat.plot_unit_cell`.
"""
crystal_structure = self.ebsd_map.phases[self.phase_id].crystal_structure
plot = Quat.plot_unit_cell(self.ref_ori, fig=fig, ax=ax, plot=plot,
crystal_structure=crystal_structure, **kwargs)
return plot
[docs] def plot_mis_ori(self, component=0, **kwargs):
"""Plot misorientation map for a given grain.
Parameters
----------
component : int, {0, 1, 2, 3}
0 gives misorientation, 1, 2, 3 gives rotation about x, y, z.
kwargs
All other arguments are passed to :func:`defdap.ebsd.plot_grain_data`.
Returns
-------
defdap.plotting.GrainPlot
"""
component = int(component)
# Set default plot parameters then update with any input
plot_params = {
'plot_colour_bar': True
}
if component == 0:
if self.mis_ori_list is None: self.build_mis_ori_list()
plot_params['clabel'] = r"Grain reference orientation " \
r"deviation (GROD) ($^\circ$)"
plot_data = np.rad2deg(2 * np.arccos(self.mis_ori_list))
elif 0 < component < 4:
if self.mis_ori_axis_list is None: self.build_mis_ori_list(calc_axis=True)
plot_params['clabel'] = r"Rotation around {:} ($^\circ$)".format(
['X', 'Y', 'Z'][component-1]
)
plot_data = np.rad2deg(np.array(self.mis_ori_axis_list)[:, component - 1])
else:
raise ValueError("Component must between 0 and 3")
plot_params.update(kwargs)
plot = self.plot_grain_data(grain_data=plot_data, **plot_params)
return plot
# define load axis as unit vector
[docs] def calc_average_schmid_factors(self, load_vector, slip_systems=None):
"""Calculate Schmid factors for grain, using average orientation.
Parameters
----------
load_vector : numpy.ndarray
Loading vector, i.e. [1, 0, 0]
slip_systems : list, optional
Slip planes to calculate Schmid factor for. Maximum for all planes
used if not set.
"""
if slip_systems is None:
slip_systems = self.phase.slip_systems
if self.ref_ori is None:
self.calc_average_ori()
# orientation of grain
grain_av_ori = self.ref_ori
# Transform the load vector into crystal coordinates
load_vector_crystal = grain_av_ori.transform_vector(load_vector)
self.average_schmid_factors = []
# flatten list of lists
# slip_systems = chain.from_iterable(slip_systems)
# Loop over groups of slip systems with same slip plane
for i, slip_system_group in enumerate(slip_systems):
self.average_schmid_factors.append([])
# Then loop over individual slip systems
for slip_system in slip_system_group:
schmidFactor = abs(np.dot(load_vector_crystal, slip_system.slip_plane) *
np.dot(load_vector_crystal, slip_system.slip_dir))
self.average_schmid_factors[i].append(schmidFactor)
return
[docs] def calc_rdr(self):
"""Calculate Relative Displacement Ratio values."""
self.rdr = []
# Loop over groups of slip systems with same slip plane
for i, slip_system_group in enumerate(self.phase.slip_systems):
self.rdr.append([])
# Then loop over individual slip systems
for slip_system in slip_system_group:
slip_dir_sample = self.ref_ori.conjugate.transform_vector(slip_system.slip_dir)
self.rdr[i].append(-slip_dir_sample[0] / slip_dir_sample[1])
@property
def slip_traces(self):
"""Returns list of slip trace angles.
Returns
-------
list
Slip trace angles based on grain orientation in calc_slip_traces.
"""
if self.slip_trace_angles is None:
self.calc_slip_traces()
return self.slip_trace_angles
[docs] def print_slip_traces(self):
"""Print a list of slip planes (with colours) and slip directions
"""
self.calc_slip_traces()
if self.average_schmid_factors is None:
raise Exception("Run 'calc_average_grain_schmid_factors' on the EBSD map first")
for ss_group, colour, sf_group, slip_trace in zip(
self.phase.slip_systems,
self.phase.slip_trace_colours,
self.average_schmid_factors,
self.slip_traces
):
print('{0}\tColour: {1}\tAngle: {2:.2f}'.format(ss_group[0].slip_plane_label, colour, slip_trace * 180 / np.pi))
for ss, sf in zip(ss_group, sf_group):
print(' {0} SF: {1:.3f}'.format(ss.slip_dir_label, sf))
[docs] def calc_slip_traces(self, slip_systems=None):
"""Calculates list of slip trace angles based on grain orientation.
Parameters
-------
slip_systems : defdap.crystal.SlipSystem, optional
"""
if slip_systems is None:
slip_systems = self.phase.slip_systems
if self.ref_ori is None:
self.calc_average_ori()
screen_plane_norm = np.array((0, 0, 1)) # in sample orientation frame
grain_av_ori = self.ref_ori # orientation of grain
screen_plane_norm_crystal = grain_av_ori.transform_vector(screen_plane_norm)
self.slip_trace_angles = []
self.slip_trace_inclinations = []
# Loop over each group of slip systems
for slip_system_group in slip_systems:
# Take slip plane from first in group
slip_plane_norm = slip_system_group[0].slip_plane
# planeLabel = slip_system_group[0].slip_plane_label
# Calculate intersection of slip plane with plane of screen
intersection_crystal = np.cross(screen_plane_norm_crystal, slip_plane_norm)
# Calculate angle between slip plane and screen plane
inclination = np.arccos(np.dot(screen_plane_norm_crystal, slip_plane_norm))
if inclination > np.pi / 2:
inclination = np.pi - inclination
# print("{} inclination: {:.1f}".format(planeLabel, inclination * 180 / np.pi))
# Transform intersection back into sample coordinates and normalise
intersection = grain_av_ori.conjugate.transform_vector(intersection_crystal)
intersection = intersection / np.sqrt(np.dot(intersection, intersection))
# Calculate trace angle. Starting vertical and proceeding
# counter clockwise
if intersection[0] > 0:
intersection *= -1
trace_angle = np.arccos(np.dot(intersection, np.array([0, 1.0, 0])))
# Append to list
self.slip_trace_angles.append(trace_angle)
self.slip_trace_inclinations.append(inclination)
[docs]class BoundarySet(object):
# boundaries : numpy.ndarray
# Map of boundaries. -1 for a boundary, 0 otherwise.
# phaseBoundaries : numpy.ndarray
# Map of phase boundaries. -1 for boundary, 0 otherwise.
def __init__(self, ebsd_map, points_x, points_y):
self.ebsd_map = ebsd_map
self.points_x = set(points_x)
self.points_y = set(points_y)
[docs] @classmethod
def from_image(cls, ebsd_map, image_x, image_y):
return cls(
ebsd_map,
zip(*image_x.transpose().nonzero()),
zip(*image_y.transpose().nonzero())
)
[docs] @classmethod
def from_boundary_segments(cls, b_segs):
points_x = []
points_y = []
for b_seg in b_segs:
points_x += b_seg.boundary_points_x
points_y += b_seg.boundary_points_y
return cls(b_segs[0].ebsdMap, points_x, points_y)
@property
def points(self):
return self.points_x.union(self.points_y)
def _image(self, points):
image = np.zeros(self.ebsd_map.shape, dtype=bool)
image[tuple(zip(*points))[::-1]] = True
return image
@property
def image_x(self):
return self._image(self.points_x)
@property
def image_y(self):
return self._image(self.points_y)
@property
def image(self):
return self._image(self.points)
@property
def lines(self):
_, _, lines = self.boundary_points_to_lines(
boundary_points_x=self.points_x,
boundary_points_y=self.points_y
)
return lines
[docs] @staticmethod
def boundary_points_to_lines(*, boundary_points_x=None,
boundary_points_y=None):
boundary_data = {}
if boundary_points_x is not None:
boundary_data['x'] = boundary_points_x
if boundary_points_y is not None:
boundary_data['y'] = boundary_points_y
if not boundary_data:
raise ValueError("No boundaries provided.")
deltas = {
'x': (0.5, -0.5, 0.5, 0.5),
'y': (-0.5, 0.5, 0.5, 0.5)
}
all_lines = []
for mode, points in boundary_data.items():
lines = []
for i, j in points:
lines.append((
(i + deltas[mode][0], j + deltas[mode][1]),
(i + deltas[mode][2], j + deltas[mode][3])
))
all_lines.append(lines)
if len(all_lines) == 2:
all_lines.append(all_lines[0] + all_lines[1])
return tuple(all_lines)
else:
return all_lines[0]
[docs]class BoundarySegment(object):
def __init__(self, ebsdMap, grain1, grain2):
self.ebsdMap = ebsdMap
self.grain1 = grain1
self.grain2 = grain2
# list of boundary points (x, y) for horizontal (X) and
# vertical (Y) boundaries
self.boundary_points_x = []
self.boundary_points_y = []
# Boolean value for each point above, True if boundary point is
# in grain1 and False if in grain2
self.boundary_point_owners_x = []
self.boundary_point_owners_y = []
def __eq__(self, right):
if type(self) is not type(right):
raise NotImplementedError()
return ((self.grain1 is right.grain1 and
self.grain2 is right.grain2) or
(self.grain1 is right.grain2 and
self.grain2 is right.grain1))
def __len__(self):
return len(self.boundary_points_x) + len(self.boundary_points_y)
[docs] def addBoundaryPoint(self, point, kind, owner_grain):
if kind == 0:
self.boundary_points_x.append(point)
self.boundary_point_owners_x.append(owner_grain is self.grain1)
elif kind == 1:
self.boundary_points_y.append(point)
self.boundary_point_owners_y.append(owner_grain is self.grain1)
else:
raise ValueError("Boundary point kind is 0 for x and 1 for y")
[docs] def boundary_point_pairs(self, kind):
"""Return pairs of points either side of the boundary. The first
point is always in grain1
"""
if kind == 0:
boundary_points = self.boundary_points_x
boundary_point_owners = self.boundary_point_owners_x
delta = (1, 0)
else:
boundary_points = self.boundary_points_y
boundary_point_owners = self.boundary_point_owners_y
delta = (0, 1)
boundary_point_pairs = []
for point, owner in zip(boundary_points, boundary_point_owners):
other_point = (point[0] + delta[0], point[1] + delta[1])
if owner:
boundary_point_pairs.append((point, other_point))
else:
boundary_point_pairs.append((other_point, point))
return boundary_point_pairs
@property
def boundary_point_pairs_x(self):
"""Return pairs of points either side of the boundary. The first
point is always in grain1
"""
return self.boundary_point_pairs(0)
@property
def boundary_point_pairs_y(self):
"""Return pairs of points either side of the boundary. The first
point is always in grain1
"""
return self.boundary_point_pairs(1)
@property
def boundary_lines(self):
"""Return line points along this boundary segment"""
_, _, lines = BoundarySet.boundary_points_to_lines(
boundary_points_x=self.boundary_points_x,
boundary_points_y=self.boundary_points_y
)
return lines
[docs] def misorientation(self):
mis_ori, minSymm = self.grain1.ref_ori.mis_ori(
self.grain2.ref_ori, self.ebsdMap.crystal_sym, return_quat=2
)
mis_ori = 2 * np.arccos(mis_ori)
mis_ori_axis = self.grain1.ref_ori.mis_ori_axis(minSymm)
# should this be a unit vector already?
mis_ori_axis /= np.sqrt(np.dot(mis_ori_axis, mis_ori_axis))
return mis_ori, mis_ori_axis
# compVector = np.array([1., 1., 1.])
# deviation = np.arccos(
# np.dot(mis_ori_axis, np.array([1., 1., 1.])) /
# (np.sqrt(np.dot(mis_ori_axis, mis_ori_axis) * np.dot(compVector,
# compVector))))
# print(deviation * 180 / np.pi)
[docs]class Linker(object):
"""Class for linking multiple EBSD maps of the same region for analysis of deformation.
Attributes
----------
ebsd_maps : list(ebsd.Map)
List of `ebsd.Map` objects that are linked.
links : list(tuple(int))
List of grain link. Each link is stored as a tuple of
grain IDs (one from each map stored in same order of maps).
plots : list(plotting.MapPlot)
List of last opened plot of each map.
"""
def __init__(self, ebsd_maps):
"""Initialise linker and set ebsd maps
Parameters
----------
ebsd_maps : list(ebsd.Map)
List of `ebsd.Map` objects that are linked.
"""
self.ebsd_maps = ebsd_maps
self.links = []
self.plots = None
[docs] def set_origin(self, **kwargs):
"""Interactive tool to set origin of each EBSD map.
Parameters
----------
kwargs
Keyword arguments passed to :func:`defdap.ebsd.Map.plot_default`
"""
self.plots = []
for ebsd_map in self.ebsd_maps:
plot = ebsd_map.plot_default(make_interactive=True, **kwargs)
plot.add_event_handler('button_press_event', self.click_set_origin)
plot.add_points([ebsd_map.origin[0]], [ebsd_map.origin[1]],
c='w', s=60, marker='x')
self.plots.append(plot)
[docs] def click_set_origin(self, event, plot):
"""Event handler for clicking to set origin of map.
Parameters
----------
event
Click event.
plot : defdap.plotting.MapPlot
Plot to capture clicks from.
"""
# check if click was on the map
if event.inaxes is not plot.ax:
return
origin = (int(event.xdata), int(event.ydata))
plot.calling_map.origin = origin
plot.add_points([origin[0]], [origin[1]], update_layer=0)
print(f"Origin set to ({origin[0]}, {origin[1]})")
[docs] def start_linking(self):
"""Start interactive grain linking process of each EBSD map.
"""
self.plots = []
for ebsd_map in self.ebsd_maps:
plot = ebsd_map.locate_grain(click_event=self.click_grain_guess)
# Add make link button to axes
plot.add_button('Make link', self.make_link,
color='0.85', hovercolor='0.95')
self.plots.append(plot)
[docs] def click_grain_guess(self, event, plot):
"""Guesses grain position in other maps, given click on one.
Parameters
----------
event
Click handler.
plot : defdap.plotting.Plot
Plot to capture clicks from.
"""
# check if click was on the map
if event.inaxes is not plot.ax:
return
curr_ebsd_map = plot.callingMap
if curr_ebsd_map is self.ebsd_maps[0]:
# clicked on 'master' map so highlight and guess grain on others
# set current grain in 'master' ebsd map
self.ebsd_maps[0].click_grain_id(event, plot, False)
# guess at grain in other maps
for ebsd_map, plot in zip(self.ebsd_maps[1:], self.plots[1:]):
# calculated position relative to set origin of the
# map, scaled from step size of maps
x0m = curr_ebsd_map.origin[0]
y0m = curr_ebsd_map.origin[1]
x0 = ebsd_map.origin[0]
y0 = ebsd_map.origin[1]
scaling = curr_ebsd_map.step_size / ebsd_map.step_size
x = int((event.xdata - x0m) * scaling + x0)
y = int((event.ydata - y0m) * scaling + y0)
grain_id = int(ebsd_map.data.grains[y, x]) - 1
grain = self[grain_id]
ebsd_map.sel_grain = grain
print(grain_id)
# update the grain highlights layer in the plot
plot.add_grain_highlights([grain_id],
alpha=ebsd_map.highlight_alpha)
else:
# clicked on other map so correct guessed selected grain
curr_ebsd_map.click_grain_id(event, plot, False)
[docs] def make_link(self, event, plot):
"""Make a link between the EBSD maps after clicking.
"""
# create empty list for link
curr_link = []
for i, ebsd_map in enumerate(self.ebsd_maps):
if ebsd_map.sel_grain is not None:
curr_link.append(ebsd_map.sel_grain.grain_id)
else:
raise Exception(f"No grain selected in map {i + 1}.")
curr_link = tuple(curr_link)
if curr_link not in self.links:
self.links.append(curr_link)
print("Link added " + str(curr_link))
[docs] def reset_links(self):
"""Reset links.
"""
self.links = []
# Analysis routines
[docs] def set_ref_ori_from_master(self):
"""Loop over each map (not first/reference) and each link.
Sets refOri of linked grains to refOri of grain in first map.
"""
for i, ebsd_map in enumerate(self.ebsd_maps[1:], start=1):
for link in self.links:
ebsd_map[link[i]].ref_ori = copy.deepcopy(
self.ebsd_maps[0][link[0]].ref_ori
)
[docs] def update_misori(self, calc_axis=False):
"""Recalculate misorientation for linked grain (not for first map)
Parameters
----------
calc_axis : bool
Calculate the misorientation axis if True.
"""
for i, ebsd_map in enumerate(self.ebsd_maps[1:], start=1):
for link in self.links:
ebsd_map[link[i]].build_mis_ori_list(calc_axis=calc_axis)