# Copyright 2021 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 defdap import plotting
from typing import Union, Tuple, List, Optional
[docs]class Quat(object):
"""Class used to define and perform operations on quaternions. These
are interpreted in the passive sense.
"""
__slots__ = ['quatCoef']
def __init__(self, *args, allow_southern: Optional[bool] = False) -> None:
"""
Construct a Quat object from 4 quat coefficients or an array of
quat coefficients.
Parameters
----------
*args
Variable length argument list.
allow_southern
if False, move quat to northern hemisphere.
"""
# construct with array of quat coefficients
if len(args) == 1:
if len(args[0]) != 4:
raise TypeError("Arrays input must have 4 elements")
self.quatCoef = np.array(args[0], dtype=float)
# construct with quat coefficients
elif len(args) == 4:
self.quatCoef = np.array(args, dtype=float)
else:
raise TypeError("Incorrect argument length. Input should be "
"an array of quat coefficients or idividual "
"quat coefficients")
# move to northern hemisphere
if not allow_southern and self.quatCoef[0] < 0:
self.quatCoef = self.quatCoef * -1
[docs] @classmethod
def fromEulerAngles(cls, ph1: float, phi: float, ph2: float) -> 'Quat':
"""Create a quat object from 3 Bunge euler angles.
Parameters
----------
ph1
First Euler angle, rotation around Z in radians.
phi
Second Euler angle, rotation around new X in radians.
ph2
Third Euler angle, rotation around new Z in radians.
Returns
-------
defdap.quat.Quat
Initialised Quat object.
"""
# calculate quat coefficients
quatCoef = np.array([
np.cos(phi / 2.0) * np.cos((ph1 + ph2) / 2.0),
-np.sin(phi / 2.0) * np.cos((ph1 - ph2) / 2.0),
-np.sin(phi / 2.0) * np.sin((ph1 - ph2) / 2.0),
-np.cos(phi / 2.0) * np.sin((ph1 + ph2) / 2.0)
], dtype=float)
# call constructor
return cls(quatCoef)
[docs] @classmethod
def fromAxisAngle(cls, axis: np.ndarray, angle: float) -> 'Quat':
"""Create a quat object from a rotation around an axis. This
creates a quaternion to represent the passive rotation (-ve axis).
Parameters
----------
axis
Axis that the rotation is applied around.
angle
Magnitude of rotation in radians.
Returns
-------
defdap.quat.Quat
Initialised Quat object.
"""
# normalise the axis vector
axis = np.array(axis)
axis = axis / np.sqrt(np.dot(axis, axis))
# calculate quat coefficients
quatCoef = np.zeros(4, dtype=float)
quatCoef[0] = np.cos(angle / 2)
quatCoef[1:4] = -np.sin(angle / 2) * axis
# call constructor
return cls(quatCoef)
[docs] def eulerAngles(self) -> np.ndarray:
"""Calculate the Euler angle representation for this rotation.
Returns
-------
eulers : numpy.ndarray, shape 3
Bunge euler angles (in radians).
References
----------
Melcher A. et al., 'Conversion of EBSD data by a quaternion
based algorithm to be used for grain structure simulations',
Technische Mechanik, 30(4)401 – 413
Rowenhorst D. et al., 'Consistent representations of and
conversions between 3D rotations',
Model. Simul. Mater. Sci. Eng., 23(8)
"""
eulers = np.empty(3, dtype=float)
q = self.quatCoef
q03 = q[0]**2 + q[3]**2
q12 = q[1]**2 + q[2]**2
chi = np.sqrt(q03 * q12)
if chi == 0 and q12 == 0:
eulers[0] = np.arctan2(-2 * q[0] * q[3], q[0]**2 - q[3]**2)
eulers[1] = 0
eulers[2] = 0
elif chi == 0 and q03 == 0:
eulers[0] = np.arctan2(2 * q[1] * q[2], q[1]**2 - q[2]**2)
eulers[1] = np.pi
eulers[2] = 0
else:
cosPh1 = (-q[0] * q[1] - q[2] * q[3]) / chi
sinPh1 = (-q[0] * q[2] + q[1] * q[3]) / chi
cosPhi = q[0]**2 + q[3]**2 - q[1]**2 - q[2]**2
sinPhi = 2 * chi
cosPh2 = (-q[0] * q[1] + q[2] * q[3]) / chi
sinPh2 = (q[1] * q[3] + q[0] * q[2]) / chi
eulers[0] = np.arctan2(sinPh1, cosPh1)
eulers[1] = np.arctan2(sinPhi, cosPhi)
eulers[2] = np.arctan2(sinPh2, cosPh2)
if eulers[0] < 0:
eulers[0] += 2 * np.pi
if eulers[2] < 0:
eulers[2] += 2 * np.pi
return eulers
[docs] def rotMatrix(self) -> np.ndarray:
"""Calculate the rotation matrix representation for this rotation.
Returns
-------
rotMatrix : numpy.ndarray, shape (3, 3)
Rotation matrix.
References
----------
Melcher A. et al., 'Conversion of EBSD data by a quaternion
based algorithm to be used for grain structure simulations',
Technische Mechanik, 30(4)401 – 413
Rowenhorst D. et al., 'Consistent representations of and
conversions between 3D rotations',
Model. Simul. Mater. Sci. Eng., 23(8)
"""
rotMatrix = np.empty((3, 3), dtype=float)
q = self.quatCoef
qbar = q[0]**2 - q[1]**2 - q[2]**2 - q[3]**2
rotMatrix[0, 0] = qbar + 2 * q[1]**2
rotMatrix[0, 1] = 2 * (q[1] * q[2] - q[0] * q[3])
rotMatrix[0, 2] = 2 * (q[1] * q[3] + q[0] * q[2])
rotMatrix[1, 0] = 2 * (q[1] * q[2] + q[0] * q[3])
rotMatrix[1, 1] = qbar + 2 * q[2]**2
rotMatrix[1, 2] = 2 * (q[2] * q[3] - q[0] * q[1])
rotMatrix[2, 0] = 2 * (q[1] * q[3] - q[0] * q[2])
rotMatrix[2, 1] = 2 * (q[2] * q[3] + q[0] * q[1])
rotMatrix[2, 2] = qbar + 2 * q[3]**2
return rotMatrix
# show components when the quat is printed
def __repr__(self) -> str:
return "[{:.4f}, {:.4f}, {:.4f}, {:.4f}]".format(*self.quatCoef)
def __str__(self) -> str:
return self.__repr__()
def __eq__(self, right: 'Quat') -> bool:
return (isinstance(right, type(self)) and
self.quatCoef.tolist() == right.quatCoef.tolist())
def __hash__(self) -> int:
return hash(tuple(self.quatCoef.tolist()))
def _plotIPF(
self,
direction: np.ndarray,
symGroup: str,
**kwargs
) -> 'plotting.PolePlot':
Quat.plotIPF([self], direction, symGroup, **kwargs)
# overload * operator for quaternion product and vector product
def __mul__(self, right: 'Quat', allow_southern: bool = False) -> 'Quat':
if isinstance(right, type(self)): # another quat
newQuatCoef = np.zeros(4, dtype=float)
newQuatCoef[0] = (
self.quatCoef[0] * right.quatCoef[0] -
np.dot(self.quatCoef[1:4], right.quatCoef[1:4])
)
newQuatCoef[1:4] = (
self.quatCoef[0] * right.quatCoef[1:4] +
right.quatCoef[0] * self.quatCoef[1:4] +
np.cross(self.quatCoef[1:4], right.quatCoef[1:4])
)
return Quat(newQuatCoef, allow_southern=allow_southern)
raise TypeError("{:} - {:}".format(type(self), type(right)))
[docs] def dot(self, right: 'Quat') -> float:
""" Calculate dot product between two quaternions.
Parameters
----------
right
Right hand quaternion.
Returns
-------
float
Dot product.
"""
if isinstance(right, type(self)):
return np.dot(self.quatCoef, right.quatCoef)
raise TypeError()
# overload + operator
def __add__(self, right: 'Quat') -> 'Quat':
if isinstance(right, type(self)):
return Quat(self.quatCoef + right.quatCoef)
raise TypeError()
# overload += operator
def __iadd__(self, right: 'Quat') -> 'Quat':
if isinstance(right, type(self)):
self.quatCoef += right.quatCoef
return self
raise TypeError()
# allow array like setting/getting of components
def __getitem__(self, key: int) -> float:
return self.quatCoef[key]
def __setitem__(self, key: int, value: float) -> None:
self.quatCoef[key] = value
[docs] def norm(self) -> float:
"""Calculate the norm of the quaternion.
Returns
-------
float
Norm of the quaternion.
"""
return np.sqrt(np.dot(self.quatCoef[0:4], self.quatCoef[0:4]))
[docs] def normalise(self) -> 'Quat':
""" Normalise the quaternion (turn it into an unit quaternion).
Returns
-------
defdap.quat.Quat
Normalised quaternion.
"""
self.quatCoef /= self.norm()
return
# also the inverse if this is a unit quaternion
@property
def conjugate(self) -> 'Quat':
"""Calculate the conjugate of the quaternion.
Returns
-------
defdap.quat.Quat
Conjugate of quaternion.
"""
return Quat(self[0], -self[1], -self[2], -self[3])
[docs] def misOri(
self,
right: 'Quat',
symGroup: str,
returnQuat: Optional[int] = 0
) -> Tuple[float, 'Quat']:
"""
Calculate misorientation angle between 2 orientations taking
into account the symmetries of the crystal structure.
Angle is 2*arccos(output).
Parameters
----------
right
Orientation to find misorientation to.
symGroup
Crystal type (cubic, hexagonal).
returnQuat
What to return: 0 for minimum misorientation, 1 for
symmetric equivalent with minimum misorientation, 2 for both.
Returns
-------
float
Minimum misorientation.
defdap.quat.Quat
Symmetric equivalent orientation with minimum misorientation.
"""
if isinstance(right, type(self)):
# looking for max of this as it is cos of misorientation angle
minMisOri = 0
# loop over symmetrically equivalent orientations
for sym in Quat.symEqv(symGroup):
quatSym = sym * right
currentMisOri = abs(self.dot(quatSym))
if currentMisOri > minMisOri: # keep if misorientation lower
minMisOri = currentMisOri
minQuatSym = quatSym
if returnQuat == 1:
return minQuatSym
elif returnQuat == 2:
return minMisOri, minQuatSym
else:
return minMisOri
raise TypeError("Input must be a quaternion.")
[docs] def misOriAxis(self, right: 'Quat') -> np.ndarray:
"""
Calculate misorientation axis between 2 orientations. This
does not consider symmetries of the crystal structure.
Parameters
----------
right : defdap.quat.Quat
Orientation to find misorientation axis to.
Returns
-------
numpy.ndarray, shape 3
Axis of misorientation.
"""
if isinstance(right, type(self)):
Dq = right * self.conjugate
Dq = Dq.quatCoef
misOriAxis = 2 * Dq[1:4] * np.arccos(Dq[0]) / np.sqrt(1 - Dq[0]**2)
return misOriAxis
raise TypeError("Input must be a quaternion.")
[docs] def plotIPF(
self,
direction: np.ndarray,
symGroup: str,
projection: Optional[str] = None,
plot: Optional['plotting.Plot'] = None,
fig: Optional['matplotlib.figure.Figure'] = None,
ax: Optional['matplotlib.axes.Axes'] = None,
plotColourBar: Optional[bool] = False,
clabel: Optional[str] = "",
makeInteractive: Optional[bool] = False,
markerColour: Optional[Union[List[str], str]] = None,
markerSize: Optional[float] = 40,
**kwargs
) -> 'plotting.PolePlot':
"""
Plot IPF of orientation, with relation to specified sample direction.
Parameters
----------
direction
Sample reference direction for IPF.
symGroup
Crystal type (cubic, hexagonal).
projection
Projection to use. Either string (stereographic or lambert)
or a function.
plot
Defdap plot to plot on.
fig
Figure to plot on, if not provided the current
active axis is used.
ax
Axis to plot on, if not provided the current
active axis is used.
makeInteractive
If true, make the plot interactive.
plotColourBar : bool
If true, plot a colour bar next to the map.
clabel : str
Label for the colour bar.
markerColour: str or list of str
Colour of markers (only used for half and half colouring,
otherwise use argument c).
markerSize
Size of markers (only used for half and half colouring,
otherwise use argument s).
kwargs
All other arguments are passed to :func:`defdap.plotting.PolePlot.addPoints`.
"""
plotParams = {'marker': '+'}
plotParams.update(kwargs)
# Works as an instance or static method on a list of Quats
if isinstance(self, Quat):
quats = [self]
else:
quats = self
alphaFund, betaFund = Quat.calcFundDirs(quats, direction, symGroup)
if plot is None:
plot = plotting.PolePlot(
"IPF", symGroup, projection=projection,
ax=ax, fig=fig, makeInteractive=makeInteractive
)
plot.addPoints(
alphaFund, betaFund,
markerColour=markerColour, markerSize=markerSize,
**plotParams
)
if plotColourBar:
plot.addColourBar(clabel)
return plot
[docs] def plotUnitCell(
self,
crystalStructure: 'defdap.crystal.CrystalStructure',
OI: Optional[bool] = True,
plot: Optional['plotting.CrystalPlot'] = None,
fig: Optional['matplotlib.figure.Figure'] = None,
ax: Optional['matplotlib.axes.Axes'] = None,
makeInteractive: Optional[bool] = False,
**kwargs
) -> 'plotting.CrystalPlot':
"""Plots a unit cell.
Parameters
----------
crystalStructure
Crystal structure.
OI
True if using oxford instruments system.
plot
Plot object to plot to.
fig
Figure to plot on, if not provided the current active axis is used.
ax
Axis to plot on, if not provided the current active axis is used.
makeInteractive
True to make the plot interactive.
kwargs
All other arguments are passed to :func:`defdap.plotting.CrystalPlot.addVerts`.
"""
# Set default plot parameters then update with any input
plotParams = {}
plotParams.update(kwargs)
# TODO: most of this should be moved to either the crystal or
# plotting module
vert = crystalStructure.vertices
faces = crystalStructure.faces
if crystalStructure.name == 'hexagonal':
szFac = 0.18
if OI:
# Add 30 degrees to phi2 for OI
eulerAngles = self.eulerAngles()
eulerAngles[2] += np.pi / 6
gg = Quat.fromEulerAngles(*eulerAngles).rotMatrix().T
else:
gg = self.rotMatrix().T
elif crystalStructure.name == 'cubic':
szFac = 0.25
gg = self.rotMatrix().T
# Rotate the lattice cell points
pts = np.matmul(gg, vert.T).T * szFac
# Plot unit cell
planes = []
for face in faces:
planes.append(pts[face, :])
if plot is None:
plot = plotting.CrystalPlot(
ax=ax, fig=fig, makeInteractive=makeInteractive
)
plot.ax.set_xlim3d(-0.15, 0.15)
plot.ax.set_ylim3d(-0.15, 0.15)
plot.ax.set_zlim3d(-0.15, 0.15)
plot.ax.view_init(azim=270, elev=90)
plot.ax._axis3don = False
plot.addVerts(planes, **plotParams)
return plot
# Static methods
[docs] @staticmethod
def createManyQuats(eulerArray: np.ndarray) -> np.ndarray:
"""Create a an array of quats from an array of Euler angles.
Parameters
----------
eulerArray
Array of Bunge Euler angles of shape 3 x n x ... x m.
Returns
-------
quats : numpy.ndarray(defdap.quat.Quat)
Array of quat objects of shape n x ... x m.
"""
ph1 = eulerArray[0]
phi = eulerArray[1]
ph2 = eulerArray[2]
oriShape = eulerArray.shape[1:]
quatComps = np.zeros((4,) + oriShape, dtype=float)
quatComps[0] = np.cos(phi / 2.0) * np.cos((ph1 + ph2) / 2.0)
quatComps[1] = -np.sin(phi / 2.0) * np.cos((ph1 - ph2) / 2.0)
quatComps[2] = -np.sin(phi / 2.0) * np.sin((ph1 - ph2) / 2.0)
quatComps[3] = -np.cos(phi / 2.0) * np.sin((ph1 + ph2) / 2.0)
quats = np.empty(oriShape, dtype=Quat)
for i, idx in enumerate(np.ndindex(oriShape)):
quats[idx] = Quat(quatComps[(slice(None),) + idx])
return quats
[docs] @staticmethod
def multiplyManyQuats(quats: List['Quat'], right: 'Quat') -> List['Quat']:
""" Multiply all quats in a list of quats, by a single quat.
Parameters
----------
quats
List of quats to be operated on.
right
Single quaternion to multiply with the list of quats.
Returns
-------
list(defdap.quat.Quat)
"""
quatArray = np.array([q.quatCoef for q in quats])
tempArray = np.zeros((len(quatArray),4), dtype=float)
tempArray[...,0] = ((quatArray[...,0] * right.quatCoef[0]) -
np.dot(quatArray[...,1:4], right.quatCoef[1:4]))
tempArray[...,1:4] = ((quatArray[...,0,None] * right.quatCoef[None,1:4]) +
(right.quatCoef[0] * quatArray[...,1:4]) +
np.cross(quatArray[...,1:4], right.quatCoef[1:4]))
return [Quat(coefs) for coefs in tempArray]
[docs] @staticmethod
def calcSymEqvs(
quats: np.ndarray,
symGroup: str,
dtype: Optional[type] = float
) -> np.ndarray:
"""Calculate all symmetrically equivalent quaternions of given quaternions.
Parameters
----------
quats : numpy.ndarray(defdap.quat.Quat)
Array of quat objects.
symGroup
Crystal type (cubic, hexagonal).
dtype
Data type used for calculation, defaults to `float`.
Returns
-------
quatComps: numpy.ndarray, shape: (numSym x 4 x numQuats)
Array containing all symmetrically equivalent quaternion components of input quaternions.
"""
syms = Quat.symEqv(symGroup)
quatComps = np.empty((len(syms), 4, len(quats)), dtype=dtype)
# store quat components in array
quatComps[0] = Quat.extract_quat_comps(quats)
# calculate symmetrical equivalents
for i, sym in enumerate(syms[1:], start=1):
# sym[i] * quat for all points (* is quaternion product)
quatComps[i, 0, :] = (
quatComps[0, 0, :] * sym[0] - quatComps[0, 1, :] * sym[1] -
quatComps[0, 2, :] * sym[2] - quatComps[0, 3, :] * sym[3])
quatComps[i, 1, :] = (
quatComps[0, 0, :] * sym[1] + quatComps[0, 1, :] * sym[0] -
quatComps[0, 2, :] * sym[3] + quatComps[0, 3, :] * sym[2])
quatComps[i, 2, :] = (
quatComps[0, 0, :] * sym[2] + quatComps[0, 2, :] * sym[0] -
quatComps[0, 3, :] * sym[1] + quatComps[0, 1, :] * sym[3])
quatComps[i, 3, :] = (
quatComps[0, 0, :] * sym[3] + quatComps[0, 3, :] * sym[0] -
quatComps[0, 1, :] * sym[2] + quatComps[0, 2, :] * sym[1])
# swap into positive hemisphere if required
quatComps[i, :, quatComps[i, 0, :] < 0] *= -1
return quatComps
[docs] @staticmethod
def calcAverageOri(
quatComps: np.ndarray
) -> 'Quat':
"""Calculate the average orientation of given quats.
Parameters
----------
quatComps : numpy.ndarray
Array containing all symmetrically equivalent quaternion components of given quaternions
(shape: numSym x 4 x numQuats), can be calculated with :func:`Quat.calcSymEqvs`.
Returns
-------
avOri : defdap.quat.Quat
Average orientation of input quaternions.
"""
avOri = np.copy(quatComps[0, :, 0])
currMisOris = np.empty(quatComps.shape[0])
for i in range(1, quatComps.shape[2]):
# calculate misorientation between current average and all
# symmetrical equivalents. Dot product of each symm quat in
# quatComps with refOri for point i
currMisOris[:] = abs(np.einsum(
"ij,j->i", quatComps[:, :, i], avOri
))
# find min misorientation with current average then add to it
maxIdx = np.argmax(currMisOris[:])
avOri += quatComps[maxIdx, :, i]
# Convert components back to a quat and normalise
avOri = Quat(avOri)
avOri.normalise()
return avOri
[docs] @staticmethod
def calcMisOri(
quatComps: np.ndarray,
refOri: 'Quat'
) -> Tuple[np.ndarray, 'Quat']:
"""Calculate the misorientation between the quaternions and a reference quaternion.
Parameters
----------
quatComps
Array containing all symmetrically equivalent quaternion components of given quaternions
(shape: numSym x 4 x numQuats), can be calculated from quats with :func:`Quat.calcSymEqvs` .
refOri
Reference orientation.
Returns
-------
minMisOris : numpy.ndarray, len numQuats
Minimum misorientation between quats and reference orientation.
minQuatComps : defdap.quat.Quat
Quaternion components describing minimum misorientation between quats and reference orientation.
"""
misOris = np.empty((quatComps.shape[0], quatComps.shape[2]))
# Dot product of each quat in quatComps with refOri
misOris[:, :] = abs(np.einsum("ijk,j->ik", quatComps, refOri.quatCoef))
maxIdxs0 = np.argmax(misOris, axis=0)
maxIdxs1 = np.arange(misOris.shape[1])
minMisOris = misOris[maxIdxs0, maxIdxs1]
minQuatComps = quatComps[maxIdxs0, :, maxIdxs1].transpose()
minMisOris[minMisOris > 1] = 1
return minMisOris, minQuatComps
[docs] @staticmethod
def polarAngles(x: np.ndarray, y: np.ndarray, z: np.ndarray):
"""Convert Cartesian coordinates to polar coordinates, for an
unit vector.
Parameters
----------
x : numpy.ndarray(float)
x coordinate.
y : numpy.ndarray(float)
y coordinate.
z : numpy.ndarray(float)
z coordinate.
Returns
-------
float, float
inclination angle and azimuthal angle (around z axis from x
in anticlockwise as per ISO).
"""
mod = np.sqrt(x**2 + y**2 + z**2)
x = x / mod
y = y / mod
z = z / mod
alpha = np.arccos(z)
beta = np.arctan2(y, x)
return alpha, beta
[docs] @staticmethod
def calcIPFcolours(
quats: np.ndarray,
direction: np.ndarray,
symGroup: str,
dtype: Optional[type] = np.float32
) -> np.ndarray:
"""
Calculate the RGB colours, based on the location of the given quats
on the fundamental region of the IPF for the sample direction specified.
Parameters
----------
quats : numpy.ndarray(defdap.quat.Quat)
Array of quat objects.
direction
Direction in sample space.
symGroup
Crystal type (cubic, hexagonal).
dtype
Data type to use for calculation.
Returns
-------
numpy.ndarray, shape (3, numQuats)
Array of rgb colours for each quat.
References
-------
Stephen Cluff (BYU) - IPF_rgbcalc.m subroutine in OpenXY
https://github.com/BYU-MicrostructureOfMaterials/OpenXY/blob/master/Code/PlotIPF.m
"""
numQuats = len(quats)
alphaFund, betaFund = Quat.calcFundDirs(
quats, direction, symGroup, dtype=dtype
)
# revert to cartesians
dirvec = np.empty((3, numQuats), dtype=dtype)
dirvec[0, :] = np.sin(alphaFund) * np.cos(betaFund)
dirvec[1, :] = np.sin(alphaFund) * np.sin(betaFund)
dirvec[2, :] = np.cos(alphaFund)
if symGroup == 'cubic':
poleDirections = np.array([[0, 0, 1],
[1, 0, 1]/np.sqrt(2),
[1, 1, 1]/np.sqrt(3)], dtype=dtype)
if symGroup == 'hexagonal':
poleDirections = np.array([[0, 0, 1],
[np.sqrt(3), 1, 0]/np.sqrt(4),
[1, 0, 0]], dtype=dtype)
rvect = np.broadcast_to(poleDirections[0].reshape((-1, 1)), (3, numQuats))
gvect = np.broadcast_to(poleDirections[1].reshape((-1, 1)), (3, numQuats))
bvect = np.broadcast_to(poleDirections[2].reshape((-1, 1)), (3, numQuats))
rgb = np.zeros((3, numQuats), dtype=dtype)
# Red Component
RDirPlane = np.cross(dirvec, rvect, axis=0)
GBplane = np.cross(bvect, gvect, axis=0)
Rintersect = np.cross(RDirPlane, GBplane, axis=0)
NORM = np.linalg.norm(Rintersect, axis=0, keepdims=True)
NORM[NORM == 0] = 1 #Prevent division by zero
Rintersect /= NORM
temp = np.arccos(np.clip(np.einsum("ij,ij->j", dirvec, Rintersect), -1, 1))
Rintersect[:, temp > (np.pi / 2)] *= -1
rgb[0, :] = np.divide(
np.arccos(np.clip(np.einsum("ij,ij->j", dirvec, Rintersect), -1, 1)),
np.arccos(np.clip(np.einsum("ij,ij->j", rvect, Rintersect), -1, 1))
)
# Green Component
GDirPlane = np.cross(dirvec, gvect, axis=0)
RBplane = np.cross(rvect, bvect, axis=0)
Gintersect = np.cross(GDirPlane, RBplane, axis=0)
NORM = np.linalg.norm(Gintersect, axis=0, keepdims=True)
NORM[NORM == 0] = 1 #Prevent division by zero
Gintersect /= NORM
temp = np.arccos(np.clip(np.einsum("ij,ij->j", dirvec, Gintersect), -1, 1))
Gintersect[:, temp > (np.pi / 2)] *= -1
rgb[1, :] = np.divide(
np.arccos(np.clip(np.einsum("ij,ij->j", dirvec, Gintersect), -1, 1)),
np.arccos(np.clip(np.einsum("ij,ij->j", gvect, Gintersect), -1, 1))
)
# Blue Component
BDirPlane = np.cross(dirvec, bvect, axis=0)
RGplane = np.cross(gvect, rvect, axis=0)
Bintersect = np.cross(BDirPlane, RGplane, axis=0)
NORM = np.linalg.norm(Bintersect, axis=0, keepdims=True)
NORM[NORM == 0] = 1 #Prevent division by zero
Bintersect /= NORM
temp = np.arccos(np.clip(np.einsum("ij,ij->j", dirvec, Bintersect), -1, 1))
Bintersect[:, temp > (np.pi / 2)] *= -1
rgb[2, :] = np.divide(
np.arccos(np.clip(np.einsum("ij,ij->j", dirvec, Bintersect), -1, 1)),
np.arccos(np.clip(np.einsum("ij,ij->j", bvect, Bintersect), -1, 1))
)
rgb /= np.amax(rgb, axis=0)
return rgb
[docs] @staticmethod
def calcFundDirs(
quats: np.ndarray,
direction: np.ndarray,
symGroup: str,
dtype: Optional[type] = float
) -> Tuple[np.ndarray, np.ndarray]:
"""
Transform the sample direction to crystal coords based on the quats
and find the ones in the fundamental sector of the IPF.
Parameters
----------
quats: array_like(defdap.quat.Quat)
Array of quat objects.
direction
Direction in sample space.
symGroup
Crystal type (cubic, hexagonal).
dtype
Data type to use for calculation.
Returns
-------
float, float
inclination angle and azimuthal angle (around z axis from x in anticlockwise).
"""
# convert direction to float array
direction = np.array(direction, dtype=dtype)
# get array of symmetry operations. shape - (numSym, 4, numQuats)
quatCompsSym = Quat.calcSymEqvs(quats, symGroup, dtype=dtype)
# array to store crystal directions for all orientations and symmetries
directionCrystal = np.empty(
(3, quatCompsSym.shape[0], quatCompsSym.shape[2]), dtype=dtype
)
# temp variables to use below
quatDotVec = (quatCompsSym[:, 1, :] * direction[0] +
quatCompsSym[:, 2, :] * direction[1] +
quatCompsSym[:, 3, :] * direction[2])
temp = (np.square(quatCompsSym[:, 0, :]) -
np.square(quatCompsSym[:, 1, :]) -
np.square(quatCompsSym[:, 2, :]) -
np.square(quatCompsSym[:, 3, :]))
# transform the pole direction to crystal coords for all
# orientations and symmetries
# (quatCompsSym * vectorQuat) * quatCompsSym.conjugate
directionCrystal[0, :, :] = (
2 * quatDotVec * quatCompsSym[:, 1, :] +
temp * direction[0] +
2 * quatCompsSym[:, 0, :] * (
quatCompsSym[:, 2, :] * direction[2] -
quatCompsSym[:, 3, :] * direction[1]
)
)
directionCrystal[1, :, :] = (
2 * quatDotVec * quatCompsSym[:, 2, :] +
temp * direction[1] +
2 * quatCompsSym[:, 0, :] * (
quatCompsSym[:, 3, :] * direction[0] -
quatCompsSym[:, 1, :] * direction[2]
)
)
directionCrystal[2, :, :] = (
2 * quatDotVec * quatCompsSym[:, 3, :] +
temp * direction[2] +
2 * quatCompsSym[:, 0, :] * (
quatCompsSym[:, 1, :] * direction[1] -
quatCompsSym[:, 2, :] * direction[0]
)
)
# normalise vectors
directionCrystal /= np.sqrt(np.einsum(
'ijk,ijk->jk', directionCrystal, directionCrystal
))
# move all vectors into north hemisphere
directionCrystal[:, directionCrystal[2, :, :] < 0] *= -1
# convert to spherical coordinates
alpha, beta = Quat.polarAngles(
directionCrystal[0], directionCrystal[1], directionCrystal[2]
)
# find the poles in the fundamental triangle
if symGroup == "cubic":
# first beta should be between 0 and 45 deg leaving 3
# symmetric equivalents per orientation
trialPoles = np.logical_and(beta >= 0, beta <= np.pi / 4)
# if less than 3 left need to expand search slighly to
# catch edge cases
if np.any(np.sum(trialPoles, axis=0) < 3):
deltaBeta = 1e-8
trialPoles = np.logical_and(beta >= -deltaBeta,
beta <= np.pi / 4 + deltaBeta)
# now of symmetric equivalents left we want the one with
# minimum alpha
min_alpha_idx = np.nanargmin(np.where(trialPoles==False, np.nan, alpha), axis=0)
betaFund = beta[min_alpha_idx, np.arange(len(min_alpha_idx))]
alphaFund = alpha[min_alpha_idx, np.arange(len(min_alpha_idx))]
elif symGroup == "hexagonal":
# first beta should be between 0 and 30 deg leaving 1
# symmetric equivalent per orientation
trialPoles = np.logical_and(beta >= 0, beta <= np.pi / 6)
# if less than 1 left need to expand search slighly to
# catch edge cases
if np.any(np.sum(trialPoles, axis=0) < 1):
deltaBeta = 1e-8
trialPoles = np.logical_and(beta >= -deltaBeta,
beta <= np.pi / 6 + deltaBeta)
# non-indexed points cause more than 1 symmetric equivalent, use this
# to pick one and filter non-indexed points later
first_idx = (trialPoles==True).argmax(axis=0)
betaFund = beta[first_idx, np.arange(len(first_idx))]
alphaFund = alpha[first_idx, np.arange(len(first_idx))]
else:
raise Exception("symGroup must be cubic or hexagonal")
return alphaFund, betaFund
[docs] @staticmethod
def symEqv(symGroup: str) -> List['Quat']:
"""Returns all symmetric equivalents for a given crystal type.
LEGACY: move to use symmetries defined in crystal structures
Parameters
----------
symGroup : str
Crystal type (cubic, hexagonal).
Returns
-------
list of defdap.quat.Quat
Symmetrically equivalent quats.
"""
# Dirty fix to stop circular dependency
from defdap.crystal import crystalStructures
try:
return crystalStructures[symGroup].symmetries
except KeyError:
# return just identity if unknown structure
return [Quat(1.0, 0.0, 0.0, 0.0)]