# --------------------------------------------------------------------------- #
# #
# Import Modules #
# #
# --------------------------------------------------------------------------- #
from __future__ import division
import warnings
import numpy as np
import scipy.spatial.distance
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pyquaternion import Quaternion
from microstructpy import _misc
from microstructpy.geometry.ellipse import Ellipse
__author__ = 'Kenneth (Kip) Hart'
# --------------------------------------------------------------------------- #
# #
# Ellipsoid Class #
# #
# --------------------------------------------------------------------------- #
[docs]class Ellipsoid(object):
"""A 3D Ellipsoid
This class contains the data and functions for a 3D ellispoid.
It is defined by its center, axes, and orientation.
If multiple keywords are given for the shape of the ellipsoid, there
is no guarantee for which keywords are used.
Args:
a (float): *(optional)* First semi-axis of ellipsoid. Default is 1.
b (float): *(optional)* Second semi-axis of ellipsoid. Default is 1.
c (float): *(optional)* Third semi-axis of ellipsoid. Default is 1.
center (list): *(optional)* The ellipsoid center.
Defaults to (0, 0, 0).
axes (list): *(optional)* List of 3 semi-axes.
Defaults to (1, 1, 1).
size (float): *(optional)* The diameter of a sphere with equal volume.
Defaults to 2.
ratio_ab (float): *(optional)* The ratio of a to b.
ratio_ac (float): *(optional)* The ratio of a to c.
ratio_bc (float): *(optional)* The ratio of b to c.
ratio_ba (float): *(optional)* The ratio of b to a.
ratio_ca (float): *(optional)* The ratio of c to a.
ratio_cb (float): *(optional)* The ratio of c to b.
rot_seq (list): *(optional)* List of rotations (deg). Each element of
the list should be an (axis, angle) tuple. The options for the
axis are: 'x', 'y', 'z', 1, 2, or 3.
For example::
rot_seq = [('x', 10), (2, -20), ('z', 85), ('x', 21)]
rot_seq_deg (list): *(optional)* Alias for ``rot_seq``, with degrees
stated explicitly.
rot_seq_rad (list): *(optional)* Same format as ``rot_seq``, except the
angles are expressed in radians.
matrix (numpy.ndarray): *(optional)* A 3x3 rotation matrix expressing
the orientation of the ellipsoid. Defaults to the identity.
position : *(optional)* Alias for ``center``.
orientation: *(optional)* Alias for ``matrix``.
"""
def __init__(self, **kwargs):
# Position
if 'center' in kwargs:
self.center = kwargs['center']
elif 'position' in kwargs:
self.center = kwargs['position']
else:
self.center = (0, 0, 0)
# Axes
self.a = None
self.b = None
self.c = None
ratio_ab = None
ratio_ac = None
ratio_bc = None
size = None
if 'a' in kwargs:
self.a = kwargs['a']
if 'b' in kwargs:
self.b = kwargs['b']
if 'c' in kwargs:
self.c = kwargs['c']
if 'ratio_ab' in kwargs:
ratio_ab = kwargs['ratio_ab']
elif 'ratio_ba' in kwargs:
ratio_ab = 1 / kwargs['ratio_ba']
if 'ratio_ac' in kwargs:
ratio_ac = kwargs['ratio_ac']
elif 'ratio_ca' in kwargs:
ratio_ac = 1 / kwargs['ratio_ca']
if 'ratio_bc' in kwargs:
ratio_bc = kwargs['ratio_bc']
elif 'ratio_cb' in kwargs:
ratio_bc = 1 / kwargs['ratio_cb']
if 'size' in kwargs:
size = kwargs['size']
elif 'volume' in kwargs:
size = 2 * np.cbrt(3 * kwargs['volume'] / (4 * np.pi))
if 'axes' in kwargs:
self.a, self.b, self.c = kwargs['axes']
if (self.a is not None) and (self.b is not None):
ratio_ab = self.a / self.b
if (self.a is not None) and (self.c is not None):
ratio_ac = self.a / self.c
if (self.b is not None) and (self.c is not None):
ratio_bc = self.b / self.c
if (ratio_ab is not None) and (ratio_bc is not None):
ratio_ac = ratio_ab * ratio_bc
if (ratio_ab is not None) and (ratio_ac is not None):
ratio_bc = ratio_ac / ratio_ab
if (ratio_ac is not None) and (ratio_bc is not None):
ratio_ab = ratio_ac / ratio_bc
for _ in range(2):
if self.a is None:
if (ratio_ab is not None) and (self.b is not None):
self.a = ratio_ab * self.b
elif (ratio_ac is not None) and (self.c is not None):
self.a = ratio_ac * self.c
if self.b is None:
if (ratio_ab is not None) and (self.a is not None):
self.b = self.a / ratio_ab
elif (ratio_bc is not None) and (self.c is not None):
self.b = ratio_bc * self.a
if self.c is None:
if (ratio_ac is not None) and (self.a is not None):
self.c = self.a / ratio_ac
elif (ratio_bc is not None) and (self.b is not None):
self.c = self.b / ratio_bc
if (self.a is not None) and (self.b is not None):
ratio_ab = self.a / self.b
if (self.a is not None) and (self.c is not None):
ratio_ac = self.a / self.c
if (self.b is not None) and (self.c is not None):
ratio_bc = self.b / self.c
if (ratio_ab is not None) and (ratio_bc is not None):
ratio_ac = ratio_ab * ratio_bc
if (ratio_ab is not None) and (ratio_ac is not None):
ratio_bc = ratio_ac / ratio_ab
if (ratio_ac is not None) and (ratio_bc is not None):
ratio_ab = ratio_ac / ratio_bc
if (size is not None):
r_eff = 0.5 * size
r3 = r_eff * r_eff * r_eff
if (self.a is not None) and (self.b is not None):
self.c = r3 / (self.a * self.b)
elif (self.b is not None) and (self.c is not None):
self.a = r3 / (self.b * self.c)
elif (self.a is not None) and (self.c is not None):
self.b = r3 / (self.a * self.c)
elif (self.a is not None) and (ratio_bc is not None):
r2 = r3 / self.a
# r2 = b * c
# r2 = c * ratio_bc * c = ratio_bc * c^2
# c = sqrt(r2 / ratio_bc) and b = ratio_bc * c
self.c = np.sqrt(r2 / ratio_bc)
self.b = ratio_bc * self.c
elif (self.b is not None) and (ratio_ac is not None):
r2 = r3 / self.b
self.c = np.sqrt(r2 / ratio_ac)
self.a = ratio_ac * self.c
elif (self.c is not None) and (ratio_ab is not None):
r2 = r3 / self.c
self.b = np.sqrt(r2 / ratio_ab)
self.a = ratio_ab * self.b
else:
# r3 = a * b * c
# r3 = a * (a / ratio_ab) * (a / ratio_ac)
# r3 * ratio_ab * ratio_ac = a^3
# a = r_eff * cbrt(ratio_ab * ratio_ac)
self.a = r_eff * np.cbrt(ratio_ab * ratio_ac)
self.b = self.a / ratio_ab
self.c = self.a / ratio_ac
if self.a is None:
self.a = 1
if self.b is None:
self.b = 1
if self.c is None:
self.c = 1
# Orientation
if 'rot_seq' in kwargs:
self.rot_seq = kwargs['rot_seq']
elif 'rot_seq_deg' in kwargs:
self.rot_seq = kwargs['rot_seq_deg']
elif 'rot_seq_rad' in kwargs:
rs_rad = kwargs['rot_seq_rad']
self.rot_seq = [(ax, 180 * ang / np.pi) for ax, ang in rs_rad]
elif ('matrix' in kwargs) or ('orientation' in kwargs):
# adapted from:
# https://www.learnopencv.com/rotation-matrix-to-euler-angles/
if 'matrix' in kwargs:
R = kwargs['matrix']
else:
R = kwargs['orientation']
sy = np.sqrt(R[0][0] * R[0][0] + R[1][0] * R[1][0])
if sy > 1e-6:
x = np.arctan2(R[2][1], R[2][2])
y = np.arctan2(-R[2][0], sy)
z = np.arctan2(R[1][0], R[0][0])
else:
x = np.arctan2(-R[1][2], R[1][1])
y = np.arctan2(-R[2][0], sy)
z = 0
axes = ('z', 'y', 'x')
angs = [180 * ang / np.pi for ang in (z, y, x)]
self.rot_seq = list(zip(axes, angs))
else:
self.rot_seq = []
# ----------------------------------------------------------------------- #
# Ellipsoid of Best Fit #
# ----------------------------------------------------------------------- #
[docs] def best_fit(self, points):
"""Find ellipsoid of best fit.
This function takes a list of 3D points and computes the ellipsoid of
best fit for the points. It uses a published algorithm to fit the
ellipsoid, then attempts to define the axes in such a way that they
most align with this ellipsoid's axes. [#turner]_
Args:
points (list): Points to fit ellipsoid
Returns:
Ellipsoid: The ellipsoid that best fits the points.
.. [#turner] Turner, D. A., Anderson, I. J., Mason, J. C., and Cox,
M. G., "An Algorithm for Fitting an Ellipsoid to Data," *National
Physical Laboratory*, 1999, The United Kingdom.
(http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.36.2773&rep=rep1&type=pdf)
""" # NOQA: E501
pts = np.array(points)
pts_mean = pts.mean(axis=0)
trans_pts = pts - pts_mean
x, y, z = trans_pts.T
L = np.zeros((len(x), 9), dtype='float')
L[:, 0] = x * x + y * y - 2 * z * z
L[:, 1] = x * x - 2 * y * y + z * z
L[:, 2] = 4 * x * y
L[:, 3] = 2 * x * z
L[:, 4] = 2 * y * z
L[:, 5] = x
L[:, 6] = y
L[:, 7] = z
L[:, 8] = 1
e = x * x + y * y + z * z
u, v, m, n, p, q, r, s, t = np.linalg.lstsq(L, e, rcond=None)[0]
axx = 1 - u - v
ayy = 1 - u + 2 * v
azz = 1 + 2 * u - v
axy = -4 * m
axz = -2 * n
ayz = -2 * p
ax = - q
ay = - r
az = - s
ac = - t
hom_mat = np.array([[axx, 0.5 * axy, 0.5 * axz, 0.5 * ax],
[0.5 * axy, ayy, 0.5 * ayz, 0.5 * ay],
[0.5 * axz, 0.5 * ayz, azz, 0.5 * az],
[0.5 * ax, 0.5 * ay, 0.5 * az, ac]])
cen = np.linalg.lstsq(-hom_mat[:3, :3], hom_mat[-1, :3], rcond=None)[0]
glbl_cen = cen + pts_mean
T = np.eye(4)
T[-1, :3] = cen
R = T.dot(hom_mat.dot(T.T))
evals, evecs = np.linalg.eigh(- R[:3, :3] / R[-1, -1])
if np.any(evals <= 0):
w_str = 'Some of the eigenvalues of the quadratic form for'
w_str += 'the ellipsoid are non-positive. Using absolute value, '
w_str += 'but fit may be poor.'
warnings.warn(w_str, RuntimeWarning)
# extract ellipsoid parameters
try:
# reorganize the axes
opt_axes = 1 / np.sqrt(np.abs(evals))
axes = np.zeros(3)
avail_axes = np.copy(opt_axes)
rearr = np.zeros((3, 3))
for ax_i, ax in enumerate(self.axes):
opt_i = np.argmin(np.abs(avail_axes - ax))
rearr[ax_i, opt_i] = 1
avail_axes[opt_i] = np.inf
axes[ax_i] = opt_axes[opt_i]
new_evecs = evecs.dot(rearr)
# settle the orientation ambiguity
ori_dot = self.matrix.T.dot(new_evecs)
vec_flip = np.sign(np.diagonal(ori_dot))
ori_matrix = new_evecs.dot(np.diag(vec_flip))
except AttributeError:
axes = 1 / np.sqrt(np.abs(evals))
ori_matrix = evecs
return type(self)(center=glbl_cen, axes=axes, matrix=ori_matrix)
# ----------------------------------------------------------------------- #
# String and Representation Functions #
# ----------------------------------------------------------------------- #
def __str__(self):
str_str = 'center: ' + str(tuple(self.center)) + '\n'
str_str += 'a: ' + str(self.a) + '\n'
str_str += 'b: ' + str(self.b) + '\n'
str_str += 'c: ' + str(self.c)
if len(self.rot_seq) > 0:
str_str += '\nrot_seq: ('
for i, (ax, ang) in enumerate(self.rot_seq):
str_str += '(' + str(ax) + ', ' + str(ang) + ')'
if i < len(self.rot_seq) - 1:
str_str += ', '
else:
str_str += ')'
return str_str
def __repr__(self):
repr_str = 'Ellipsoid('
repr_str += 'center=' + repr(tuple(self.center))
repr_str += ', axes=' + repr(tuple([self.a, self.b, self.c]))
rot_str = repr(tuple([tuple(r) for r in self.rot_seq]))
repr_str += ', rot_seq=' + rot_str
repr_str += ')'
return repr_str
# ----------------------------------------------------------------------- #
# Size and Orientation Getters #
# ----------------------------------------------------------------------- #
@property
def size(self):
"""float: diameter of equivalent volume sphere"""
return 2 * np.cbrt(self.a * self.b * self.c)
@property
def axes(self):
"""tuple: the 3 semi-axes of the ellipsoid"""
return self.a, self.b, self.c
@property
def ratio_ab(self):
"""float: ratio of x-axis length to y-axis length"""
return self.a / self.b
@property
def ratio_ba(self):
"""float: ratio of y-axis length to x-axis length"""
return self.b / self.a
@property
def ratio_ac(self):
"""float: ratio of x-axis length to z-axis length"""
return self.a / self.c
@property
def ratio_ca(self):
"""float: ratio of z-axis length to x-axis length"""
return self.c / self.a
@property
def ratio_bc(self):
"""float: ratio of y-axis length to z-axis length"""
return self.b / self.c
@property
def ratio_cb(self):
"""float: ratio of z-axis length to y-axis length"""
return self.c / self.b
@property
def rot_seq_deg(self):
"""list: rotation sequence, with angles in degrees"""
return self.rot_seq
@property
def rot_seq_rad(self):
"""list: rotation sequence, with angles in radiands"""
return [(ax, np.pi * ang / 180) for ax, ang in self.rot_seq]
@property
def matrix(self):
"""numpy.ndarray: A 3x3 rotation matrix"""
ax_dict = {1: 0, 2: 1, 3: 2, 'x': 0, 'y': 1, 'z': 2}
q = Quaternion()
for ax, ang in self.rot_seq_deg:
vec = np.eye(3)[ax_dict[ax]]
q *= Quaternion(axis=vec, degrees=ang)
return q.rotation_matrix
@property
def orientation(self):
"""numpy.ndarray: A 3x3 rotation matrix"""
return self.matrix
# ----------------------------------------------------------------------- #
# Quadratic Form Matrix #
# ----------------------------------------------------------------------- #
@property
def matrix_quadform(self):
"""numpy.ndarray: Matrix of the quadratic form"""
R = np.array(self.orientation)
scl_mat = np.diag(1 / (np.array(self.axes) * np.array(self.axes)))
return R.dot(scl_mat).dot(R.T)
@property
def matrix_quadeq(self):
"""numpy.ndarray: Matrix of the quadratic equation"""
A33 = self.matrix_quadform
grad_vec = - A33.dot(self.center)
cen_vec = np.array(self.center).reshape(-1, 1)
const = cen_vec.T.dot(A33.dot(cen_vec))[0, 0] - 1
quad_mat = np.zeros((4, 4))
quad_mat[:3, :3] = A33
quad_mat[-1, :3] = grad_vec.T
quad_mat[:3, -1] = grad_vec
quad_mat[-1, -1] = const
return quad_mat
@property
def coefficients(self):
"""tuple: coeffificients of equation,
:math:`(A, B, C, D, E, F, G, H, K, L)` in
:math:`Ax^2 + Bxy + Cxz + Dy^2 + Eyz + Fz^2 + Gx + Hy + Kz + L = 0`
"""
quad_eq = self.matrix_quadeq
A = quad_eq[0, 0]
B = 2 * quad_eq[0, 1]
C = 2 * quad_eq[0, 2]
D = quad_eq[1, 1]
E = 2 * quad_eq[1, 2]
F = quad_eq[2, 2]
G = 2 * quad_eq[3, 0]
H = 2 * quad_eq[3, 1]
K = 2 * quad_eq[3, 2]
L = quad_eq[3, 3]
return A, B, C, D, E, F, G, H, K, L
# ----------------------------------------------------------------------- #
# Number of Dimensions #
# ----------------------------------------------------------------------- #
@property
def n_dim(self):
"""int: number of dimensions, 3"""
return 3
# ----------------------------------------------------------------------- #
# Volume Property #
# ----------------------------------------------------------------------- #
@property
def volume(self):
""" float: volume of ellipsoid, :math:`V = \\frac{4}{3}\\pi a b c`"""
return 4 * np.pi * self.a * self.b * self.c / 3
[docs] @classmethod
def volume_expectation(cls, **kwargs):
r"""Expected value of volume.
This function computes the expected value for the volume of an
ellipsoid. The keyword arguments are the same as the input parameters
for the class, :class:`microstructpy.geometry.Ellipsoid`. The
values for these keywords can be either constants or distributions from
the SciPy :mod:`scipy.stats` module.
The expected value is computed by the following formula:
.. math::
\mathbb{E}[V] &= \mathbb{E}[\frac{4}{3}\pi A B C] \\
&= \frac{4}{3}\pi \mathbb{E}[A] \mathbb{E}[B] \mathbb{E}[C] \\
&= \frac{4}{3}\pi \mu_A \mu_B \mu_C
If the ellisoid is specified by size and aspect ratios, then the
expected volume is computed by:
.. math::
\mathbb{E}[V] &= \mathbb{E}[\frac{\pi}{6} S^3] \\
&= \frac{\pi}{6} (\mu_S^3 + 3 \mu_S \sigma_S^2 + \gamma_{1, S} \sigma_S^3)
If the ellipsoid is specified using a combination of semi-axes and
aspect ratios, then the expected volume is the mean of 1000 random
samples:
.. math::
\mathbb{E}[V] \approx \frac{1}{n} \sum_{i=1}^n V_i
where :math:`n=1000`.
Args:
**kwargs: Keyword arguments, see :class:`microstructpy.geometry.Ellipsoid`.
Returns:
float: Expected value of the volume of the sphere.
""" # NOQA: E501
# Check for size distribution
if 'size' in kwargs:
s_dist = kwargs['size']
if type(s_dist) in (float, int):
return 0.5 * np.pi * s_dist * s_dist * s_dist / 3
return 0.5 * np.pi * s_dist.moment(3) / 3
if 'volume' in kwargs:
v_dist = kwargs['volume']
try:
v_exp = v_dist.moment(1)
except AttributeError:
v_exp = v_dist
return v_exp
# check for a, b, and c distribution
try:
exp_vol = 4 * np.pi / 3
for kw in ('a', 'b', 'c'):
dist = kwargs[kw]
try:
mu = dist.moment(1)
except AttributeError:
mu = dist
exp_vol *= mu
return exp_vol
except KeyError:
pass
# Use Monte Carlo to determine expected volume
n_trials = 1000
kws = set(kwargs.keys()) - set(_misc.ori_kws)
total_vol = 0
for i in range(n_trials):
params = {}
for kw in kws:
try:
params[kw] = kwargs[kw].rvs()
except AttributeError:
params[kw] = kwargs[kw]
total_vol += Ellipsoid(**params).volume
avg_vol = total_vol / n_trials
return avg_vol
# ----------------------------------------------------------------------- #
# Bounding Circles #
# ----------------------------------------------------------------------- #
@property
def bound_max(self):
"""tuple: maximum bounding sphere, (x, y, z, r)"""
r = max(self.a, self.b, self.c)
return tuple(list(self.center) + [r])
@property
def bound_min(self):
"""tuple: minimum interior sphere, (x, y, z, r)"""
r = min(self.a, self.b, self.c)
return tuple(list(self.center) + [r])
# ----------------------------------------------------------------------- #
# Sphere Approximation of Ellipsoid #
# ----------------------------------------------------------------------- #
[docs] def approximate(self, x1=None):
"""Approximate Ellipsoid with Spheres
This function approximates the ellipsoid by a set of spheres.
It does so by approximating the x-z and y-z elliptical cross sections
with circles, then scaling those circles and promoting them to spheres.
See the documentation for
:meth:`microstructpy.geometry.Ellipse.approximate` for more details.
Args:
x1 (float): *(optional)* Center position of the first sphere.
Default is 0.75x the minimum semi-axis.
Returns:
numpy.ndarray: An Nx4 list of the (x, y, z, r) data of the spheres
that approximate the ellipsoid.
"""
if (self.a == self.b) and (self.a == self.c):
return np.array([self.bound_max])
if x1 is None:
x1 = 0.25 * min(self.axes)
# Perform approximation such that a > b > c
if (self.a >= self.b) and (self.b >= self.c):
a = self.a
b = self.b
c = self.c
inds = [0, 1, 2]
elif (self.a >= self.c) and (self.c >= self.b):
a = self.a
b = self.c
c = self.b
inds = [0, 2, 1]
elif (self.b >= self.a) and (self.a >= self.c):
a = self.b
b = self.a
c = self.c
inds = [1, 0, 2]
elif (self.c >= self.a) and (self.a >= self.b):
a = self.c
b = self.a
c = self.b
inds = [0, 2, 1]
else:
a = self.c
b = self.b
c = self.a
inds = [2, 1, 0]
# Prolate Ellipsoid
if np.isclose(b, c):
circs_xz = Ellipse(a=a, b=c).approximate()
spheres = np.insert(circs_xz, 1, 0, axis=1)[:, np.append(inds, 3)]
cens = np.array(self.center) + spheres[:, :-1].dot(self.matrix.T)
spheres[:, :-1] = cens
return spheres
# Points in First Octant
surface_pts = _ellipse_pts(a, b, c)
plane_mask = np.isclose(surface_pts[:, -1], 0)
# Circle centers
circs_xz = Ellipse(a=a, b=c).approximate(x1)
x_grid = np.append(circs_xz[circs_xz[:, 0] >= 0, 0][:-1], a)
circs_yz = Ellipse(a=b, b=c).approximate(x1)
y_grid = np.append(circs_yz[circs_yz[:, 0] >= 0, 0][:-1], b)
xx, yy = np.meshgrid(x_grid / a, y_grid / b)
# Conformal mapping
xx_mesh = a * xx * np.sqrt(1 - 0.5 * yy * yy)
yy_mesh = b * yy * np.sqrt(1 - 0.5 * xx * xx)
xx_inner = xx_mesh[:-1, :-1].flatten()
yy_inner = yy_mesh[:-1, :-1].flatten()
pts_inner = np.array([xx_inner, yy_inner, 0 * xx_inner]).T
# Interior spheres
inner_dists = scipy.spatial.distance.cdist(pts_inner, surface_pts)
nearest_ind = np.argmin(inner_dists, axis=1)
nearest_dist = inner_dists[np.arange(len(nearest_ind)), nearest_ind]
in_plane = plane_mask[nearest_ind]
cens_inner = pts_inner[~in_plane]
rads_inner = nearest_dist[~in_plane]
# Exterior spheres
xx_outer = np.concatenate((xx_mesh[:, -1], xx_mesh[-1, :-1]))
yy_outer = np.concatenate((yy_mesh[:, -1], yy_mesh[-1, :-1]))
u_outer = np.arctan2(yy_outer / b, xx_outer / a)
n_x = b * np.cos(u_outer)
n_y = a * np.sin(u_outer)
norm_vec = np.array([n_x, n_y]).T
norm_vec /= np.linalg.norm(norm_vec, axis=1).reshape(-1, 1)
rads_outer = c * c * np.sqrt(n_x * n_x + n_y * n_y) / (a * b)
rel_pos_outer = rads_outer.reshape(-1, 1) * norm_vec
cens_outer = np.array([xx_outer, yy_outer]).T - rel_pos_outer
cens_outer = np.hstack((cens_outer, 0 * xx_outer.reshape(-1, 1)))
# First octant spheres
cens = np.vstack((cens_inner, cens_outer))
rads = np.concatenate((rads_inner, rads_outer))
spheres = np.hstack((cens, rads.reshape(-1, 1)))
# All octants
for dim in range(3):
refl_spheres = np.copy(spheres)
mask = spheres[:, dim] > 0
refl_spheres[mask, dim] *= -1
spheres = np.vstack((spheres, refl_spheres[mask]))
spheres = spheres[:, np.append(inds, 3)]
cens = np.array(self.center) + spheres[:, :-1].dot(self.matrix.T)
spheres[:, :-1] = cens
return spheres
# ----------------------------------------------------------------------- #
# Plot Function #
# ----------------------------------------------------------------------- #
[docs] def plot(self, **kwargs):
"""Plot the ellipsoid.
This function uses the :meth:`mpl_toolkits.mplot3d.axes3d.Axes3D.plot_surface`
method to add an ellipsoid to the current axes. The keyword arguments
are passes through to the plot_surface function.
Args:
**kwargs (dict): Keyword arguments for matplotlib.
""" # NOQA: E501
if plt.gcf().axes:
ax = plt.gca()
else:
ax = plt.gcf().add_subplot(projection=Axes3D.name)
u = np.linspace(0, 2 * np.pi, 11)
cv = np.linspace(-1, 1, 12)
uu, cvv = np.meshgrid(u, cv)
svv = np.sin(np.arccos(cvv))
grid_shape = uu.shape
a = self.a
b = self.b
c = self.c
xp = a * np.cos(uu) * svv
yp = b * np.sin(uu) * svv
zp = c * cvv
pts = np.array([xp.reshape(1, -1).flatten(),
yp.reshape(1, -1).flatten(),
zp.reshape(1, -1).flatten()])
rot_pts = self.orientation.dot(pts)
xr = rot_pts[0].reshape(grid_shape)
yr = rot_pts[1].reshape(grid_shape)
zr = rot_pts[2].reshape(grid_shape)
xc, yc, zc = self.center
xx = xc + xr
yy = yc + yr
zz = zc + zr
mod_kwargs = {}
for key, val in kwargs.items():
if key == 'facecolors' and type(val) != list:
mod_kwargs['color'] = val
else:
mod_kwargs[key] = val
ax.plot_surface(xx, yy, zz, **mod_kwargs)
# ----------------------------------------------------------------------- #
# Limits #
# ----------------------------------------------------------------------- #
@property
def limits(self):
"""list: List of (lower, upper) bounds for the bounding box"""
if np.all(np.isclose(self.matrix, np.eye(3))):
ax = np.array(self.axes)
cen = np.array(self.center)
return [(x - r, x + r) for x, r in zip(cen, ax)]
n = 4
u = np.linspace(0, 2 * np.pi, 1 + 4 * n)
cv = np.linspace(-1, 1, 1 + 2 * n)
uu, cvv = np.meshgrid(u, cv)
svv = np.sin(np.arccos(cvv))
xp = self.a * np.cos(uu) * svv
yp = self.b * np.sin(uu) * svv
zp = self.c * cvv
pts = np.array([xp.flatten(), yp.flatten(), zp.flatten()])
r_pts = self.matrix.dot(pts)
lbs = r_pts.min(axis=-1) + np.array(self.center)
ubs = r_pts.max(axis=-1) + np.array(self.center)
return list(zip(lbs, ubs))
@property
def sample_limits(self):
"""list: List of (lower, upper) bounds for the sampling region"""
return self.limits
# ----------------------------------------------------------------------- #
# Within Test #
# ----------------------------------------------------------------------- #
[docs] def within(self, points):
"""Test if points are within ellipsoid.
This function tests whether a point or set of points are within the
ellipsoid. For the set of points, a list of booleans is returned to
indicate which points are within the ellipsoid.
Args:
points (list or numpy.ndarray): Point or list of points.
Returns:
bool or numpy.ndarray: Set to True for points in geometry.
"""
pts = np.array(points)
single_pt = pts.ndim == 1
if single_pt:
pts = pts.reshape(1, -1)
rel_pos = pts - np.array(self.center)
rot_pos = rel_pos.dot(self.orientation)
scl_pos = rot_pos / np.array(self.axes).reshape(1, -1)
sq_dist = np.sum(scl_pos * scl_pos, axis=-1)
mask = sq_dist <= 1
if single_pt:
return mask[0]
else:
return mask
# ----------------------------------------------------------------------- #
# Reflect #
# ----------------------------------------------------------------------- #
[docs] def reflect(self, points):
"""Reflect points across surface.
This function reflects a point or set of points across the surface
of the ellipsoid. Points at the center of the ellipsoid are not
reflected.
Args:
points (list or numpy.ndarray): Points to reflect.
Returns:
numpy.ndarray: Reflected points.
"""
pts = np.array(points)
single_pt = pts.ndim == 1
if single_pt:
pts = pts.reshape(1, -1)
rel_pos = pts - np.array(self.center)
rot_pos = rel_pos.dot(self.orientation)
scl_pos = rot_pos / np.array(self.axes).reshape(1, -1)
dist = np.sqrt(np.sum(scl_pos * scl_pos, axis=-1))
mask = dist > 0
if not np.any(mask):
return np.array([])
new_dist = 2 - dist[mask]
scl = new_dist / dist[mask]
new_scl_pos = scl_pos[mask] * scl
new_rel_pos = new_scl_pos.dot(self.orientation.T)
new_pos = new_rel_pos + np.array(self.center)
if single_pt:
return new_pos[0]
else:
return new_pos
def _ellipse_arc(a, b, n):
horiz_x = a * np.linspace(1, -1, n)
verti_y = np.linspace(0, 2 * b, n)
t_denom = (a - horiz_x) * verti_y + 4 * a * b
t_numer = 4 * a * b
t = t_numer / t_denom
x_arc = (2 * t - 1) * a
y_arc = t * verti_y
return np.array([x_arc, y_arc]).T
def _ellipse_pts(a, b, c, n=81):
arc_xy = _ellipse_arc(a, b, n)
arc_xz = _ellipse_arc(a, c, n)
arc_yz = _ellipse_arc(b, c, n)
ii_xy, ii_z = np.meshgrid(np.arange(n - 1), np.arange(n - 1))
theta_xy = np.arctan2(arc_xy[:, 1] / b, arc_xy[:, 0] / a)
tt_xy = theta_xy[ii_xy]
lat_xz = np.arctan2(arc_xz[:, 1] / c, arc_xz[:, 0] / a)
lat_yz = np.arctan2(arc_yz[:, 1] / c, arc_yz[:, 0] / b)
ll_f_xz = lat_xz[ii_z] * np.cos(tt_xy)
ll_f_yz = lat_yz[ii_z] * np.sin(tt_xy)
latlat = np.sqrt(ll_f_xz * ll_f_xz + ll_f_yz * ll_f_yz)
xx = arc_xy[:, 0][ii_xy] * np.cos(latlat)
yy = arc_xy[:, 1][ii_xy] * np.cos(latlat)
zz = c * np.sin(latlat)
pts = np.array([xx.flatten(), yy.flatten(), zz.flatten()]).T
pts = np.append(pts, [[0, 0, c]], axis=0)
return pts